A global spatially contiguous solar-induced fluorescence (CSIF) dataset using neural networks

,

Obtaining a spatiotemporal continuous photosynthetic carbon fixation or gross primary production (GPP) dataset is crucial to food security, ecosystem service and health evaluation, and global carbon cycle studies (Beer et al., 2010).However, this is not possible without remote sensing data, since in situ carbon flux measurements, such as FLUXNET (Baldocchi et al., 2001), are usually costly and have limited spatial and temporal coverage (Schimel et al., 2015).Many remote-sensing-based productivity efficiency models (PEMs) have been built, but the model structure and parameterizations differ from each other and the performance of most models is not satisfactory in terms of simulated interannual variability and trends (Anav et al., 2015;Chen et al., 2017).Müller (1874) found that the chlorophyll fluorescence (Chl F) from a dilute chlorophyll solution was much stronger than the Chl F from a green leaf, suggesting that an alternative energy pathway exists for leaves in vivo.In the 1980s, scientists found that plant photosynthesis and heat dissipation are two alternatives to quench the excited chlorophyll, and there is a close linkage between Chl F and carbon assimilation rate (Genty et al., 1989;Krause and Weis, 1991).Leaf-level photosynthesis (A leaf ) and fluorescence (Chl F) share the same source of energy originating from photosynthetically active radiation (PAR) absorbed by chlorophyll (APAR chl ), which can be written using a light-use efficiency (LUE) approach (Monteith, 1972): where ϕ F and ϕ P represent the efficiencies for Chl F emission and photochemistry, respectively.fPAR chl , being different from the conventional definition of fraction of photosynthetically active radiation absorption, only considers the fractions absorbed by chlorophyll pigments where the photosynthesis and fluorescence originate (Zhang et al., 2018c).However, Chl F measurements have been mostly conducted at the leaf level, using pulse amplitude modulation (PAM) fluorometers (Porcar-Castell et al., 2008;Roháček and Barták, 1999).In this case, the measured Chl F intensity is not induced by the Sun but by the modulated light source.Although the absolute value of the Chl F intensity does not directly link to A leaf , it can still be used to calculate the fluorescence yield and investigate the reaction mechanism of the energy partitioning during the light reaction, and to calculate the quantum yield for photochemistry or as a tool to detect plant reactions under stress (Adams and Demmig-Adams, 2004;Flexas et al., 2002).
The successful retrieval of solar-induced (steady-state) chlorophyll fluorescence (SIF) from satellites has made it possible for vegetation photosynthetic activities to be observed at the global scale (Frankenberg et al., 2011;Guanter et al., 2012;Joiner et al., 2011Joiner et al., , 2013)).Satellite SIF can be expressed as a function similar to the Chl F at the leaf level but with extra terms considering the radiative transfer within the canopy and through the atmosphere (Joiner et al., 2014): where the satellite-retrieved SIF (SIF sat ), fluorescence yield (Θ F ), f esc and τ atm are all functions of the wavelength (λ); in addition, f esc and τ atm are also affected by sun-sensor geometry characterized by Sun zenith angle (SZA; θ s ), view zenith angle (θ v ) and relative azimuth angle (ϕ).f esc is a factor describing how much SIF emitted by the chloroplast leaves the canopy, and τ atm is a function of atmospheric optical depth, which indicates how much SIF that leaves the canopy top passes through the atmosphere before it is captured by the satellite sensors.It should be noted that the fraction of PAR for fluorescence (fPAR F ) may have a different activation spectrum than that for photosynthesis (fPAR chl ), but this difference is ignored here for simplicity.Although additional factors come into play during this process, satelliteretrieved SIF shows high consistency with GPP using both model simulations and ground-based measurements from eddy covariance (EC) flux towers, at least at the monthly timescale (Guanter et al., 2014;Li et al., 2018a;Zhang et al., 2016c, b).In addition, recent studies suggest that the GPP-SIF relationship is consistent across biome types (Sun et al., 2017).This finding, if valid across all biomes, would greatly benefit the usage of SIF for model benchmarking (Luo et al., 2012) and global GPP estimation.
However, several issues hinder exploring the relationship between SIF and in situ GPP estimates.Since the SIF signal is very small and sensors used to retrieve SIF were not initially built to estimate SIF, the satellite-retrieved SIF usually has a large footprint and large uncertainties in individual retrievals (Frankenberg et al., 2014;Joiner et al., 2013Joiner et al., , 2016)).For instance, the SIF retrieval from the Global Ozone Monitoring Experiment-2 (GOME-2) has a footprint of 40 km×40 km or larger, and the SIF from the Greenhouse gases Observing SATellite (GOSAT) has a circular footprint with 10.5 km in diameter.Direct comparison between the satellite-retrieved SIF signal and GPP estimates from EC flux tower sites thus faces the problem of spatial inconsistency except in areas of large homogenous landscape, e.g., the US Midwest cropland (Zhang et al., 2014) or boreal evergreen forests (Walther et al., 2016).However, corn (C4 pathway) and soybean (C3 pathway) in SIF footprints have different electron use efficiencies (Guan et al., 2016), which should affect the relationship between SIF and GPP.The low precision of SIF measurements also leads to a need for averaging multiple pixels either in space or time before being used.SIF retrieved from the Orbiting Carbon Observatory-2 (OCO-2) satellite partially solved this issue with a much smaller footprint size (1.3 km×2.25 km), higher signal-to-noise ratio compared to GOSAT (relatively higher SIF retrieval accuracy) and much larger numbers of observations per day (Frankenberg et al., 2014;Sun et al., 2018).However, due to the sparse sampling strategy and long revisit cycle, the OCO-2 SIF data have large gaps between nearby swaths, and the average sampling frequency for each flux tower site is only 3.21 year −1 during 2015-2016 (Lu et al., 2018).In addition, OCO-2 is often aggregated to a monthly dataset at relatively coarse spatial resolution, typically at 1∘×1∘, which limits its application in small regions.Although several statistical methods have been proposed to downscale satellite observations to finer spatial-temporal resolutions (Tadić et al., 2015(Tadić et al., , 2017)), considering the large land surface heterogeneity and wide gaps between OCO-2 swaths (∼100 km), it could be challenging to apply these methods to OCO-2 SIF.
A high spatiotemporal resolution SIF dataset is needed to improve our understanding of the relationship between SIF and GPP and provide accurate GPP estimates at the global scale.As discussed previously, the satelliteobserved SIF contains signals from APAR chl , fluorescence yield, and canopy and atmospheric attenuation.APAR chl is considered to be the first-order approximation of SIF as it exhibits high correlation with SIF at the canopy scale (Du et al., 2017;Rossini et al., 2016;Verrelst et al., 2015;Zhang et al., 2018c).Previous studies have shown that fPAR chl can be inversely estimated using the surface reflectances and radiative transfer models (Zhang et al., 2005(Zhang et al., , 2016a)).The canopy structure information that affects the SIF reabsorption within canopy is also embedded in the near-infrared reflectance (Badgley et al., 2017;Knyazikhin et al., 2013;Yang and van der Tol, 2018).Many previous studies have shown high correlation between SIF and vegetation indices (VIs), especially VIs related to the chlorophyll concentration (Frankenberg et al., 2011;Guanter et al., 2012).Therefore, broadband surface reflectances may have the potential to be used to estimate vegetation information and reconstruct global SIF (Duveiller and Cescatti, 2016;Gentine and Alemohammad, 2018a).However, physical models that can predict SIF (e.g., the Soil Canopy Observation, Photochemistry and Energy fluxes, SCOPE; van der Tol et al., 2009) often require many parameters, making it difficult to use reflectance and modeling to predict SIF at a larger scale.
Neural networks (NNs), together with many other machine learning algorithms, have been used with remote sensing datasets in the Earth sciences, especially for carbon and water fluxes estimation (Alemohammad et al., 2017;Jung et al., 2011;Tramontana et al., 2016), land cover mapping (Kussul et al., 2017;Zhu et al., 2017), soil moisture retrievals and downscaling (Alemohammad et al., 2018;Kolassa et al., 2018) or to bypass parameterization (Gentine et al., 2018).These studies mostly attempted to link the satellite signals with limited in situ observation or model simulations for model training, while taking advantage of the large amount of data in remote sensing observations; they applied the trained algorithm to generate a regional or global dataset.Reconstructing SIF from surface reflectance, on the other hand, uses no in situ observations but faces more problems related to the satellite data quality assurance.The SIF-reflectance relationship is complicated, and the NN benefits from the fact that an explicit physical and radiative transfer relationship is not required.
In this study, we aim to generate a global contiguous SIF (CSIF) product based on the SIF retrievals from OCO-2 and surface reflectances from Moderate-resolution Imaging Spectroradiometer (MODIS) aboard the Terra and Aqua satellites.The CSIF dataset aims to fill the spatial gaps between the OCO-2 swaths and temporal gaps due to the long revisit cycle of OCO-2.Specifically, we first trained and validated the NN using the satelliteobserved instantaneous SIF under clear-sky conditions so that the relationship is not affected by cloud-related artifacts.We further generated two SIF products, namely the clear-sky instantaneous SIF (CSIF clear-inst ) and the all-sky daily SIF (CSIF all-daily ).The spatiotemporal variations of these CSIF products were analyzed and compared with SIF from OCO-2 and three other GOME-2 SIF datasets.Finally, we showed two applications of CSIF datasets: (1) monitoring drought impact using CSIF clear-inst and OCO-2 SIF; (2) evaluating the GPP-SIF relationship by comparing CSIF with GPP estimates from 40 flux tower sites.
2 Materials and methods

