Sun-Induced Chlorophyll Fluorescence III: Benchmarking Retrieval Methods and Sensor Characteristics for Proximal Sensing

Abstract

The interest of the scientific community on the remote observation of sun-induced chlorophyll fluorescence (SIF) has increased in the recent years. In this context, hyperspectral ground measurements play a crucial role in the calibration and validation of future satellite missions. For this reason, the European cooperation in science and technology (COST) Action ES1309 OPTIMISE has compiled three papers on instrument characterization, measurement setups and protocols, and retrieval methods (current paper). This study is divided in two sections; first, we evaluated the uncertainties in SIF retrieval methods (e.g., Fraunhofer line depth (FLD) approaches and spectral fitting method (SFM)) for a combination of off-the-shelf commercial spectrometers. Secondly, we evaluated how an erroneous implementation of the retrieval methods increases the uncertainty in the estimated SIF values. Results show that the SFM approach applied to high-resolution spectra provided the most reliable SIF retrieval with a relative error (RE) ≤6% and <5% for F₆₈₇ and F₇₆₀, respectively. Furthermore, although the SFM was the least affected by an inaccurate definition of the absorption spectral window (RE = 5%) and/or interpolation strategy (RE = 15–30%), we observed a sensitivity of the SIF retrieval for the simulated training data underlying the SFM model implementation.

1. Introduction

Plant photosynthesis is the primary process in terrestrial ecosystems. Accurate estimates of photosynthesis and its dynamics are pivotal to understand complex feedbacks and exchange interactions in the land–atmosphere system [1,2]. The assessment of photosynthesis also presents a key challenge of the remote sensing (RS) community. Sun-induced chlorophyll fluorescence (SIF) is considered the most direct RS signal to track photosynthetic activity and its dynamics at leaf, canopy, ecosystem, or even global scale [3,4]. Therefore, an accurate retrieval of SIF in crucial to understand photosynthesis and its dynamics.

SIF is emitted by chlorophyll-a (Chl_a), whose intensity depends on the incoming radiation and Chl_a concentration. The full chlorophyll fluorescence spectrum covers the wavelength range from 650 up to 800 nm, from the red (F_R) to the far-red (F_FR) range of the spectrum [5] (Figure 1). The weak SIF signal is convolved with the vegetation’s reflected radiance. This radiance is typically 10 orders of magnitude greater than SIF, making the decoupling of both a significant challenge. In order to retrieve SIF from spectroradiometric measurements, absorption features in the solar or Earth’s atmosphere can be exploited. In particular, the solar Fraunhofer lines—Fe (758.8 nm) and KI (770.1 nm) [6]—or the Earth’s two O₂ absorption features—O₂B (687 nm) and O₂A (760 nm) bands—are typically used due to their spectral proximity to the peaks of the chlorophyll SIF emission spectrum [7]. SIF retrieval using Fraunhofer lines requires high spectral resolution (SR) (e.g., << 0.1 nm), high radiometric resolution, and very high signal-to-noise ratio (SNR) [8]. Most proximal sensing instruments do not satisfy these requirements. Hence, this paper will focus on the retrieval of SIF at the oxygen absorption features, which do not require such a high-performance instrument.

Meroni et al. [9] proposed that the retrieval methods used to quantify SIF can be divided into two major categories: Reflectance-based and radiance-based approaches. Particularly, radiance-based retrieval schemes using on the Fraunhofer line depth (FLD) principle, e.g., standard FLD (sFLD [10]), three bands FLD (3FLD [11]), improved FLD (iFLD [12]), and the spectral fitting method (SFM [13]) are the most frequently used methods for the proximal sensing of SIF (i.e., based on a bibliographic survey in Web of Science—core collection—during the last year 21 papers were published using the keyword sFLD, 22 using 3FLD, 22 using iFLD, and 23 using SFM) (FLD methods and SFM are described more in Section 2).

Many recent publications focus on the retrieval and interpretation of SIF from satellite platforms like GOSAT, GOME-2, and OCO-2 [6,14,15]. The first approach proposed for GOSAT took advantage of its very high spectral resolution (0.025 nm) to evaluate the in-filling by SIF of single solar Fraunhofer lines [6,16,17]. The same approach is now being used for OCO-2 [15,18]. This type of retrieval was later extended to wider fitting windows also including Earth’s atmospheric absorption features in the case of GOME-2, SCIAMACHY, and more recently, TROPOMI, which have a coarser spectral resolution (0.5 nm) than GOSAT. Data-driven methods are used to empirically model atmospheric radiative transfer effects (e.g., References [19,20,21]). On the other hand, with the recent development of new in situ-based SIF measuring systems (e.g., FloX, JB Hyperspectral, Dusseldorf, Germany, PhotoSpec [8], Piccolo [22,23], and HyScreen [24]), the number of publications investigating proximal sensing approaches of SIF steadily increases [7,25,26,27,28,29] and will possibly grow in the near future. Moreover, hyperspectral ground measurements play a crucial role in the calibration and validation of recent (e.g., Sentinel-2) and future (e.g., FLEX) satellite missions and can be used to upscale ground signals to the satellite.

