Systematic Orbital Geometry-Dependent Variations in Satellite Solar-Induced Fluorescence (SIF) Retrievals

While solar-induced fluorescence (SIF) shows promise as a remotely-sensed measurement directly related to photosynthesis, interpretation and validation of satellite-based SIF retrievals remains a challenge. SIF is influenced by the fraction of absorbed photosynthetically-active radiation at the canopy level that depends upon illumination geometry as well as the escape of SIF through the canopy that depends upon the viewing geometry. Several approaches to estimate the effects of sun-sensor geometry on satellite-based SIF have been proposed, and some have been implemented, most relying upon satellite reflectance measurements and/or other ancillary data sets. These approaches, designed to ultimately estimate intrinsic or physiological components of SIF related to photosynthesis, have not generally been applied globally to satellite measurements. Here, we examine in detail how SIF and related reflectance-based indices from wide swath polar orbiting satellites in low Earth orbit vary systematically due to the host satellite orbital characteristics. We compare SIF and reflectance-based parameters from the Global Ozone Mapping Experiment 2 (GOME-2) on the MetOp-B platform and from the TROPOspheric Monitoring Instrument (TROPOMI) on the Sentinel 5 Precursor satellite with a focus on high northern latitudes in summer where observations at similar geometries and local times occur. We show that GOME-2 and TROPOMI SIF observations agree nearly to within estimated uncertainties when they are compared at similar observing geometries. We show that the cross-track dependence of SIF normalized by PAR and related reflectance-based indices are highly correlated for dense canopies, but diverge substantially as the vegetation within a field-of-view becomes more sparse. This has implications for approaches that utilize reflectance measurements to help account for SIF geometrical dependences in satellite measurements. To further help interpret the GOME-2 and TROPOMI SIF observations, we simulated cross-track dependences of PAR normalized SIF and reflectance-based indices with the one dimensional Soil-Canopy Observation Photosynthesis and Energy fluxes (SCOPE) canopy radiative transfer model at sun–satellite geometries that occur across the wide swaths of these instruments and examine the geometrical dependencies of the various components (e.g., fraction of absorbed PAR, SIF yield, and escape of SIF from the canopy) of the observed SIF signal. The simulations show that most of the cross-track variations in SIF result from the escape of SIF through the scattering canopy and not the illumination.