OCO-2 solar-induced chlorophyll fluorescence dataset
The 8100r OCO-2 SIF data between September 2014 and December 2017 were used for NN training and evaluation (Frankenberg, 2015;Frankenberg et al., 2014;Sun et al., 2018).The daily sounding-based SIF retrievals at 757 nm were first aggregated to 0.05 ∘ (around 5.6 km×5.6 km at the Equator), consistent with MODIS Climate Modeling Grid (CMG) resolution.The reasons for using this resolution include the following: (1) it is directly comparable (of the same order of magnitude) to the OCO-2 SIF footprint size (around 1.3 km×2.25 km) and the samples within each grid cell can be more evenly distributed and thus more representative of the grid cell SIF values than using much coarser 1∘×1∘ or 2∘×2∘ grids; (2) by averaging multiple observations, the uncertainty in the SIF signal can be approximately reduced by a factor of n−−√ (n is the number of observations within this grid cell), assuming independent estimates and homogeneous SIF value within each grid cell (Frankenberg et al., 2014).During this aggregation, we only used cloud-free observations indicated by the OCO-2 cloud flag.For each 0.05 ∘ grid cell, the SIF value was only calculated when it contained more than five cloud-free SIF soundings.Although several studies have shown that SIF at different wavelengths has different sensitivity to stress and leaf and canopy reabsorption (Porcar-Castell et al., 2014;Rossini et al., 2015Rossini et al., , 2016)), we only use SIF at 757 nm since it showed superior performance to SIF at 771 nm in predicting GPP (Li et al., 2018a).The years 2015 and 2016 were used for training and 2014 and 2017 were used for validation.Altogether, 2 947 819 SIF grid cells passed quality check during 2014-2017.Figure 1 shows the spatial distribution of the SIF grid cells used for training and validation (test).It should be noted that the OCO-2 satellite started obtaining data from September 2014 and experienced some malfunctioning during August and September in 2017, causing lower coverage for validation samples in boreal regions.
In addition to these cloud-free observations, we also calculated the all-sky SIF at 0.05 ∘ resolution.All SIF retrievals that passed the suggested quality checks (documented in detailed by Sun et al., 2018) were used for the aggregation.The aggregated all-sky instantaneous SIF retrievals were converted to daily values based on the solar zenith angle (Zhang et al., 2018a).We used this dataset to validate the all-sky daily SIF (CSIF all-daily ) (see Sect. 2.5).In both cloud-free and all-sky aggregations, only observations from the nadir mode were used since glint mode tends to underestimate SIF (Sun et al., 2018).