Since SIF only represents a fraction of the radiance measured by the spectrometer, insufficient characterization of the sensor, inadequate measurement protocols, or an incorrect implementation of the retrieval method will lead to an erroneous estimation and interpretation of the SIF signal. For this reason, many researchers have worked within the framework of a European cooperation in science and technology (COST) action, “Innovative optical tools for proximal sensing of ecophysiological processes” (OPTIMISE, ES1309; https://optimise.dcs.aber.ac.uk/) over the last four years to address these challenges. The OPTIMISE community has compiled three papers: One on instrument characterization [30], one on measurement setups and protocols [31], and this paper on retrieval methods to make gathered information available. This compilation of papers aims to summarize the state-of-the-art for high quality SIF measurements and bring up-to-date the review papers of Meroni et al. [9] and Damm et al. [32].

This study aims to provide a guiding document to evaluate possible SIF retrieval accuracies for F_R and F_FR considering the characteristics of novel instrumentation (e.g., SR, spectral sampling interval (SSI), signal-to-noise ratio (SNR)) in combination with frequently applied retrieval schemes (e.g., sFLD, 3FLD, iFLD, and SFM). Further, this document aims to outline critical steps in the implementation of the retrieval schemes. The insights provided need to be considered in the total uncertainty budget of SIF retrieval. However, it is important to note that other factors such as canopy structure and atmospheric effects also need to be taken into account when estimating SIF uncertainties (cf. Reference [31]).

The paper is structured as follows: In Section 2 we give an overview of the most commonly used retrieval methods and outline their advantages and disadvantages. Section 3 describes sensitivity analysis evaluating how the retrieval accuracy depends on the retrieval methods and the dataset used to evaluate them. Section 4 presents and discusses the results of this work and, finally, the main findings are summarized in Section 5.

2. Background on Frequently Used SIF Retrieval Methods for Proximal Sensing

Vegetation fluorescence under solar illumination conditions adds a weak signal to the reflectance solar radiation. In the red and far-red spectral region, the up-welling radiance (L^↑) is the result of two contributions: The reflected radiance $\left(L_R=R·E^{\downarrow }/\pi \right)$ from the vegetation and fluorescence,

$$L^{\uparrow }\left(\lambda \right)=\frac{R\left(\lambda \right)·E^{\downarrow }\left(\lambda \right)}{\pi }+F\left(\lambda \right)$$

(1)

where λ is the wavelength, R is the actual reflectance, E^↓ is the down-welling irradiance incident to the surface, and F the top of the canopy SIF radiance in the direction of observation. According to the international system of units, spectral $L^{\uparrow }$ is defined as the radiant flux emitted, reflected, transmitted, or received by a surface, per unit projected area per unit solid angle per wavelength (e.g., mW m⁻²sr⁻¹nm⁻¹). $L^{\uparrow }$ is a directional quantity. On the other hand, spectral E^↓ is defined as the radiant flux received by a surface per unit area per wavelength (e.g., mW m⁻²nm⁻¹). E^↓ is not a directional quantity. Therefore, to compute the reflected radiance, E^↓ needs to be divided by $\pi$. Finally, the reflectance obtained from L^↑ and E^↓ measurements and containing the contribution of F is defined as apparent reflectance (R_app) [9] (see Figure 2),

$$R_{app}\left(\lambda \right)=\frac{\pi ·L^{\uparrow }\left(\lambda \right)}{E^{\downarrow }\left(\lambda \right)}=R\left(\lambda \right)+\frac{\pi ·F\left(\lambda \right)}{E^{\downarrow }\left(\lambda \right)}.$$

(2)

Most retrieval algorithms used to estimate SIF from ground measurements are based on the FLD principle. Conceptually, FLD approaches exploit the different relative contribution of F to the L^↑ and the E^↓ spectra inside and outside of an absorption feature (Figure 3).

FLD approaches exploit the ratio between these up- and down-welling fluxes in the form of the apparent reflectance factor (R_app). The method aims to separate F from L^↑ estimating, therefore, the actual reflectance (R), the ratio between reflected radiance and E^↓, assuming a smooth R profile in the absorption region. Importantly, an inaccurate estimation R has a high impact on the retrieval of F [12]. Therefore, a precise and accurate function to interpolate R inside the absorption line is needed. This function should take into account the following factors: (1) The appropriate definition of both shoulders of the used absorption feature, (2) the selection of measurements outside the used absorption feature as end points for the interpolation, and (3) the selection of an appropriate interpolation method. The choice of these three factors differentiates each FLD-based retrieval scheme. In contrast, the SFM method decouples F and R using spectral curve fitting. The method relies on parameterized mathematical functions representing R and F within narrow spectral windows centered at the O₂ absorption features.

2.1. FLD-Based SIF Retrieval Methods

Retrieval methods based on the FLD principle take advantage of the reduced E^↓ in the O₂B and O₂A absorption features reaching the surface, which increases the relative contribution of F to L^↑ (Figure 3). Such retrieval allows estimating F in a specific spectral band, 687 nm (F₆₈₇) at O₂B and 760 nm (F₇₆₀) at O₂A, which correspond to the smallest E^↓ and L^↑ value in the oxygen absorption windows (Figure 1). While the O₂B absorption window is close to the maximum peak of fluorescence at 685 nm (maxF₆₈₅), the O₂A window is shifted towards longer wavelengths. Hence, F₇₆₀ represents only around 70% of the radiated F signal at the F_FR peak (maxF_740, Figure 1).

Since the seminal review provided by Meroni et al. [9], various new SIF retrieval methods were proposed by the scientific community, for instance, the singular vector decomposition (SVD) [17,33], the nFLD approach [34], and the peak height method [35]. However, as presented in the introduction, sFLD, 3FLD, and iFLD are still the most frequently applied methods for proximal sensing of SIF. Therefore, only these approaches are included for further discussion and their main assumptions, strengths, and shortcomings are described.

2.1.1. sFLD

The Fraunhofer line depth (FLD) principle [10] relies on measurements of E^↓ and L^↑ inside and outside an absorption window (λ_in and λ_out) such as solar Fraunhofer lines or Earth’s atmosphere O₂ absorption features.

$$\begin{gathered} L^{\uparrow }\left({\lambda }_{in}\right)=\frac{R\left({\lambda }_{in}\right)·E^{\downarrow }\left({\lambda }_{in}\right)}{\pi }+F\left({\lambda }_{in}\right) \\ {L}^{\uparrow }\left({\lambda }_{out}\right)=\frac{R\left({\lambda }_{out}\right)·E^{\downarrow }\left({\lambda }_{out}\right)}{\pi }+F\left({\lambda }_{out}\right) \end{gathered}$$

(3)

Fluorescence is calculated from a linear system of equations by comparing the measured signal inside and outside an absorption feature,

$$\begin{gathered} R=\frac{L^{\uparrow }\left({\lambda }_{out}\right)-L^{\uparrow }\left({\lambda }_{in}\right)}{E^{\downarrow }\left({\lambda }_{out}\right)-E^{\downarrow }\left({\lambda }_{in}\right)}·\pi \\ F=\frac{E^{\downarrow }\left({\lambda }_{out}\right)·L^{\uparrow }\left({\lambda }_{in}\right)-L^{\uparrow }\left({\lambda }_{out}\right)·E^{\downarrow }\left({\lambda }_{in}\right)}{E^{\downarrow }\left({\lambda }_{out}\right)-E^{\downarrow }\left({\lambda }_{in}\right)} \end{gathered}$$

(4)

A fundamental assumption is that F and R remain constant over the exploited absorption feature. The original FLD (sFLD) approach is relatively simple and requires the measurement of E^↓ and L^↑ inside the absorption window and only one measurement at the shoulder of the absorption feature. The main limitation of this approach is the frequent violation of the base assumption—that F and R are spectrally constant over the used absorption feature [12,36,37].

2.1.2. 3FLD

In order to overcome the limitations of the sFLD approach, Maier et al. [11] proposed the use of three bands to solve Equation (4). The single reference band (λ_out) of sFLD is replaced by a linear interpolation of two bands on the left and right shoulders of the absorption feature. The advantage of this modification is that it considers a linear variation of R and F over the absorption window, which is a good approximation when the red-edge shoulder is displaced towards short wavelengths. Even though this approach is theoretically more advanced than the sFLD, non-linear variations of R and F often result in inaccurate estimates in particular for the O₂B band.

2.1.3. iFLD

The iFLD method uses two correction factors (α_R and α_F) to account for the non-linear variation of R and F inside and outside the absorption feature and to account for using the apparent reflectance instead of actual reflectance in the measurement. The method proposed by Alonso et al. [12] makes use of the full E^↓ and L^↑ spectral information in the region around the absorption to estimate the α_R factor. F is expressed as:

$$F\left({\lambda }_{in}\right)=\frac{{\alpha }_R·E^{\downarrow }\left({\lambda }_{out}\right)·L^{\uparrow }\left({\lambda }_{in}\right)-L^{\uparrow }\left({\lambda }_{out}\right)·E^{\downarrow }\left({\lambda }_{in}\right)}{{\alpha }_R·E^{\downarrow }\left({\lambda }_{out}\right)-{\alpha }_F·E^{\downarrow }\left({\lambda }_{in}\right)}$$

(5)

In order to calculate the value of α_R and α_F correction factors, it is necessary to know, in advance, F and the actual reflectance. Since this is not possible due to the variable nature of both parameters, Alonso et al. [12] proposed the use of R_app to calculate the correction factors for reflectance (${\widehat{\mathsf{\alpha }}}_{\mathrm{R}}$) and F (${\widehat{\mathsf{\alpha }}}_{\mathrm{F}}$):

$${\widehat{\alpha }}_R=\frac{R_{app}\left({\lambda }_{out}\right)}{{\tilde{R}}_{app}\left({\lambda }_{in}\right)}$$

(6)

$${\widehat{\alpha }}_F\approx \frac{E^{\downarrow }\left({\lambda }_{out}\right)}{{\tilde{E}}^{\downarrow }\left({\lambda }_{in}\right)}·{\widehat{\alpha }}_R$$

(7)

where $R_{\mathrm{app}}\left({\mathsf{\lambda }}_{\mathrm{out}}\right)$ is the apparent reflectance measured outside the absorption feature, and ${\tilde{R}}_{\mathrm{app}}\left({\mathsf{\lambda }}_{\mathrm{in}}\right)$ is the apparent reflectance inside the absorption feature λ_in. ${\tilde{R}}_{\mathrm{app}}\left({\mathsf{\lambda }}_{\mathrm{in}}\right)$ is obtained from the non-linear interpolation of the apparent reflectance using the continuous reflectance spectrum at the left and right shoulders. Analogously, ${\tilde{E}}^{\downarrow }\left({\mathsf{\lambda }}_{\mathrm{in}}\right)$ is obtained by interpolation of the irradiance in order to predict $E^{\downarrow }$ unaffected by atmospheric absorption. Replacing α_R and α_F by the apparent coefficients (${\widehat{\alpha }}_R,{\widehat{\alpha }}_F$), F can be expressed as:

$$F^{*}\left({\lambda }_{in}\right)=\frac{{\widehat{\alpha }}_R·E^{\downarrow }\left({\lambda }_{out}\right)·L^{\uparrow }\left({\lambda }_{in}\right)-L^{\uparrow }\left({\lambda }_{out}\right)·E^{\downarrow }\left({\lambda }_{in}\right)}{{\widehat{\alpha }}_R·E^{\downarrow }\left({\lambda }_{out}\right)-{\widehat{\alpha }}_F·E^{\downarrow }\left({\lambda }_{in}\right)}.$$

(8)

2.1.4. Implementation of the FLD-Based Methods

In summary, these FLD-based retrieval schemes can be distinguished by the approximation of R and F over the spectral extent of the absorption feature used, while the sFLD requires a single measurement of E^↓ and L^↑ inside and outside the absorption feature. The 3FLD requires one measurement of E^↓ and L^↑ inside the absorption feature and two outside, at the left and right shoulders. The iFLD makes use of the complete E^↓ and L^↑ spectral information to compute ${\widehat{\alpha }}_R$ and ${\widehat{\alpha }}_F$ correction factors. Importantly, E^↓ and L^↑ must be defined in the same units. In the case of measuring E^↓ and L^↑ pointing a fiber optics towards a Spectralon ^® panel (Labsphere, Inc., North Sutton, NH, USA) and vegetation target, respectively, both measurements are directional, i.e., are dependent of the observation direction. Therefore, both are defined as the radiant flux reflected by a surface, per unit projected area per unit solid angle per wavelength (e.g., mW m⁻²sr¹nm¹). However, if E^↓ is measured with a fiber optic attached to a cosine diffuser pointing the sky, E^↓ is not directional, i.e., E^↓ is integrated over the entire hemisphere. It is defined as the radiant flux received by a surface per unit area per wavelength (e.g., mW m⁻²nm⁻¹). However, provided that the cosine receptor responds according to Lambert’s cosine law, the total luminous flux (global sky and sun irradiance) $F{\downarrow }_{tot}=\pi \times E{\downarrow }_{max}$, being E^↓_max the peak luminous intensity at nadir incidence of the global irradiance. It is possible to divide E^↓ by $\pi$ to define E^↓ and L^↑ with the same units (e.g., mW m²sr¹nm¹).

The first, and a critical, step is to define the wavelength range covering the full O₂B and O₂A absorption windows. This is particularly relevant when the spectral resolution is sub-optimal and spectral bands not affected by absorption are not easily identifiable. In our study, we determined that the spectral ranges affected by O₂ absorption are: 686–697 nm (O₂B) and 759–770 nm (O₂A) (see Appendix C for further discussion). We consider that outside these ranges, the influence of O₂ absorption is minimal. Once the absorption region is identified, the spectral band representing maximum absorption can be easily identified by searching for the smallest E^↓ value in the respective wavelength range (Appendix A–sFLD, 3FLD, and iFLD red dot). It is recommended not to use a fixed wavelength position to define the band representing the minimum absorption. If the spectroradiometer used to measure is not accurately spectrally calibrated or there is a spectral shift due to variable measurement temperature (cf. Reference [24]), the use of a fixed wavelength position will introduce an error in the F retrieved. In addition, in case of using a spectroradiometer with low SNR, it is recommendable to define $L^{\uparrow }\left({\lambda }_{in}\right)$ and $E^{\downarrow }\left({\lambda }_{in}\right)$ as the average value of $L^{\uparrow }\left({\lambda }_{in}\right)$ and $E^{\downarrow }\left({\lambda }_{in}\right)$ of adjacent wavelengths.

The next step is to determine the spectral lines defining the left and right shoulders of the absorption feature. This can be achieved by searching for E^↓ local maxima at both shoulders. According to the method assumptions, one or several wavelengths have to be used to compute E^↓ and L^↑ outside the absorption feature. Again, it is not recommended to use a fix wavelength to define the left and right shoulders position and by computing the average values of $L^{\uparrow }\left({\lambda }_{out}\right)$ and $E^{\downarrow }\left({\lambda }_{out}\right)$ using adjacent wavelengths it is possible to increase the SNR.

The interpolation method, and its implementation, is a key step in the most advanced retrieval algorithms. The main purpose of the interpolation is to eliminate the effect of the absorption feature (Earth or solar atmosphere) by constructing the L^↑, E^↓, and/or R_app spectra shape between the absorption features shoulders (Figure 4, green line). Through the interpolation, a reference band inside the absorption feature can be computed, mitigating the effect of having a different reflectance inside and outside of the absorption feature, which has been the weak point in the earlier retrieval algorithms, such as sFLD. In the case of the E^↓, interpolation is easier, since the solar spectrum is very well characterized. The reflected radiance is more complex, since it is modulated by the reflectance and the fluorescence spectral shapes, which are complex, especially in the red region around the O₂B absorption. Therefore, it requires more flexible interpolating functions capable to fit its curvature (see Figure 4B,D).

For the interpolation, it is necessary to determine the known data points from which construct the new data points (Figure 4, green dots). Properly selecting those points is the most critical step for an accurate interpolation. Similar to previous steps, this can be achieved by searching for E^↓ local maxima at both shoulders of the absorption feature (green dots and lines in Figure 4A,C). Then, those bands can be used to interpolate L^↑ and/or R_app. Importantly, the small Earth’s atmospheric absorptions (e.g., water) present in the left and right shoulders of the oxygen absorption feature must be avoided. Otherwise, if the points inside those absorptions are included as key points for the interpolation, the resulting curve will be shifted downwards, which can introduce an error of a similar or larger magnitude than the fluorescence itself (blue dots and lines in Figure 4A,C). This problem is great for spectrometers of finer spectral resolution, since they resolve deeper fine features. Instruments with coarser spectral resolution are less sensitive, since they cannot resolve this level of detail.

The presence of noise also needs to be considered as it distorts the spectrum, and this might, in turn, distort the interpolation. Therefore, the use of spectrometers with very high SNR is recommended. R_app is particularly affected by noise, since noise will be amplified by dividing L^↑ by E^↓ (Figure 4B,D). In this case, it is possible to mitigate the effect of noise by two approaches: One would be to avoid using spline or any other method that obliges the resulting curve to pass precisely through the key points used for interpolating; the other is to use noise-reduction algorithms to mitigate noise, pixel binning being the simplest one, or by applying a low-pass filter. However, care should be exercised to ensure that spectral features are not smoothed out when denoising the spectrum.

Appendix C provides additional detail on wavelength setting, spectral window intervals, and interpolation method used in this study for the sFLD, 3FLD, and iFLD retrieval methods. We would like to highlight that these settings worked for this specific dataset (cf. Section 3) and should not be seen as a standard setting. It is important to adapt the herein-described setting to the characteristics of the used instruments and observed vegetation [38].

2.2. Spectral Fitting Method-Based SIF Retrieval

2.2.1. Background of the Spectral Fitting Method

The SFM is a technique to decouple SIF and reflectance from high spectral resolution radiance observations. The method relies on general mathematical functions representing canopy R and F within narrow spectral windows centered at oxygen absorption features. As previously defined, in the case of ground-based spectrometer, radiance detected by sensors includes contributions from canopy L_R and emitted F. The parameters of the functions employed to represent F and R are thus optimized by non-linear least square optimization process, comparing instrument observations with radiance computed accordingly with Equation (1). In this way, the a-priori F and reflectance functions can be spectrally decoupled on the base of their contribution at the different spectral lines.

The use of Spectral Fitting for proximal sensing measurements is easier compared with airborne or satellite measurement because the surface irradiance is directly measured, therefore Equation (1) can be solved without additional information (i.e., indirectly estimating top of canopy irradiance). In the past, several functions were proposed to fit the SIF and reflectance spectral behavior within narrow spectral windows as suggested in References [9,13,37,39]. However, SIF is generally represented as peak-like functions such as Gaussian, Lorenzian, or Voigt; whereas canopy reflectance is usually described as second- or third-order polynomial or more complex piecewise cubic splines.

The major advantage of SFM is the exploitation of a large number of spectral bands because the mathematical system, which describes L^↑ at different wavelengths, is better conditioned and, consequently, estimations are less affected by instrument noise. The SFM method is therefore particularly useful for spectroradiometers with very high spectral resolution that suffer from higher noise levels. On the other hand, this approach has a number of limitations that are due to its complex setup and the overall computational time. Recently, the SFM approach has been further extended to a unique and broad spectral window that covers the entire SIF emission region (650−800 nm) to estimate the full SIF emission spectrum consistently [13].

2.2.2. SFM Retrieval Method Implementation

As previously described, this method relies on general mathematical functions representing the R and F spectrum within a narrow spectral window. Before going into the description of the SFM implementation, the mathematical functions used to model R and F are defined. Since R (not distorted by F) and F spectral emission are not available from field measurements, first, a set of representative R and F spectra are modeled using a canopy reflectance-SIF radiative transfer model. Then, the accuracy of different functions used to represent R and F reference spectra are tested using the reference dataset. For instance, Cogliati et al. [13] used the FluoSpec radiative transfer model (RTM) to build a reference dataset simulating the spectral characteristics of the FLEX spectrometer FLORIS. In their study, different functions were tested to model R (i.e., linear, quadratic, cubic spline) and F (i.e., linear, quadratic, cubic, Gaussian, Lorentzian, and Voigt). From this analysis, they concluded that a cubic spline function and a Voigt function are the best functions to model R and F, respectively. On the other hand, Julitta et al. [40], also using FluoSpec to build a reference dataset but simulating the spectral response of Ocean Optic’s QE Pro spectrometer, suggest to use a cubic spline function to model R and a Gaussian function to model F.

In the study presented here, the Gaussian function performed better than the Voigt function in modeling F (and R² = 0.96 and R² < 0.5 for the Gaussian and Voigt functions, respectively). Cogliati et al. [13] propose that different sensor spectral responses (e.g., SR, SSI, and SNR) may lead to different functions to best represent R and F. This may explain why, for proximal sensing spectrometers, a Gaussian function was the most suitable algorithm to model F [40] instead of a Voigt function as suggested by Cogliati et al. [13] for FLORIS, FLEX space-borne sensor.

Therefore, in this study, R was approximated using a cubic spline function and a Gaussian function was used to model F,

$$F=a·e^{-\frac{{\left(\left(\lambda -{\lambda }_0\right)-\left(c-{\lambda }_0\right)\right)}^2}{2·b^2}}$$

(9)

where $a$ stands for the heights of the red (O₂B) or far-red (O₂A) F curve’s peak,$\lambda$ represents the wavelengths of the defined F and R fitting window interval, ${\lambda }_0$ is the first wavelength of the defined F and R fitting window interval, $c$ is the central wavelength of F peaks, and $b$ controls the width of the red and far-red F spectrum (Appendix C).

To retrieve F with the SFM, high-resolution E^↓ and L^↑ spectral measurements are needed. Similar to FLD retrieval methods, E^↓ and L^↑ must be defined with the same units (cf. Section 2.1.4). For this study, the F and R fitting spectral windows are defined between 750–780 nm (O₂A) and 680–698 nm (O₂B) (Appendix C). The first step to retrieve SIF using the SFM is to compute R_app and interpolate within the defined absorption range (O₂A: 759–770 nm and O₂B: 686–697 nm) (Appendix C). The output parameters of the interpolated R_app are used in the optimization cost function as the R ‘first guess’ parameters. Secondly, by using the iFLD method (Section 2.1), F at each O₂ absorption feature is estimated. The retrieved F is used in the optimization cost function as a first-guess parameter for the height of the red (O₂B) and far-red (O₂A) F peak curve (a parameter of the Gaussian function used in this study). Third, the first guess as well as the lower and upper boundaries for the R and F functions need to be defined. Finally, the parameters of the functions employed to model F and R are optimized by non-linear least square optimization process, comparing instrument observations with radiance computed by Equation (1). In order to estimate the accuracy of the retrieval method, it is recommended to compare the modeled and reference L^↑. For a good implementation of the SFM, the differences between modeled and measured L^↑ must be minimal.

In Appendix C, a diagram that sets out the most relevant steps when implementing the SFM retrieval method is presented. Appendix C then summarizes the input parameters, wavelength settings, wavelength intervals, and R and F model functions used in this study. Here, again, we would like to highlight that these settings specific for this dataset (cf. Section 3). These should not be considered as standard settings.

3. Assessment of SIF Retrieval Uncertainties-Sensitivity Analysis

The aim of this analysis is to assess possible SIF retrieval uncertainties for the combination of consolidated SIF retrieval approaches presented using off-the-shelf spectrometers currently used for proximal sensing of SIF. A summary of considered sensors and their characteristics is presented in

Table 1. Characteristics of the typical off-the-shelf spectrometers for the measurement of sun-induced fluorescence: Spectral range, spectral sampling interval (SSI), optical resolution measured as full width at half maximum (FWHM), and signal-to-noise ratio (SNR). * OceanOptics spectroradiometers can be built with these or other characteristics by demand.

Spectrometer	Range (nm)	SSI (nm)	FWHM (nm)	SNR
ASD FieldSpec	350−1000	1.4	3	4000
OceanOptics MAYA *	650−803	0.08	0.44	450:1
OceanOptics HR4000 *	650−840	0.05	0.28	300:1
OceanOptics QEPRO *	651−803	0.13	0.38	1100:1

and Figure 5.

For the data simulation, we employed a combination of the RMT FluorSAIL3 and MODTRAN-5 [41] to simulate radiance (L_TOC) and top-of-canopy (TOC) irradiance (E_TOC), respectively. FluorSAIL3, an extended version of the physically based FluorSAIL model [42], provides simulations of surface reflectance and emitted SIF with a spectral resolution of 1 nm. The coupling with the atmosphere was realized using the four-stream theory [43], while atmospheric functions were calculated with MODTRAN-5 [41].

Table 2. Parameters and variables used to calculate the reference radiance signatures for incoming and reflected light using FluorSAIL3.

Parameter	Unit	Values
FluorMODleaf
Internal structure parameter N	-	1.5
Chlorophyll ab	µg cm−2	20, 80
Leaf water	cm	0.025
Dry matter	g cm−2	0.01
Fluorescence efficiency factor		0.02, 0.04
Temperature	°C	20
FluorSAIL
LAI	-	1.4
LIDF	-	erectophile, planophile
Hot-spot parameter	-	0.1
MODTRAN 5
Correlated-K option	-	yes
DISORT number of streams	-	8
Molecular band model resolution	cm−1	0.1
Atmospheric profile	-	midlatitude summer
Aerosol model	-	rural
Visibility	km	23
Surface height	m	200
Water vapor	g cm−2	2.65
CO2	ppm	385
Solar zenith angle	deg	30
Viewing zenith angle	deg	0

and

Table 3. Structural and physiological canopy properties of the 16 vegetation canopies investigated. Varied parameters are the leaf angle distribution (LAD) (erectophile (erec.) or planophile (Plan.), the leaf area index (LAI, m² m⁻²), the chlorophyll content (µg cm⁻²) and the Fs efficiency factor.

Chlorophyll Content	Fs-Efficiency Factor	LAI	LAD	Canopy Type (Example)
20	0.02	1	Plan.	Sparse young unstressed wheat
80	0.02	1	Plan.	Sparse old unstressed wheat
20	0.04	1	Plan.	Sparse young stressed wheat
80	0.04	1	Plan.	Sparse old stressed wheat
20	0.02	4	Plan.	Dense senescent unstressed wheat
80	0.02	4	Plan.	Dense mid old unstressed wheat
20	0.04	4	Plan.	Dense senescent stressed wheat
80	0.04	4	Plan.	Dense mid old stressed wheat
20	0.02	1	Erec.	Sparse young unstressed bean
80	0.02	1	Erec.	Sparse old unstressed bean
20	0.04	1	Erec.	Sparse young stressed bean
80	0.04	1	Erec.	Sparse old stressed bean
20	0.02	4	Erec.	Dense senescent unstressed bean
80	0.02	4	Erec.	Dense mid old unstressed bean
20	0.04	4	Erec.	Dense senescent stressed bean
80	0.04	4	Erec.	Dense mid old stressed bean

list the simulation parameters and variables used in this study representing 16 typical vegetation targets, observed under a cloudless sky. We deliberately considered vegetation types as defined in Damm et al. [32] to relate our results with previous investigations and allow progress in the field of SIF proximal sensing to be tracked. The advantage of using simulated data is that SIF emission and surface reflectance are precisely known (Figure 6). E_TOC^↓ and L_TOC^↑ were simulated assuming Lambertian surface reflectance and neglecting adjacency effects following

$${L_{\mathrm{TOC}}}^{\uparrow }=\frac{<E_s^{\mathrm{o}}>cos{\theta }_s}{\pi }\lceil \frac{\left(\langle {\tau }_{\mathrm{ss}}\rangle +\langle {\tau }_{\mathrm{sd}}\rangle \right)·R}{1-R\langle {\rho }_{\mathrm{dd}}\rangle }\rceil +F$$

(10)

$${E_{\mathrm{TOC}}}^{\downarrow }=\frac{<E_{\mathrm{s}}^{\mathrm{o}}>cos{\theta }_{si}}{\pi }\lceil \frac{\left(\langle {\tau }_{\mathrm{ss}}\rangle +\langle {\tau }_{\mathrm{sd}}\rangle \right)}{1-R\langle {\rho }_{\mathrm{dd}}\rangle }\rceil ,$$

(11)

where $E_s^o$ represents the combined direct and diffuse top-of-atmosphere irradiance flux for a given sun-incident angle ${\theta }_{\mathrm{s}}$. ${\tau }_{\mathrm{ss}}$ is the atmospheric transmittance for direct downwelling irradiance, ${\tau }_{\mathrm{sd}}$ is the atmospheric transmittance for diffuse downwelling irradiance, and ${\rho }_{dd}$ is the spherical albedo. Angle brackets indicate functions that need to be convolved individually to avoid artefacts stemming from this processing step.

To evaluate the impact of sensor parameters on the SIF retrieval accuracy, simulated E_TOC^↓ and L_TOC^↑ were spectrally resampled using a two-step approach that includes the convolution of the input signal using a Gaussian function and the sampling of the convolved signal considering predefined sensor characteristics, following the approach described in Damm et al. [32].

We then added random Gaussian distributed noise with a mean value of zero and a standard deviation equal to the ratio between signal level (e.g., L_TOC^↑) and SNR to E_TOC^↓ and L_TOC^↑ as:

$$L_{\mathrm{TOC}}^{\mathrm{noise},\uparrow }={L_{\mathrm{TOC}}}^{\uparrow }+\mathcal{N}\left(0,\frac{{L_{\mathrm{TOC}}}^{\uparrow }}{SNR}\right)$$

(12)

$$E_{\mathrm{TOC}}^{\mathrm{noise},\downarrow }={E_{\mathrm{TOC}}}^{\downarrow }+\mathcal{N}\left(0,\frac{{E_{\mathrm{TOC}}}^{\downarrow }}{SNR}\right)$$

(13)

where $\mathcal{N}$ is a function providing a random value according to a normal distribution (μ,σ) defined by its mean (μ) and standard deviation (σ). For the SNR, we identified representative values per sensor from literature and assumed the noise contribution constant across wavelengths (cf. Table 2).

Finally, retrieved SIF signals were compared to the modeled reference SIF at the wavelength representing the maximum absorption (e.g., lowest L_TOC^↑). We then calculated several quality statistics including, the total relative error (RE), the coefficient of determination (R²), and the root mean square error (RMSE).

4. Results and Discussion

4.1. Impact of Sensor Specification on SIF Retrieval Methods

The following Section describes the impact of the different sensor specifications (i.e., SR, SSI, and SNR) on the precision, accuracy, and relative error of four different retrieval methods: sFLD, 3FLD, iFLD, and SFM. First, we compare the reference (modeled) and retrieved F₆₈₇ and F₇₆₀ of the 16 simulated canopy cases without noise (Figure 7 and

Table 4. Overview of the total relative error (RE), the coefficient of determination (R²), and the root mean square error (RMSE) of reference (modeled) and estimated F₇₆₀ and F₆₈₇ for 16 simulated canopy signatures and four different sensors (ASD, QE Pro, MAYA, and HR4000). Bold numbers represent the best performing retrieval method.

	F760				F687
	RE
	sFLD	3FLD	iFLD	SFM	sFLD	3FLD	iFLD	SFM
ASD	234.5	31.9	11.8	---	370.1	101.9	41.2	---
QE Pro	26.2	7.7	4.7	4.5	56.0	50.8	13.8	6.2
MAYA	33.4	14.4	7.0	4.9	62.9	50.8	10.4	7.2
HR4000	40.9	25.4	9.6	4.8	66.5	34.5	9.7	5.9
	R2
ASD	0.44	0.91	0.91	---	0.15	0.00	0.26	---
QE Pro	0.93	0.98	0.98	0.98	0.52	0.22	0.88	0.90
MAYA	0.93	0.98	0.98	0.98	0.48	0.13	0.88	0.90
HR4000	0.94	0.98	0.98	0.99	0.54	0.40	0.90	0.91
	RMSE
ASD	2.48	0.53	0.18	---	8.80	4.99	0.84	---
QE Pro	0.31	0.13	0.09	0.09	1.39	1.58	0.36	0.37
MAYA	0.32	0.13	0.08	0.08	1.55	2.04	0.37	0.38
HR4000	0.30	0.12	0.08	0.08	1.31	1.06	0.34	0.36

) to assess the effect of SR and SSI. Afterwards, the relative impact of SNR on F₆₈₇ and F₇₆₀ retrieval is quantified by adding the corresponding noise (Figure 8 and Table 4). The coarse SR and SSI of the ASD only allows assessing the performance of sFLD, 3FLD and iFLD method, while for sensors with higher SR and SSI (e.g., MAYA, QE Pro and HR4000), all four approaches were tested.

Considering retrieval in the O₂A band using high-resolution instruments, we obtain nearly identical results for F₇₆₀ estimates using 3FLD, iFLD, and SFM methods where the regression lines of F₇₆₀ (Figure 7E–H) and R (Figure 9E–H) are close to the 1:1 line (R² ≈ 0.98 and RMSE ≈ 0.08 mW m⁻² sr⁻¹ nm⁻¹). In more detail, the SFM approach has a tendency to underestimate R, while 3FLD, iFLD, and SFM overestimate F₇₆₀ for higher F values. On the other hand, sFLD and 3FLD applied to low-resolution sensors constantly overestimated F₇₆₀. However, when using the iFLD, an R² = 0.91 and RMSE = 0.18 mW m⁻² sr⁻¹ nm⁻¹ was found. The sFLD performs worst in estimating F₇₆₀ with a significant overestimation for all sensors (R² = 0.93–0.94 and RMSE = 0.3–0.32 mW m⁻² sr⁻¹ nm⁻¹). As described by Damm et al. [32], the strong overestimation of F₇₆₀ by the sFLD method is caused by a violation of the assumption that R and F are spectrally constant over the O₂A band. The 3FLD and iFLD methods assume a linear and non-linear relationship respectively between R and SIF, which results in a better prediction of R and F₇₆₀.

The retrieval of F₆₈₇ at the O₂B band is more challenging and was found to be unfeasible with low-resolution instruments such as the ASD (Figure 7A–D). Higher resolution sensors MAYA, QE Pro, and HR4000, however, all perform reasonably well for F₆₈₇ retrieval (e.g., for iFLD and SFM a R² between 0.88–0.91 and an RMSE between 0.34–0.38 mW m⁻² nm⁻¹ sr⁻¹ was observed). Due to the non-linear relationship of R and F across the O₂B band, both sFLD and 3FLD show a significant larger error compared to iFLD and SFM (R² ≤ 0.5 and RMSE = 1–8 mW m⁻² sr⁻¹ nm⁻¹). By using a more complex interpolation method, iFLD and SFM can better estimate R at the O₂B (R² = 1 and RMSE = 0.001 mW m⁻² sr⁻¹ nm⁻¹), a fact that substantially impacts F retrieval.

When retrieving F₆₈₇, the iFLD shows a tendency to underestimate high F₆₈₇ values, while low F₆₈₇ values are overestimated by both iFLD and SFM. This result does not agree with results presented by Jullita et al. [40] and Zhang et al. [44]. In their study, the same SFM implementation [40] was used as in the present study, but they found a 1:1 relation between retrieved and reference F₆₈₇. We hypothesize that this could be explained by differences on how the F₆₈₇ and R reference data were simulated. In the case of Zhang et al. [44], FluoSpec was used to generate the reference dataset. Jullita et al. [42] also defined the R and F₆₈₇ functions implemented in their R code using FluoSpec reference data. In our study, we used FluorSAIL3 to generate the reference dataset (to be consistent with the study by Damm et al. [32]) but used the SFM fitted for FluoSpec. It is known that F₆₈₇ modeled by FluorSAIL3 presents a taller and narrower red SIF peak compared to FluoSpec. Thus, most probably the shape of spectra of F₆₈₇ and R simulated with FluorSAIL3 may not be captured with the FluoSpec-trained SFM approach. To ensure that the results presented were not due to incorrect implementation of the SFM approach, we retrieved F₆₈₇ and F₇₆₀ with the original R code from Julitta et al. [40] available on GitHub platform at https://github.com/tommasojulitta and we obtained the same results as with our SFM approach. We suggest the systematic assessment of the effect in the near future to ensure that the SFM is flexible enough to cope with possible dynamics of SIF.

Figure 8 shows the RE of retrieved F obtained when including noise. In spite of its high SNR, the ASD shows the highest RE for F₇₆₀ and F₆₈₇, where sFLD performs worst (RE = 234.5 mW m⁻² nm⁻¹ sr⁻¹, 370.1%) and iFLD best (RE = 11.8 mW m⁻² nm⁻¹ sr⁻¹, 41.2%). The high-resolution sensors QEPro, MAYA, and HR4000 provide better F estimates with lower RE for SFM (F₇₆₀ = 4.5–4.9%, F₆₈₇ = 5.9–7.2%) compared to the iFLD method (F₇₆₀ = 4.7–9.3%, F₆₈₇ = 9.7–13.4%). From these three sensors, the QE Pro produces lowest RE for F₇₆₀ and the HR4000 lowest RE for F₆₈₇. The SFM furthermore shows the lowest error range and lowest mean RE for F₇₆₀ and F₆₈₇ (6.2–6.5% and 25.3–27.8%).

Importantly, when retrieving SIF with the sFLD or 3FLD, the relative error was 2–6 (QE Pro-F₇₆₀) and 6–11 (HR4000–F₆₈₇) times higher than with the SFM. Moreover, when retrieving SIF using the ASD simulated data, the error was 2–10 times higher than when using spectrometers with a finer spectral resolution, suggesting the need for care when interpreting SIF results retrieved with the sFLD and 3FLD and/or with spectrometers with SSI ≥ 1.4 nm and FWHM ≥ 3 nm.

In summary, the results presented in this study corroborate Damm et al. [32] results, where SNR was the most important parameter affecting the retrieval accuracy, followed by SR and SSI. Nevertheless, since in the study by Damm et al. [32], the spectrometers used present a relative high SNR, the accuracy and precision of retrieved F₇₆₀ and F₆₈₇ mainly depend on the SR and SSI. A finer spectral resolution generally reduces retrieval uncertainties and allows the use of the more sophisticated iFLD and SFM methods.

4.2. Uncertainties Caused by the Setup of Retrieval Methods

This section describes the impact of an incorrect implementation of the FLD and SFM methods for the SIF retrieval. For this purpose, retrieved and modeled F signals of a representative vegetation target were compared. The selected vegetation target represent the spectrum of a planophile canopy with a leaf area index of 1, a chlorophyll content of 80 μg/cm², and SIF efficiency factor of 2 (Table 3). The major sources of error, when retrieving SIF, can be attributed to:

(1)

Definition of the wavelength interval affected by O₂B and O₂A absorption according to instrumental spectral characteristics (FLD-based and SFM methods).

(2)

Selection of wavelengths at the shoulder and inside the absorption feature (FLD-based methods).

(3)

Interpolation strategy to estimate F and R at wavelengths affected by O₂ absorption (FLD-based and SFM method).

(4)

Definition of the b parameter (SFM method).

Error source 1: The F retrieval errors related to the definition of the wavelength representing the full oxygen absorption feature are shown in Figure 10. For this experiment, the left and right shoulder limits were simultaneously changed. For both retrievals, F₆₈₇ (Figure 10A,B) and F₇₆₀ (Figure 10C), the relative error of FLD-based methods, increases with an increased absorption range. The 3FLD was the most affected retrieval method (RE-F₆₈₇ 100–1200% and RE-F₇₆₀ 10–20%), followed by sFLD (RE-F₆₈₇ 10–70% and RE-F₇₆₀ 30–70%), and iFLD (RE-F₆₈₇ 2–70% and RE- F₇₆₀ 1–30%), which led to an over- and/or underestimation of E^↓ and L^↑ outside the absorption feature (Appendix D, green and blue crosses). The SFM was not affected by changes in the absorption range. For the FLD-based methods, the oxygen absorption range affects the definition of the shoulder wavelength position bands (error source 2) as well as the interpolation to model E^↓, L^↑, and R (error source 3). This is the reason an incorrect characterization of the oxygen absorption window has a significant impact in the retrieval of SIF when using FLD-based method. For instance, for the 3FLD method, an uncorrected definition of the O₂B absorption range (Figure 10A, RE-F₆₈₇ = 1200%) led to an underestimation of E^↓ and to an overestimation of L^↑ outside the absorption feature (Appendix D, blue crosses), which in the end results in an underestimation of the retrieved fluorescence (F = −0.1049 mW m⁻² nm⁻¹).

Error source 2: Retrieval errors related to definition of the shoulder wavelength position is shown in Figure 11. For this experiment, the absorption range was fixed between 686-697 nm (O₂B, Appendix C) and 759–770 nm (O₂A, Appendix C), and systematically changed the shoulder wavelength position and computed corresponding F retrieval errors (i.e., RE). First, the left shoulder wavelength was set to a fixed position and the right shoulder wavelength changed (Figure 11, O₂B from 695–710 nm and O₂A from 770–780 nm). Second, the left shoulder wavelength was modified by fixing the right one at a specific band (Figure 11 O₂B from 670–685 nm and O₂A from 750–760 nm). In general, the further away the shoulder wavelength position is from the absorption range feature, the greater is the error. For F₆₈₇ (Figure 11A,B), the 3FLD was most sensitive to changes in the shoulder wavelengths position (RE-F₆₈₇ 5–2500%), followed by sFLD (RE-F₆₈₇ 20–200%), and iFLD (RE-F₆₈₇ 5–35%). For F₇₆₀ (Figure 11C,D), the sFLD was most sensitive to changes of the shoulder position (RE- F₇₆₀ 200–800%), followed by 3FLD (RE- F₇₆₀ 10–20%), and iFLD (RE- F₇₆₀ 5–10%) method. SIF retrieval errors by the sFLD were two times higher when using the right shoulder, which indicates that the use of reference bands at the right shoulder should be avoided. The iFLD method is less sensitive to the shoulder wavelength position because, in contrast to sFLD and 3FLD methods, it uses the full spectrum to interpolate E^↓ and R in the absence of the oxygen absorption feature. Another source of error can occur when selecting the local maximum on both shoulders (Figure 11, green dots). The incorrect selection of these points can result in a RE of 60% for F₆₈₇ and RE 35% for F₇₆₀ (Figure 11, blue dots).

Error source 3: The interpolation strategy is of great importance, as can be seen from Figure 12. The distance of the selected shoulder points defines the interpolation range used for the SIF retrieval inside the absorption feature. The interpolation range was modified as follows: (i) Changing the left and right shoulder positions simultaneously (Figure 12A,C); (ii) only the right shoulder position (Figure 12B, O₂B 700–730 nm and Figure 12D; O₂A 780–820 nm); and (iii) only the left shoulder position (Figure 12B, O₂B 650–680 nm and Figure 12D, O₂A 720–750 nm). Figure 12 shows that the RE increase with increasing interpolation range resulting in values of 120% at F₆₈₇ and 30% at F₇₆₀ for the iFLD and a RE of 15% at F₆₈₇ and 30% at F₇₆₀ for the SFM. In this experiment, for both retrieval methods, a cubic spline function was used to interpolate R in the absorption feature. A splines function will go exactly through all the key points provided to perform the interpolation. This decreases the interpolation accuracy at the absorption feature when the key points selected are located some distance away from the absorption region. The SFM yields a lower RE than the iFLD because the parameters of the R interpolation are used as a first guess in a non-linear least square optimization function. Therefore, a small error in the R parameters first guess is corrected in the optimization process (Appendix C).

Error source 4 only relates to the SFM approach. For a good implementation of this method, the crucial issue is to define a set of functions that are able to model the broad range of R and F dynamics found under natural conditions. As introduced in Section 2 in this study, and after testing a different set of functions to model R and F, a cubic spline function to model R and a Gaussian function to model F were implemented [40]. Hence, for this specific implementation of the algorithm, the parameters that could be considered for the sensitivity analysis are (i) absorption feature range, (ii) F and R fitting interval, (iii) R interpolation, and iv) the Gaussian function parameter $b$, which controls the width of the red and far-red SIF modeled peaks. We will now focus on the $b$ parameter characterization. Figure 13 shows the impact on the retrieved SIF when the $b$ factor first guess changes from 0 to 100. The F retrieval is mostly sensitive to low $b$ values, for instance the RE varies between 10–310% (F₆₈₇) and 5–20% (F₇₆₀) for $b$ values between 1 to 3. When $b$ > 3, the RE decreases to 6% (F₆₈₇) and 1% (F₇₆₀). Similar to the R parameters, the value assigned to the b factor is used as a first estimate in a non-linear least square optimization function. A small error in the b factor first estimate is corrected in the optimization process (Appendix B).

5. Conclusions

For the first time, the retrieval errors of F₆₈₇ and F₇₆₀ have been systematically analyzed, taking into account two known error sources, the sensor configuration (i.e., SR, SSI, or SNR), and the retrieval method applied. Our extended analysis reveals further sources of error caused by the definition and implementation of retrieval approaches, including the selection of appropriate wavelengths representing the absorption feature and related shoulder values.

We identified the HR4000 and QE Pro spectrometers in combination with the SFM approach to provide the most reliable SIF retrieval with a relative error ≤ 6% and < 5% for F₆₈₇ and F₇₆₀, respectively, followed by iFLD with a relative error for F₆₈₇ ≤ 10% and F₇₆₀ < 5%. Regarding the retrieval method implementation, again, the SFM was the least affected by the incorrect definition of the oxygen absorption spectral window, interpolation strategy, or model parameters characterization. Although the SFM approach was found to yield highest SIF retrieval accuracies, we also observed a sensitivity of the SIF retrieval to the simulated training data underlying the SFM model implementation. In contrast to iFLD approach, the SFM is less sensitive to the instrument noise but more sensitive to R and F modeling. We suggest further investigation of this sensitivity is required.

We would like to encourage the scientific community to validate their implementation of the different retrieval methods by generating graphs of the most common error sources in the retrieval of SIF described in this paper (cf. Section 4.2). Furthermore, we recommend measuring fluorescence (e.g., grass) and non-fluorescence targets (e.g., soil) as part of the validation exercise (Dr. Andreas Burkart, personal communication). In the case of retrieving negative and/or extremely high SIF values in a fluorescence target (cf. Meroni et al. [9], Section 5 for SIF range variation) and SIF values different to zero in a non-fluorescence target, a revision of the implemented retrieval method must be performed.

Finally, we would like to emphasize the importance of instrument characterization [30], measurement protocol [31], and implementation of the retrieval methods (this paper). We propose that the insights provided in the three review papers can assist the scientific community to achieve higher accuracies and replicability of their SIF measurements to significantly contribute to advancing vegetation research.