Introduction
Solar-induced fluorescence (SIF) is produced by the photosynthetic machinery of terrestrial vegetation and is emitted across a range of wavelengths spanning red to far-red (near-infrared, NIR) wavelengths. It has been measured with a number of aircraft-and ground-based instruments (e.g., [1,2] and references therein). SIF has also been mapped globally with several different satellite sensors over a range of spatial and temporal resolutions (see [2] for a review). Canopy-level SIF in the far-red spectral region, where most current satellite-based SIF observations are reported (wavelengths > ∼735 nm), displays a highly linear relationship with respect to gross primary production (GPP), the amount of carbon dioxide taken up by vegetation during photosynthesis, at weekly or longer time-scales (e.g., [3][4][5][6][7][8][9][10]); this result is also supported by modeling studies [11,12]. Therefore, satellite-based SIF data sets are of interest for carbon cycle and terrestrial ecosystem research and have been used in a number of studies (e.g., [7,[13][14][15][16][17]). In contrast to traditional broad-band reflectance measurements, SIF has shown an excellent sensitivity to highly productive regions such as the US corn belt [3] and to vegetation in regions affected by snow cover [18][19][20][21].
SIF can be expressed in a congruent form of the light-use efficiency model [22,23] as SIF(λ em , SZA, VZA, RAA) = e(λ em , SZA, VZA, RAA) Φ F (SZA, λ em ) where λ em is the SIF emission wavelength (peak emission at 740 nm for far-red spectral emission feature), e is the fractional amount of leaf-level emitted fluorescence that escapes the canopy at λ em in the direction of the observer, Φ F is the fluorescence yield, fPAR chl is the fraction of photosynthetically-active radiation for fluorescence (PAR F ) absorbed by chlorophyll, SZA is the solar zenith angle, VZA is the view zenith angle, and RAA is the relative azimuth angle between the sun and satellite. Another angle of interest is the phase angle,γ, the angle at a given point between the sun and satellite: The minimum amount of shadowing occurs at zero phase angle, also known as the hot spot. Equation (1) is sometimes further simplified by combining e and Φ F into a single canopy-level yield quantity for a given sun-sensor geometry. fPAR chl and e are determined by canopy structure and leaf biochemistry, while physiological effects, including responses to various types of stress, are embedded in Φ F (e.g., [24]). Yang et al. [25] defined the product of fPAR and e as Γ rt that they termed the radiative transfer factor.
At present, data producers typically provide an adjustment factor to compute daily-average SIF from a single time of day measurement with the assumption that Φ F , fPAR chl , and e are constant throughout the day [50]. Daily-averaged SIF from different satellite instruments has been directly compared (e.g., [51]). While this accounts for differences in the top-of-atmosphere incoming irradiance, it does not account for structurally-dependent aspects of the illumination and their dependence on the ratio of direct to diffuse radiation. It also does not account for the fact that only a fraction of the total emitted SIF will be observed in a given direction at the top of canopy and that this fraction depends upon the observation geometry relative to the illumination.
Several approaches have been proposed and/or implemented to more fully account for sun-satellite geometrical dependencies of satellite SIF measurements (e.g., [25,30,[37][38][39][42][43][44][45]52]). These methods rely on various ancillary data including reflectances along with theoretical and/or machine learning constructs (see e.g., [25] for a review). They have been tested using simulated measurements and observations, showing promise for dealing with angular variations in global satellite SIF data sets. Various limitations of these methods, such as for sparsely vegetated scenes where there is a substantial contribution from soil reflectance, have also been discussed (e.g., [25,44,46] and references therein).
Here, we examine the sun-sensor geometry dependence of far-red SIF, near infrared (NIR) reflectance, and two reflectance-based indices that have been proposed as proxies for Γ rt with large swath satellite sensors: the Global Ozone Monitoring Experiment 2 (GOME-2) and the TROPOspheric Monitoring Instrument (TROPOMI). We consider the range of local solar times (and thus solar zenith angles) observed systematically across the large swaths of GOME-2 and TROPOMI. Although the host satellites of GOME-2 and TROPOMI have equator crossing times that differ by hours, their orbital attributes allow SIF retrievals from these sensors to be directly compared at the same local overpass times at high latitudes in the northern hemisphere during summer, sometimes with similar viewing geometry. Using a one dimensional (1D) canopy radiative transfer model, we examine the geometry dependences of SIF, its basic components in Equation (1), as well as NIR reflectance and related indices including systematic satellite cross-track variations.

GOME-2 SIF
There are now three GOME-2 instruments flying on the operational European Meteorological Satellites (EUMETSAT MetOp) series satellites [53]. Here we used data from the MetOp-B GOME-2 (GOME-2B), launched in 2012, that operates in the nominal wide swath mode concurrently with TROPOMI. GOME-2B uses a scanning mirror to measure across its swath that has a width of 1920 km. An individual pixel in its forward scan across the swath is 40 km × 80 km at nadir and the backward scans are approximately three times larger. GOME-2 has four spectral bands ranges. Band 4 measures solar radiance and Earthshine radiance in the 593-790 nm range with 1024 spectral elements and a full-width at half maximum (FWHM) bandpass per spectral element of 0.48 nm. This band was used to retrieve SIF.
SIF was retrieved with GOME-2 using variants of a principal component analysis (PCA) technique [54][55][56]. Similar PCA approaches have been used in SIF retrievals with other satellite and ground-based instruments (e.g., [35,57,58]). Here, we used the GOME-2B SIF from the version 27 (v27) GOME-2 NASA dataset [54,59] produced with updates as noted below applied to the base unadjusted retrievals provided in that data set. This version uses a spectral fitting window covering 734-758 nm that has relatively weak water vapor absorption. Along with SIF, reflectances are provided at 670 and 780 nm. We used only the highest quality GOME-2B retrievals that passed all quality control checks described in Joiner et al. [59]. GOME-2B data are limited to solar zenith angle (SZA), < 75 • and local solar times earlier than 13:00.
One change from v27 is that a new bias adjustment scheme is implemented with a machine learning neural network (NN) approach to remove documented biases [55]. These biases are thought to primarily result from temperature variations that occur in the instrument along its orbit that may lead to changes in the instrument response function or other effects that that mimic the spectral response of SIF. This is applied to the uncorrected data from the v27 data set. Five inputs to the regression approach of Joiner et al. [59] are used for the NN; the inputs are cosines of the view zenith angle (VZA), SZA, and relative azimuth angle (RAA) as well as latitude, and reflectance at 780 nm. The single output is the bias, defined as the retrieved far-red SIF over ocean for all pixels, both clear-sky and cloudy, with the assumption that SIF for the ocean pixels should be zero (an assumption justified by satellite measurements). This bias is then subtracted from the terrestrial SIF retrievals. The basic NN structure is a three layer feed forward network that was found by trial and error to produce reproduce the basic features of the bias. The two hidden layers in the NN architecture each contain eight nodes or two times the number of inputs. Additional details of the NN structure are given in Appendix A.
Training was done separately for each day of GOME-2B observations to estimate the bias adjustment. There were of the order of 600,000 samples available for training on a typical day. Previous versions of the regression bias adjustment scheme showed small discontinuities at latitudes boundaries where the scheme was applied; these were removed with the NN approach. The monthly mean bias adjustment for June 2018 (see Figure A1) shows a significant north to south gradient of more than 0.5 mW m −2 nm −1 sr −1 as well as a dependence on surface reflectance. This latitudinal dependence is thought to result from temperature changes in the instrument along its orbit that may affect the instrument response function. Note that many of the predictors of the bias adjustment (reflectance and the sun-satellite angles) overlap with those that we will relate to SIF measurements below. This bias adjustment was thus critical for proper interpretation of the GOME-2B data that follows. Because the bias adjustment uses data over ocean, where SIF was expected to be zero, it was completely independent of SIF emissions over vegetated land.
An effective cloud fraction, f c , was estimated by inverting where I meas , I clr , and I cld were the measured, clear-sky, and cloudy sky radiances at a particular wavelength (here we use 680 nm). Within the context of the mixed Lambertian-equivalent reflectivity (MLER) model and the independent pixel approximation (IPA), clouds were considered as opaque with a reflectivity of 0.8 (e.g., [60]), and f c is an effective rather than a geometrical cloud fraction. The MLER concept was shown to reproduce the correct amount of absorption or Rayleigh scattering for pixels with either optically thin and/or broken clouds and has been used in numerous trace-gas and cloud algorithms (e.g., [61,62]). At our wavelengths of interest, aerosols are not highly absorbing and scatter light similar to cloud particles. As such, they were treated like clouds within the MLER framework.
The calculation of f c has also been updated from v27. In v27, a climatology was used for I clr , and the f c calculation did not include atmospheric scattering or surface bidirectional effects. In this work, we used the collocated MODIS MCD43D reflectance parameters [63][64][65] to account for both of these. Specifically, we trained an NN to predict I clr at 680 nm. The NN has the same basic architecture as that applied to the bias adjustment (see Appendix A) but with a different number of inputs and hidden-layer nodes. The six inputs are the three bidirectional reflectance distribution function (BRDF) parameters for MODIS band 1 (620-670 nm) as well as the cosines of the SZA, VZA, and RAA. The number of nodes in each of the hidden layers was twelve or two times the number of inputs.
The training used ten days (June 20-30, 2018, 43605 GOME-2 samples) of global collocated MCD43D and GOME-2B data. Only GOME-2B pixels with 670 nm reflectance < 0.6 were considered for training and of those reflectance data, only the 10% with the closest agreement to the band 1 reflectances were used. This was found to provide a sufficient amount of data for training. I meas from the v27 data set was used to compute f c . We use a strict filter of f c < 0.1 in our analysis involving SIF and reflectance data from GOME-2.

TROPOMI SIF
TROPOMI was launched on the European Space Agency (ESA) Sentinel-5 precursor (S5P) satellite in October 2017. TROPOMI measures solar irradiance and Earthshine radiance across a 2600 km swath with 448 spatial elements [66]. SIF is retrieved from TROPOMI pixels (3.5 km × 7 km at nadir) using a subset of the NIR band 6 that covers wavelengths 725-755 nm in 497 spectral elements with a FWHM of 0.38 nm. The TROPOMI SIF retrieval is similar to that of GOME-2 v27 in that it uses a PCA approach. The fitting window for the TROPOMI SIF retrieval is 743-758 nm [67]; this simplifies the retrieval as this wavelength range was nearly free of absorption by atmospheric molecules and thus absorption of the SIF signal does not need to be accurately accounted for. Quality control was applied as in Köhler et al. [67].
The TROPOMI retrievals currently include cloud information derived from collocated Visible Infrared Imaging Radiometer Suite (VIIRS) data on the Suomi National Polar-orbiting Partnership (NPP) flying in tandem with the S5P satellite. A cloud fraction for a TROPOMI pixel was defined by weighting the cloud fractions of all the VIIRS pixels contained within the TROPOMI pixel as described in Köhler et al. [67]. Therefore, the TROPOMI cloud fraction was not the same as f c computed above for GOME-2B. R NIR is not explicitly provided in the TROPOMI data set, but reflectance across the SIF retrieval fitting window can be computed using the absolute and relative SIF values provided in the data set as R 740 nm =SIF/(relative SIF)× (100×π)/(cos(SZA) ×F), where F is the solar irradiance at 740 nm and relative SIF is defined as SIF relative to its reflectance background, defined as a percent. Table 1 summarizes the main observing characteristics of the two SIF-capable instruments considered here. Both the instruments and retrievals differ in several respects. Below, we provide to our knowledge the most detailed comparison of these two data sets to date by accounting for differences in their orbital geometries.  Figure 1 shows SZA, VZA, the phase angle, and local solar times (LST) of the GOME-2B and TROPOMI observations over single orbits on 28 June 2018. The inclination angles (near 98.7 • ) and nodes of the host satellites impact the local overpass times along their orbits. With GOME-2B in a descending orbit (daylight on descending node), LSTs are generally later than the 09:30 equator crossing time in the northern hemisphere (NH); for TROPOMI in its ascending node, LSTs are generally earlier than its 13:30 equator crossing time in the NH. The reverse situation occurs in the southern hemisphere.

GOME-2B and TROPOMI Instrumental and Orbital Characteristics
The LSTs corresponding to individual ground footprints vary across the wide swaths of these instruments. This leads to an interesting scenario for SIF measurements; in the NH boreal summer, particularly at middle to high latitudes, there is coverage spanning a wide range of LSTs throughout the morning and early afternoon with the two instruments albeit with systematically varying sun-sensor geometries. There are measurements at the same local time with the two instruments at high latitudes in the NH, sometimes with similar viewing geometries and other times not.

Data Analysis Approach
We calculated the zonal means of the native GOME-2 or TROPOMI SIF retrievals and related parameters for 10 • latitude, ten minute bins, for different plant functional types (PFTs). We did this for the month of June 2018 around the peak of the northern hemisphere growing season. To compute the zonal averages, we first averaged the data monthly into 0.5 • latitude × 0.5 • longitude grid boxes then aggregated the gridded data into the appropriate latitude, time, and PFT bins. Only data with cloud fractions <0.1 were used. We repeated our analysis using all data with cloud fractions <0.3 and did not find qualitatively different results although correlations between SIF/PAR and reflectances were somewhat lower on average. At the highest high northern latitudes examined here in summer months (60-70 • ), GOME-2B and TROPOMI provided nearly continuous temporal coverage from 09:00 to 13:00. We computed the percentages of PFTs in each of our 0.5 • gridboxes using 2017 land cover types on the climate modeling grid (CMG, 0.05 • grid) derived from MODIS observations (MCD12C1 product) [68]. Only gridboxes that have a majority PFT with coverage >50% were used in our analyses.
The two approaches proposed to serve as proxies for Γ rt that we evaluate are 1) the far-red fluorescence correction vegetation index (FCVI) defined as the difference between the NIR directional reflectance (R NIR ) and the broadband visible (VIS) directional reflectance (R vis ) proposed by Yang et al. [25], and 2) the NIR V , defined as the product of R NIR and the normalized difference vegetation index (NDVI = (R NIR − R red )/(R NIR + R red )), i.e., NIR V = NDVI × R NIR , as proposed by Zeng et al. [44]. As the GOME-2 SIF data set includes R NIR and R red , we were able to readily compute NIR V collocated with SIF. We computed the DVI (DVI= R NIR − R red ) as a proxy for FCVI.

