Spaceborne Sun-Induced Vegetation Fluorescence Time Series from 2007 to 2015 Evaluated with Australian Flux Tower Measurements

Abstract

A global, monthly averaged time series of Sun-induced Fluorescence (SiF), spanning January 2007 to June 2015, was derived from Metop-A Global Ozone Monitoring Experiment 2 (GOME-2) spectral measurements. Far-red SiF was retrieved using the filling-in of deep solar Fraunhofer lines and atmospheric absorption bands based on the general methodology described by Joiner et al, AMT, 2013. A Principal Component (PC) analysis of spectra over non-vegetated areas was performed to describe the effects of atmospheric absorption. Our implementation (SiF KNMI) is an independent algorithm and differs from the latest implementation of Joiner et al, AMT, 2013 (SiF NASA, v26), because we used desert reference areas for determining PCs (as opposed to cloudy ocean and some desert) and a wider fit window that covers water vapour and oxygen absorption bands (as opposed to only Fraunhofer lines). As a consequence, more PCs were needed (35 as opposed to 12). The two time series (SiF KNMI and SiF NASA, v26) correlate well (overall R of 0.78) except for tropical rain forests. Sensitivity experiments suggest the strong impact of the water vapour absorption band on retrieved SiF values. Furthermore, we evaluated the SiF time series with Gross Primary Productivity (GPP) derived from twelve flux towers in Australia. Correlations for individual towers range from 0.37 to 0.84. They are particularly high for managed biome types. In the de-seasonalized Australian SiF time series, the break of the Millennium Drought during local summer of 2010/2011 is clearly observed.

1. Introduction

The Gross Primary Production (GPP) of the terrestrial biosphere is a key quantity in the understanding of the global carbon cycle. GPP is the amount of atmospheric carbon fixed through the process of photosynthesis by living biomass and represents the largest gross flux of CO₂ between the atmosphere and the Earth surface [1]. The interactions between climate, GPP and atmospheric CO₂ result in a strong “carbon cycle—climate coupling” [2]. The magnitude of the coupling is a major source of uncertainty in the understanding of future carbon assimilation by crops (impacting the food production), carbon storage in forests (representing a large sink for anthropogenic CO₂), and climate forcing by CO₂. Estimation of GPP over time at the global scale is therefore of great societal and scientific interest.

To date, monitoring of GPP has not been possible at scales beyond that of a single agricultural field or natural ecosystem [3]. At those scales, networks of eddy-covariance towers [4] provide a platform to measure Net Ecosystem Exchange (NEE) of carbon at high temporal resolution. GPP can be estimated from these measurements using various statistical techniques, with estimates differing by up to 10% [5]. Upscaling of eddy-covariance measurements using observation-derived carbon-climate relations and satellite data has been shown to be one of the best methods to estimate GPP at larger scales [6].

Numerical modelling to estimate GPP may employ either a land surface model [7] or approaches that rely on remotely-sensed vegetation indices (VI) combined with mechanistic descriptions of plant photosynthesis [8,9,10,11]. A typical VI is the Normalized Difference Vegetation index (NDVI) [12,13,14], which is related to the fraction of Absorbed Photosynthetically Active Radiation (fAPAR) and the Leaf Area Index (LAI). These products are inputs in light-use-efficiency approaches [15] or even in a full Farquhar photosynthesis method [16] that drive several carbon cycle models (e.g., C-Fix [17], and SibCASA [18]). VIs are linked to chlorophyll content and the temporal resolution of spaceborne composite VI products is typically eight days or more. Consequently, short term environmental changes affecting plant productivity (e.g., water availability, temperature, nutrient deficiency, diseases) may not be well accounted for; this can lead to biases in the carbon uptake capacity of vegetation [19,20]. Hence the use of “greenness indicators” tends to return the potential carbon uptake by plants rather than the actual uptake [21,22]. It is therefore not surprising that the available methods’ global GPP estimates have a large range (from 100 to 180 PgC/year [6,23]). At smaller scales, such as individual ecosystems or plant species, the uncertainty in the level of photosynthesis is likely to be larger, because of differing drought responses, management practices, soil conditions, and nitrogen availability [24].

Photosynthetic activity is tightly related to plant fluorescence by plants in the red and near-infrared wavelength range (see below). Sun-induced fluorescence (SiF) is the re-emission of solar radiation absorbed by leaf chlorophyll at longer wavelengths. It has a distinct spectral shape with one peak at 730–740 nm (far-red fluorescence) and another at 685–690 nm (red fluorescence). This emission varies smoothly with wavelength [25]. SiF can be measured in the field or by airborne remote sensing instruments (e.g., [26,27]; see also [28] and references therein). Since approximately 1%–2% of the absorbed solar energy by vegetation is released as fluorescence (against typically about 19% released as heat and about 80% used for photosynthesis [29]), capturing the SiF requires sensors with good signal-to-noise ratios (and fine spectral resolutions, see below). At the same time, SiF may potentially interfere with retrieval of other atmospheric properties [30,31].

It has recently been shown that a number of satellite spectrometers can be used to retrieve SiF (particularly the far-red peak). These include the Greenhouse gases Observing SATellite (GOSAT) [32,33,34,35], the Global Ozone Monitoring Experiment 2 (GOME-2) [21,36], the SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY (SCIAMACHY) [33,37], and the Orbiting Carbon Observatory 2 (OCO-2) [38]. Future instruments capable of measuring SiF are the TROPOspheric Monitoring Instrument (TROPOMI) [39] and the dedicated FLuorescence EXplorer mission (FLEX) [40]. From the current instruments, particularly the GOME-2 instruments profit from a good spatial coverage and the availability of long historical records: GOME-2 on Metop-A has been fully operational since 2007 and GOME-2 on Metop-B since 2013. Retrieving reliable SiF signals from spaceborne remote sensing instruments has been an active field of research in the past decennium.