MODIS reflectance dataset (MCD43C4 V006)
We used the 0.05 ∘ daily nadir bidirectional reflectance distribution adjusted reflectance (NBAR) product from MODIS (MCD43C4 V006) during 2000-2017 as input variables for the NN.The NBAR product computed the reflectance at a nadir viewing angle for each pixel at local solar noon.Compared to MOD09 or MYD09 surface reflectance product, it removed the angle effects and therefore should be more stable and consistent (Schaaf et al., 2002).This dataset was processed in two different ways for training and prediction.For the training process, following Gentine and Alemohammad (2018a), we extracted the reflectance from the first four bands of MODIS (centered at 645, 858, 469 and 555 nm, respectively) for the corresponding pixels and days when the cloud-free SIF observations were obtained.It should be noted that although the MCD43C4 is generated for each day and can match the daily SIF observations, the MCD43C4 NBAR uses 16 days worth of inputs and so that the reflectances includes the information on other days than the day of interest.However, we consider this to have limited effects since (1) the vegetation growth/changes are continuous in time, (2) the NBAR product uses 16-day data but also emphasizes the specific day of interest (Schaaf, 2018).These four bands were selected because the visible and near-infrared bands included most of the vegetation information and drives the variation of SIF (Verrelst et al., 2015).We also tested using all seven bands with/without the meteorological variables (temperature and vapor pressure deficit, obtained from the OCO-2 SIF lite files) to train the NN, but the improvements in training and validation were very minor (R 2 increased by less than 0.01; data not shown), and thus we decided not to use it.Since SIF is very sensitive to the incoming solar radiation, using cloud-free training samples can minimize the uncertainty of using cosine of the solar zenith angle as the proxy of incoming PAR.It should be noted that the training dataset may contain snow-affected samples, but these were not removed to get a more realistic prediction of SIF during winter.
For prediction, we first aggregated the daily reflectance to 4 days.The 4-day temporal resolution is selected to reach a balance among application requirements, information redundancy and dataset sizes.During this process, we used a gap-filling and smoothing algorithm to reconstruct the surface reflectance for the four bands.The detailed description of the gapfilling algorithm can be found in Zhang et al. (2017a).In this study, we slightly modified the algorithm by not applying the best index slope extraction (BISE) algorithm and Savitzky-Golay (SG) filter.The reconstructed 4-day 0.05 ∘ reflectance together with other datasets allowed us to predict SIF at 4-day 0.05 ∘ resolution during 2000-2017.Since this processing does not involve any extra information and only uses the reflectance observations from the successful model inversion, it should be comparable to the reflectance used for NN training.

