Radiative Transfer Modeling of Phytoplankton Fluorescence Quenching Processes

Abstract

We report the first radiative transfer model that is able to simulate phytoplankton fluorescence with both photochemical and non-photochemical quenching included. The fluorescence source term in the inelastic radiative transfer equation is proportional to both the quantum yield and scalar irradiance at excitation wavelengths. The photochemical and nonphotochemical quenching processes change the quantum yield based on the photosynthetic active radiation. A sensitivity study was performed to demonstrate the dependence of the fluorescence signal on chlorophyll a concentration, aerosol optical depths and solar zenith angles. This work enables us to better model the phytoplankton fluorescence, which can be used in the design of new space-based sensors that can provide sufficient sensitivity to detect the phytoplankton fluorescence signal. It could also lead to more accurate remote sensing algorithms for the study of phytoplankton physiology.

1. Introduction

Solar induced fluorescence of chlorophyll a in phytoplankton is an important source of information on phytoplankton biomass variation [1,2]. Chlorophyll a concentration is a primary factor that influences the fluorescence signal [3], however, the signal is also strongly impacted by phytoplankton health (or physiology) as well as environmental and physical factors such as species composition, pigment packaging, nonphotochemical quenching (NPQ), and the availability of light and nutrients (notably biologically available iron) in the ocean surface layer [3,4,5]. A variety of algorithms were developed to study phytoplankton physiological state from satellite observations of chlorophyll a fluorescence [3,5,6,7,8,9].

To understand the fluorescence signal emanating from the ocean, it is important to have a rigorous radiative transfer model that is capable of simulating fluorescence strength in response to different physiological and environmental factors. In recent years, several approximate analytical models [10,11] and rigorous radiative transfer models [12,13,14,15,16,17,18,19] were developed to predict inelastic scattering processes in ocean waters, which included the methods of Monte Carlo [12,13], invariant embedding [14], matrix operator [15], successive order of scattering [16,17], and discrete ordinates [18,19]. The main inelastic scattering mechanisms considered in these works were Raman scattering, fluorescence by dissolved organic matter (FDOM), and fluorescence by chlorophyll a phytoplankton [14]. There remains, however, a lack of radiative transfer models that can account for the photo-compensation mechanisms [8]. This paper reports the first radiative transfer model that can simulate the impacts of both photochemical and nonphotochemical quenching effects on the phytoplankton fluorescence signal within the ocean and that is observed with space-based sensors. The primary inputs of the model are the inherent optical properties of the atmospheric and oceanic components and the physiological parameters which fully describe the fluorescence quantum yield as a function of spectral scalar irradiance in ocean waters.

Our model can be used to perform systematic studies of fluorescence in response to environmental factors and ocean water inherent optical properties, which in turn leads to better and more accurate remote sensing algorithms. The paper is organized in the following way. Section 2 outlines the theoretical implementation of photo-compensation in a radiative transfer model; Section 3 presents representative simulations on fluorescence in response to different atmosphere and oceanic conditions; Section 4 discusses a few limitations and some potential applications of our model; Section 5 summarizes the main conclusion.

2. Methods

Our radiative transfer model is based on the successive order of scattering (SOS) method for coupled atmosphere and ocean systems [20,21]. Recently inelastic scattering in ocean waters was integrated in the SOS code [16,17,22]. This paper uses the theoretical formulation in [17] for the implementation of inelastic scattering in radiative transfer. A key parameter which determines the fluorescence strength is the quantum yield ${\mathsf{\Phi }}_C$, which typically ranges between 0.01 and 0.05, and depends on phytoplankton species, physiological state, and the inherent optical properties of the water column [1,3,8,14]. We used a constant number for the quantum yield in Zhai et al. [17] to demonstrate the sensitivity of fluorescence on environmental parameters, which is applicable for a wide range of waters [6].

To account for the fluorescence quenching processes in the radiative transfer, we adopted here the model for quantum yield developed by Morrison et al. [1,8]:

$${\mathsf{\Phi }}_C=q_I{e}^{-\mathrm{IPAR}/E_T}\left({\mathsf{\Phi }}_{Cmin}{e}^{-\mathrm{IPAR}/E_k}+{\mathsf{\Phi }}_{Cmax}(1-e^{-\mathrm{IPAR}/E_k})\right),$$

(1)

where the parts in and outside of the parentheses represent photochemical and nonphotochemical quenching, respectively; $q_I$ is a parameter between 0 and 1 which stands for thermal dissipation of excess excitation energy through reaction center quenching. Alternatively, $q_I$ can be viewed as the relative probability of a photon absorbed by pigment molecules being delivered to a photosystem II reaction center [1,8]. IPAR is the instantaneous photosynthetically available radiation, which was obtained from the scalar irradiance $E_o$ via [23]:

$$\begin{array}{cc}\hfill \mathrm{IPAR}=& {\int }_{{\lambda }_1}^{{\lambda }_2}{E}_o\left(\lambda \right)\frac{\lambda }{hc}d\lambda ,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill E_o\left(\lambda \right)=& {\int }_{4\pi }I(\mathsf{\Omega },\lambda )d\mathsf{\Omega },\hfill \end{array}$$

(3)

where the symbol $\mathsf{\Omega }$ represents the solid angle; and I is the radiance. In the IPAR definition, ${\lambda }_1=400$ nm and ${\lambda }_2=700$ nm are normally used for the lower and higher wavelength limits. In this work, we adopted ${\lambda }_1=370$ nm and ${\lambda }_2=690$ nm to be consistent with the inelastic fluorescence source term limits (see Equation (4) below) [14,17]. The different limits are not expected to have significant impacts to our analysis and conclusions in this paper.

$E_k$ and $E_T$ are the light saturation parameters for photochemical and energy-dependent quenching, respectively; ${\mathsf{\Phi }}_{Cmin}$ and ${\mathsf{\Phi }}_{Cmax}$ are the minimum and maximum quantum yields, respectively. In our radiative transfer model, $q_I$, $E_T$, $E_k$, ${\mathsf{\Phi }}_{Cmin}$, and ${\mathsf{\Phi }}_{Cmax}$ are all input parameters, which may vary considerably depending on phytoplankton species and physiological states. $E_o$ and IPAR are calculated by the radiative transfer solver, given the solar spectral irradiance and inherent optical properties of the atmosphere and oceanic constituents.

First we run the radiative transfer solver through the whole range of excitation wavelengths to determine the scalar irradiance $E_o$ at different depths. IPAR is obtained through Equation (2). Then the quantum yield is calculated from Equation (1). The fluorescence scattering coefficient $b_C$ is determined by [14,17]:

$$b_C(z,\lambda ,{\lambda }_e)=a_p{\mathsf{\Phi }}_C\frac{1}{\sqrt{2\pi }{\sigma }_C}exp\left[-\frac{{(\lambda -{\lambda }_{C,0})}^2}{2{\sigma }_C^2}\right]\frac{{\lambda }_e}{\lambda },370\mathrm{nm}<{\lambda }_e<690\mathrm{nm}$$

(4)

where z denotes the vertical position; $\lambda$ and ${\lambda }_e$ are the emission and excitation wavelengths, respectively; $a_p$ is the phytoplankton absorption coefficient; ${\lambda }_{C,0}=685$ nm is the wavelength of maximum emission; and ${\sigma }_C=$ 10.6 nm is the Gaussian standard deviation [14]. The source function $S_C$ at emission wavelength for fluorescence is [14,17]:

$$S_C(z,\lambda )=\frac{1}{4\pi }{\int }_0^{\infty }{b}_C(z,\lambda ,{\lambda }_e)·E_o(z,{\lambda }_e)d{\lambda }_e,$$

(5)

which assumes that the fluorescence emission is isotropic [24].

The plane-parallel inelastic vector radiative transfer equation without thermal emission is:

$$\mu \frac{d\mathbf{L}(z,\mu ,\varphi ,\lambda )}{dz}=-c(z,\lambda )\mathbf{L}(z,\mu ,\varphi ,\lambda )+\mathbf{S}(z,{\mu }_0,\mu ,\varphi ,\lambda )+{\mathbf{S}}_i(z,{\mu }_0,\mu ,\varphi ,\lambda ),$$

(6)

where ${\mu }_0=cos{\theta }_0$; $\mu =cos\theta$; ${\theta }_0$ and $\theta$ are the solar and viewing zenith angles, respectively; $\varphi$ is the viewing azimuth; $\mathbf{L}$ = [I, Q, U, V]${}^{\mathrm{T}}$; I, Q, U, V are the Stokes parameters; and superscript T stands for transpose; $\mathbf{S}$ and ${\mathbf{S}}_i$ are the elastic and inelastic source terms, respectively; c is the extinction coefficient, also called beam attenuation coefficient in ocean optics . The inelastic source term is the summation of three mechanisms:

$${\mathbf{S}}_i={\mathbf{S}}_R+{\mathbf{S}}_Y+{\mathbf{S}}_C,$$

(7)