Early work on satellite remote sensing of vegetation fluorescence investigated the possibility of using the filling-in of strong Fraunhofer lines and specifically the filling-in of the dark and spectrally wide O₂ A (~760 nm) and O₂ B (~687 nm) line absorption bands (e.g., [28,41,42]). The proposed methods typically require accurate modelling of atmospheric radiative transfer in the absorption bands. A subsequent series of studies showed that filling-in of deep solar Fraunhofer lines can also be used to detect vegetation fluorescence from space [32,33,35,37,43], which in some respects simplifies the retrieval problem. The first maps of global terrestrial fluorescence using filling-in of Fraunhofer lines were achieved with the high-spectral-resolution interferometer aboard GOSAT [32,34,35]. GOSAT’s coarse temporal and spatial resolution, however, requires considerable averaging of the individual fluorescence retrievals into spatiotemporal composites to produce global maps, which may introduce sampling or representativeness errors.

It was shown by Joiner et al. [21,33] that the retrieval approaches developed for GOSAT could also be applied to the imaging spectrometers SCIAMACHY and GOME-2. This improves the spatial representativeness because these instruments achieve global coverage within a few days. In addition, the statistical or data-driven retrieval approach developed by Guanter et al. [35] for micro-windows centred at isolated Fraunhofer lines has been extended by Joiner et al. [21] to much broader spectral windows covering also the O₂ A and water vapour absorption bands. In this method, the modelling of atmospheric transmittance is approximated using a Principal Component (PC) analysis of a set of reference spectra. A large sample of observations over fluorescence-free scenes is needed for generating a representative set of PCs. The selection of these fluorescence-free scenes, the statistical analyses of atmospheric transmittances and the choice of the fitting window can however strongly impact the retrieved fluorescence signal (e.g., [21,36]). For example, optimization of PC selection for the fluorescence retrieval is addressed by Köhler et al. [36].

The physiological relationships between SiF and GPP have been described in the literature (e.g., [44,45,46,47,48,49,50,51,52]). SiF is one of many quenching mechanisms that plants use to release excess energy. The amount of excess energy depends on the efficiency with which absorbed light is used for photosynthesis, which in turn depends on rate-limiting factors such as the atmospheric CO₂ concentration, canopy temperature, and the availability of water and nutrients. In high light conditions, studies have shown that SiF positively correlates with photosynthesis and relates to both absorbed PAR and photosynthetic efficiency [52]. In addition, we mention that each of the two SiF peaks is an expression of two related but individual photosystems.

Despite the complexity of the underlying mechanisms, satellite-retrieved far-red SiF has been directly compared with GPP from data-driven and global dynamic vegetation models as well as tower-based flux estimates [34,53,54]. Significantly positive, biome-specific correlations were found between absolute SIF and absolute GPP values with the highest correlations reported for managed vegetation systems [53,54]. Note that the satellites are in sun-synchronous orbits, which means that the local (solar) overpass time is nearly constant and lighting conditions are similar across the globe.

The goal of this paper is twofold. Firstly, we further explore the applicability of the data-driven retrieval approaches developed for broad spectral windows by Joiner et al. [21] (see also [36]). We set up our own algorithm for retrieval of far-red SiF from GOME-2A and follow the general method of Joiner et al. [21], while the details of the implementation differ in a number of aspects. Rather than performing detailed retrieval simulations, we can now make a comparison between the SiF time series produced by the two algorithm implementations and provide scientific verification of the data-driven method. Secondly, we further investigate the relationship between satellite-retrieved far-red SiF and flux tower GPP. To this end, we produce a SiF time series retrieved from GOME-2A spanning more than eight years. This time series is then compared with that of GPP derived from flux tower measurements. Here, we focus on flux tower measurements from Australia (OzFlux). These sites have the advantage of being representative of relatively large areas due to the relative homogeneity of ecosystems in Australia.