Machine learning algorithms
A feed-forward NN is a number of computational nodes (called neurons) structured in a single or multi-layer architecture.Each neuron is connected with all neurons in the previous layer and next layer.The neuron values are calculated using an activation function with a pre-activated value, i.e., the weighted sum of all neurons in previous layer plus biases.The training of the NN attempts to optimize these weights and biases so that the differences between the output variable in the training data and NN prediction are minimized.In this study, we used Tensorflow (https://www.tensorflow.org,last access: 27 September 2018) and built feed-forward networks with one to three layers and two to nine neurons for each layer.After training models with data from 2015 and 2016, we validated the models using the test dataset from the years 2014 and 2017.We then picked the one with best performance and simplest structure for SIF prediction.The rectified linear unit (ReLU) was used as the activation function since it has shown better performance in our application, and the cost function used is the root mean square error (RMSE).We used 50 epochs with a batch size of 1024.Before training, each variable was normalized by its mean and standardized deviation.Since the NN is not deep and there is no sign of overfitting, we did not use any regularization methods during the training.

Reconstructing the clear-sky instantaneous SIF and daily SIF
During the NN training process, we only used the SIF and reflectance data in clear-sky conditions, and therefore cos(SZA) was used as a proxy of the incoming photosynthetically active radiation at top of canopy.In the prediction process, we also used the calculated cos(SZA) based on the satellite overpass local solar time and latitude.Since we did not consider the cloud and aerosol attenuation of the PAR, this product was referred to as the "clear-sky instantaneous SIF (CSIF clear-inst )".
In addition to the clear-sky instantaneous SIF, we also calculated two daily SIF data by assuming that the incoming solar radiation is the only factor that drives the diurnal cycle (Zhang et al., 2018a).All-sky daily SIF (CSIF all-daily ) can be calculated using the clear-sky top-of-canopy radiation (PAR clear-inst ) and the daily average radiation from the Breathing Earth System Simulator (BESS) (Ryu et al., 2018): where PAR clear−inst was calculated following previous studies that only considered atmospheric scattering (see Appendix A1).Clear-sky daily SIF (CSIF clear-daily ) assumes no cloud throughout the day and can be calculated by multiplying CSIF clear-inst with a daily correction factor (γ) (Zhang et al., 2018a): γ is calculated as the ratio between the cos(SZA) during the satellite overpass and the daily averaged cos(SZA).
2.5 GOME-2 SIF (SIF GOME-2 ), reconstructed SIF from GOME-2 (RSIF GOME-2 ) and SIF * datasets In this study, we also used the GOME-2 SIF (SIF GOME-2 ), reconstructed SIF from GOME-2 (RSIF GOME-2 ) using machine learning and the SIF * dataset in comparison with our contiguous SIF from OCO-2.The GOME-2 SIF V27 was retrieved using a principle component analysis algorithm in the wavelength range 734-758 nm (Joiner et al., 2013(Joiner et al., , 2016)).The V27 version, compared to the widely used V26, provides daily correction factor and improved bias correction and calibration (https://avdc.gsfc.nasa.gov/,last access: 27 September 2018).The level-3 monthly 0.5 ∘ daily average SIF was used to compare with CSIF all-daily .RSIF GOME-2 (Gentine and Alemohammad, 2018a) uses a similar machine learning technique approach to CSIF but the training is based on the biweekly gridded SIF product from GOME-2 and the 8-day MYD09A1 reflectance dataset.Both clear-sky and cloudy-sky SIF are used for NN training.This dataset has a spatial resolution of 0.05 ∘ and 8-day temporal resolution.Both RSIF GOME-2 and CSIF all-daily were aggregated to the 0.5 ∘ and semi-monthly to facilitate the comparison.
The SIF * dataset (Duveiller and Cescatti, 2016) applies a statistical method and calibrates a model that links monthly 0.5 ∘ SIF to the normalized difference vegetation index (NDVI), evapotranspiration (ET) and land surface temperature (LST) dataset for each moving window.The model and its spatiotemporally varied parameters were then applied to finer resolution dataset (NDVI, ET, LST) with a weighted average to generate SIF at 0.05 ∘ resolution.In this study, we used the 0.5 ∘ monthly SIF * dataset during 2007-2013 to compare with CSIF.

Comparing CSIF with GPP at flux tower sites
We further compared the CSIF dataset to GPP estimates from the tier 1 FLUXNET2015 datasets (http://fluxnet.fluxdata.org,last access: 27 September 2018) to investigate the SIF-GPP relationship.Since the CSIF dataset is continuous in space and time, it provides many more samples pairs compared to the original OCO-2 SIF data (Lu et al., 2018).However, because of the landscape heterogeneity and inconsistency between the flux tower footprint and CSIF pixel size, a rigorous site selection is needed.We took the vegetation growth condition into consideration during this process: (1) the annual average, minimum, maximum and seasonal variability (represented by standard deviation) of NDVI (from MOD13Q1 C6) for the target pixel (where the flux tower is located, 250 m by 250 m) need to be similar (within 20 % difference or 0.05 NDVI) to the neighboring (5 km by 5 km) area; (2) the maximum NDVI value for target pixel and neighboring area needs to be greater than 0.2 (not barren).The daily GPP estimates, estimated using nighttime method (Reichstein et al., 2005), were averaged and aggregated into 4-day values to compare with CSIF.The 4-day GPP based on more than 80 % of half-hourly valid (not gap-filled) net ecosystem exchange was retained.Only sites that have at least 92 valid observations (1 year) were used.Only 40 out of 166 sites passed these criteria and were grouped into different biome types (Table S1).In addition to CSIF all-daily , we also calculated CSIF clear-daily and CSIF site which used flux-tower-observed radiation instead of in Eq. ( 4).

NN training and validation
The NN with one layer and five neurons generally predicts the OCO-2 SIF during the training with a coefficient of determination (R 2 ) around 0.8 and an RMSE of 0.18 mW m −2 nm −1 sr −1 (Fig. 2).The model also performs well in the validation (R 2 =0.79,RMSE = 0.18) and does not show effects of overfitting.Using a variety of layer (one to three) and neuron (two to nine) combinations, we found that one layer with five neurons exhibited slightly higher model performance during the validation compared with a more complex NN (Fig. A1 in Appendix).Therefore, we chose to use the four-band reflectances to feed the one-layer-five-neuron NN to generate the contiguous SIF for 2000 to 2017 when MCD43C4 NBAR dataset is available.
We also investigated the bias of our prediction among different biome types in Fig. 3.For 9 out of 14 biome types, the differences between the CSIF clear-inst and the satellite-retrieved SIF are less than 10 %, and most of the biases were within 5 %.Wetlands and urban ecosystem show a 15 % bias compared to the satellite-retrieved SIF, which may be caused by the water or built-up contamination on the reflectance signal and the relatively small sample numbers.For savannas and grassland, the changes in fluorescence yield due to seasonal drought may be important, which cannot be considered in the NN based on reflectances only.Over croplands, CSIF exhibits a 12 % underestimation.The croplands usually have high nitrogen/chlorophyll concentration that may not be fully captured by the four broadband reflectances (Wu et al., 2008).Because we did not build biome-specific NNs for the training, we do not expect biome-specific (especially needleleaf vs. broadleaf) relationships between SIF and reflectance.Interestingly, we still reproduced SIF with very high accuracy regardless of the plant function traits (PFTs), i.e., leaf types and canopy characteristics (leaf clumping, etc.).This suggests that the escape factor and long-term changes in mean fluorescence yield might be correctly accounted for by the NN across PFTs, through the information available in the reflectances only.However, it should be noted that this does not suggest that the NN and reflectances can fully replicate the fluorescence yield variations due to short-term variations caused by stresses.
We also compared the time series of predicted CSIF and OCO-2 SIF for 12 typical biome types (Fig. 4).The predicted CSIF accurately captures the seasonal and interannual variation for most biome types, while the standard deviation for each DOY is usually smaller than OCO-2 SIF.This may suggest that the uncertainty of SIF is smaller in CSIF dataset.For some ecosystems, e.g., DBF, MF and CRO, CSIF shows slight underestimation during the peak growing season.
When comparing the daily average SIF from satellite retrievals with the predicted all-sky daily CSIF (CSIF all-daily ) dataset (Fig. 5), the predicted SIF exhibits ∼7 % underestimation, with an R 2 of 0.71 and a RMSE of 0.08 mW m −2 nm −1 sr −1 .The clear-sky daily CSIF (CSIF clear-daily ) shows ∼11 % overestimation, with a slightly higher R 2 and lower RMSE.Considering the uncertainty in SIF retrievals and the inconsistency in time of the comparison (satellite SIF was based on instantaneous PAR at the time of satellite overpass and converted to daily values assuming the atmospheric condition did not change within a day; predicted CSIF was based on 4-day average PAR), the all-sky daily CSIF performs reasonably well.

Spatial-temporal variation of the global 0.05 ∘ SIF datasets
Using the trained NN with the gap-filled reflectance datasets, we produced two global CSIF datasets at 4-day temporal and 0.05 ∘ spatial resolution.
Figure 6 shows the spatial patterns of the 90th percentile for each pixel and the annual average for both clear-sky instantaneous CSIF (CSIF clear-inst ) and the all-sky daily average CSIF (CSIF all-daily ).For the 90th percentile, CSIF clear-inst exhibits hotspots in the tropical rainforest, south Asia and the North American corn belt, consistent with regions with high peak productivity (Guanter et al., 2014); CSIF all-daily shows similar spatial patterns but with relatively lower values in the tropical forest, due to the persistent cloud coverage.For the annual average SIF, tropical forests exceed temperate cropland and show very high values for instantaneous clear-sky SIF.In all conditions, African tropical forests exhibit lower values than Amazonian and Southeast Asian tropical forests.
We further investigated the seasonal and interannual variations of the all-sky daily SIF across the latitudes.The tropical regions show continuous high SIF values across seasons, and the northern mid-to high-latitude regions also exhibit recurrent high values during the Northern Hemisphere summers (Fig. 7a).Near 40 ∘ S, a hot spot is present in austral summer, with high interannual variability.Low SIF values can be found in dry years (2006)(2007)(2009)(2010), while high values were observed in wet or normal years (2010)(2011)(2012)(2013)(2014)(2015).The global average SIF also displays a strong seasonality coinciding with the Northern Hemisphere growing season (Fig. 7b).For the annual total SIF values, a statistically significant increasing trend (Mann-Kendall test, p<0.0001) is found with around 0.39 % increase per year.The year 2015 exhibited a low anomaly after detrending, which may be caused by the El Niño events (Fig. 7c).
The spatial pattern of the trend in CSIF all-daily is displayed in Fig. 8.An increasing trend dominates Europe, southeast Asia and south Amazon.A decreasing trend is mostly found in east Brazil, east Africa and some areas of inland Eurasia.The histogram also shows a positive shift with a magnitude (0.00027 mW nm −1 sr −1 yr −1 ) similar to the average global trend in Fig. 7c.The spatial pattern of CSIF all-daily is very similar to the trend pattern of MODIS enhanced vegetation index (EVI) (C6) (Zhang et al., 2017b), but the south Brazilian Amazon forest shows a more positive trend than that of EVI.

Comparison between SIF from GOME-2 and CSIF
We then compared the CSIF datasets with the reconstructed SIF (RSIF) and SIF * based on coarser-scale and all-sky GOME-2.Although these datasets were trained based on different satellites, the relationship between CSIF and RSIF or CSIF and SIF * is consistent across most regions across the globe (Fig. 9).The R 2 values are generally high (> 0.8) for most regions except over tropical rainforests, barren regions in western US, northwestern China, and northern Canada and Russia.The low R 2 values are mostly due to the relatively low variability in the temporal domain in the tropics but are also indicative of regions strongly polluted by cloud cover in which CSIF might have a competitive advantage, as the training OCO-2 data better observe the surface due to smaller footprint and with higher signal-to-noise ratio.The regression slopes are higher for regions with persistent cloud cover (e.g., tropical forest).In the time series comparison (Fig. 9e-p), all three SIF datasets show similar seasonal patterns, while GOME-2-based RSIF and SIF * generally show higher values than CSIF.In addition, RSIF exhibits larger fluctuation during the non-growing season for some sites, which may be caused by snow contamination.
We further compared the CSIF all-daily with GOME-2 daily average SIF (Fig. 10).
In general, the correlation is much lower as compared with RSIF for most regions.For regions with high variability in temporal domain, the CSIF all-daily still shows high R 2 values with respect to GOME-2 SIF.The regression slopes exhibit smaller variation except for the Amazonian tropical rainforests, southeast Asia and barren regions in the Sahara, western US, northwestern China, central Australia and the Andes mountains in South America.In general, considering the various uncertainties and different satellite overpass times, sensors used and retrieval algorithms, CSIF all-daily well captured the GOME-2 SIF variations both in space and time.In addition, since GOME-2 SIF in most Argentina is affected by the South Atlantic Anomaly (SAA), the coefficient of determination values are also lower as compared with Fig. 9.

Using CSIF for drought monitoring
Since the CSIF dataset only uses broadband reflectances, it should not contain the SIF yield information.Compared to the SIF retrieved from OCO-2, the difference can be mostly attributed to the SIF yield .Therefore, the difference or ratio between SIF OCO-2 and CSIF can reflect the environmental stress on SIF yield . Figure 11 shows the difference between instantaneous clear-day OCO-2 SIF and CSIF clear-inst .Except for Fig. 11c, the difference mostly captures the physiological limitation of drought on energy partitioning after being absorbed by chlorophyll.The spatial extent of drought is also wellcaptured by the difference, where the most severe drought-impacted places also exhibited the largest decline (e.g., Namibia, Botswana, Zimbabwe in Fig. 11a, northeast Amazon in Fig. 11b and southern Spain, southernmost France, central Italy, Croatia and Bosnia and Herzegovina).The drought impact on California is less pronounced, possibly because of the irrigation systems and sparse sampling points.
We further focused on the 2015 European drought to compare the drought response of CSIF and two vegetation indices (NDVI and EVI).Because the OCO-2 samples were not collected at the same swath for each DOY, a large fluctuation can be found in OCO-2 SIF and on the CSIF (which are using the same pixels for a fair comparison) (Fig. 12a-d).However, when calculating the ratio between CSIF and OCO-2 SIF, its variation can be mostly attributed to the variation in SIF yield , which can quantify the drought stress on plant physiology.In all three regions, the ratio between OCO-2 SIF and CSIF experienced a decrease during the drought period, but the signal is only obvious after applying a smoothing filter.The two vegetation indices, NDVI and EVI, on the other hand, show a reduced response in Spain and Italy, perhaps due to the plants' adaption or very short drought duration.

GPP-CSIF relationship across biome types
With this contiguous SIF all-daily dataset, we finally evaluated the GPP-CSIF relationship using GPP estimates from 40 flux tower sites from FLUXNET tier 1 dataset.The regression slope between GPP and CSIF (a GPP∕CSIF ) spreads across sites with a regression slope ranging from 11.91 to 68.59 (g C m −2 day −1 /mW m −2 nm −1 sr −1 ) for CSIF all-daily , 11.61 to 72.10 (g C m −2 day −1 /mW m −2 nm −1 sr −1 ) for CSIF site and 11.37 to 62.75 (g C m −2 day −1 /mW m −2 nm −1 sr −1 ) for CSIF clear-daily .The R 2 value for each individual site ranges from 0.01 to 0.93 with a median value of 0.64, 0.62 and 0.69 for all-daily, site and clear-daily CSIF, respectively.The RMSE is 1.67 g C m −2 day −1 on average.
Although the CSIF-GPP relationship varies across 40 sites, when lumping all observations within each biome type, the variation is smaller (c.v.= 0.16, rhombus in Fig. 13c, f, i).Specifically, ENF exhibited a significant larger a GPP∕CSIF (two-tailed Student's t test, p=0.036), which is caused by a stronger canopy reabsorption/scattering of SIF.OSH only has one site and also showed very high value.If both biomes are eliminated, the a GPP∕CSIF for the other biomes exhibited smaller variation (c.v.= 0.08).
The CSIF-GPP relationship not only varies across biomes but also varies within each biome type, especially for evergreen needleleaf forest (ENF, nine sites), grassland (GRA, eight sites) and wetland (WET, two sites) (Fig. 13c, f).
For CSIF all-daily , the average within-biome variation of a GPP∕CSIF (c.v.= 0.26±0.08) is comparable to cross-site variations (c.v.= 0.34) but larger than the crossbiome variations (c.v.= 0.16, using the biome-specific CSIF-GPP factor).A similar pattern can be found using CSIF site or CSIF clear-daily .

Information in contiguous SIF produced by machine learning
Vegetation photosynthetic activity has variations in several respects controlled by vegetation type, phenology, coverage and interactions with the environment.These variations can be expressed in the spatial, seasonal, diurnal and/or interannual domains (Zhang et al., 2018a).Machine learning algorithms try to minimize the differences between the predicted SIF and the satellite-observed SIF.For OCO-2 SIF and the MODIS reflectance used for NN training, the variance in the spatial and seasonal domains is largest.Therefore, the NN generally predicts SIF well in these two domains.The interannual variations (i.e., the variations caused by year-to-year anomalies, e.g., due to drought) typically have much smaller variance and are more difficult to capture.This is why some machine learning products fail to reproduce interannual variability accurately (Jung et al., 2011).Using additional variables that are sensitive to this interannual anomaly in the model training can improve the model performance (Alemohammad et al., 2017;Gentine and Alemohammad, 2018b;Tramontana et al., 2016).
In this study, since the variations in SIF yield are relatively small (Lee et al., 2015) and cannot be detected by broadband surface reflectances, the SIF yield information may not be reproduced by our CSIF data.Because the environmental limitation on SIF yield may be complicated (may not be a linear combination of temperature, vapor pressure deficit (VPD) or surface reflectance in the shortwave infrared) and biome specific (van der Tol et al., 2014), inclusion of other environmental variables and reflectances in shortwave bands during NN training did not greatly increase the SIF prediction accuracy.It should also be noted that SIF yield is relatively stable when no strong environmental limitation is present (Zhang et al., 2018c).Therefore, the CSIF product should be considered as a good proxy of OCO-2 SIF.
The satellite-retrieved SIF has a relatively large uncertainty for each individual sounding, typically ranging between 0.3 and 0.5 mW m −2 nm −1 sr −1 (Frankenberg et al., 2014).Previous site-level studies usually use SIF averaged over a large buffered area (Li et al., 2018a;Verma et al., 2017) to reduce the uncertainty.Assuming the uncertainty is unbiased and has a Gaussian distribution, machine learning algorithms are designed to reproduce SIF with lower uncertainty.Compared with previous studies that use light-use efficiency models to downscale SIF to higher resolution (Duveiller and Cescatti, 2016), this study does not rely on multiple modeled input (evapotranspiration, for example) that may introduce additional uncertainties.
Considering there is no significant trend (−0.02 % yr −1 , p>0.1) in BESS PAR (Ryu et al., 2018), this increase is likely caused by the greening of the Earth (Zhang et al., 2017b;Zhu et al., 2016) as captured in the MODIS reflectance data.This increasing trend is also within the range of most Earth system models' predictions (Anav et al., 2015).We also observed a more pronounced increasing trend in the southern Amazon than when using MODIS EVI (Zhang et al., 2017b).This may suggest that CSIF is less likely to suffer from high biomass saturation than optical vegetation indices and can more effectively detect changes in tropical rainforests or over high leaf area regions such as croplands.

The use of satellite SIF for drought monitoring
Drought can be categorized into different stages.At an early stage, when plants sense water deficit in the soil and higher vapor pressure deficit in the atmosphere, they reduce water loss through stomatal closure.This, in turn, also reduces the CO 2 exchange from stomatal closure and inhibits photosynthesis.The quantum yield for heat dissipation will increase accompanied by a decrease in quantum yield for photochemical quenching and fluorescence (Genty et al., 1989;Porcar-Castell et al., 2014).This should allow satellites to potentially capture this decrease in the SIF signal (especially during the mid-afternoon when stress is more pronounced) as an indicator of vegetation stress.In the second stage, with prolonged dry conditions, plants will recycle the nitrogen in the leaves as represented by a decrease of the greenness (chlorophyll content) of leaves.In the third stage, if the drought continues, leaf senescence and vegetation mortality may follow.SIF can potentially detect changes during all those drought stages, whereas broadband reflectance-based indices (NDVI, EVI) should only see the second and third stages.
Previous drought monitoring studies have mostly used vegetation indices (VIs) as a indictor of drought stress (Ji and Peters, 2003;Zhang et al., 2013).However, vegetation indices can only respond to drought changes in the plants' optical properties (mostly during the second and third stages).For most plants, there might be a tipping point where plants will not recover from drought-induced xylem cavitation (Urli et al., 2013).Since most VIs (e.g., NDVI, EVI) are most sensitive to the canopy changes, drought monitoring based on VIs may not be useful for drought mitigation and agricultural irrigation management.SIF retrievals from satellite, compared to optical reflectance signals, carry the information not only about the PAR absorption by chlorophyll but also about the drought stress on plant physiology.Although previous studies used satellite-based SIF datasets for post-drought impact assessment (Lee et al., 2013;Yoshida et al., 2015;Sun et al., 2015;Wang et al., 2016), these studies did not separate the contribution of decreased APAR chl or deceased SIF yield .A more recent study compared the SIF and VIs in India during heat stress (Song et al., 2018) and found that SIF is more sensitive to heat stress than VIs.Similarly, since NDVI and EVI cannot well capture the change in chlorophyll concentration, heat stress on APAR chl and SIF yield cannot be fully separated.This study developed a new method to compare the difference between SIF signals and the reflectances, which can be applied for early drought warning at global scale.
Although daily OCO-2 data have large gaps between swaths, combining several days of observation can provide enough spatial coverage considering the spatial extent for most drought events.The spatial coverage issues could be further improved using geostatistical-based methods (Tadić et al., 2017), but this may need further investigation.Compared to other meteorological drought indices, this drought monitoring technique uses only near real-time data and avoids the interannual anomalies caused by other factors (land cover change, crop rotation, etc.).The MCD43C4 dataset uses 16 days of inputs for the model inversion, and although this may lead to temporal inconsistencies for the comparison between CSIF and OCO-2 SIF, it may have limited effect due to the higher data quality during drought because of the reduced cloud coverage.

Cross-biome and within-biome GPP-CSIF relationship
In contrast to Sun et al. (2017), we found a large variation of GPP-CSIF relationship across sites.Compared to previous studies, our study gave higher a GPP∕CSIF estimates, probably due to a much higher a GPP∕CSIF value for evergreen needleleaf forest (10 out of 40 sites are ENF) (Table S1 in the Supplement) and slight underestimation of CSIF all-daily dataset.This higher a GPP∕CSIF value for ENF was also suggested by the comparison between OCO-2 SIF and FLUXCOM GPP dataset (Sun et al., 2018) and other comparisons using GOSAT SIF (Guanter et al., 2012).Consistent with Li et al. (2018b), we also found small cross-biome variation of the GPP-SIF relationship.However, a large within-biome variation of a GPP∕CSIF is also found, which contributes to a large proportion of the observed cross-site variations rather than the cross-biome variation.Compared to studies that use OCO-SIF within a large buffering area (e.g., 40 km diameter circle in Verma et al., 2017), we made the comparison over a much smaller area and much higher temporal frequency.
There are several explanations for the observed site-specific GPP-SIF relationship.
(1) Leaf morphology may directly affect the reabsorption and scattering of SIF that leaves the foliage (Atherton et al., 2017); however, this factor is not considered in current SIF modeling (van der Tol et al., 2009;Verrelst et al., 2015) and will directly affect the model simulation of the GPP-SIF relationship at the ecosystem scale (Verrelst et al., 2016;Zhang et al., 2016c).( 2) Vegetation canopy characteristics also affect the reabsorption and scattering of SIF before leaving the canopy (Romero et al., 2018;Yang and van der Tol, 2018).
(3) Atmospheric condition may attenuate and bias satellite SIF retrievals to some extent, but this effect is assumed to be small unless thick clouds are present (Frankenberg and Berry, 2017).( 4) SIF and GPP likely have different sensitivities to environmental stresses (Flexas et al., 2002); therefore, ecosystems with frequent environmental stresses (e.g., drought) during the growing season tend to have relatively lower GPP-to-SIF ratio.( 5) Since light saturations have less effect on SIF than GPP (Damm et al., 2015;Zhang et al., 2016c), the growing-season averaged light intensity (affected by latitude, average cloud coverage), vegetation canopy structure and leaf characteristics that relate to the light saturation will also affect GPP-SIF relationship.For example, the evergreen needleleaf forests have much higher specific leaf area and usually lower Sun zenith angle, making them less prone to light saturation.These factors may vary not only across biomes but also across sites.Therefore, within one biome type, the GPP-SIF relationship can also be different.
It is also noteworthy that clear-sky daily SIF exhibited stronger correlation with GPP (Fig. 13); a possible explanation would be that the light-use efficiency increases with diffused radiation, which partly compensates for the decrease in incoming PAR when clouds are present (Gu et al., 2002;Turner et al., 2006).Because the satellite SIF retrieval algorithm discarded observations that were affected by thick clouds (Sun et al., 2018), the SIF retrievals from OCO-2 are more positively biased than the actual SIF emission of the plants.However, during periods when thick clouds are present, the LUE also increases and so does the GPP ∕ SIF ratio.The positive SIF retrieval biases compensated the increase in the GPP ∕ SIF ratio and therefore contributed to a stronger correlation between satellite-retrieved SIF (rather than the actual SIF emission) and GPP.

Uncertainties and caveats
Although our CSIF clear-inst showed good performance as supported by the comparison with the clear-sky instantaneous SIF retrievals from OCO-2, the CSIF all-daily exhibits a slight underestimation.A possible explanation is that most SIF retrievals during overcast conditions did not pass the quality checks, such that OCO-2 SIF are more likely obtained during clear-sky conditions.This is supported by the fact that if we compare OCO-2 SIF with clear-daily SIF, the R 2 is even higher (Fig. 6).
The canopy structure and sun-sensor geometry were not explicitly considered in our modeling and only implicitly embedded in the machine learning retrieval.Several recent studies suggest that canopy structure will affect the PAR absorption and re-absorption of SIF before leaving the canopy (f esc in Eq. 3) (Knyazikhin et al., 2013;Liu et al., 2016;Yang and van der Tol, 2018) and further affect the GPP-SIF relationship (He et al., 2017;Migliavacca et al., 2017;Zhang et al., 2016c).However, most of these studies made assumptions requiring either a dense canopy or non-reflecting soil and thus cannot be easily applied at the global scale.In addition, OCO-2 SIF data used in this study are from nadir observations, while both the MODIS and GOME-2 sensors acquire images both at nadir and near nadir.Such discrepancy in observation angles may induce bidirectional effects.
Since CSIF is trained based on the satellite-observed SIF instead of the canopy SIF emission, and as previously discussed, it did not consider the atmospheric attenuation of SIF signal in the presence of clouds.The CSIF values are expected to be closer to the canopy SIF emission than the satellite-observed SIF at the top of atmosphere.
The BESS PAR 4-day dataset has high overall accuracy (RRMSE of 15.2 %) and very little bias (1.4 %).For different climate zones, the uncertainties are typically under 20 %.These uncertainties do not affect the CSIF clear-inst data but will propagate to CSIF all-daily .

Conclusions
In this study, using the surface reflectance from the MODIS instrument and a NN algorithm, we developed two spatially contiguous and high-temporalresolution SIF datasets (CSIF).These two SIF products not only show high accuracy when validated against the satellite-retrieved OCO-2 SIF but also exhibit reasonably high consistency with both reconstructed and satelliteretrieved GOME-2 SIF.CSIF all-daily exhibits an increasing trend globally during 2001-2016, which is attributed to the Earth greening and not to changes in PAR.Since the CSIF dataset includes most information on PAR absorption of chlorophyll, the difference between OCO-2 SIF and CSIF mostly contains the information on physiological stress on fluorescence yield.This indicator is found to be more effective for early drought warning than vegetation indices.By comparing CSIF all-daily with GPP estimates across 40 EC flux tower sites, the GPP-SIF relationship is found to vary across sites, and a large proportion of this comes from within-biome variation.However, this finding still requires further examination using SIF from both new satellites instruments (e.g., TROPOMI) and ground-based measurements.The high-resolution CSIF dataset can be further used for regional to global carbon and water flux analysis.