where the subscripts R, Y, and C stand for Raman scattering, FDOM, and chlorophyll fluorescense, respectively. The fluorescence source matrix ${\mathbf{S}}_C$ is assumed to be independent of polarization so that ${\mathbf{S}}_C={[S_C,0,0,0]}^{\mathrm{T}}$, where $S_C$ is given by Equation (5). The source function ${\mathbf{S}}_R$ and ${\mathbf{S}}_Y$ are the same as in [17]. With both elastic and inelastic source terms are known, the radiative transfer solver is used at emission wavelengths to find the total radiation field.

3. Results

We used the model to simulate the fluorescence signals for sensors located at both the top of the atmosphere (TOA) and the top of the ocean (TOO). The simulation was for a hypothetical hyperspectral sensor similar to the ocean color instrument onboard NASA’s next generation satellite: the Plankton, Aerosol, Cloud and ocean Ecosystem mission (PACE) [25]. The light field at the emission wavelength depends on the inherent optical properties at both the emission and excitation spectral ranges. Therefore the inherent optical properties of aerosols and hydrosols used in the simulation need to be spectrally consistent and realistic.

It was assumed that the inherent optical properties of the ocean are described using three components: pure sea water, colored dissolved organic matter (CDOM), and particles with properties covariant with the pigment chlorophyll a concentration [Chla], which are referred as the [Chla] covariant particles hereafter. Both pure sea water and the [Chla] covariant particles scatter and absorb light, while CDOM is assumed to absorb only. The absorption of pure sea water in Zhai et al. [17] was based on Pope and Fry [26]. In this work, we updated the absorption coefficient of pure sea water using a new dataset by Lee et al. [27]. The scattering coefficient of pure sea water is from [14,28]. The absorption coefficient of CDOM decays exponentially as wavelength increases with a decay constant of 0.018 nm${}^{-1}$ [17,29,30]. The absorption coefficient for the [Chla] covariant particles follows the parameterization in Bricaud et al. [31], which represents an average behavior of phytoplankton absorption with pigment packaging included. The extinction coefficient of the [Chla] covariant particles is the same as in Zhai et al. [17] (see Equation (16)), which is based on Voss et al. [32] and the International Ocean Color Coordinating Group (IOCCG) report [30]. The scattering coefficient is modeled as the difference between the extinction and absorption coefficients. The phase function of the ocean water is a weighted sum of the pure ocean water [14] and the [Chla] covariant particles, which is determined by the backscattering ratio (see Equations (18) and (19) in Zhai et al. [17]). The ocean water scattering Mueller matrix is the average ocean water measurement by Voss and Fry [33].

The atmosphere is a mixture of aerosols and molecules with molecular number density determined by the 1976 US standard atmosphere [34]. The depolarization of 0.0284 is used to in the calculation of the Rayleigh scattering matrix [35]. The aerosol scattering matrix follows the maritime aerosol with a humidity of 80% [36], and the aerosol vertical distribution is the average height distribution in [37]. The aerosol and Rayleigh scattering matrix are internally mixed by using their scattering coefficients as weighting factors [38]. Gas absorption has been incorporated, which includes ozone, oxygen, water vapor, nitrogen dioxide, methane, and carbon dioxide. First, the atmospheric radiative transfer simulator (ARTS) [39] and the HITRAN2012 database [40] are used to build a hyperspectral lookup tables for gas absorption coefficients of water vapor, oxygen, carbon dioxide, and methane. Absorptions by ozone and nitrogen dioxide are included separately based on the measurements in [41] and [42]. respectively. A hypothetical Gaussian instrument line shape (ILS) function is assumed with full width at half maximum of 5 nm. Within each spectral interval of 5 nm, the radiance is assumed to be an exponential function of the gas absorption optical depth so that we can solve a few wavelengths to fit the hyper-spectral variation due to gas absorption, which is then convolved with the ILS to simulate the response of a satellite sensor. This philosophy is similar to the double–k method [43]. The solar irradiance spectra from [44] is used as incident source of the radiative transfer system.

To simulate the quenching processes, the parameters in Equation (1) need to be assigned. We used $q_I=0.35$, ${\mathsf{\Phi }}_{Cmin}=0.03$, ${\mathsf{\Phi }}_{Cmax}=0.09$, $E_k=55$ μmol quanta m${}^{-2}$s${}^{-1}$, $E_T=634$ μmol quanta m${}^{-2}$s${}^{-1}$. These values are consistent with the findings from [8], though we recognized they are variable for different phytoplankton populations. In addition to the chlorophyll fluorescence, the Raman scattering and FDOM were also included in the simulation using the same scheme as [17]. This includes chlorophyll a concentration [Chla] = 0, 0.03, 0.1, 1, 10 mg/m${}^3$; aerosol optical depth at 550 nm ${\tau }_a$ = 0, 0.1, 0.2, 0.5; and solar zenith angle ${\theta }_s$ = 0, 30, 60, 78${}^{\circ }$. It is more convenient to study the radiance field in terms of reflectance, which is defined as:

$$\rho (z,\theta ,\varphi )=\pi ·\frac{I(z,\theta ,\varphi )}{E_d\left(z\right)}$$

(8)

where I is the radiance and $E_d$ is the downwelling irradiance.

Equation (1) builds a connection between the quantum yield and IPAR. It is informative to show how IPAR changes with depths in water and how the quantum yield responds to this variation. Figure 1 shows IPAR and the quantum yield as a function of depth measured from the water surface. For all [Chla] values, IPAR is approximately 1500 μmol quanta m${}^{-2}$s${}^{-1}$ at the surface, as IPAR is dominated by the atmospheric condition, which is kept the same in these simulation cases. As the depth increases, IPAR decreases faster for larger [Chla] values. As a consequence, the quantum yield starts with the same value at the surface, and increases to a maximum at different levels, and decreases again after that, for different [Chla] values. The vertical location of the peak quantum yield values is deeper for lower [Chla] values. The fluorescence signal at the ocean surface and TOA is essentially a vertical integration of the quantum yield, scalar irradiance, and particle absorption coefficients as shown in Equations (4)–(6). A higher location of the quantum yield peak indicates a stronger fluorescence signal, even though the minimum and maximum values of quantum yield are generally the same as suggested by Figure 1b.

Figure 2a shows the TOA reflectance ${\rho }_{\mathrm{TOA}}$ viewing at nadir ($\theta =0^{\circ }$ and $\varphi =0^{\circ }$) as a function of wavelength from 640 to 750 nm. Gas absorption features due to oxygen B band centered around 690 nm and water vapor centered around 650 nm and 720 nm are clearly seen. The fluorescence signal centered at 685 nm increases as chlorophyll a concentration increases. Figure 2b shows the relative change of reflectances ${\rho }_{\mathrm{TOA},\left[\mathrm{Chla}\right]}$ with respect to ${\rho }_{\mathrm{TOA},0}$, where the subscript 0 stands for [Chla] = 0 mg/m${}^3$. The peak value of relative change varies between 0.2% to 3.8%, which can be detected by ocean color sensors. The requirements for advanced ocean radiometers specify a signal to noise ratio of 1000 in the visible [45], which could theoretically detect signal variations at a level of 0.1%.

The fluorescence line height (FLH) algorithm developed for the NASA Moderate Resolution Imaging Specroradiometer (MODIS) calculates the baseline signals at 678 nm from measurements at 667 nm and 748 nm [46,47]. The algorithm was a compromise of available bands that sometimes lead to negative FLH [46]. In future satellite missions that measure hyperspectral radiances, fluorescence line height algorithms can be designed using optimal combination of wavelengths. In Figure 2a, we used the TOA reflectance at 640 nm, 710 nm, and 745 nm to obtain a second order polynomial fitting of the baseline signal level without fluorescence for the [Chla] = 10 mg/m${}^3$ case. The gas absorptions are minimal (column transmittance >97%) at these wavelengths. The peak of the fluorescence are near 680 nm because the signal at 685 nm is attenuated by oxygen B band absorption. The fluorescence line height (FLH) at 680 nm for [Chla] = 10 mg/m${}^3$, i.e., the difference between the observed measurement and the baseline, is shown in Figure 2a. The FLH algorithm proposed here can be applied to hyperspectral radiometric sensors onboard future satellites, such as the PACE mission [25].

The signal at TOA is dominated by the atmospheric scattering contribution such that the relative strength of chlorophyll a fluorescence is weak. This is not the case, however, for sensors in ocean waters. Figure 3 shows the nadir viewing reflectance at the top of the ocean (TOO) just below ocean surface as a function of wavelength. The reflectance at 680 nm varies from 0.0005 to 0.0035 for different [Chla] values. The spectral range of the fluorescence signals is approximately between 660 nm and 720 nm, which is primarily determined by ${\lambda }_{C,0}=$ 685 nm and ${\sigma }_C=$ 10.6 nm used in the simulation.