Simulations with a 1D Canopy Radiative Transfer Model
To help interpret the satellite SIF observations, in particular the sun-satellite geometry dependences, we examined the results of SIF and reflectance simulations, including R NIR , NIR V , and DVI, from the Soil-Canopy Observation Photosynthesis and Energy fluxes (SCOPE) model version 1.70 [34]. SCOPE simulates the full radiative transfer of fluorescence from the photosystem to the top of canopy as well as reflected sunlight through the canopy. We provide a list of SCOPE input parameters for two different scenarios used in our simulations in Appendix B: (1) a C4 corn (maize) setup with a 2 m canopy height; and (2) a C3 30 m canopy designed to mimic a deeper forested type canopy. These are referred to as the C4 and C3 setups. We ran simulations for two values of leaf area index (LAI) for each setup. For the C4 case, we performed simulations for LAI = 3 and 1, and for the C3 case, we used LAI = 5 and LAI = 1. Most other parameters were similar for the two scenarios.
We note that SCOPE is a 1D model, meaning that it models vegetation as homogeneous, i.e., LAI is distributed evenly throughout the pixel with no clumping. Therefore, it does not fully capture the complexity of radiative transfer for complex and varied canopies that exists within large satellite footprints. Nevertheless, SCOPE is a useful tool for modeling general sun-satellite geometrical dependences of both SIF, as well as its individual components from Equation (1), and reflectance for a given canopy structure. It can therefore be used for a qualitative comparison with satellite measurements.
Note that the top-of-canopy irradiance in our SCOPE runs was specified as a single value regardless of the output at different solar zenith angles. That input irradiance value corresponds to a single atmospheric realization (determined by an external atmospheric radiative transfer calculation using the MODTRAN code) at an SZA of 45 • . To adjust the incoming irradiance to appropriate conditions at other SZAs, the SCOPE SIF outputs were multiplied by a factor of cos(SZA)/cos(45 • ). With further SCOPE simulations (1) adjusting the input irradiance to the actual magnitude of solar irradiance to fully account for the effects of irradiance on Φ F and (2) using separate external radiative transfer calculations at the appropriate SZA, we determined that the simple irradiance magnitude adjustment to SCOPE outputs was appropriate (see Appendix B for details).