The paper is structured as follows. Section 2 describes the retrieval setup that we used for producing the SiF time series from GOME-2A measurements. Section 3 provides a brief description of the resulting time series. Section 4 then compares retrieved SiF values with GPP derived from Australian flux tower measurements. An anomaly in the SiF series for Australia in early 2011 is discussed. Section 5 compares the SiF time series from our implementation (“SiF KNMI”) with the latest implementation of Joiner et al. [21] (“SiF NASA (v26)”; http://avdc.gsfc.nasa.gov/)). Section 6 presents sensitivity experiments with real GOME-2A spectra. A final conclusion is provided in Section 7.

2. SiF Retrieval Method and Setup

As stated, we follow the data-driven retrieval approach developed by Joiner et al. [21] for moderate spectral resolution instruments such as GOME-2. We performed far-red SiF retrievals for GOME-2 on Metop-A and produced a time series covering the time frame from the start of the mission in January 2007 to June 2015. We compare our GOME-2A SiF time series with that from version 26 produced by Joiner et al. [21]. A number of differences exist between the two retrieval implementations. We briefly summarize the methodology with a focus on implementation differences. For the three main differences, sensitivity studies are performed (see Section 6).

2.1. The GOME-2 Instrument

The GOME-2 spectrometer aboard the Metop-A and Metop-B satellites [55] has a spectral resolution in the near-infrared channels of ~0.5 nm and a sampling interval of ~0.2 nm. The satellites are in polar sun-synchronous orbits with local equator crossing times of the ascending node around 9:30 a.m. The instrument is scanning in nadir and its scanning mirror has a variable angular speed such that footprints have about equal size. GOME-2 aboard Metop-A had a nadir footprint of 80 km × 40 km until July 2013 at which point the swath width was reduced and footprint sizes became 40 km × 40 km. Owing to the large footprints. A single pixel often represents a mixture of different surface types. Global coverage for the wider swath is achieved in almost a day.

2.2. Forward Model

The expression for the top-of-atmosphere reflectance of an atmosphere bounded by an isotropically reflecting ground surface (e.g., [56]) is extended to include a fluorescence term. Here, we assume that the upwelling fluorescence radiance field at the surface (I_F) is isotropic as well. We further simplify this expression by assuming that atmospheric scattering in the near-infrared wavelength region is small (see also Frankenberg et al. [31]). The monochromatic reflectance (R) can then be written as

$$R(\lambda ;\mu ,{\mu }_0)\approx A_s(\lambda )t(\lambda ;\mu )t(\lambda ;{\mu }_0)+\frac{\pi I_F(\lambda )t(\lambda ;\mu )}{{\mu }_0{E}_0(\lambda )}$$

(1)

where μ and μ₀ are the cosines of the viewing zenith angle and solar zenith angle, respectively; A_s is the surface reflectance; and E₀ is the solar irradiance at the top of the atmosphere. The transmission factors t(μ) and t(μ₀) here include only atmospheric extinction (no diffuse transmission). Any strong spectral variation of the top-of-atmosphere reflectance is due to these quantities as well as to (the in-filling of) the solar irradiance. The surface albedo A_s and fluorescence radiance I_F are spectrally smooth. The transmission factors can be combined and rewritten in terms of the two-way atmospheric slant absorption optical thickness ${\tau }_s^{\downarrow \uparrow }$, so that Equation (1) becomes

$$R(\lambda ;\mu ,{\mu }_0)\approx A_s(\lambda )e^{-{\tau }_s^{\downarrow \uparrow }(\lambda )}+\frac{\pi I_F(\lambda )}{{\mu }_0{E}_0(\lambda )}{e}^{-\frac{{\mu }^{-1}}{{\mu }^{-1}+{\mu }_0^{-1}}{\tau }_s^{\downarrow \uparrow }(\lambda )}$$

(2)

In our retrieval window, which runs from 712 nm to 783 nm, the absorption optical thickness comprises absorption by water vapour up to about 740 nm and by oxygen between about 760 and 770 nm. The retrieval window used by Joiner et al. [21] for data set version 26 runs from 734 nm to 758 nm, which excludes the oxygen A band but still contains some water vapour absorption lines. This is the first major difference between the two retrieval implementations.

Next, the monochromatic reflectance R(λ) in Equation (2) is replaced by the measured reflectance R(λ_i) from GOME-2. The spectral structure of the two-way slant absorption optical thickness is modelled using PCs from an analysis of measured reflectance spectra over fluorescence-free regions (see below). Convolution with the instrument’s spectral response function and any other instrument effects causing spectral structures are implicitly described by the PCs. In the actual fluorescence retrieval, the surface albedo is described as a fourth-order polynomial in wavelength (n = 4). The spectral shape of the fluorescence radiance across the retrieval window, which roughly corresponds to the second peak of the SiF spectrum, is parameterised with a Gaussian function located (μ_F) at 737 nm and having a sigma-width (σ_F) of 34 nm following Zarco-Tejada et al. [57].

Finally, we replace the monochromatic solar irradiance E₀(λ) by the high-resolution solar irradiance reference spectrum from Chance and Kurucz [58] convolved with the GOME-2 instrument response function and scaled to the actual Earth-Sun distance (denoted as ${\overline{E}}_0({\lambda }_i)$). We use this spectrum instead of the actually measured irradiance spectrum E₀(λ_i) because of radiometric degradation of the GOME-2 instrument, which largely cancels out in the calculation of the reflectance but would show up as a confounding trend in the SiF time series (for example, the SiF in 2013 would have been underestimated by more than 10% compared to 2007 if the Level-1b irradiance had been used). The use of ${\overline{E}}_0({\lambda }_i)$ instead of E₀(λ_i) is another difference compared with Joiner et al. [21]; in SiF NASA (v26) a (different) correction for radiometric degradation applied to the measured irradiance is now implemented [59].

The forward model for the fluorescence retrieval is then completed and written as

$$R({\lambda }_i;\mu ,{\mu }_0)\approx \left({\displaystyle \sum _{j=0}^n}{a}_j{\lambda }_i^j\right)e^{-{\sum }_{k=1}^m{b}_k{f}_k({\lambda }_i)}+\frac{\pi I_{F,0}{e}^{-\frac{1}{2}{(\frac{{\lambda }_i-{\mu }_F}{{\sigma }_F})}^2}}{{\mu }_0{\overline{E}}_0({\lambda }_i)}{e}^{-\frac{{\mu }^{-1}}{{\mu }^{-1}+{\mu }_0^{-1}}{\sum }_{k=1}^m{b}_k{f}_k({\lambda }_i)}$$

(3)

In our baseline retrieval setup, the number (m) of PCs (f_k) is 35. Thus, the total number of fit parameters is 5 + 35 + 1 or 41 (a_j + b_k + I_F,0). Extensive retrieval simulations were performed by Joiner et al. [21], which indicated that the accuracy and precision of the retrieved fluorescence peak value for a fit window running from 747 nm to 780 nm hardly changes when increasing the amount of PC from 25 to 35. In this paper, we will show results from real data experiments further investigating the effect of the number of PCs, which confirm the conclusion that at least 25 components are needed. In our real data experiments, we observed, however, a small improvement when using 35 compared to 25 components. Therefore this number has been chosen for the baseline retrieval setup. In SiF NASA (v26), m was reduced to twelve, because the fit window did not include the oxygen A band. The number of PCs is the second major difference between the two retrieval implementations. We remark that selection of the PCs can be further optimized by, for example, setting a lower threshold on the PC’s explained variance [39] or using more sophisticated statistical criteria [36]. Finally, we note that by making additional simplifications, the forward model can be linearized to speed up the calculations [36,39].

2.3. PC Analysis of Atmospheric Absorption Optical Thickness

The assumption is made that the same set of empirical basis functions can be used to describe the atmospheric absorption optical thickness both over fluorescent target scenes and over non-fluorescent reference scenes. In the baseline retrieval setup, we construct a set of non-fluorescent reference spectra from a year of observations over the Saharan desert. The PCs for the two-way slant absorption optical thickness are determined from this set of reference spectra in the following way.

When there is no fluorescence, the second term in Equation (3) cancels. A number of sub-windows are selected in which the absorption optical thickness becomes very small (712–713 nm, 748–757 nm and 775–783 nm) so that the top-of-atmosphere reflectance approaches the surface reflectance. We fit a second-order polynomial to the measured reflectance in these sub-windows to determine the surface albedo across the entire retrieval window (the surface reflectivity for non-vegetated surfaces shows much less spectral variation across the retrieval window, so that a polynomial of second-order can be used here). The two-way slant absorption optical thickness for a particular ground pixel is then given by

$${\tau }_s^{\downarrow \uparrow }(\lambda )=-\text{ln}\left(\frac{R^{ref}({\lambda }_i;\mu ,{\mu }_0)}{{\sum }_{j=0}^2{c}_j{\lambda }_i^j}\right)$$

(4)

Finally, PCs are determined from an analysis of slant absorption optical thicknesses for this set of non-fluorescent reference scenes. The set of reference spectra is constructed from observations over the Saharan desert (here: the latitude–longitude box defined by corners 16N, 8W and 30N, 29E). GOME-2A spectra are included only if they are completely free of vegetation according to the 1 km resolution USGS Global Land Cover Characterization database (version 2; https://lta.cr.usgs.gov/GLCC) and if the pixel’s effective cloud fraction in the near-infrared is below 0.4. The effective cloud fraction is taken from the FRESCO product [60,61] as included in the GOME-2 radiance product. We may make different choices for the region of vegetation-free reference areas. Initial experiments indicated that the Saharan desert works well in our retrieval. In Section 6.3 we return to the question of PC reference area.

For SiF retrievals in a particular month, the set of reference spectra is collected from the twelve months preceding and including this month. For example, SiF retrievals for the month of June 2015 use PCs determined from reference spectra in the time frame of July 2014 to June 2015. For SiF retrievals in 2007, which is the first year of the GOME-2A time series, reference spectra are taken from that year.

The number of reference spectra is very large (in the order of 30,000 spectra) and calculating the PCs using the covariance method is problematic. Instead, we calculate PCs sequentially in order of decreasing explained variance using the Non-linear Iterative Partial Least-Squares algorithm [62,63].

The selection of reference spectra is a third difference with SiF NASA (v26), for which observations mainly over cloudy ocean, and to some extent also desert and snow covered land are used, but for a single day only.

2.4. Monthly SiF Maps

SiF is retrieved for GOME-2A pixels that have solar zenith angles below 70° and for which the FRESCO effective cloud fraction is below 0.4. This means that if FRESCO is in snow/ice mode, pixels are excluded from the retrieval. Fluorescence retrievals are rejected whenever their root-mean-square residual is larger than 1%. Finally, retrieved SiF values are averaged into 0.5° by 0.5° latitude–longitude grid boxes for every month. An overview of the algorithm flow is given in Figure 1. In SiF NASA (v26), a first bias correction to SiF maps is made based on retrieved SiF values where zero SiF is expected ([59]; see also [37]).

3. Seasonal Global Maps of Far-Red SiF Derived from GOME-2A

Figure 2 shows several monthly mean maps that reveal the expected spatial patterns of vegetation activity (four months in 2013). For example, in the temperate and subtropical regions of the Northern Hemisphere, such as the eastern US, western Europe or Southeast China, SiF values start increasing from March onwards, peak in June and decrease again to their minimum towards November. In the more continental climate zones of eastern Europe, Russia and northeastern US, the SiF seasonal cycle is somewhat shorter: SiF values increase from May onward, have a clear peak in June and reach background levels again in October. In the subarctic climate zones of the northern parts of Eurasia and North America, the SiF cycle shortens even further. In the lower latitudes, significant SiF values occur in all months, but there is a latitudinal shift in peak values during the year. In the Southern Hemisphere, SiF values in southern Australia are high in September, October and November.

Apart from these large-scale biome patterns, very detailed spatially and temporally heterogeneous SiF patterns can also be observed from GOME-2A when evaluating over specific regions. Especially in arid regions with narrow strips of agricultural land such as in North Africa and the Fertile Crescent (Figure 3) the possibility of satellite monitoring of vegetation activity is demonstrated. The agricultural land use type shows substantially higher SiF values compared to its arid background with the highest values in March, April and May. From Fall on, retrieved values over the Fertile Crescent fade to background levels of the desert. Note that the surface albedo of desert areas is comparable to the surface albedo of vegetated land.

Finally, we calculated monthly mean SiF values for the PC reference area and for open ocean grid boxes. In both cases, no significant SiF is expected and these values can serve as an error estimate for the SiF time series. The overall mean SiF bias and its standard deviation for the PC reference area are 0.03 and 0.06 mW/sr/m²/nm, respectively. The overall mean SiF bias and its standard deviation for open ocean are 0.17 and 0.05 mW/sr/m²/nm, respectively. Small but systematically negative values over shallow or murky waters (e.g., lakes, river deltas) are observed, which are perhaps related to the interaction of visible light with suspended matter in the water.

4. Evaluating Satellite-Retrieved Far-Red SiF with Data on Vegetation Activity

4.1. Comparison with NDVI

To date, no suitable independent measurements of plant fluorescence exist for proper validation of SiF retrieved from satellite sensors. However, local fluorescence measurement techniques are rapidly advancing and show their value in studying the ecophysiology of photosynthesis at the leaf level and in situ ([52] and references therein). During the last decade, SiF has been successfully measured from tower [26,64,65,66], aircraft [27,67,68] and satellite platforms [32,34,35,41]. For evaluation purposes, vegetation indices from satellite measurements (such as the NDVI) can be used as first proxies for SiF from spaceborne measurements due to their similar spatio-temporal patterns; they can therefore be used as a qualitative indicator of vegetation activity. Another evaluation method is to compare the spatial SiF patterns derived from different satellites based on different methodologies [21,34,53]. We note that the absolute fluorescence values will depend on the assumed spectral shape of the emissions [21,52] and on other algorithm assumptions.

Figure 4 illustrates the good qualitative agreement between the spatial patterns of NDVI and SiF over North America for July 2013 suggesting that the results of our retrieval setup indeed indicate the presence of green vegetation. Over the eastern US, the northwestern US coastal area and the southern part of Mexico high NDVI values are measured, while over the western US and northern Mexico low NDVI values are observed in line with the general patterns of SiF values.

However, there are also clear differences between the spatial patterns of NDVI and SiF, especially in areas covered with tundra, taiga and rainforest vegetation. While in September SiF is already small over the tundra/taiga, NDVI is still high. Rainforest areas show slightly lower SiF radiance values compared to other areas with similar NDVI values (not shown). SiF values of evergreen forests also tend to show some minor seasonality, which in turn is hardly observed from NDVI data. It is expected that SiF measurements provide additional information that is not covered by vegetation indices due to signal saturation, bidirectional reflectance distribution function (BRDF) effects [69], and because materials such as soil, wood, and dead biomass also absorb PAR but do not contribute to photosynthesis [70]. The correlation between NDVI and SiF at the global scale is moderately strong (R = 0.78, a value similar to correlations given by Guanter et al. [53]), but large NDVI values can be observed at low SiF levels. In the next paragraph, we therefore compare SiF with GPP values derived from flux towers.

4.2. Comparing Far-Red SiF with GPP Derived from Australian Flux Tower Observations

Here, we compare SiF retrieved from GOME-2A based on the presented methodology and parameter settings with GPP derived from the regional Australian flux tower network (OzFlux—www.ozflux.org.au) [71,72]. A full description of the network can be found in Beringer et al. [72]. We used 12 Australian flux towers (

Table 1. Overview of the Australian flux tower sites from the OzFlux network (www.ozflux.org.au) used in this study. The first seven are the northern sites, and the last five are the southern sites.

Site Name	Latitude (S)/Longitude (E)	Vegetation
Sturt Plains	17°9′2′′, 133°21′1′′	Grassland
Daly Uncleared	14°09′33′′, 131°23′17′′	Woodland savannah
Adelaide River	13°4′37′′, 131°7′4′′	Woodland savannah
Daly Pasture	14°03′48′′, 131°19′05′′	Tropical pasture
Dry River	15°15′32′′, 132°22′14′′	Woodland savannah
Fogg Dam	12°32′42′′, 131°18′26′′	Flooded wetland
Howard Springs	12°29′42′′, 131°09′00′′	Open-forest savannah
Riggs	36°39′00′′, 145°34′34′′	Dryland agriculture (pasture)
Wallaby	37°25′34′′, 145°11′14′′	Eucalyptus forest
Whroo	36°40′23′′, 145°01′34′′	Box woodland
Wombat	37°25′20′′, 144°05′40′′	Secondary forest
Yanco Jaxa	34°59′16′′, 146°17′27′′	Dryland agriculture (pasture)

) available for the period January 2007 to June 2015. The tower sites cover different biomes in northern and southern Australia, from agriculture to woodlands and from pastures to savannah [73,74]. The spatial footprint of flux tower observations is much smaller than the spatial sampling of a satellite instrument such as GOME-2, and this limits the spatial representativeness. However, by selecting flux towers in areas with homogenous land cover, the vegetation type in the flux tower’s footprint will agree better with the dominant vegetation type in the satellite pixel. Since large homogeneous biomes are more common in Australia, we focus our comparison on Australian flux towers.

At the flux tower sites, meteorological parameters, carbon, water and energy fluxes were measured and GPP was derived. For each flux tower we selected the 0.5° by 0.5° grid cell with centre coordinates closest to the tower site. We then averaged the tower-derived GPP at the overpass time of GOME-2 (time window of 2 h centred at 9:30 local time) to monthly values. In this way, we only take GPP into account for the same lighting conditions in which SiF is retrieved (e.g., instantaneous GPP). For each tower site, the derived time series of GPP and the corresponding SiF series, either from the KNMI or the NASA (v26) product, are shown in Figure 5. Subsequently, we have conducted a simple linear regression of monthly GPP data on monthly satellite-retrieved SiF values (

Table 2. Comparison of Sun-induced Fluorescence (SiF) time series from KNMI and NASA (v26) with Gross Primary Productivity (GPP) derived from Australian flux tower sites for the period January 2007 to June 2015. For each site the slope, intercept and correlation between the time series are given: the linear relation reads as GPP = intercept + slope × SiF.

	GOME-2 SIF KNMI vs. Flux Tower GPP			GOME-2 SIF NASA (v26) vs. Flux Tower GPP
	Slope	Intercept	R	Slope	Intercept	R
Sturt Plains	11.5	−3.0	0.739	13.4	−2.0	0.776
Daly Uncleared	7.9	2.5	0.699	8.6	3.4	0.637
Adelaide River	17.5	−1.7	0.835	12.6	3.0	0.691
Daly Pasture	20.0	−6.0	0.831	24.0	−4.4	0.871
Dry River	5.4	4.0	0.622	8.4	3.4	0.821
Fogg Dam	12.1	−0.2	0.763	7.0	3.6	0.554
Howard Springs	12.2	5.1	0.554	15.0	6.4	0.646
Riggs	8.2	−1.0	0.627	11.8	−1.0	0.781
Wallaby	3.7	4.0	0.265	7.6	3.0	0.366
Whroo	2.5	5.0	0.371	2.2	5.8	0.299
Wombat	4.4	8.6	0.397	6.8	8.7	0.433
Yanco Jaxa	4.4	0.0	0.570	6.5	0.3	0.800
All	6.8	2.8	0.456	9.7	2.9	0.510

). In Figure 6 we show two scatterplots for a flux site with a high and a low correlation between SiF and GPP with slopes, intercepts, and correlation values given in Table 2.

Despite large differences in spatial scales, SiF is related to the flux tower GPP data with R-values up to 0.84–0.87, although some tower sites show small correlations (Table 2). The SiF KNMI-GPP slopes and correlations are generally larger for agricultural biomes as opposed to more natural biomes such as savannah and woodland. Earlier reported observations confirm the higher correlation between SiF and GPP estimated at towers located near agricultural sites [53]. The slopes of the GPP-SiF relation for grass type biomes are very similar for European grasslands and US cropland [53] except for tropical grasslands. Interestingly, biomes with lower SiF-GPP slopes have also lower correlation and in general this happens at higher LAI values, which should be investigated further. Analysing differences in the SiF-GPP relationships for different biomes requires advanced canopy and SiF modelling using radiative transfer models such as SCOPE [44,75]. By applying models, scale-related issues can be identified and separated from biophysical factors.

Overall, both KNMI and NASA (v26) SiF products perform similarly, sometimes with either higher or lower correlations for each product. Slopes of the SiF-GPP regression are smaller for the KNMI product as compared with NASA (v26). This may be due to the fact that in general the KNMI SiF is higher as compared with the NASA (v26) product (see Figure 8). The vegetation onsets (Spring) and offsets (Fall) from the SiF time series broadly follow the GPP cycle.

4.3. GOME-2A Retrieved SiF Time Series Analysis over Australian Flux Tower Sites: A Case Study

The GOME-2A retrieved SiF (January 2007 to June 2015) over Australia covers a period with interesting climatological contrasts. The 2007–2010 time range is part of the Australian decade-long “Millennium Drought”, the most severe drought period since instrumental records began in the 1990s [76]. This drought ended in the Australian summer of 2010/2011 due to the strong 2010–2012 La Niña phenomenon [77]. La Niña is typically associated with increased rainfall in northern and eastern Australia (both years experienced rainfall well above the long-term average of 465 mm–703 mm in 2010 and 708 mm in 2011). During La Niña, winter and spring daytime temperatures are below average across southern Australia, but are above average in the northern parts, while summer daytime temperatures are below average across most of Australia [77].

The relatively large timespan and spatial coverage of the GOME-2A or any other satellite-retrieved SiF series, increases the ability to study anomalies and trends in vegetation dynamics over areas with limited surface data. In addition, GPP data derived from flux tower measurements may suffer from data gaps as not all the sites provide continuous data records (see for example Figure 5). SiF data may be used as a proxy for GPP (as shown in Section 4.1, Table 2 and Figure 5) to examine trends in the carbon cycle.

Across all the Australian flux tower sites in this study, a strong positive anomaly in the de-seasonalized SiF time series (computed by subtracting from each monthly value the eight-year average for that specific month and dividing by that value to obtain percentages) was observed in 2010/2011 (blue area in Figure 7). In contrast, strong negative anomalies were observed at the start and especially at the end of the de-seasonalized time series (red areas in Figure 7). This resulted in an increase (+2.7% to +17.2% year⁻¹; on average 10.3% year⁻¹) of de-seasonalized SiF from January 2007 to April 2011. This increase is higher over southern sites (Riggs, Wallaby, Whroo, Wombat, Yanco Jaxa: +15.7% to +17.2% year⁻¹) compared to the northern ones (Sturt Plains, Daly Uncleared, Adelaide River, Daly Pasture, Dry River, Fogg Dam, Howard Springs: +2.7% to +7.5% year⁻¹). For the period May 2011 to June 2015 a decrease is observed (+0.7% to −10.3% year⁻¹; on average −6.0% year⁻¹). This decrease is stronger over the southern sites (−5.9% to −10.3% year⁻¹) as compared with the northern sites (+0.7% to −9.5% year⁻¹). We attribute the increasing SiF values until 2010/2011 and the decreasing SiF values afterwards to the meteorological extremes occurring during that period: drought at the start of the time series indicated by the negative anomalies, the severe rainfall in 2010 and 2011 connected to La Niña indicated by the positive anomalies, and again the drier period at the end of the time series [77,78,79,80]. The small differences in responses (in magnitude, not sign) between northern and southern sites can perhaps be explained by the somewhat different impact of La Niña on precipitation and temperature in both parts of Australia [77].

This short time series analysis illustrates the strong potential of SiF data to monitor vegetation activity in relation with meteorological anomalies [81,82]. Future research should further investigate the (complex) relation between satellite-retrieved solar-induced fluorescence and derived parameters, such as the fluorescence yield, and photosynthetic activity.

5. Comparison with SiF Retrieved from GOME-2A by Joiner et al. (2013)

We now compare our global SiF time series with SiF NASA (v26). In Section 2 we discussed our retrieval setup and indicated where it differed from the implementation in SiF NASA (v26), the most current NASA version at this time. The main differences are the selection of the reference scenes for determination of the PCs (desert areas, yearly vs. mainly cloudy ocean, daily), the fit window (712–783 nm vs. 734–758 nm) and the number of PCs (35 vs. 12). Both time series are spatiotemporal composites of retrieved SiF values averaged per month and in 0.5° by 0.5° grid boxes. Figure 8 shows a scatter plot of the two SiF data sets where a single point represents a monthly grid box average from the global time series between January 2007 and June 2015.

Only grid boxes that are over land and for which the monthly average in each time series is calculated from at least three successful SiF retrievals are taken into account. The correlation between the two datasets is 0.78 and the RMS difference is 0.35 mW/sr/m²/nm. Slightly different effective cloud fractions are used to filter the retrievals that go into the global SiF maps (0.4 in our implementation vs. 0.3 in SiF NASA, v26), but this has a negligible effect on the comparison. We find that filtering with cloud fractions other than 0.4 hardly changes the correlation. This is in agreement with Köhler et al. [36], who show that filtering with higher cloud fractions decreases SiF values but leaves temporal patterns mostly unaffected.

A global map of correlations is depicted in Figure 9. The number of data points per grid box is high because of the long time series (up to 102 months per grid box). Firstly, we see that correlations of SiF retrieved over the ocean are quite small but slightly positive and do not exhibit any spatial patterns. This may be due to noise having similar effects on the retrieval solution for the two implementations (the fit windows partly overlap). The absence of spatial patterns in correlations over the oceans argues in favour of this explanation. Similarly, we see small but slightly positive correlations over land free of vegetation (e.g., the Saharan desert). Secondly, correlations between the two time series are in general substantial over vegetated areas. However, there is a clear dependency on the vegetation type. In particular, we see that correlations for tropical rainforests drop to the same background levels observed for non-vegetated areas (see, for example, the rainforests in Central and South America, West and Central Africa, and Indonesia). We have checked and excluded the possibility of a confounding effect of a reduced number of valid monthly SiF averages because of increased cloud contamination for these areas. Either the two retrieval implementations respond differently to SiF over tropical rainforests or the natural variability of SiF for these regions is too small (compared to the noise error) to reveal a correlation in the first place.

6. Sensitivity Analyses for the Retrieval of GOME-2A SiF Using PCs

In this section, we discuss a number of experiments investigating the three most important retrieval sensitivities. Rather than performing closed-loop retrieval simulations, we do retrieval experiments with real GOME-2A spectra. We investigate the effect of the number of PCs, the fit window, and the PC reference area on retrieval and use the four months of January, April, July and October 2013 as test months for reasons of computational time.

6.1. Number of PCs

In addition to the baseline retrieval, which uses 35 PCs, we performed retrievals for the four test months using 15, 25 or 45 PCs while keeping all other settings the same. The results are summarised in

Table 3. Various statistical parameters describing the difference between retrieved SiF values for the months of January, April, July and October 2013 when using 15, 25 or 45 PCs as compared to 35 PCs (baseline). Only grid boxes over land are taken into account. Parameters are the root-mean-square difference (RMS), the mean difference and the standard deviation (STD) of the difference, the correlation coefficient (R), the slope and intercept of the regression line, and finally the total number of valid monthly grid box pairs (N) from which differences are calculated. The difference is defined as SiF(test) − SiF(baseline) and the regression line as SiF(baseline) = intercept + slope × SiF(test). All correlations are highly significant. The unit of SiF values is mW/sr/m²/nm, as before.

Number of PCs	Month	RMS	Mean	STD	R	Slope	Intercept	N
15	January	0.85	−0.11	0.85	0.67	0.39	0.44	28,844
	April	0.80	0.54	0.59	0.83	0.39	0.00	37,706
	July	1.52	1.00	1.14	0.85	0.29	0.16	44,326
	October	0.74	0.32	0.67	0.79	0.40	0.13	40,786
	All	1.06	0.49	0.94	0.74	0.32	0.19	151,662
25	January	0.20	0.09	0.18	0.99	0.82	0.04	28,844
	April	0.25	0.19	0.16	0.97	0.76	−0.06	37,706
	July	0.90	0.69	0.58	0.94	0.48	0.00	44,328
	October	0.51	0.34	0.39	0.95	0.58	−0.01	40,786
	All	0.57	0.36	0.45	0.92	0.56	0.03	151,664
45	January	0.15	−0.06	0.13	0.99	1.13	−0.01	28,844
	April	0.18	0.09	0.15	0.97	0.78	0.01	37,706
	July	0.17	0.11	0.13	0.99	0.84	0.01	44,328
	October	0.22	0.12	0.18	0.98	0.77	0.01	40,786
	All	0.18	0.07	0.17	0.97	0.86	0.01	151,664

and presented in terms of differences with respect to the baseline of 35 PCs. The comparison includes only grid boxes over land.

In a qualitative sense, we find that global maps when using only 15 PCs show patches with negative values as well as large patches with highly positive values. The patches with negative SiF values more or less disappear when using 25 PCs or more and SiF patterns correspond much better to actual vegetation patterns. SiF values tend to be lowest for the baseline retrieval with 35 PCs (also when compared to the retrieval with 45 PCs).

A sufficient number of PCs is needed to capture the spectral variability of atmospheric transmission. When this number has been reached, one would in principle not expect the retrieval outcome to change much when more PCs are added. The correlation between SiF values for 25 and 35 PCs is about as high as between SiF values for 45 and 35 PCs: this correlation is well above 0.9. However, we see that differences between SiF values for 25 and 35 PCs, in terms of RMS, absolute value of the mean difference and the standard deviation, are higher. The mean difference between SiF values for 45 and 35 PCs has reduced to 0.07 mW/sr/m²/nm and it seems justified to not further increase the number of PCs from 35 to 45 as the two conditions already resemble each other and we do not want to use more PCs than needed.

When more PCs are used there is the potential danger of overfitting or noise fitting. This would show up in the retrieval results as retrieved SiF values becoming noisier. The comparison with flux tower GPP is important here. Comparing retrieved SiF values for the months of January, April, July and October 2013 for the different PC numbers (15, 25, 35, or 45) with flux tower-derived GPP values, we find that the correlation between GPP and SiF increases with five percentage points when using 25 PCs instead of 15 PCs, with 28 percentage points when using 35 PCs instead of 15 PCs but with only 20 percentage points when using 45 PCs instead of 15 PCs. Thus, there is a slight decrease of the correlation when going from 35 to 45 PCs. These results indicate that adding more PCs increases the ability to better retrieve SiF, but that there is indeed an optimum amount at which adding more PCs will deteriorate the retrieval (see also [21]).

6.2. Spectral Window

In a second experiment, we test the effect of the spectral window. We now run the retrieval for the four months using spectral windows 734–758 nm, 712–758 nm and 734–783 nm in addition to the baseline window running from 712 to 783 nm. All other settings are the same as in the baseline retrieval. The four windows differ among each other with respect to the atmospheric absorption bands they cover. The 734–758-nm window contains mainly Fraunhofer lines, the 712–758-nm window covers a water vapour absorption band and the 734–783-nm window covers the oxygen A band.

Table 4. Similar to Table 3 but for tests with different fit windows. Spectral windows 734–758 nm, 712–758 nm or 734–783 nm are compared to 712–783 nm (baseline). The windows contain Fraunhofer lines (FH), water vapour and/or oxygen absorption, as indicated in the table. All correlations are highly significant.

Spectral Window	Month	RMS	Mean	STD	R	Slope	Intercept	N
734–758 nm (FH)	January	1.16	−0.11	1.15	0.45	0.23	0.53	28,805
	April	0.41	−0.08	0.40	0.74	0.51	0.21	37,677
	July	0.43	−0.09	0.42	0.72	0.67	0.27	44,297
	October	1.58	−0.19	1.57	0.22	0.07	0.41	40,592
	All	1.01	−0.12	1.00	0.42	0.21	0.43	151,371
712–758 nm (H2O, FH)	January	0.41	0.02	0.41	0.91	0.67	0.20	28,813
	April	0.57	0.25	0.51	0.85	0.43	0.09	37,700
	July	0.46	0.21	0.41	0.93	0.59	0.15	44,325
	October	0.33	0.04	0.33	0.95	0.62	0.13	40,761
	All	0.45	0.14	0.43	0.90	0.58	0.13	151,599
734–783 nm (FH, O2 A)	January	0.46	−0.28	0.37	0.88	1.48	0.10	28,844
	April	0.21	−0.03	0.21	0.86	1.06	0.01	37,706
	July	0.40	−0.30	0.26	0.90	1.34	0.18	44,328
	October	0.28	−0.14	0.24	0.90	1.25	0.07	40,786
	All	0.34	−0.19	0.29	0.87	1.31	0.09	151,664

summarizes the results and reports differences in retrieved SiF values with respect to the baseline fit window of 712–783 nm.

Overall, SiF maps are most noisy for the shortest fit window, which is understandable as there are simply less spectral points and because there is no filling-in of atmospheric absorption bands that can provide information. Differences typically get smaller and correlations get better when the spectral window resembles the baseline window more. Mean differences are most meaningful here: Importantly, we see that when the H₂O band is mostly excluded SiF values generally decrease and when the O₂ A band is excluded SiF values typically increase. We remark that the shorter fit window containing only Fraunhofer lines is used in SiF NASA (v26), but that the retrieval setup also differs in a number of other aspects from the experiment reported here.

6.3. PC Reference Area

Finally, we have investigated the effect of the reference area from which spectra are selected for construction of the PCs. For the four test months in 2013, we have also used cloudy ocean and areas covered with snow or ice as reference areas. For selection of cloudy ocean pixels, we follow the approach by Joiner et al. [21] and select pixels over ocean with a minimum cloud fraction of 0.9 but also with a minimum cloud pressure of 850 hPa to avoid shielding of lower parts of an atmospheric column by high clouds. Pixels are taken from two large predefined latitude–longitude boxes over the North and the South Atlantic Ocean. Pixels over snow- or ice-covered areas are selected from a predefined latitude–longitude box over Greenland and the Canadian Arctic Archipelago.

The retrieval using PCs determined from spectra over snow/ice reference areas performs worst (

Table 5. Similar to Table 3 but for tests with different reference spectra. Reference spectra for the PC analysis from cloudy ocean or snow/ice-covered areas are compared to desert areas (baseline). All correlations are highly significant.

Reference Area	Month	RMS	Mean	STD	R	Slope	Intercept	N
Cloudy ocean	January	0.33	−0.19	0.27	0.91	1.07	0.15	28,844
	April	0.37	0.24	0.28	0.95	0.61	−0.01	37,706
	July	0.58	0.41	0.41	0.93	0.58	0.03	44,328
	October	0.27	0.11	0.25	0.89	0.82	−0.01	40,786
	All	0.42	0.17	0.38	0.86	0.65	0.07	151,664
Snow/ice	January	1.06	−0.68	0.82	0.06	0.08	0.65	28,844
	April	0.79	−0.42	0.67	0.09	0.06	0.35	37,706
	July	0.84	−0.46	0.71	0.55	0.36	0.58	44,328
	October	0.55	−0.21	0.51	0.52	0.49	0.32	40,786
	All	0.81	−0.42	0.69	0.35	0.29	0.49	151,664

). Global SiF maps show very large patches of negative SiF values, particularly in a large band around the equator. SiF maps using cloudy ocean or desert as PC reference areas, however, are largely consistent with each other: the overall correlation R is as high as 0.86. As mentioned before, the baseline retrieval shows slightly positive SiF values over the oceans, which indeed reduce when cloudy ocean is used as reference area. However, these positive values do not completely disappear, i.e. a positive bias is present for the cloudy ocean reference case as well. There is a tendency for SiF values using the cloudy ocean reference to be larger than those using the desert reference.

6.4. Conclusion of the Sensitivity Experiments

Taking the results from the three experiments together, we can now better understand the comparison with our SiF values versus NASA (v26) set as discussed in the previous section. In the implementation of SiF NASA (v26), the O₂ A band is excluded, less PCs are used, and cloudy ocean spectra are selected to determine PCs for atmospheric transmission. Tested in isolation in the experiments discussed here, these three differences would all increase retrieved SiF values. However, the water vapour absorption band is also excluded in SiF NASA (v26) and this change would mean a decrease of retrieved SiF values as compared with our implementation. In the direct comparison of the SiF KNMI and SiF NASA (v26) data sets where all effects are combined, we see that SiF NASA values are smaller that SiF KNMI values (smaller slope in Figure 8). This suggests that the water vapour band has a strong impact on the retrieval. Further research may indicate whether this explains the low correlations over tropical rainforests.

7. General Conclusions

In this paper, we describe the satellite-based retrieval of a sun-induced far-red vegetation fluorescence derived from GOME-2A spectra. We followed the general retrieval approach described in Joiner et al. [21] and implemented our own retrieval setup that differs in a number of important aspects from version 26 of the NASA data set. In particular, we have used desert only as reference areas for determining the PCs to describe atmospheric absorption optical thickness. Furthermore, reference spectra were collected from a year of data. We used a much wider fit window (712–783 nm) that also covers water vapour absorption up to about 740 nm and absorption in the oxygen A band. Finally, we needed a larger number of PCs because atmospheric transmission shows more spectral variation across our fit window.

The two time series correlate well with an overall R of 0.78. However, we also observe that over tropical rain forests correlations drop to background levels. Subsequent sensitivity experiments investigating the effect of the number of PCs, the fit window, and the reference area for PC determination, suggest the strong effect on the retrieval outcome of whether or not the water vapour absorption band is included in the fit window. Future research may indicate whether this or insufficient natural variability (compared to noise levels) explains the drop in correlations for tropical rain forests.

Furthermore, we have compared the SiF time series with Gross Primary Productivity derived from flux tower measurements in Australia. The correlation between SiF and site GPP ranged from 0.37 to 0.84. Managed biome types typically have higher correlations than natural types. Finally, by using the de-seasonalized SiF time series for Australia, the break of the Millennium Drought during the local summer 2010/2011 could be detected. Vegetation activity increased from January 2007 to April 2011 due to the transition from droughts to increased rainfall in the summer 2010/2011. Decreases in SiF were seen after April 2011 and may be due to less rainfall and higher temperatures towards 2015. Analysing SiF time series illustrates the strong potential of SiF data to monitor vegetation activity in relation with meteorological anomalies.