The quantum yield is strongly affected by the atmospheric conditions through their influence on IPAR. Figure 4a shows the nadir viewing reflectance at TOO for different aerosol optical depths at the reference wavelength of 550 nm. Figure 4b shows the relative signal variation with respect to ${\tau }_a$ = 0, which corresponds to a Rayleigh scattering only atmosphere. The fluorescence signal at TOO increases as the aerosol optical depth increases, which is somewhat counterintuitive. A detailed analysis revealed that this is a combined effect of a larger quantum yield value due to the decreased IPAR in Equation (1) and a smaller $E_o$ term in Equation (5). The effect of larger quantum yield dominates the effect so that the fluorescence signal is larger for smaller IPAR values. This phenomenon has been observed by the Geostationary Ocean Color Imager [5]. The relative difference is from 6% to 10 % for optical depth ranging form ${\tau }_a$ = 0.1 to 0.5.

The solar zenith angle is a primary factor which impacts the irradiance field in the ocean. Figure 5a shows the nadir viewing reflectance at TOO for different ${\theta }_s$ values. The relative difference with respect to the case of ${\theta }_s=0^{\circ }$ is shown in Figure 5b. As the solar zenith angle increases, IPAR in ocean water decreases and hence the fluorescence signal increases for the same reason indicated by Figure 4b [5]. The fluorescence signal has variations from 10% to 60% when solar zenith increases from 0 to 78${}^{\circ }$.

4. Discussion

Our treatment of photochemical and nonphotochemical quenching processes in chlorophyll fluorescence relies on the fidelity of Equation (1), a model developed by Morrison et al. [1,8]. The quantum yield quenching model (Equation (1)) was used in satellite remote sensing of quantum yield which reveals striking correlation with soluble iron deposition in ocean [3]. The same concept was also applied to geostationary satellite observation to study daily to seasonal dynamics in phytoplankton photophysiology [5]. Both Behrenfeld et al. [3] and O’Malley et al. [5] used a semianalytical fluorescence model, which involved several approximations. The associate uncertainty of these simplifications has to be evaluated with rigorous radiative transfer simulation. Our newly developed model serves well on this purpose. It can also be used to generate synthetic dataset for testing future satellite algorithms which retrieves chlorophyll and/or phytoplankton physiology parameters

In vertically inhomogeneous water columns, the vertical location of phytoplankton has significant impacts on the water leaving signal and consequently on retrieved [Chla] values [48]. The effects of vertical distribution of phytoplankton to the water leaving fluorescence signal were not studied in this paper but will be the subject for a future study. The fluorescence signal variation at TOA is much suppressed by the dominating effect of atmosphere contribution to the radiation field. Nevertheless, the FLH has been detected and retrieved by MODIS, though nonphysical negative FLH data have been obtained in complex waters. This can be improved by using hyperspectral measurements made by the PACE ocean color instrument. One possible way to achieve this is to derive the fluorescence line height using a combination of 640 nm, 680 nm, 710 nm, and 745 nm.

The radiative transfer model can also be used to study the impacts of system configurations to the retrieval algorithms, for example, vertical distribution of phytoplankton blooms and atmospheric particles like aerosols and clouds. Recently a new trend in ocean color remote sensing is to inverse system parameters using nonlinear least squares fitting procedures, in which radiative transfer models are called dynamically to fit radiometric measurements [49,50,51]. Our model can be used for fitting fluorescence signals by varying phytoplankton physiology parameters in the same mathematically frame. Moreover, our radiative transfer model can be used to explore the information content regarding desired retrieval parameters [52].

The physiology parameters $q_I$, $E_T$, and $E_k$ in Equation (1) depend on phytoplankton species and environmental factors including temperature, nutritions, etc. [8]. There is a lack of understanding in the literature on how these parameters varies with different environmental factors, which remains a future research direction in the marine biological community. In addition, we used the “average” bio-optical models of the inherent optical properties outlined in Section 3. It has been demonstrated that large natural variations exist in these relations [53]. To resolve the natural variations, we could use specific models in which regional and seasonal variability is included [54].

5. Conclusions

We report an exact radiative transfer solution that includes both photochemical and nonphotochemical quenching effects. Major gas absorptions in the atmosphere are considered in order to accurately simulate the scalar irradiance in ocean waters, which in turn regulates the quantum yield due to the quenching processes. Moreover, the scalar irradiance in ocean waters changes the inelastic source terms for fluorescence directly. The actual fluorescence signal depends on the compound effects of these two effects. The utility of this model is illustrated by using it to simulate the effect of photochemical and nonphotochemical quenching on radiance measured at TOA for an upcoming satellite and at TOO for an in situ sensor. Besides the obvious dependence on the chlorophyll a concentration, the fluorescence signal at TOO also varies for different aerosol optical depths and solar zenith angles. The relative variation of fluorescence ranges from 5% to 60% under different illumination conditions.