GOME-2 and TROPOMI Satellite-Based SIF and Reflectance Dependence on Illumination and View Geometry
Figures 2-7 show SIF normalized by incoming PAR (plotted as SIF/cos(SZA) and denoted as SIF/PAR), R NIR , and satellite and solar angle information as a function of LST (or alternatively VZA) for late morning through early afternoon averaged over gridboxes dominated by different plant functional types for different latitude bins in June 2018. DVI and NIR V are shown only for GOME-2 as R red was not available in the TROPOMI SIF data set. The standard error of the mean, i.e., standard deviation of all gridbox values (σ) divided by √ N for each binned zonal average is shown by the vertical error bar and represents an uncertainty due to random errors. A conservative empirical estimate of the systematic error (0.1 mW/m 2 /nm/sr for SIF/PAR) is shown with the shading. Correlations between the binned SIF/PAR averages and the corresponding R NIR and available reflectance-based indices are also provided, separately for GOME-2 and TROPOMI. Correlations will similarly be computed using SCOPE 1D simulations below and then may be compared with those from satellite observations. Individual satellite SIF observations are dominated by retrieval noise and this makes it difficult to interpret results at the pixel level; this is why we have conducted our analysis on zonal averages. The correlations are computed for the binned averages and contain only a limited number of points per PFT and latitude range. We note that even though we average hundreds of individual pixels per bin, one or two noisy or questionable binned data points may substantially degrade the correlations. Error bars indicate the standard error, hatching shows additional estimated uncertainties due to systematic instrument and algorithmic effects. Also shown for reference are the cosines of SZA and phase (γ) angles. Correlation coefficients (r) for reflectance (near infrared), DVI, and NIR V with respect to SIF/PAR are provided for GOME-2. Only NIR reflectance is shown for comparison with TROPOMI SIF/PAR as reflectance at a red wavelength needed to compute DVI and NIR V was not available in the TROPOMI data set. SZA: solar zenith angle.     Particularly at the higher northern latitudes, the top-of-atmosphere solar irradiance, that is proportional to cos(SZA) and indicated in the figures, does not change dramatically over the hours shown (varies by ∼20% from 09:00-13:30). The systematic variation in VZA as a function of LST is also apparent owing to the fact that LST varies across the satellite swath. The cos(γ) variation across the swath is also shown with the highest values for GOME-2 (lowest amount of shadowing) on the western part of its swath and the highest values on eastern swath portion for TROPOMI.
For GOME-2, the highest values of SIF/PAR, R NIR , NIR V , and DVI generally occur on the early morning (west) edge of the swath that has the highest cos(γ) values as well as the highest VZAs and SZAs (lowest values of cosine SZA). For TROPOMI, the highest SIF/PAR values are generally seen on the later afternoon (east) edge of the swath, with the highest cos(γ) values. For the areas dominated by woody savannas, mixed forests, and deciduous broadleaf forests in the range 30-40 • N where cos(γ) is highest, there is a peak in SIF/PAR as well as R NIR on the right hand side of the swath for TROPOMI observations, corresponding to a maximum in cos(γ), approaching the hot spot where cos(γ) = 1. Particularly for TROPOMI, SIF/PAR at the higher northern latitudes over grasslands shown in Figure 3 and also for croplands shown in Figure 4, a "U"-shaped curve emerges that roughly corresponds with the VZA.
R NIR derived from GOME-2 and TROPOMI (as well as NIR V and DVI from GOME-2) shows similar though not identical cross-track patterns as SIF/PAR; highest values for GOME-2 on the west part of the swath corresponding to the lowest phase angles and a peak for woody savannas at the middle latitudes also corresponding to minima in phase angle. Note especially for the large GOME-2 footprints that have perhaps more residual cloud contamination than TROPOMI, R NIR and reflectance-based indices would more affected by cloud contamination than SIF/PAR [69].
An interesting phenomenon occurs near 11:30 LST for the highest latitude bin (60-70 • N); observations are made by both satellite instruments at similar VZAs and γ angles. Under these conditions, SIF/PAR values from TROPOMI and GOME-2B agree nearly to within their estimated uncertainties for most PFTs. Note that this is not the case for woody savannas (Figure 2) where GOME-2 SIF values are higher at 60-70 • N latitude than at 50-60 • N in contrast with TROPOMI.
The purpose of the binning process in our analysis was to reduce the effects of instrumental noise, which are relatively large for individual SIF retrievals. We note that vegetation is inhomogeneous within each satellite footprint, gridbox, and PFT bin. In addition, vegetation phenology and structure within a given PFT bin may vary substantially across a given latitude bin [49]. We conducted our analysis by additionally binning by longitude, computing averages separately for three different longitude bins of 60 • . This had an impact of increasing noise for GOME-2 across the swath while leading to more smooth behavior of SIF/PAR across LST for TROPOMI. GOME-2 results tended to look a bit noisier (more random bumps when plotting SIF as a function of LST) likely owing to the fact that with fewer pixels GOME-2 binned averaged are affected by retrieval noise. In contrast, TROPOMI binned averages looked more smooth as a function of LST, presumably owing to decreased inhomogeneity; with many more pixels, TROPOMI binned averages are less impacted by retrieval noise. The overall magnitudes of SIF/PAR varied for some PFTs in different longitude bins, but the general variations with respect to satellite swath geometry were similar across all longitude bins.  (1) vary with sun-satellite geometry that aid in the interpretation of the GOME-2 and TROPOMI SIF swath dependences. With SCOPE we were able to vary certain angles (e.g., VZA) while holding others (e.g., SZA) constant. We were not able to reproduce these analyses with GOME-2 and TROPOMI owing to the cross-track viewing geometry that varies systematically with SZA. Figure 8 shows SCOPE simulations of various radiation-and canopy-related parameters in Equation (1) that vary only as a function of SZA for the C4 canopy simulation with LAI = 3. APAR variations with SZA are dominated by PAR rather than fPAR.  Figure 8b shows the wavelength-integrated canopy SIF yield, defined as observed total SIF divided by APAR. A SIF emission efficiency factor, defined as the hemispherically averaged SIF at 740 nm divided by APAR chl , is also shown in Figure 8b. Both quantities are normalized by their maximum values. These show a very similar SZA dependence with the highest values at high SZAs. However, there is little diurnal variation in the simulated canopy SIF yield or emission efficiency through a large range of irradiances (SZA < ∼50 • ). Figure 9 shows for the same setup a SIF directional escape factor at 740 nm (i.e., in the direction of the observer), defined as TOC SIF at 740 nm divided by SIF emitted by all photosystems at 740 nm normalized by π. This quantity is shown as a function of SZA at a constant VZA (40 • ) and as a function VZA for a constant SZA (40 • ) for a range of RAA from 0 • to 180 • . The hemispherically averaged SIF divided by the SIF emitted by all photosystems is also shown with heavy black lines and only varies with SZA. The hot spot is visible at angles of 40 • for the curves at or near the principal plane (e.g., the red solid line for the principal plane). The SIF directional escape variation with SZA is relatively small except in the principal plane (red solid curve) near the hot spot. The directional escape factor varies more with VZA, particularly near the hot spot.      Figure 11a shows simulated TOC SIF normalized by PAR (to a value of incident PAR equivalent to that at 45 • SZA) as a function of SZA. PAR normalized SIF generally increases with SZA in contrast with the unnormalized SIF. This dependence is a result of the increasing fPAR and SIF yield with SZA as shown in Figure 10. Figure 11b shows how much of the total PAR normalized TOC SIF is coming from sunlit leaves, shaded leaves, and SIF scattered within the canopy. Here, we see that only a small fraction comes from shaded leaves and of the remainder, it is split similarly between scattered SIF and SIF from sunlit leaves, with an expected peak in SIF from sunlit leaves at the hotspot and a corresponding drop to zero from shaded leaves at the hotspot.   Figure 13 shows how SIF/PAR and its components vary with LAI for the C4 setup. For the LAI = 1 case, fPAR chl is lower but increases more rapidly with SZA, producing a larger SZA dependence. The canopy yield is higher overall for the LAI = 1 case but has a similar SZA dependence as LAI = 3. SIF/PAR has a larger dependence on VZA and SZA for the LAI = 1 case. The larger VZA dependence for LAI = 1 is due to the effect of seeing more leaves and less ground at high VZA that is more important for lower LAI values. The magnitude of SIF from the sunlit leaves is similar in both cases but there is more scattered SIF from the LAI = 3 canopy. This leads to higher SIF for the LAI = 3 case in general as may be expected, though not at the highest SZAs and VZAs. Similar differences between the deeper canopy C3 LAI = 5 and LAI = 1 cases are shown in Figure A3.  (b) wavelength integrated SIF yield; decomposition of total normalized SIF into that from sunlit leaves, shaded leaves, and from SIF scattered within the canopy as a function of (c) solar zenith angle and (d) view zenith angle. Figure 14 shows SCOPE simulations of TOC SIF, SIF/PAR, R NIR , NIR V , and DVI as a function of VZA for SZAs and phase angles corresponding to typical swath sun-satellite geometries for GOME-2 and TROPOMI; VZAs with negative values correspond to the west sides of the satellite swaths. SCOPE calculations computed at 10 SZAs, 10 VZAs, and 19 RAAs were interpolated to the average geometries shown in Figure 3 for grasslands in the 30-40 • latitude band corresponding to June 2018. TOC SIF is shown with and without normalization for incoming PAR. The PAR normalization produces cross-track SIF variations that are more similar to that of reflectance and reflectance-based indices for GOME-2, but overall PAR normalization does not much change the cross-track SIF variation. Similar cross-track behavior of SIF/PAR and reflectance-based indices are shown for the C4 case for geometries corresponding to the 50-60 • latitude band in Figure A4 and for the deeper canopy C3 LAI = 5 and LAI = 1 cases in Figure A5.   Figure 14. SCOPE simulations of SIF, SIF/PAR, R NIR , NIR V , and DVI for a maize setup (C4) to mimic typical sun-sensor geometry across a swath (shown as a function of view zenith angle but at each angle for appropriate values of solar zenith angle and phase angle corresponding to the geometry for 30-40 • N in June for (a,c,e) GOME-2B (left) and (b,d,f) TROPOMI (right) for LAI values of 3 (top), 1 (middle), and 0.5 (bottom). SIF/PAR shows SIF normalized to a constant solar zenith angle representative of the center of the swath. Correlations (r) between R NIR , NIR V , and DVI with respect to SIF/PAR are also displayed.

Comparison of GOME-2 and TROPOMI SIF
That the GOME-2 and TROPOMI SIF values overlap (to nearly within estimated uncertainties) when viewing conditions are similar (at the highest latitudes) provides confidence in the consistency of the relative calibration and retrievals from the two satellite-based instruments. It should be noted that the satellites do not sample the Earth in the same way; TROPOMI has greater sampling in general owing to smaller pixels that tend to have less cloud contamination. Also, TROPOMI, with its wider swath can provide greater sampling at the highest latitudes. It should also be noted that NIR reflectances from GOME-2 and TROPOMI are not reported at the same wavelength; GOME-2 reflectance is reported near 780 nm while TROPOMI reflectance is computed across its fitting window near 740 nm. This may contribute some of the differences in behaviors of SIF and reflectance between the two satellites.

Geometrical Dependencies of SIF and Reflectance
The SCOPE simulated variation in SIF from nadir to a VZA of −45 • is substantial at ∼40% for the LAI = 3 case for the GOME-2 geometries and less so for the simulated TROPOMI swath and for the LAI = 1 and LAI = 0.5 cases as seen in Figure 14. The LAI = 1 and LAI = 0.5 cases show a distinctly more symmetric bowl-like cross-track dependence similar to that seen in the TROPOMI data for grasslands, particularly at the higher latitudes (above 40 • N). From these simulations, LAI is shown to be an important factor in determining the cross-track SIF dependence in the satellite data. We note here again that the 1D SCOPE model may not capture all of the complexity in a 3D canopy as observed from satellite and the current version does not allow for inhomogeneity in the vertical structure.
Substantial cross-track asymmetry in both observed SIF/PAR and reflectance is shown for the deeper canopy forest PFTs that tend to have higher SIF/PAR and reflectances related to higher LAI values. For some of these canopies (forested) with particularly high SIF/PAR and reflectance, we can also see a well defined hot spot in both SIF and reflectance. This is particularly the case at the middle latitudes as phase angles approach zero. For TROPOMI, the hot spot occurs on the right (East) part of the swath and produces a peak in the cross-track SIF and reflectance (see for example, Figures 2 and 7 for woody savanna and deciduous broadleaf, respectively, at latitudes between 30 • N and 50 • N). For GOME-2, the hotspot occurs near the western most edge of the swath, leading to a high degree of asymmetry across its swath for these PFTs and latitudes. This is generally captured in the SCOPE simulations (Figure 14a,c).
A bowl-like cross-track pattern that corresponds roughly to VZA with less asymmetry is observed by GOME-2 and TROPOMI for croplands, grasslands, and evergreen needleleaf forests. This variation is particularly apparent at the higher latitudes where SIF values tend to be lower in general, corresponding to lower values of LAI. This is explained by scattering processes affecting both SIF and reflected light within the canopy (i.e., mutual shadowing) and by the bidirectional gap fraction [37,46,70] and is captured in the SCOPE simulations (Figure 14a,c). SCOPE simulations show that SIF/PAR, R NIR , NIR V , and DVI should all share a similar cross-track dependence for satellite instruments, particularly at higher values of LAI (1 < LAI < 5). However, SCOPE simulations also show that the SIF/PAR and R NIR cross-track variations begin to diverge at lower LAI values (0.5). For example, correlations between SIF/PAR and R NIR are near to one for LAI = 3 in Figure 14 but drop below 0.9 for the LAI = 0.5 case. The similarities between SIF/PAR and R NIR were explained by (e.g., Yang and van der Tol [39]) through the relationship between canopy scattering of SIF and reflectance (see also e.g., [43,44]). This relationship may be straight-forward under certain conditions that are more likely to hold for far-red SIF examined here than for SIF at red wavelengths. Derivations of straight-forward SIF-reflectance relationships are generally developed under the condition of "black soil". In reality, soil (as well as other non-green components of the canopy such as woody branches) has a direct contribution to the observed reflectance through the single scattering effects, but does not directly contribute to SIF as long as the soil is not fluorescing (e.g., Yang and van der Tol [39], Colombo et al. [49]). The derived straight-forward relationship between canopy scattering of SIF and reflectance therefore will start to break down with increasing gap fraction. SCOPE simulations of NIR V and DVI show higher correlations with SIF/PAR, consistent with results of Zeng et al. [44] and Yang et al. [25].
This aspect of the SCOPE simulations is generally consistent with the observed correspondence (or lack thereof) between SIF/PAR, R NIR , NIR V , and DVI across the tracks of GOME-2 and TROPOMI. GOME-2 and TROPOMI measurements show similar cross-track behavior of SIF/PAR, R NIR , NIR V , and DVI for the deciduous broadleaf forests in Figure 7 and croplands with DVI values > ∼0.12 (typical for latitudes > 30N) in Figure 4. In contrast, TROPOMI and GOME-2 observations show more mixed results for cross-track correlations between SIF/PAR and reflectance-based parameters for grasslands and other PFTs that display lower values of SIF/PAR and DVI.
There are generally larger differences between SIF and reflectance for GOME2 and TROPOMI measurements than in SCOPE simulations. In addition, neither DVI nor NIR V has consistently higher correlations with SIF/PAR, particularly at lower values of DVI corresponding to lower LAI canopies. There are several possible explanations. We note that SCOPE simulates a homogeneous canopy, i.e., LAI is distributed evenly within a simulated scene. Therefore, the contribution of the directly scattered light from the soil surface to the reflectance is typically underestimated for sparse vegetation or scenes with larger gap fractions. While SCOPE would produce a larger divergence between SIF and reflectance at even lower LAI values (<0.5), the gap fraction would still likely be underestimated for inhomogeneous vegetation. It should also be noted that reflectances are more sensitive than SIF to the effects of clouds [69], and this could lead to some of the differences between the cross-track variations of SIF/PAR and reflectance-based indices. In addition, there may still be residual biases in the SIF/PAR retrievals as the signal to noise ratios are fairly low. Three dimensional effects that are not captured by SCOPE may play a role. Finally, we note that there could be variations in SIF yield as a function of local solar time (illumination). SCOPE simulations generally do not support that the large cross-track differences between SIF/PAR and reflectances for lower LAI scenes would result from SIF yield variations across the track.
Differences in the geometrical dependences of observed SIF/PAR and reflectance-based indices appear to be substantial in many cases corresponding to low DVI or LAI values. This is a concern for methods that rely upon reflectance or reflectance-based indices to normalize geometrical variations in SIF observations from wide swath instruments to estimate a more intrinsic property of SIF, such as SIF emitted at the photosystem level. While some approaches attempt to address the issue of soil reflectance (e.g., [43][44][45]), it is not clear to what extent and under which conditions these methodologies applied to satellite data can provide the needed accuracies for various applications. Yang et al. [25] found that FCVI, NIR V , and FCVI × NDVI (an attempt to improve upon FCVI to better account for soil effects) all showed errors in approximating Γ rt in the most challenging cases of sparse canopies with low LAI with SCOPE simulations. They suggest that errors of less than 30% could be obtained for cases of FCVI > 0.18. Our results are generally consistent with these conclusions. Further work is needed to evaluate these methods with satellite data under a variety of conditions including at different spatio-temporal scales. One method of evaluating such approaches will be to examine the resulting cross-track dependences of the derived intrinsic SIF components (for example normalizing SIF by PAR and fPAR). Averaging over and across large regions, as done here, is needed to mitigate the effect of retrieval noise, which is large even for the highest quality satellite SIF data currently available.
Ground-based measurements over various canopies at different geometries also show different bidirectional distributions of SIF including asymmetric bowl-like patterns for erectophile and spherical canopies and a more dome-like shape with a peak at the hotspot for phanophile canopy types [37,46]. These measurements showed varied amounts of correlation between SIF and NIR reflectance. Our satellite measurements for high LAI canopies are generally consistent with these observations. Biriukova et al. [46] also performed SCOPE simulations and showed that the bidirectionality of SIF is mainly controlled by the leaf inclination distribution function (LIDF); the magnitude and shape of the hotspot is controlled by the LAI and the ratio of leaf width to canopy height, with low ratios and high LAI leading to a stronger and more pronounced hotspot, as is present in our satellite measurements of deciduous broadleaf and mixed forests.
In summary, our analysis of GOME-2 and TROPOMI data show that for very dense canopies such as deciduous broadleaf forests and mature croplands, R NIR , NIR V , and DVI all show a very similar cross-track patterns as SIF/PAR. We are not able to discern if any of these reflectance-based indices provides a superior approximation of Γ rt . For sparse and even some moderately-vegetated canopies, the cross-track behavior of all of these parameters diverges substantially from that of SIF/PAR, in contrast to the SCOPE simulations. It is unclear how well the reflectance-based indices will work for scenarios in between these two extremes and if a single threshold value derived from the reflectances will apply generally. The fluorescence yield (Φ F in Equation (1)) should be PAR dependent and different for sunlit and shaded leaves. As Φ F also varies with stress, the SIF yield derived for a particular geometry may vary from the canopy SIF yield. This may contribute to the differences between cross-track SIF and reflectance observations and 1D simulations and needs to be considered for interpretation of directional SIF.

Implications for the Use of Time Series of Gridded Averages from Large Swath Sensors
Many studies have utilized time series of gridbox SIF averages from wide-swath sensors including GOME-2 to investigate interannual variations such as the effects of drought (e.g., [13,71,72]). Our analysis of GOME-2 and TROPOMI SIF data indicates that some of the day-to-day SIF variability for a given location arises from daily changes in the illumination and viewing geometry as the orbits of host satellites do not follow a common ground track every day. The use of time-series of gridded SIF data from such sensors is valid so long as there is adequate and consistent sampling over the swath throughout the record of interest. This is generally the case for monthly data from GOME-2; a typical 0.5 • gridbox over a month encompasses about 10-20 samples. Higher spatial and/or temporal resolution averages would be appropriate for TROPOMI.
One change in sampling of a long record SIF data set occurs with the GOME-2A (on MetOp-A); its observing mode changed from the nominal large swath to a reduced swath with smaller pixels in July 2013. To assess the impact of the change in sampling, we computed gridded SIF monthly means with GOME-2B using all swath positions and using only those with VZA < 35 • similar to the small swath mode GOME-2A. The GOME-2B swath is the same width as that of GOME-2A in its earlier observation mode. The bias between the two samples is very small (restricted VZA averages lower by 0.02 mW m −2 nm −1 sr −1 , see Figure A6). This provides some confidence in time series analyses of monthly mean gridded data with GOME-2A (biases due to instrument degradation not withstanding).
Zhang et al. [40] pointed out that caution should be applied when using SIF from the various observation modes of the Orbiting Carbon Observatory 2 (OCO-2) that have different samplings of VZAs (not necessarily symmetric about the nadir view). Our results with TROPOMI span the range of observing conditions of OCO-2 and agree with the findings of Zhang et al. [40] that the different geometries sampled by the three observation modes (nadir, glint, and target) should produce systematic differences in SIF. The general caution on the use of different observing modes for scientific studies, supported by our results, also applies to OCO-3 as well as GOSAT; the latter underwent a change in its cross-track swath sampling during the mission.

Conclusions
We have shown that for low Earth orbit (LEO) sensors measuring far-red SIF, variations of sun-satellite geometry across the satellite swath are also related to local solar time of the observation, making it difficult to disentangle physiological effects, including those related to diurnal variations in light-use efficiency, from sun-satellite geometry effects. We show that there are substantial systematic variations in far-red SIF across the satellite swaths of GOME-2 and TROPOMI related to sun-sensor geometry and that these variations differ according to canopy properties, in particular the LAI. The variations are roughly consistent with those predicted using the SCOPE canopy radiative transfer model as well as previous ground-based measurements. SIF tends to increase with decreasing phase angle and increasing VZA and SZA. Cross track variations in SIF/PAR are highly correlated with R NIR , NIR V , and DVI for dense canopies such as deciduous broadleaf forests near the summer solstice and crops in a mature state as predicted by SCOPE simulations. This is encouraging as these parameters have been proposed as a means to account for the angular dependencies of SIF in order to extract the embedded physiological information. However, the observed differences in cross-track dependences of SIF are not well captured with any of the reflectance-based indices for sparsely vegetated canopies as well as some moderately vegetated canopies. We note that the physiological status of the vegetation is unknown, and this may affect the cross track correlations shown here as well. More testing of these methods needs to be applied to satellite-as well as ground-based data under a variety of conditions in order to develop reliable thresholds for which the proposed approaches would be applicable. We recommend that both NIR and red wavelength reflectances be provided with satellite SIF data sets if possible for further studies. We find that GOME-2 and TROPOMI far-red SIF values agree to nearly within their estimated uncertainties where their sun-satellite geometries are similar (high northern latitudes in summer). We also note that it is imperative to continue to identify biases and improve bias correction schemes within the SIF retrieval process.
SIF geometrical dependencies will also be important for SIF-capable sensors in geostationary Earth orbit (GEO). Several such sensors are planned to be launched over the next decade including the NASA Earth Ventures Tropospheric Emissions: Monitoring Pollution (TEMPO) and the Geostationary Carbon Observatory (GeoCarb) as well as the ESA Sentinel 4 ultraviolet visible near-infrared (UVN) sounder [73]. These sensors will face a different scenario as compared with those in LEO. For GEO sensors, VZAs at a given point on the ground will be fixed, while the SZA will change throughout the day. A significant portion of the GEO coverage will have moderate to high VZAs (latitudes higher than ∼30 • ) [74] where the SIF dependence on SZA and phase angle will be substantial. There will also be a need to connect SIF measurements from the different GEO instruments as well as the GEO measurements from a single sensor at different locations and thus with different VZAs. SIF measurements from LEO will be helpful in this respect provided that the sun-sensor geometry can be accurately accounted for. We used an adaptive moment estimation optimizer to minimize the error function with a learning rate of 0.1 and default values of β 1 and β 2 , parameters used by the algorithm (0.9 and 0.999, respectively). The main advantage of the neural network for bias correction is that it eliminated small discontinuities along latitude lines used as boundaries in the regression approach of Joiner et al. [59]. Figure A1 shows averaged SIF (gridded) bias adjustments computed for June 2018. The bias is seen to have a strong latitudinal gradient. Bias adjustments are larger on average over land in the northern hemisphere owing to higher average reflectances, particularly for Greenland.    Figure A2 shows three different methods of computing a solar zenith angle dependence within SCOPE. In the first simulation (blue line), the magnitude of the irradiance remains the same for all simulations, while the incident angle changes. The adjustment to the appropriate magnitude of irradiance was applied to the SCOPE-modelled SIF afterwards. In this scenario, SCOPE runs with an SZA-independent irradiance, and thus variations of Φ F due to irradiance were not accounted for by the model. This is what was done in the simulations shown in the main part of the paper. We tested the accuracy of this approach with two other simulations. The pink curve shows results when the magnitude of the irradiance was adjusted to the appropriate SZA before input into SCOPE. In this scenario, the model also simulates variations in Φ F due to SZA. This somewhat reduces the sensitivity to SZA, because the model predicts a lower Φ F when irradiance is high. However, the difference is not large, and this also suggests that there is not much signal of diurnal variation in Φ F in canopy level SIF. The black line shows SCOPE results using separate MODTRAN simulations (for a given aerosol loading, etc.) for each SZA as input to SCOPE instead of correcting the magnitude (with cosine of the SZA). This takes into account the fact that the contribution of diffuse radiation also changes with SZA. This has an even smaller effect than whether a correction is made to either inputs or outputs.  Figure A2. Three different approaches to account for the SIF solar zenith angle dependence in SCOPE showing very little differences in the approaches (see text for more details). In these simulations, the viewing geometry was nadir. Figures A3 and A5 show SCOPE simulation results for the C3 setup, contrasting results for LAI = 1 and LAI = 5 cases. These figures have counterparts for the C4 setup in the main part of the paper and show similar differences between high and lower LAI cases. For example, the hotspot effect is smaller for the LAI = 1 case, but the VZA and SZA dependences are generally larger as shown in Figure A3, similar to those shown in Figure 13. The C3 cross-track dependence, shown in Figure A5 displays less asymmetry for the LAI = 1 case, similar to what was shown for C4 in Figure 14. (b) wavelength integrated SIF yield; decomposition of total normalized SIF into that from sunlit leaves, shaded leaves, and from SIF scattered within the canopy as a function of (c) solar zenith angle and (d) view zenith angle.      Figure A5. Similar to Figure 14 but for the C3 setup; SCOPE simulations of SIF and reflectance for the C3 setup to mimic typical sun-sensor geometry across a swath (shown as a function of view zenith angle but at each angle for appropriate values of solar zenith angle and phase angle corresponding to the geometry for 30-40 • N in June for LAI=5: (a) GOME-2B (left) and (b) TROPOMI (right); and LAI=1: (c) GOME-2B (left) and (d) TROPOMI (right). SIF PAR normalized shows SIF normalized to a constant solar zenith angle representative of the center of the swath.

Appendix C. Monthly Means from GOME-2B with Different Samplings of VZA
To check whether the change in operations of GOME-2A from nominal to small swath mode in July 2013 produces a bias in the time series, we examined one month of GOME-2B data where we gridded the values using all swath positions versus only those pixels with VZA < 35 • . The results are shown in Figure A6. The overall bias between the two different sampling schemes is very small. Therefore, the change in GOME-2 operations is not expected to cause substantial discontinuities in gridded SIF time series. None-the-less, the change in observation mode, leading to changes in swath width and viewing angles, should be mentioned when using GOME-2A time series SIF data for scientific studies.