Canopy structure explains the relationship between photosynthesis and sun-induced chlorophyll fluorescence in crops

1 Department of Landscape Architecture and Rural Systems Engineering, Seoul National University, South Korea 2 Department of Global Ecology, Carnegie Institution for Science, Stanford, USA 3 International Institute for Earth System Sciences, Jiangsu Provincial Key Laboratory of Geographic Information Science and Technology, Nanjing University, Nanjing, China 4 Dynamic Meteorology Laboratory, Ecole Polytechnique, Palaiseau, France 5 National Center for AgroMeteorology, Seoul, South Korea


Introduction
Sun-induced chlorophyll fluorescence (SIF) is increasingly used as a proxy for gross primary productivity (GPP) at large scales (Mohammed et al., 2019;Ryu et al., 2019). SIF is an optical signal emitted in the spectral range 650-850 nm from chlorophyll a molecules in vegetation (Lichtenthaler et al., 1986). The far-red part of SIF at about 760 nm, on which our study strongly concentrates, has a well-documented, empirical relationship to estimate GPP at both the site Yang et al., 2015;, regional (Guanter et al., 2014;Yao Zhang et al., 2016), and global scales (Frankenberg et al., 2011;Guanter et al., 2012;Joiner et al., 2011). The precise reason for the SIF-GPP relationship at the canopy scale, however, lacks a clear mechanistic explanation. This is mostly due to an insufficient understanding of the relative contributions of leaf physiological and canopy structure effects to SIF and how the physiological and structural components of SIF relate to photosynthetic light use efficiency. Therefore, it is helpful to revisit the basic processes and equations relevant for SIF and GPP as the basis for more closely examining the possible mechanisms that might underlie their strong empirical correspondence.
SIF and GPP differ in one important respect with regard to their fundamental processes. GPP is related to leaf-level gas exchange processes and therefore, observed top-of-canopy GPP is simply the cumulated GPP of all leaves, as gases that might temporarily accumulate in the canopy eventually will diffuse out of it. Far-red SIF, however, is an optical signal in the near-infrared spectral range, where light is strongly scattered by leaves allowing only a certain fraction to escape the canopy (Knyazikhin et al., 2013;Yang and van der Tol, 2018;Zeng et al., 2019). Therefore, the top-of-canopy SIF as observed from tower, airborne or satellite platforms is not simply the cumulative signal of SIF emitted by all leaves but contains an extra term quantifying the effect of canopy scattering. While reabsorption of SIF photons is a strong effect for red SIF around 690 nm, this effect can be neglected for the far-red SIF at 760 nm (Yang and van der Tol, 2018), which we exclusively consider in this study.
Both SIF and GPP can be understood conceptually in the light use efficiency framework originally introduced for net primary productivity (Monteith, 1972;Monteith and Moss, 1977). Thus, for GPP we have where GPP is defined as the product of the photosynthetically active radiation (PAR) absorbed by the canopy (APAR) and the photosynthetic light use efficiency of the canopy (LUE P ). Similarly, for SIF as observed above the canopy (SIF obs ), we have where Φ F is the physiological SIF emission yield of the whole canopy and f esc is the fraction of all SIF photons, emitted from all leaves, that escape from the canopy (Guanter et al., 2014;Zeng et al., 2019).
When comparing the basic equations for GPP and SIF, Φ F corresponds to the LUE P term and f esc is the extra term that quantifies the effects of canopy scattering.
Several experimental and modelling studies demonstrated that APAR is the dominant factor in the variability of SIF obs and the relationship of SIF obs to GPP at sub-daily temporal resolution (Li et al., 2020;Miao et al., 2018;Wieneke et al., 2018;.
However, the separate contributions of Φ F and f esc have not been quantified so far to the best of our knowledge.
Modeling studies from as early as 2012 recognized that f esc plays a significant role in controlling the amount of SIF observed at top-of-canopy (Fournier et al., 2012). In terms of functional dependencies, f esc has been shown to respond strongly to changes in both leaf area index (LAI) (Fournier et al., 2012;Yang and van der Tol, 2018) and leaf angle distribution (Du et al., 2017;Migliavacca et al., 2017;Zeng et al., 2019). More generally, it follows from the results of Zeng et al.
(2019) that any canopy structure parameter that influences the near-infrared reflectance can have a considerable effect on f esc . In particular, apart from LAI and leaf angle distribution, the clumping index is expected to play an important role for f esc (Zeng et al., 2019).
Despite the advances in our theoretical understanding of f esc , little work has gone into understanding what, if any, role f esc has in explaining the SIF-GPP relationship. Instead, the effects of f esc for relating SIF to GPP have largely been ignored in the SIF literature (e.g. Guanter et al., 2014;Wieneke et al., 2016;Yang et al., 2017Yang et al., , 2015. However, a number of more recent studies have reported direct evidence of distorting effects of f esc on the SIF obs -APAR relationship (Du et al., 2017; and indirect evidence of how f esc affects the SIF obs -GPP relationships using both process-based modelling and observations (Migliavacca et al., 2017). While it is clear from the literature and Eqn. 2 that f esc partially obscures the Φ F signal, there is no decisive conclusion so far whether f esc is helpful or detrimental for GPP estimation. Two apparently opposing views on this question are presented in more detail in the following paragraphs.
Several studies have argued that canopy structure effects on SIF obs should be corrected for the purpose of optimal GPP estimation based on the assumption that Φ F has a positive relationship to LUE P (Du et al., 2017;Yang and van der Tol, 2018). In terms of quantitative evidence, however, such reasoning has been largely based on the improvement of the APAR-SIF relationship when accounting for canopy structure (Du et al., 2017;. However, in considering only the SIF-APAR relationship, these studies are insufficient for evaluating the possible role that f esc might have in explaining the SIF-GPP relationship. In particular, as GPP estimation has been the ultimate goal of most SIF research, GPP needs to be explicitly considered in the analysis. Logically, an improvement in the APAR-SIF relationships could very well go together with a degradation of the SIF-GPP relationship; the better the relationship of SIF to APAR, the smaller the variation of the efficiency term SIF/APAR but LUE P = GPP/APAR actually represents a special case of an efficiency term with considerable variation. Moreover, we are unaware of any experimental or modelling studies that explicitly considered the relationship between Φ F and LUE P at the canopy scale. Instead, several studies have evaluated the relationship between LUE P and the joint influence of canopy physiology (Φ F ) and canopy structure (f esc ), as obtained by dividing canopy escaping SIF obs by APAR (Li et al., 2020;Eqn. 2;Miao et al., 2018;Wieneke et al., 2018Wieneke et al., , 2016Yang et al., 2015). Wieneke et al. (2018) and Li et al. (2020) attempted to empirically correct for canopy structure effects in the apparent SIF yield by dividing with the MTVI2 vegetation index that was shown to be related to LAI.
However, such an approach not only assumes a linear relationship between LAI and f esc but also neglects the potentially strong impact of temporally varying leaf angle distribution and clumping.
Interestingly, the reasoning that Φ F underlies the SIF-GPP relationship is not well supported by several strands within the existing SIF literature. At the leaf scale, Φ F is known both to vary less than LUE P and to be nonlinearly related to LUE P in light response curves (Gu et al., 2019;van der Tol et al., 2014). Similarly, in situ measurements have shown low seasonal correlation of Φ F and LUE P (Goulas et al., 2017). At the canopy level, these findings are supported by several experimental studies which have shown higher correlation of SIF obs to APAR compared with SIF obs to GPP Wieneke et al., 2018;.
In contrast to the Φ F -based reasoning above highlighting the role of leaf physiology, (Badgley et al., , 2017 have argued that it is precisely the canopy structure that explains the SIF obs -GPP relationship. Their reasoning is based on the NIR reflectance of vegetation (NIR V ), a multi-spectral reflectance-based measurement that is strongly linear with both SIF obs and GPP at large spatial and long temporal scales. However, the NIR V -GPP relationship has so far not been tested with groundlevel spectral observations at the shorter time scales that include diurnal variations of both APAR and LUE P . Nor have previous studies of NIR V gone beyond documenting the empirical NIR V -GPP relationship. The strong connection between the NIR reflectance and f esc Yang and van der Tol, 2018;Zeng et al., 2019) hints at potential links between f esc and LUE P which could be investigated more directly.
For the sake of clarity and simplicity, the two opposing views presented above can be condensed into the following two hypotheses: (Hphys) dominance of leaf physiology: Φ F is correlated to LUE P , while f esc is not and, therefore, the product of APAR and Φ F outperforms SIF obs for GPP estimation.
(Hstruc) dominance of canopy structure: f esc is correlated to LUE P , while Φ F is not and, therefore, the product of APAR and f esc outperforms SIF obs for GPP estimation.
Our goal in this study is to test which of these two hypotheses is better supported by observations.
To do this, we applied the f esc estimation approach of Zeng et al. (2019) to three independently collected datasets in rice, wheat and corn, which consist of continuous measurements covering the whole crop growing season. As the above hypotheses represent the two extreme cases, it is clear that the truth might actually lie somewhere in between the two. Furthermore, the results are expected to depend on the time scale as canopy structure changes that affect f esc , which mostly occur at the seasonal time scale, while physiological changes that might affect Φ F occur over the course of individual days as well as over the growing season . Our analysis therefore considers these different time scales

Theoretical framework: decomposing far-red SIF into its mechanistic parts
The basic equation for canopy-level SIF obs in terms of its separate mechanistic components was given in the introduction (Eqn. 2). For convenience in later sections, we define LUE F = SIF obs /APAR as the apparent light use efficiency of SIF in analogy to the definition of photosynthetic light use efficiency, LUE P = GPP/APAR. The escape fraction of whole canopy far-red SIF emissions can be estimated from NIR v and fPAR following the approach outlined by Zeng et al. (2019): where NIR V is the product of NIR reflectance and the NDVI (Badgley et al., 2017). Zeng et al. (2019) demonstrated the good performance of Eqn. 3 with comprehensive radiative transfer simulations as well as supporting evidence from satellite observations. In particular, they showed that f esc derived from NIR V performs well even for sparsely vegetated canopies and is minimally affected by changes in soil brightness.
When fPAR, PAR and NIR V observations are available, in addition to SIF obs , we can combine Eqns. 2 and 3 to calculate the following two SIF-related variables. First, the product of APAR and Φ F , which represents the total emitted SIF and includes signals emitted from the leaf adaxial and abaxial parts, and is therefore a mixture of the canopy radiation absorption and the physiological emission efficiency (Fig. 1). Second, the product of APAR and f esc , which represents a combined canopy structure and radiation factor and includes both absorption and scattering aspects of the canopy structure ( Fig. 1). While f esc is a pure canopy structure property for far-red SIF, APAR is, in principle, not only affected by canopy structure but also by leaf pigments. dfs However, the impact of leaf pigments on APAR was found to be marginal and strongly dominated by canopy structure characteristics such as leaf angle distribution, clumping and LAI (Luo et al., 2019). Only when distinguishing total vs. green APAR would we expect to see an impact of leaf pigments on certain parts of our analysis, which is discussed in detail in section 4.4. We first conducted analyses focusing on Φ F and f esc to directly examine the roles of canopy scattering and physiological processes. In a second step, we investigated the relationships of SIF obs , APAR × Φ F , and APAR × f esc to GPP. This second step is important as it quantifies the contributions of Φ F or f esc combined with APAR to the estimation of GPP (Fig. 1). All relevant SIF-related variables as well as Φ F and f esc can be derived from Eqns. 2 and 3 as shown in Table 1. In particular, it is important that APAR × f esc is estimated using only NIR V and PAR without any contribution from SIF obs .

In-situ data sets
Three in-situ data sets were combined for the analysis. The data sets differ not only in terms of crop type but also in terms of geographical location, instruments used, observation geometry, and retrieval method used to estimate SIF obs . An overview is given in Table 3 and more detailed descriptions are provided in the following subsections. Half-hourly data from a single growing season for each crop were used. All sites had in situ observations of APAR and high spectral resolution instruments for SIF retrieval.  Goulas et al. (2017). For the hemispheric field-of-view configuration in the rice paddy, the footprint diameter is given as range from the cumulative 50% to 80% of the total footprint. The whole row for wheat is shaded in grey for easier visual distinction between the different rows.

Rice
The rice paddy site is located in Cheorwon, South Korea (38.2013°N, 127.2506°E) and is part of the national eddy flux network, KoFlux (Huang et al., 2018). Oryza sativa L. ssp. Japonica is grown in rows each year there, with a maximum canopy height of 70 cm and maximum LAI of 6 in the 2016 growing season that lasted from April 29 to September 3. Measurements instruments are operated by Seoul National University and the National Center for Agro Meteorology. An eddy covariance system consisting of a three-dimensional sonic anemometer (Model CSAT3, Campbell Scientific Inc., Logan, UT, USA) and a closed-path infrared gas analyzer (Model LI-7200, LI-COR Inc., Lincoln, NE, USA) was used to measure CO2 fluxes. Net CO2 flux partitioning into gross primary production (GPP) and ecosystem respiration was conducted according to the night time-based method (Reichstein et al., 2005). SIF was monitored with a very high spectral resolution instrument (full width at half maximum, FWHM = 0.17 nm, QEpro, Ocean Optics, Dunedin, FL, USA) and a fiber switch to measure up-and downwelling irradiances in sequence (K. . The singular vector decomposition (SVD) method (Guanter et al., 2013)

Wheat
The wheat field is located close to Avignon in southeastern France ( The site is also part of CarboEurope (Dolman et al., 2006) which provided the eddy flux data. The eddy covariance system consists of a three-dimensional sonic anemometer (Model 81000, R.M. Young Company, Traverse City, MI, USA) and an open path infrared gas analyzer (Model LI-7500, LI-COR Inc.).
Flux partitioning was done based on day time values (Kowalski et al., 2003). SIF and canopy reflectance was observed with the TriFLEX instrument that consists of separate spectrometers for measurements of SIF and VIS-NIR hyperspectral reflectance (Daumard et al., 2010). The spectrometer for SIF observations has a high spectral resolution (FWHM = 0.5 nm) while the broadband reflectance is observed at lower resolution (FWHM = 2 nm). The TriFLEX instrument was mounted on a crane at 21 m above ground, was operated at nadir observation angle and has a narrow field of view resulting in a footprint diameter of about 2 m (Goulas et al., 2017). In contrast to the system in the rice paddy, downwelling solar irradiance was observed using a white reference panel. SIF was retrieved from high frequency observations (about 1 Hz) using a modified Fraunhofer Line Depth (FLD) method that is described in detail in the corresponding references (Daumard et al., 2010;Goulas et al., 2017) and then averaged over 30 minute periods. fPAR was continuously monitored using a network of ten quantum sensors. The data was acquired in the year 2010. More details on the site and the methods used can be found in Daumard et al. (2010) and Goulas et al. (2017).

Corn
The corn field is located in Henan province, in central China (34.5199° N, 115.5916° E) and is also used for growing winter wheat in a rotation system. Corn (Zea mays ssp. Mays) was sown in rows in late June and harvested in early October and had maximum canopy height of 2.5 m and maximum LAI of about 3.5 (Li et al., 2020). The measurement instruments are operated by Nanjing University and the Farmland Irrigation Institute of the Chinese Academy of Agricultural Sciences. An eddy covariance system consisting of a three-dimensional sonic anemometer (WindMaster Pro, Gill Instruments Limited, Hampshire, UK) and a closed path infrared gas analyzer (LI-7500RS, LI-COR Inc., Lincoln, NE, USA) and was continuously operated. Eddy flux partitioning was done according to the night timebased method of Reichstein et al. (2005). A similar system as in the rice paddy based on the FluoSpec2 design (X. Yang et al., 2015) was used for SIF retrievals (Ocean Optics QEpro spectrometer with FWHM = 0.17 nm). SIF was retrieved from the observations taken every two minutes with a spectral fitting method (SFM) (Meroni et al., 2010;Meroni and Colombo, 2006), before averaging the observations over 30 minute intervals. Integration times were optimized before each measurement. In addition to the high spectral resolution instrument, a lower spectral resolution (FWHM = 1.1 nm) spectrometer covering the VIS-NIR range was used to continuously monitor canopy reflectance changes (HR 2000+, Ocean Optics, Dunedin, FL, USA). The latter was also operated with a shutter system to switch between observations of down-and upwelling radiation. In contrast to the rice paddy observations, the field of view of both spectrometer systems for upwelling observations was 25° as bare fibers were used. As the observations were located about 10 m above the canopy, the measurement footprint circle had a diameter of about 4.4 m (Li et al., 2020

Data processing and statistical analysis
For all three crop datasets, only time steps that had measurement values for APAR, SIF obs at 760 nm, NIR V and GPP between 8 am and 4 pm were selected in order to ensure direct comparability of the correlation values for relationships to GPP. For rice and corn, most of the growing season satisfied these criteria, while the wheat site only had data after the green-up phase (supplementary figure Fig.   S1). All datasets had some gaps and gap-filled GPP data was not used in any of the presented analyses.
For all analyses, we further restricted the data by setting a threshold for fPAR above 0.45 in order to exclude the noisy data in the early growing stage. This affected primarily the rice dataset and was only strictly necessary for analyses of LUE P , LUE F , f esc and Φ F but nevertheless also done for the corresponding variables that include APAR to be consistent. In addition to the fPAR threshold, we applied the Hampel outlier filter (Hampel, 1974) to detect and remove several severe outliers in the variables normalized for APAR. The Hampel filter uses the median absolute deviation and an estimate of the standard deviation within the moving window as robust distance metrics. We applied the filter with a conservative threshold parameter value of four and a window length of 12 days. We used the hampel function from the R package pracma (Borchers, 2018) and slightly modified it to accept data with gaps.
For the diurnal correlation and daily mean value calculation, only days with a minimum of five valid data points were selected. Combined with the fPAR threshold, this resulted in 62, 49 and 57 days of data for rice, wheat and corn, respectively (Table 3).
f esc is reported as unitless quantity, which means that in case of the nadir viewing geometry of wheat and corn, it is equal to the π times directional ratio of SIF obs to APAR × Φ F . This is a similar procedure as chosen for the reflectance factor, that uses a Lambertian surface as reference (Schaepman-Strub et al., 2006). This can also be understood in terms of NIR reflectance in Eqn. 3, which is used for the estimation of f esc .
SIF obs of the rice dataset was converted from irradiance to radiance units for better comparability with the other two datasets. This was done by a simple division by π, which can only account for the overall magnitude difference but not the more subtle differences in viewing geometry related to asymmetric angular patterns that differ from a Lambertian surface. For a discussion of these effects, see section 4.4. NIR V was calculated from the lower spectral resolution VIS-NIR spectrometers by averaging over the range 800-850 nm for the NIR band and over 600-650 nm for the red band ( Table   1). The impact of shifting the spectral range of the red band towards the red-edge (e.g. 620-670 nm, 650-700 nm) were evaluated but found to have negligible effects.
All correlation analyses were done based on the Pearson correlation coefficient of which the statistical significance was evaluated with a two-sided t-test at a confidence level of 95%. The squared Pearson correlation coefficient (R 2 ) is equal to the coefficient of determination in the special case of linear regression with one explanatory variable and an intercept term (Devore, 2011;Kasuya, 2019).
Therefore, e.g. the GPP estimation performance of different predictor variables in such linear regression models can be evaluated based on R 2 . Statistical analysis was done using the programming language R (R Core Team, 2012).

Temporal patterns
We found clear differences in both the seasonal patterns and the degree of diurnal variability of f esc , Φ F , LUE F and LUE P (Figs. 2, S2). fPAR , f esc and LUE F all appeared to be mainly characterized by considerable seasonal variation (Figs. 2 and S2). LUE P and Φ F , however, showed considerable diurnal variability. In case of LUE P , the diurnal variation was superposed to the seasonal changes, while for Φ F it seemed to be the main source of variation.
Despite some apparent general patterns for each variable, we found considerable differences between crops. Patterns of fPAR increase during green-up and decrease during senescence were similar for rice and corn while for wheat, where the data for the green-up phase were not available, in around 0.5-0.6 for all crops and minimum values were in the range of 0-0.1. As already mentioned above, Φ F was seasonally rather constant but wheat showed a clear decrease late in the senescence stage that coincided with a decrease in chlorophyll content (Goulas et al., 2017). For corn, a slight decrease followed by an increase late in the season were observed, while for rice, slight decreases appeared during both the late green-up and senescence. LUE F overall had similar seasonal patterns as f esc but also showed some smaller seasonal patterns similar to Φ F in addition to more diurnal variability from the latter (results not shown). The seasonal patterns for LUE P were partly masked by strong diurnal variations but decreases during senescence could be observed for all three crops ( Fig.   2j-l). While the decreasing patterns of LUE P for wheat and corn were similar and showed relatively steep slopes, the decrease in rice occurred at an earlier stage and after that LUE P was rather stable.
Increases in LUE P during green-up were only clearly evident in the case of rice and to a lesser degree for wheat. The figure corresponding to Fig. 2 in terms of daily mean values highlighting only the seasonal component of variation and including LUE F is shown in the supplementary material (Fig. S1).

Relationships to LUE P
To more directly test our two hypotheses, we built on the qualitative analysis of temporal patterns shown in Fig. 2, by quantifying the correlations of Φ F and f esc to LUE P . We found that in the relationships to LUE P , the following pattern of increasing correlation was consistent for all crops and for both half-hourly and daily mean values at the seasonal time scale: Φ F < LUE F < f esc (Fig. 3 and Fig.   B1 in the appendix).
For half-hourly values, i.e. seasonal+diurnal variability, Φ F had no correlation to LUE P for all three crop types. f esc , in contrast, had significant and moderate correlations, with squared Pearson correlation values (R 2 ) of 0.28, 0.28 and 0.44 for rice, wheat and corn, respectively. The corresponding values for LUE F were generally on the order of half the value of f esc (Fig. B1). Overall, LUE F was most strongly related to Φ F but wheat was an exception where the correlation of f esc was stronger (Fig. S2).
f esc and Φ F were almost uncorrelated although wheat and corn showed significant but low values (R 2 of 0.14 and 0.06, respectively). For the sake of completeness, we included full correlation tables between all four efficiency terms in the supplementary material (Fig. S2).
When considering the correlations of f esc and Φ F to LUE P for the seasonal variability based on daily mean values, we found that f esc outperformed Φ F in all cases but with considerably larger differences for the seasonal variability (Figs. 3, B1). While significant correlations between Φ F and LUE P did not increase much for the seasonal variability compared to the seasonal+diurnal variability (a large increase for rice was not significant), the corresponding f esc -LUE P correlations strongly increased up to 0.6 for wheat and corn and 0.4 for rice and were highly significant (Fig. 3).
The correlations to LUE P for diurnal variability only showed different patterns compared to the seasonal variability for Φ F and f esc . For Φ F , the correlations based on diurnal variability exceeded the R 2 values that included seasonal variability for all three crop types, but were all below 0.2 and considerably higher for wheat than for rice and corn (Fig. B1a). For f esc , the diurnal correlations were consistently lower than those for seasonal+diurnal variability and seasonal variability, although the differences were smaller for rice (Fig. B1c). The diurnal variability case was the only one where LUE F was not consistently intermediate between Φ F and f esc . In particular, for wheat, LUE F was slightly more strongly related to LUE P than either Φ F or f esc separately and for rice, LUE F had a slightly lower correlation to LUE P than Φ F (Fig.   B1b).

Relationships to GPP and APAR
The first part of the results focused on the relationships of Φ F and f esc to LUE P . However, the ultimate goal of most SIF-related research is GPP estimation. Therefore, we investigate the implications of the relationships of Φ F and f esc to LUE P for GPP estimation by comparing the relationships of APAR multiplied with Φ F , f esc or their product to GPP. APAR itself is used as a reference case both as predictor of GPP and as target variable for estimation.
There were clear differences in the GPP estimation performance of SIF obs , APAR × f esc , and APAR × Φ F , with APAR × f esc having the strongest correlation with GPP ( Figs. 4 and 5a) while the roles of rice and corn were exchanged for the differences between APAR × Φ F and APAR × f esc (Fig. S3a). While APAR × f esc showed a slight tendency to saturate at high values in the case of rice and corn, we did not observe this pattern for wheat (Fig. 4d,h,l). Furthermore, the SIF obs -GPP relationship showed clear signs of saturation for corn and wheat but not for rice ( Fig. 4c,g,k). For rice and corn, the differences in GPP estimation performance of different predictor variables were mainly caused by deviations from the ideal case of perfect correlation for values in the middle and high part of the GPP range ( Fig. 4a-d,i-l). For wheat, however, the deviations from the ideal case was more notable also in the low to middle parts of the GPP range ( Fig. 4e-h). Relaxing the fPAR threshold criterion (see section 2.3) did not change the relative GPP estimation performance of APAR × Φ F , APAR × f esc , and SIF obs and only had marginal effects on the R 2 values for wheat and corn (results not shown). However, as rice had a prolonged green-up phase resulting in a large number of low GPP values, the performances of predictors for the full dataset without applying the fPAR threshold were considerably higher for all variables (0.1≤ΔR 2 ≤0.18), in particular had an R 2 value of 0.84 (Fig. S2).
For GPP estimation, the performance rank of APAR compared to SIF-related variables differed for wheat, where APAR performed worst, while the order APAR × Φ F < APAR < SIF obs in terms of R 2 held for rice and corn (Figs. 4 and 5a). The differences between APAR and SIF obs in terms of absolute increase in R 2 were smallest for corn and largest for wheat. An overview of the time series for all GPP predictor variables and all crop types is given in Fig. A1 in the appendix. Figure 4: Overview of relationships between gross primary productivity (GPP) and SIF -related variables at half-hourly temporal resolution for the whole growing season (as far as data was available). Absorbed photosynthetically active radiation (APAR), and the two other SIF-related variables. GPP (µmol m -2 s -1 ) is compared with APAR (µmol m -2 s -1 ) in panels a, e, i, the structure and radiation factor ( × f esc , units: µmol m -2 s -1 ) in panels b, f, j, the total emitted SIF ( × Φ F ,units: mW m -2 nm -1 sr -1 ) in panels c, g, k, and with observed SIF (SIF obs , units: mW m -2 nm -1 s -1 ) in panels d, h, l. Relationships are shown on the basis of half-hourly data covering one growing season for each crop. Significant squared Pearson correlation values (R 2 ) are shown for reference. All data are shown as partially transparent, filled circles in order to visualize the density of points.
When considering relationships to APAR rather than GPP, both APAR × f esc and APAR × Φ F had a stronger correlation to APAR than SIF obs (Fig. 5b and Fig. C1 in the appendix). For APAR × Φ F , the improvement was largest for wheat with differences in R 2 reaching 0.2, intermediate for corn and smallest for rice (Figs. 5b, S3b). For APAR × f esc , the improvement in relationships to APAR was on a similar level for rice and wheat and somewhat larger for corn, reaching differences in R 2 of 0.2-0.3 with the highest values for the diurnal variability (Figs. 5b, S3d). Overall and for half-hourly values, the relationship between SIF obs and APAR was stronger in rice, than in corn and wheat (R 2 of 0.72, 0.65, 0.62, respectively; Fig. C1) but the differences were smaller for the seasonal variability. A more detailed analysis of the relationships between APAR and either SIF obs , APAR × f esc or APAR × Φ F in terms of patterns in the half-hourly scatter plots is presented in Appendix C.
Overall, the patterns in correlations to GPP and APAR were qualitatively very similar for rice and corn, but differed for wheat (Fig. 5a,b).

Hypothesis evaluation
In this study, we comprehensively investigated the relationships between SIF and GPP by decomposing observed SIF (SIF obs ) into its canopy structure and radiation component, APAR × f esc , and its total emission component, APAR × Φ F , which contains the physiological emission yield term, Φ F (Fig. 1, Table 1). This decomposition allowed us to directly examine the mechanistic basis of the SIF -GPP relationship beyond the overall dominant factor APAR by studying the underlying relationships of f esc and Φ F to LUE P . Relationships of of the SIF-related variables to APAR instead of GPP were also studied for comparison.
Overall, we found strong support for the hypothesis Hstruc, which states that canopy structure, as captured by the escape fraction (f esc ), underlies the observed SIF obs -GPP relationship beyond the dominant and common driver APAR (Fig. 6). Moreover, we found that GPP estimation based on APAR × f esc performed considerably better than or equally well as SIF obs (see Section 4.2 for a more detailed discussion and the link to NIR V ). These findings held consistently at seasonal time scales for the rice and corn datasets and to a lesser degree for the wheat dataset. Even though for wheat data Hstruc was not strictly fulfilled, the correlation of f esc to LUE P was considerably stronger than the correlation of Φ F to LUE P (Fig. B1a,c). We did not find any results directly supporting the hypothesis Hphys, which states that Φ F carries the more relevant information for LUE P estimation and that f esc is a distorting factor. However, the results based on the diurnal variability, especially for wheat, indicated that Hstruc may not hold, in the strictest sense, at the diurnal time scale. This appears to be consistent with the findings of Dechant et al. (2019) who reported larger contributions of SIF at the diurnal than the seasonal time scales for process-based simulations based on the rice paddy dataset.
The correlations to LUE P for diurnal variability showed different patterns compared to the seasonal variability for Φ F and f esc . For Φ F , the correlations based on diurnal variability exceeded the R 2 values that included seasonal (Fig. B1b). For f esc , the diurnal correlations were consistently and considerably lower than those for seasonal+diurnal variability and seasonal variability, although the differences were considerably smaller for rice than for wheat and corn (Fig. B1).
Our study made, to some degree, use of a novel framework of comparative analysis of the GPP vs. APAR estimation performance based on the three predictor variables SIF obs , APAR × Φ F , and APAR × f esc of which a detailed version is presented in the supplementary material (Table S1, Fig. S3).
The motivation for this was to ensure logical consistency of results as well as a plausibility check for a reasonable performance of our f esc estimation approach that has so far not yet been applied to measured data (Zeng et al., 2019). The necessary condition for a reasonable f esc estimation given a considerable seasonal variability is that dividing SIF obs by f esc estimate results in improved APAR estimation performance (Table S1). This was strictly fulfilled for all three crop types and for all variability cases including diurnal and seasonal variations (Fig. S3). Also, the results of rice and corn strictly fulfilled the strong form of Hstruc, while the result of wheat only fulfilled a weaker form of Hstruc (see Table 1).
To ensure the generality of our results for different sky conditions, we examined the influence of the diffuse PAR fraction on our results by analyzing clear sky and cloudy sky data separately (Fig. S5).
These results confirmed our main finding that f esc was much more related to than Φ F . Only for the daily mean wheat data, Φ F showed a much stronger correlation to LUE P for cloudy conditions than for sunny and the combined data, but this was due to a confounding effect of cloudy sky data and the senescence phase, which dominates the correlations to LUE P . More detailed analyses of the effects of the fraction of diffuse light on the SIF-GPP relationship can be found in the original publications of each dataset, i.e. K. , Goulas et al. (2017), and Li et al. (2020).
Our analysis strongly focused on the separate results for each dataset. The motivation for this was to examine the temporal relationships as also observed in a given remote sensing pixel over time. As long as it is not covered by evergreen vegetation, each such vegetation pixel, and each of our datasets, shows considerable variation in the canopy structure parameters such as LAI, leaf angle distribution and clumping, over the growing season. Nevertheless, due to overall differences in e.g. leaf angle distribution between the crop types, differences could be expected in the results for an analysis of all data pooled together. We tested this for the case of only the two C3 crop types and all three crop types by applying a scaling factor to obtain comparable magnitudes of C3 and C4 GPP. We found consistent results compared to the separate analyses of each dataset (Fig. S7). However, when combining datasets from several ecosystems with stronger differences in canopy structure, the results could differ, as dividing by f esc can account for some of the structural differences. This is related to the discussion on the spatial case in section 4.5, where a more detailed reasoning is presented. Figure 6: Schematic overview of the observed relationships between SIF-related terms and either absorbed photosynthetically active radiation (APAR) or gross primary productivity (GPP). The relationships of canopy-level observed SIF (SIF obs ) and its two components, total emitted SIF ( × Φ F = SIF obs /f esc ) and the canopy structure and radiation factor ( × f esc =SIF obs /Φ F ), are shown. In addition to these terms that all contain APAR as main driver, also the corresponding ratio terms that are normalized for APAR are indicated in grey color on top of the flux variables. LUE P is the photosynthetic light use efficiency, LUE F the apparent light use efficiency of SIF obs (LUE F = SIF obs /APAR), f esc is the escape fraction and Φ F is the physiological SIF yield. The line widths between flux terms represent strength of linear correlations (not to scale) based on the results of half-hourly time series at the seasonal time scale and include aspects of seasonal slope stability.
Our study was entirely limited to far-red SIF and therefore we cannot draw any direct conclusions on the analogous results for red SIF. While red SIF data was available for the wheat dataset, estimating f esc for red SIF is more challenging due to the strong reabsorption effects (Yang and van der Tol, 2018).
Findings from process-based simulations for the rice paddy dataset did not indicate stronger performance of total emitted red SIF compared to total emitted far-red SIF .

f esc , × f esc and the important role of NIR V
We found strong evidence of the considerable seasonal dynamics in f esc that corresponded with seasonal variations in LUE P (Figs. 2 and 3). The strong variation of f esc is consistent with previous studies using both process-based simulations and in situ observations (Du et al., 2017;Fournier et al., 2012;Yang and van der Tol, 2018;Zeng et al., 2019) and contradicts studies assuming constant f esc (e.g. Guanter et al., 2014). In particular, the increase in f esc of rice and wheat in the early growing season coinciding with an increase in LAI is consistent with the simulation results of Fournier et al. (2012) and Yang and van der Tol (2018), as well as our own process-based simulations for the rice paddy (results not shown). The overall decreasing pattern for corn, in contrast, cannot be explained by the variation in LAI alone and is likely caused by variations in leaf angle distribution and clumping.
While we do not have detailed measurements of such canopy structure parameters for corn, we observed seasonal changes in leaf angle distribution in rice (results not shown) and seasonal patterns in clumping index for rice were also reported in the literature (Fang et al., 2014). Therefore, the seasonal patterns of f esc appear consistent with the relevant literature and satisfied the plausibility test of improved correlation of SIF obs /f esc to APAR for all three datasets (section 4.1). The diurnal variation of f esc (Fig. 2) was considerably smaller than the seasonal component and can, in general, be caused by interactions of canopy structure with changing solar zenith angle as well as by leaf movements (Goulas et al., 2017).
As fPAR was rather constant for the part of the growing season we studied, it is clear from Eqn. 3 that the variations in f esc were dominated by variations in NIRV. The strong decrease in wheat despite constant total fPAR during senescence might be caused by an increase in the fraction of yellow leaves.
This interpretation is supported by the drastic change in the red to far-red SIF ratio reported in Goulas et al. (2017), which indicates a decrease in chlorophyll content. This also implies a larger difference between green and total fPAR, which could lead to an underestimation in f esc (section 4.4). Apart from pigment contents, changes in leaf angles might partly explain the observed pattern.
In terms of the absolute magnitude of f esc , our results for the rice paddy seem to suffer from an overall underestimation as for hemispheric view, values above 0.5 are expected from the theory. For the wheat and corn datasets, the absolute magnitude is less straightforward to interpret due to the narrow angle viewing geometry, which creates a stronger dependence on leaf angle distribution.
While there might be temporally constant biases in the overall magnitude of f esc , this does not affect the results of our correlation analyses as they are insensitive to a global scaling factor.
We believe our study is the first one that explicitly investigated the relationships of f esc to canopy LUE P , which makes a comparison with other literature somewhat indirect. Nevertheless, a strong link to previous studies on the NIR V -GPP relationship (Badgley et al., , 2017 can be established in the following way. APAR × f esc is calculated as the simple product of NIR V and PAR (Table 1). Therefore, our results can be considered an extension of the NIR V -GPP relationship first identified by Badgley et al. (2017) to the sub-daily scale. Indeed, we could demonstrate that the GPP estimation performance of the product of NIR V times PAR (NIR V P) converges to NIR V at time scales of about two weeks (Fig.  S5). It is noteworthy that in contrast to other studies that simply multiplied NIR V by PAR or short wave radiation without a solid mechanistic basis (e.g. Joiner et al., 2018), we found that NIR V P quite naturally emerges from Eqns. 3 and 5 as an estimate for APAR × f esc (Table 2).
Apart from this way of linking our results to NIR V , it is instructive to consider Eqn. 3 for reinterpreting NIR V in terms of its mechanistic meaning. More specifically, Eqn. 3 can be rearranged to read NIR V = fPAR × f esc , implying that NIR V integrates light absorption and f esc . This is an intriguing consequence of the findings of Zeng et al. (2019) given that the relationship of NIR V to GPP was originally only empirically motivated. Our results add further nuance by establishing an empirical relationship between f esc and LUE P .
Another relevant aspect of our findings is that the APAR × f esc -SIF obs relationship was more consistent than the APAR × Φ F -SIF obs relationship (Fig. S6). This is due to the strong seasonality of f esc and overall seasonal stability of Φ F (Fig. 2). As f esc is the slope in the SIF obs -APAR × Φ F relationship (assuming zero intercept), the seasonal variation of f esc translates into seasonally varying slopes in the SIF obs -APAR × Φ F relationship, similar to the the APAR × f esc -APAR relationship (Appendix C). The SIF obs -APAR × f esc relationship, in contrast has a seasonally stable slope. The relative invariance of Φ F when compared to f esc found across all three datasets helps explain the basis of the globally consistent relationship between SIF obs and NIR V identified by Badgley, Field and Berry (2017) at the monthly time scale.

Φ F and × Φ F
While we found that variability of Φ F , measured in terms of the coefficient of variation (CV), were similar to that of f esc (Fig. B1), the variability of Φ F did not correlate with LUE P for rice and corn. In fact, for these two crops, the considerably high seasonal CV of Φ F appeared to be mainly caused by fluctuations around the mean value without a clear seasonal trend (Fig. S1). This was not the case for wheat, however, which showed more stable values and a clearly decreasing trend coinciding with senescence as well as a weak level of correlation (Figs. 2, 3 and S2). We think that this difference between datasets might be partly explained by the better signal quality of SIF obs in the wheat dataset, which is likely related to the much higher acquisition frequency (Table 3; see Section 4.4 for a more detailed discussion). Nevertheless, our results seem consistent with known, leaf-level patterns of small variation of Φ F at diurnal (Gu et al., 2019) and seasonal time scales (Goulas et al., 2017).
Moreover, Goulas et al. (2017) only found a weak correlation between actively measured Φ F and LUE P at the leaf level, which is consistent with our canopy-level findings (Fig. 3).
We found consistently better APAR estimation performance of APAR × Φ F compared to SIF obs (Figs. 5b and S3b).  previously reported similar results based on site data and airborne images. It is with a certain irony that we conclude from our analyses that, APAR × Φ F , the physiological component of SIF obs , appears to be better suited for estimating APAR as opposed to estimating GPP (Figs. 5, S3). This conclusion is in stark contrast to the widely held hopes of the SIF community to use the inherent physiological link between chlorophyll fluorescence and photosynthesis for improved GPP estimation (see, for example, Porcar-Castell et al., 2014). Although Φ F and APAR × Φ F should theoretically help explain short-term variations in LUE P and GPP due to changes in non-photochemical quenching, such effects were not observed in the data we analyzed.
Such possibilities, as well as the prospect of using APAR × Φ F as the basis for consistently estimating APAR at the global scale (e.g., Ryu et al., 2018Ryu et al., , 2019, deserve continued attention in future research. Logically speaking, however, using APAR × Φ F for APAR estimation seems to be a somewhat circular approach as fPAR is used in the calculation of f esc that is needed to estimate APAR × Φ F (Table 1).
Approaches to circumvent this circularity, e.g. iterative methods starting from existing fPAR products, would therefore be needed. Also, given available fPAR data, the approach could be used to estimate PAR from satellite data.

Main sources of uncertainties and potential effects on our results
We are aware that despite our efforts in data collection and analysis, there are several limitations potentially affecting our results, which are discussed below.
It is a well-known problem that the radiometric footprint of the instruments used for SIF obs retrieval and NIR V observations is typically much smaller than the eddy covariance footprint, even for a hemispheric viewing geometry (Marcolla and Cescatti, 2017). Therefore, the remote measurements typically do not fully capture the effects of spatial heterogeneity in canopy-level fluxes (Gamon, 2015;Marcolla and Cescatti, 2017). While crop fields tend to be more homogeneous than natural ecosystems, this factor is still expected to affect our results to some degree although it is hard to quantify.
There are differences between observations in narrow angle nadir (wheat and corn) vs.
hemispheric (rice) viewing geometry that go beyond the size of the measurement footprint. In particular, these effects are related to the asymmetric angular patterns in canopy SIF emission, which are similar to the angular patterns in NIR reflectance . How these differences affect the relationships between SIF obs and GPP requires further investigation.
The retrieval of SIF obs introduces additional uncertainties (noise and bias), largely arising from the difficulty of retrieving the small SIF signal (~1% of the absorbed light) from the much stronger background of reflected sunlight (Damm et al., 2011;Frankenberg and Berry, 2018;Meroni et al., 2009). In particular, these uncertainties quite directly propagate into uncertainties in our estimates of Φ F , as Φ F was calculated as the 'residual' of SIF obs and APAR × f esc , assuming that APAR × f esc is much less affected by measurement uncertainties. The reason for this assumption is that only direct observations of RED and NIR reflectance and PAR were used to calculate APAR × f esc (Table 1). A similar reasoning as that for Φ F is valid for APAR × Φ F as the latter is the product of Φ F times APAR.
Apart from different choices of SIF obs retrieval algorithms in interaction with the spectral and radiometric characteristics of the instruments and processing (Table 3, Section 2.2), differences in uncertainty in SIF obs between datasets also arises from considerable differences in the number of spectra collected per minute. In particular, the wheat site typically collected about 60-120 times the number of spectra measured at the rice and corn sites (Table 3). Such a high acquisition frequency translates into stronger noise reduction in SIF retrievals via temporal averaging.
We think the uncertainties related to SIF retrieval do not affect our main conclusions in their essence. First, the strong evidence for Hstruc is not affected by this aspect at all, as APAR × f esc and f esc are not based on SIF obs retrievals but on NIR V (Table 1). Second, although the lack of evidence supporting Hphys appears to be partly caused by noise issues, even the wheat dataset, which apparently had the best SIF obs quality, did not show indications of better performance of physiological over structural information.
It is also worth noting that we used field observations of fPAR tot rather than fPAR green that only captures the light absorbed by chlorophyll. Previous studies have suggested to use the fraction of green LAI to total LAI to estimate fPAR green from fPAR tot (Gitelson et al., 2018;Gitelson and Gamon, 2015). Such data were only available for the rice dataset, however, preventing a consistent application of the fPAR green approach to all datasets. Nevertheless, we think that using in situ measured fPAR has advantages over using an indirect estimate based on canopy reflectance due to the uncertainties in the latter approach. Using fPAR tot rather than fPAR green is expected have an impact mainly during the senescence phase at the end of the growing season and logically should affect the estimates of f esc , LUE P , and APAR × Φ F (Table 1). However, our interpretation of the results in terms of the strong performance of APAR × f esc for GPP estimation is not affected as APAR × f esc was estimated as NIR V x PAR.

Implications for large scale GPP estimation
Ultimately, the main motivation to study SIF-GPP relationships is to improve the remote and largescale estimation of GPP . We found that the structural component of SIF (APAR × f esc ) showed strong linear relationships to GPP that were as good as or even better than those based on SIF obs However, our results strongly focused on crops, the site level and temporal variations.
Large scale applications using airborne or satellite SIF retrievals, however, tend to capture variations across space more than variations in time. Even when longer time series of SIF are considered, the number of pixels across space is typically much larger than the number of observation time steps. Therefore, we need to be cautious in making incorrect generalizations from the temporal case to the spatial case. Two different spatial cases should be distinguished. First, spatial variability within a given ecosystem type. This case is relevant for airborne remote sensing and could be tested in an approach that combines hyperspectral reflectance -based estimation of leaf and canopy parameters, SIF retrieval, and process-based modelling of GPP and APAR (Tagliabue et al., 2019).
Second, spatial variability within and between ecosystems, which is relevant for global-scale satellite applications. In particular, several global scale studies have shown differences in the slope of SIF obs -GPP relationships between ecosystems (Guanter et al., 2014;Sun et al., 2018;. Since the main origin of the slope differences for a given photosynthetic pathway is strongly suspected to be f esc (Zeng et al., 2019), APAR × f esc should show essentially the same tendencies as SIF obs . Therefore, using f esc to normalize the differences in slope between SIF obs and GPP in different ecosystems, is expected to improve spatial SIF-GPP relationships. Judging from our results, however, the improvement in spatial patterns goes together with a degradation of temporal relationships. In contrast to our site-level findings, Qiu et al. (2019) reported consistent and partly considerable improvements in the temporal SIF -GPP relationships when normalizing for f esc based on global process-based simulations. However, small to moderate improvements and partly even degradation of results were found when applying the f esc normalization to OCO-2 SIF retrievals and evaluating the temporal correlation to FLUXCOM GPP (Qiu et al., 2019).
For APAR × f esc , the approach of using f esc to correct for slope differences between SIF obs and GPP is 1) not possible as NIR V is already used to calculate APAR × f esc and 2) not desirable as f esc contains the LUE P -relevant information in addition to APAR. It seems, therefore, that for APAR × f esc -based GPP estimation, an ecosystem-dependent slope has to be applied as was partly done in Badgley et al. (2019). This is particularly relevant for evergreen needleleaf forests (ENF) that have a much lower NIR reflectance despite rather high GPP during the growing period, as can be inferred also from previous studies investigating SIF obs -GPP relationships (Sun et al., 2018;Yao Zhang et al., 2018a).
Based on our results, different slopes need to be applied for C3 and C4 species no matter if SIF obs or APAR × f esc is used (Fig. S12), which indicates that ecosystem-or pathway-dependent slopes will be required to some degree for global studies in any case. Further research is needed to determine if SIF obs considerably outperforms APAR × f esc for GPP estimation in ENF ecosystems where canopy structure is stable but photosynthesis and Φ F are downregulated in winter Porcar-Castell, 2011;Raczka et al., 2019).
A point of practical interest is that the product of NDVI times the upwelling NIR radiance (NIR V R) shows very strong correlation to NIR V P (Fig. D1 in the appendix) and is therefore an alternative that only needs data from the RED and NIR bands. The results of NIR V R and NIR V P showed similar levels of performance for GPP estimation (Fig. D2). As NIR V R does not require separate PAR data, it might have advantages over NIR V P for satellite-based applications at large scales where PAR data are not readily available despite recent progress .
Apart from the different impacts of temporal and spatial variation, there is another aspect to consider for satellite applications. In contrast to ground-based observations, non-geostationary satellites take snapshots at a given moment in time, not daily mean values. While the expected effects of this were previously simulated and directly examined in several studies Yao Zhang et al., 2018b), we found overall consistent patterns compared to using half-hourly values with only minor differences in terms of the performance ranking (Fig. S8). We therefore conclude that our findings are particularly relevant for satellite-based applications when considering the temporal (per-pixel) variability. Future SIF retrievals from geostationary missions such as TEMPO (Zoogman et al., 2017) or GeoCarb (Moore III et al., 2018) could provide the basis for large scale analyses similar to our site-level approach.

Conclusion
We mechanistically decomposed canopy-level SIF obs /APAR at three crop sites into the physiological emission yield, Φ F , and the canopy escape fraction, f esc , and examined their relationships to the photosynthetic light use efficiency, LUE P . We found that f esc had a considerable seasonal correlation to LUE P across the three sites, while Φ F was essentially uncorrelated and did not show any clear seasonal patterns. Even at the diurnal scale, Φ F did not outperform f esc in terms of correlation to LUE P except in the case of the wheat site. Consistent with these results, we found that the canopy structure and radiation component of SIF, APAR × f esc , was a better estimator of GPP than the total emission component of SIF, APAR × Φ F . In fact, APAR × f esc even outperformed SIF obs for GPP estimation for two of the three sites and had comparable performance to SIF obs for the remaining crop. APAR × Φ F , in contrast, improved the relationship to APAR compared to SIF obs , but had a considerably weaker performance for GPP estimation than SIF obs . Importantly, the canopy structure and radiation component of SIF can be observed on the basis of the near-infrared reflectance of vegetation NIR V with multispectral instruments without the need for SIF retrieval or explicit fPAR or even PAR information. Without independent PAR observations, the upwelling near-infrared radiance of vegetation can be used directly as estimate of APAR × f esc .
Our findings focus on the temporal relationships at the seasonal time scale and highlight the importance of the canopy structure effects on the SIF obs signal, as well as APAR × f esc as effective GPP proxy in crops. Apart from providing further evidence of the practical usefulness of NIR V in general and the NIR V -based escape fraction formula in particular, we presented a comprehensive framework of analyzing the separate contributions of Φ F and f esc that is expected to stimulate future research.
Appendix A. Temporal patterns of APAR, GPP, SIFobs, × f esc , and × Φ F Figure A1: Overview of the temporal patterns of the main variables related to GPP. Half-hourly data for the same data selection as in Fig. 4 are shown. Note that the y-axis scales are not the same for the different crop types in order to facilitate the visual comparison of temporal patterns. Units are as in Fig. 4.

Appendix C. Relationships between APAR and SIF
As in Fig. 4 only the R 2 results for the relationships between SIF obs , APAR × f esc and APAR × Φ F and GPP are shown, we include the detailed relationships between APAR on the one hand and SIF obs , APAR × f esc and APAR × Φ F on the other hand here (Fig. C1). Figure C1: Illustration of the effect of seasonal changes in f esc on the relationship between × f esc and APAR for the wheat dataset. As f esc is the ratio of × f esc and APAR, its seasonal changes (panel c) directly correspond to seasonal changes of the slope in the relationship.
It is instructive to study the patterns in the scatterplots in Fig. C1, as there are partly relatively small differences in R 2 between APAR × Φ F and APAR × f esc but clear patterns of seasonally varying slopes in the APAR × f esc -APAR relationships. This can be well illustrated with the example of the wheat dataset (Fig. C2). Assuming a zero intercept, the slope in the APAR × f esc -APAR relationship is simply the ratio of APAR /(APAR × f esc ), which is 1/f esc . As f esc is showing considerable seasonal variation, this is reflected also in the slope. An illustration of this situation is shown in Fig. C2. The reason why the APAR × Φ F -APAR relationship does not show such varying slope patterns is that Φ F did not show clear seasonal trends (Fig. 2od-f). Figure C2: Illustration of the effect of seasonal changes in f esc ) on the relationship between × f esc and APAR for the wheat dataset. As f esc is the ratio of × f esc and APAR, its seasonal changes (panel c) directly correspond to seasonal changes of the slope in the relationship (panel d).
Appendix D. Comparison of NIR V R and NIR V P for practical purposes NIR V P = NIR V × PAR requires PAR information that cannot be obtained directly from multispectral sensors designed for NDVI and NIR V observations. Although eddy covariance towers are typically equipped with separate quantum sensors for PAR measurements, some sites without eddy covariance instruments might not have such sensors. Therefore, using the radiance version NIR V R=NIR V × (downwelling nir Radiance) = NDVI × (upwelling nir radiance) could be useful in certain circumstances ("nir" was used to denote "near-infrared" here, as "NIR" was defined as nir reflectance, see Table 2). For example, when using airborne or satellite observations, the radiance data from the sensors in the atmospheric window around the O2A band does not require atmospheric correction and therefore has a clear practical advantage Zeng et al., 2019). For all three crop datasets, we found very strong correlations of NIR V R and NIR V P when using the full half-hourly time series (Fig. D1). The reason underlying the strong NIR V R and NIR V P correlations is the strong correlation of incoming PAR and downwelling near-infrared irradiance. While it is known from atmospheric radiative transfer theory that the relationship between the downwelling near-infrared radiance and incoming PAR is affected by the diffuse fraction of PAR, this effects seems negligible at the seasonal time scale. As can be expected from the very high linear correlation shown in Fig. D1, NIR V R and NIR V P have essentially the same performance for GPP estimation (Fig. D2).

Appendix E. Supplementary data
Supplementary data to this article can be found online at [insert DOI link]