Biases in estimated vegetation indices from observations under cloudy conditions

Abstract

Field observations of vegetation indices (VIs) are derived from ratios of spectral reflectance data that are collected by drones and aircraft, providing higher spatial resolution than satellites. These reflectance data require periodic reference measurements over calibrated reflectance panels under cloud-free conditions. However, the reference measurements are partly performed in cloudy situations with the effect that wavelength-dependent scattering and absorption of solar radiation by clouds affects the subsequently derived VIs. This paper quantifies these effects using combined atmosphere-vegetation radiative transfer (RT) simulations. We study the general case when VIs are obtained from reflectance ratios of two wavelengths, and for the special cases of the normalized difference vegetation index (NDVI), the normalized difference water index (NDWI), and the enhanced vegetation index (EVI). For the general case of two-band VIs the lowest sensitivity to cloud changes was found for wavelength combinations below 1400 nm and outside the water vapor absorption bands. The NDVI was almost insensitive to changes in cloud conditions, while greater biases were identified for the NDWI. The EVI was most susceptible to cloud changes, with biases of 0.2 in the selected example. This lead to biases in the estimated leaf area index of 0.9. Biophysical properties derived from EVI, such as gross primary product, are also affected with variations of up to ±2 gCm-2d-1 in the selected cases.

1 Introduction

Numerous vegetation indices (VIs) have been developed to remotely characterize vegetation. They are based on ratios of vegetation-reflected radiation in wavelength regions that are sensitive to chlorophyll and liquid water absorption, thus indirectly providing information on plant vitality, productivity, and photosynthetic activity, and to monitor anthropogenic degradation . Established VIs are, for example, the normalized difference vegetation index (NDVI), the normalized difference water index (NDWI), and the enhanced vegetation index (EVI) . They were originally developed for satellite observations and customized to spectral bands of satellites such as Landsat , the Moderate Resolution Imaging Spectroradiometer MODIS; , or Sentinel-2 . In the course of the time, the available measurement platforms diversified, additionally including aircraft and unmanned aerial vehicles UAVs, .

All remotely sensed VIs rely on the spectral reflectance $\mathit{\rho }\left(\mathit{\lambda }\right)$. In the general form, $\mathit{\rho }\left(\mathit{\lambda }\right)$ is defined as the ratio of the surface-reflected radiation flux to the incident radiation flux . The angular distribution of the reflected radiation is characterized by the bidirectional reflectance distribution function (BRDF), which describes the surface-specific scattering of radiant energy . The magnitude of the surface-reflected radiation depends on factors such as topography and varying illumination and atmospheric conditions, including solar zenith angle, aerosols, and clouds . For brevity and simplicity, we restrict the definition of $\mathit{\rho }\left(\mathit{\lambda }\right)$ by assuming idealized conditions of purely diffuse illumination and Lambertian surface reflection, under which reflectance becomes independent of viewing geometry and simplifies to the following:

1$$\mathit{\rho }\left(\mathit{\lambda }\right)={\displaystyle \frac{I^{\uparrow }\left(\mathit{\lambda }\right)\cdot \mathit{\pi }\cdot \mathrm{sr}}{F^{\downarrow }(\mathit{\theta },\mathit{\lambda })}},$$ where $\mathit{\lambda }$ symbolizes the wavelength, $\mathrm{sr}$ the unit steradian, $I^{\uparrow }\left(\mathit{\lambda }\right)$ the spectral, upward (reflected) radiance in units of W m⁻² nm⁻¹ sr⁻¹, and $F^{\downarrow }(\mathit{\theta },\mathit{\lambda })$ the downward spectral irradiance in units of W m⁻² nm⁻¹, which also depends on the solar zenith angle $\mathit{\theta }$. The cloud optical thickness $\mathit{\tau }\left(\mathit{\lambda }\right)$ is a measure of the extinction of radiation for a vertical path through the cloud, serving as vertical coordinate. The spectral reflectance ranges between 0 and 1, except under broken cloud conditions where radiation enhancement occurs, leading to reflectance values exceeding unity . In neighboring regions, reflectance is lower. On average, this results in energy conservation. Radiation that has never been scattered represents the direct radiation component $F_{\mathrm{dir}}^{\downarrow }(\mathit{\theta },\mathit{\lambda })$, while radiation that has undergone at least one scattering event contributes to diffuse irradiance $F_{\mathrm{dif}}^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })$. Thus, the total downward irradiance $F^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })$ is composed of: 2$$F^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })=F_{\mathrm{dir}}^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })+F_{\mathrm{dif}}^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda }).$$ The related direct fraction $f_{\mathrm{dir}}(\mathit{\theta },\mathit{\tau },\mathit{\lambda })$ is defined by: 3$$f_{\mathrm{dir}}(\mathit{\theta },\mathit{\tau },\mathit{\lambda })={\displaystyle \frac{F_{\mathrm{dir}}^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })}{F^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })}}.$$

Remote sensing techniques to estimate VIs from satellite measurements require cloud-free conditions. Cloud masks are used to flag cloud-contaminated pixels, which are excluded from the VI calculation. The remaining pixels are assumed to be unbiased by clouds. Figure A schematically illustrates the measurement principle. However, for all types of airborne measurements from drones or aircraft that are maybe performed below clouds, see Fig. B and C, the downward spectral irradiance $F^{\downarrow }\left(\mathit{\lambda }\right)$ is altered by the presences of clouds. The impact of clouds depends on $\mathit{\tau }$, droplet/ice crystal size, and wavelength . The contribution of scattering and absorption is spectrally dependent, with scattering most pronounced in the visible near–infrared (VNIR, 0.3–1.0 $\mathrm{\mu }$m) part of the solar spectrum and absorption in the shortwave infrared (SWIR, 1.0–2.5 $\mathrm{\mu }$m) wavelength range. The spectral dependency leads to a spectral slope in $F^{\downarrow }\left(\mathit{\lambda }\right)$, $I^{\uparrow }\left(\mathit{\lambda }\right)$, and $\mathit{\rho }\left(\mathit{\lambda }\right)$ . Consequently, changes in cloud properties and illumination lead to biases in estimated VIs.

Figure 1

(A) Schematic illustration of measured spectral upward radiance $I^{\uparrow }$ under cloud-free conditions. (B) and (C) Schematic illustration of below cloud measurements that can be subject to a change in cloud conditions from ${\mathit{\tau }}_{\mathrm{cal}}$ (during reflectance panel (RP) calibration) to cloud conditions ${\mathit{\tau }}_{\mathrm{mea}}$ (during actual measurements). The spectral total downward irradiance $F^{\downarrow }$ and spectral upward radiance $I^{\uparrow }$ are affected by the cloud conditions. The reflectance ${\mathit{\rho }}_{\mathrm{RP}}$ of the RP panel is known from its specifications and ${\mathit{\rho }}_{\mathrm{mea}}$ is the reflectance of the vegetation. The transfer functions ${\mathit{\kappa }}_{\mathrm{cal}}$ and ${\mathit{\kappa }}_{\mathrm{mea}}$ connect measured counts and $\mathit{\rho }$.

[Figure omitted. See PDF]

To obtain $\mathit{\rho }\left(\mathit{\lambda }\right)$ for VI retrievals during field studies, calibrated $I^{\uparrow }(\mathit{\tau },\mathit{\lambda })$ and $F^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })$ are required, which are sometimes measured with dedicated, calibrated sensors . Alternatively, $F^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })$ is simulated using radiative transfer (RT) models. In either case, absolute measurements are complicated and simulations introduce uncertainties. For practical applications, absolute measurements are often avoided by performing relative measurements, using reflectance panels (RPs) characterized by a well-known reflectance ${\mathit{\rho }}_{\mathrm{RP}}\left(\mathit{\lambda }\right)$ .

Using a RP allows a transfer calibration. Digital spectrometers register digital counts $S$, which are related to ${\mathit{\rho }}^{\uparrow }(\mathit{\tau },\mathit{\lambda })$ by a calibration factor $\mathcal{C}\left(\mathit{\lambda }\right)$. Including this calibration factor in Eq. () yields: 4$$\mathit{\rho }\left(\mathit{\lambda }\right)={\displaystyle \frac{S^{\uparrow }(\mathit{\tau },\mathit{\lambda })\cdot {\mathcal{C}}^{\uparrow }\left(\mathit{\lambda }\right)\cdot \mathit{\pi }\cdot \mathrm{sr}}{F^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })}},$$ where $S^{\uparrow }(\mathit{\tau },\mathit{\lambda })$, $I^{\uparrow }(\mathit{\tau },\mathit{\lambda })$, and $F^{\downarrow }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })$ below clouds are determined by $\mathit{\tau }$. Equation () can be simplified to: 5$$\mathit{\rho }\left(\mathit{\lambda }\right)=S^{\uparrow }(\mathit{\tau },\mathit{\lambda })\cdot \mathit{\kappa }(\mathit{\theta },\mathit{\tau },\mathit{\lambda }),$$ with the transfer function $\mathit{\kappa }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })$ given by: 6$$\mathit{\kappa }(\mathit{\lambda },\mathit{\tau },\mathit{\theta })={\displaystyle \frac{{\mathcal{C}}^{\uparrow }\left(\mathit{\lambda }\right)\cdot \mathit{\pi }\cdot \mathrm{sr}}{F^{\downarrow }(\mathit{\tau },\mathit{\lambda })}}.$$ In a first step during the relative measurements, the RP is overflown and the signal is taken for calibration (Fig. B), where $\mathit{\rho }\left(\mathit{\lambda }\right)$ in Eq. () is set to the well-defined ${\mathit{\rho }}_{\mathrm{RP}}\left(\mathit{\lambda }\right)$ and $S^{\uparrow }\left(\mathit{\lambda }\right)$ is the number of counts registered by the radiance sensor above the RP. The transfer function ${\mathit{\kappa }}_{\mathrm{cal}}(\mathit{\theta },{\mathit{\tau }}_{\mathrm{cal}},\mathit{\lambda })$ determined during the calibration procedure now represents the relationship between the recorded digital counts and the surface reflectance under the illumination conditions at that time. During actual measurements over vegetation (Fig. C), the reflectance ${\mathit{\rho }}_{\mathrm{mea}}\left(\mathit{\lambda }\right)$ of the vegetated surface is determined applying $\mathit{\kappa }(\mathit{\theta },{\mathit{\tau }}_{\mathrm{cal}},\mathit{\lambda })$: 7$${\mathit{\rho }}_{\mathrm{mea}}\left(\mathit{\lambda }\right)=S^{\uparrow }\left(\mathit{\lambda }\right)\cdot {\mathit{\kappa }}_{\mathrm{cal}}(\mathit{\theta },{\mathit{\tau }}_{\mathrm{cal}},\mathit{\lambda }).$$ Although RP calibration avoids the determination of the individual factors included in ${\mathit{\kappa }}_{\mathrm{cal}}(\mathit{\theta },{\mathit{\tau }}_{\mathrm{cal}},\mathit{\lambda })$, the dependence of $\mathit{\kappa }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })$ on $\mathit{\theta }$ and $\mathit{\tau }$, requires frequent RP overflights to account for illumination changes and to obtain updated transfer functions.

presented example time series of $F^{\downarrow }$ from a field campaign and compared them with cloud-free simulations. They determined single-point enhancements of $F^{\downarrow }$ up to 50 % due to cloud-side enhancements. Short-term fluctuations of $F^{\downarrow }$ have been found on time scales down to 300 s , 60 s , 20 s , and also below 1 s . Consequently, $\mathit{\kappa }(\mathit{\theta },\mathit{\tau },\mathit{\lambda })$ and thus ${\mathit{\rho }}_{\mathrm{mea}}\left(\mathit{\lambda }\right)$ vary on the same time scales. Therefore, even regular sequences of RP measurements seem insufficient to capture fluctuations caused by broken clouds, either by the clouds themselves or by their shadows, which may lead to potential errors in the retrieved VIs . Even when $F^{\downarrow }\left(\mathit{\lambda }\right)$ is measured by a dedicated sensor or frequent measurements of the RP are performed, the sole presence of clouds causes a change in ${\mathit{\rho }}_{\mathrm{mea}}\left(\mathit{\lambda }\right)$, since clouds change the ratio of direct and diffuse radiation, which determines how radiation is reflected by a non-isotropic surface .

Here we present coupled simulations using the atmospheric RT model “library for Radiative transfer” libRadtran and the vegetation RT model “Soil Canopy Observation of Photosynthesis and Energy fluxes” SCOPE2.0. The iterative model coupling allows to account for mutual influences between cloud and canopy. The coupled simulations are used to systematically investigate the impact of different cloud conditions on general two-band VIs, and in particular on NDVI, NDWI, and EVI. The aim is to quantify the impact of clouds on below-cloud airborne observations that are based on relative measurements using RP calibration thus providing approximate associated errors in estimated VI. The VIs are calculated using the wavelength bands of the multispectral instrument (MSI) onboard the Sentinel-2 satellite. This introduction is followed by a brief overview of the main terminology, the applied RT models, and a short review of the model coupling in Sect. . This is followed by the results in Sect. , where the effect of clouds and changing cloud conditions on estimated VIs are presented. Section  provides a summary of the results and outlines the limitations of the simulations.

2 Vegetation indices and radiative transfer simulations

2.1 Definition of vegetation indices

Vegetation indices are based on ratios of $\mathit{\rho }\left(\mathit{\lambda }\right)$ at several wavelengths (mostly pairs). The exact center wavelength and width used in the calculation of VIs depend on the characteristics of the observing instrument. In the case of the Sentinel-2 satellites, several wavelength combinations are suitable for VI retrievals . A subset of the Sentinel-2 wavelength bands are listed in Table . Unless otherwise noted, the Sentinel-2 wavelength bands were used to calculate VIs throughout this paper by applying a boxcar-like spectral response function. The Sentinel-2 wavelength bands are considered representative of other measurement platforms and sensors. UAVs may be equipped with spectral sensors that cover similar wavelength bands, or in the case of multispectral sensors, spectral integration is performed over similar wavelength bands .

Table 1

Sentinel-2 wavelength bands from the multispectral instrument (MSI) following .

Band	Center wavelength	Spectral width
number	[nm]	[nm]
B2	490	$\pm$33
B3	560	$\pm$35
B4	665	$\pm$15
B8	842	$\pm$56
B8a	865	$\pm$10
B9	945	$\pm$10
B10	1375	$\pm$15
B11	1610	$\pm$45
B12	2190	$\pm$180

Multiple two-band VIs with the index value $\mathit{\gamma }$ exist that follow the general form of a band transformation: 8$$\mathit{\gamma }={\displaystyle \frac{\mathit{\rho }\left({\mathit{\lambda }}_{\mathrm{1}}\right)-\mathit{\rho }\left({\mathit{\lambda }}_{\mathrm{2}}\right)}{\mathit{\rho }\left({\mathit{\lambda }}_{\mathrm{1}}\right)+\mathit{\rho }\left({\mathit{\lambda }}_{\mathrm{2}}\right)}},$$ where ${\mathit{\lambda }}_{\mathrm{1}}$ and ${\mathit{\lambda }}_{\mathrm{2}}$ represent a pair of individual wavelengths or narrow bands with band centers between 400 and 2400 nm wavelengths.

An example of a two-band VI is the normalized difference water index (NDWI). Two versions of the NDWI exist. The NDWI proposed by uses a wavelength combination of 980 and 1240 nm; it is subsequently labeled with ${\mathrm{NDWI}}_{\mathrm{1240}}$. provided the alternative ${\mathrm{NDWI}}_{\mathrm{1640}}$, which was designed to be less sensitive to saturation with respect to the water content in plant matter compared to ${\mathrm{NDWI}}_{\mathrm{1240}}$. Another commonly used VI that follows Eq. () is the normalized difference vegetation index NDVI; , with $\mathit{\rho }\left({\mathit{\lambda }}_{\mathrm{1}}\right)$ from B4 and $\mathit{\rho }\left({\mathit{\lambda }}_{\mathrm{2}}\right)$ from B8 from the MSI Sentinel-2 channels. By construction, the NDVI ranges between $-$1 and 1. An NDVI close to unity indicates vital and productive vegetation, while stressed vegetation results in NDVI values around 0.2. Measurements above bare soil will give a NDVI around 0, while measurements over water will yield negative NDVI. Even though the NDVI proofed to be useful it also showed sensitivity to surface brightness, and scattering and absorption by atmospheric constituents . The NDVI also becomes saturated and insensitive to vegetation structure when the canopies become dense, i.e., high leaf density .

These issues led to the development of the enhanced vegetation index EVI; , which was derived to overcome the shortcomings of the NDVI . The EVI attempts to account for aerosol loading, surface reflectance, and other factors by including spectral information from the “blue” region (${\mathit{\rho }}_{\mathrm{B}\mathrm{2}}$) of the VNIR spectrum. The EVI is calculated for MSI Sentinel-2 bands by: 9$$\mathrm{EVI}=G\cdot {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{B}\mathrm{8}}-{\mathit{\rho }}_{\mathrm{B}\mathrm{4}}}{{\mathit{\rho }}_{\mathrm{B}\mathrm{8}}+C_{\mathrm{1}}\cdot {\mathit{\rho }}_{\mathrm{B}\mathrm{4}}-C_{\mathrm{2}}\cdot {\mathit{\rho }}_{\mathrm{B}\mathrm{2}}+L}},$$ with the scaling factors $G=\mathrm{2.5}$, $C_{\mathrm{1}}=\mathrm{6}$, $C_{\mathrm{2}}=\mathrm{7.5}$, and $L=\mathrm{1}$ .

2.2 Radiative transfer simulations

Equation  shows that the reflectance of a surface is determined by several factors, including $\mathit{\theta }$, the observation geometry, as well as the direct and diffuse components of $F^{\downarrow }\left(\mathit{\lambda }\right)$ that are determined by $\mathit{\tau }$ (see Appendix ). Furthermore, radiation interactions may occur between the surface and the cloud, which can be accounted for by iterative coupling of atmosphere and vegetation RT models . In the present paper we use the same model coupling setup introduced and described by .

2.2.1 Atmospheric radiative transfer model libRadtran

The atmospheric RT above the canopy was simulated with the library for Radiative transfer libRadtran, . The RT equation was calculated with the one-dimensional solver “Discrete-Ordinate-Method Radiative Transfer” DISORT, using 12 streams. Clouds were assumed to be homogeneous. A low-level and mid-level warm stratus or altostratus were included by liquid water clouds between 3 and 3.5 km altitude with a fixed droplet effective radius of 10 $\mathrm{\mu }$m . A high-level cirrostratus was included by ice water clouds between 11 and 11.5 km altitude with a fixed ice effective radius of 85 $\mathrm{\mu }$m . Aggregated ice particles with moderate surface roughness were assumed to represent mature ice crystals . The ice particle parametrization after was applied. The liquid water and ice water path were scaled so that a specific $\mathit{\tau }(\mathit{\lambda }=\mathrm{550}\mathrm{nm})$ at 550 nm wavelength was reached, with all other wavelengths being scaled considering the wavelength dependence of $\mathit{\tau }$. For simplicity, $\mathit{\tau }(\mathit{\lambda }=\mathrm{550}\mathrm{nm})$ is referred to as $\mathit{\tau }$ in the following. The incoming spectral irradiance at the top of atmosphere was represented by the solar reference spectrum provided by . Molecular absorption was considered using the parameterization of the “medium” resolution of . A default aerosol distribution after was applied, representing rural type aerosol in the boundary layer, back-ground aerosol above 2 km during spring-summer, and a visibility of 50 km. Atmospheric profiles of air temperature, humidity, and gas concentrations were taken from the mid-latitude summer profile “afglms” . Absorption by water vapor and other atmospheric trace gases was accounted for in the simulations . The iteration process was first started by running libRadtran with an initial guess for the surface albedo. The “mixed-forest” albedo was taken from the IGBP database . After one iteration cycle, the surface albedo determined during the iterative model coupling process was used . The RT simulations of $F^{\uparrow }\left(\mathit{\lambda }\right)$, $F^{\downarrow }\left(\mathit{\lambda }\right)$, and $I^{\uparrow }\left(\mathit{\lambda }\right)$ spanned a wavelength range from 0.4 to 2.4 $\mathrm{\mu }$m. The output of libRadtran was specified to be at an altitude of 40 m above the canopy to be representative for UAV measurements. The radiances were simulated for a nadir sensor geometry.

2.2.2 Vegetation radiative transfer model SCOPE2.0

The solar RT within vegetation was simulated with the Soil Canopy Observation of Photosynthesis and Energy fluxes version 2 SCOPE2.0, . The surface albedo was determined using the Brightness-Shape-Moisture (BSM) model . Irradiance and radiance simulations in the solar part of the spectrum have been performed for a wavelength range from 0.4 to 2.4 $\mathrm{\mu }$m. The angular dependence of $I^{\uparrow }\left(\mathit{\lambda }\right)$ is considered for by the actual illumination and observation geometries, the direct and diffuse $F^{\downarrow }\left(\mathit{\lambda }\right)$, and the parameterized intrinsic reflectivity within SCOPE2.0. An important parameter to describe the radiative properties of a canopy is the leaf area index LAI, , typically ranging between 0 and 12. It provides a measure of the accumulated, one-sided area of leaves per unit of ground area given in units of m² m⁻². A constant LAI of 3 m² m⁻² was used in the simulations, corresponding to the typical LAI of temperate vegetation types such as pine forests or various grasslands. The leaf angle distribution (LAD) is another important parameter controlling the RT in a canopy . It is a measure of the overall orientation of the leaf ensemble of a tree. Therefore, it influences the area of a leaf illuminated by the Sun with respect to the total one-sided leaf area . The simulations were primarily considered for the general case of a spherical LAD, where all leaf angles have a similar probability . Further simulations have been performed to also include the erectophile and the planophile LADs. The static LAD values obviously ignore that short term changes in weather and illumination conditions affect the LAD , which may likewise have significant effects on VIs , but is nevertheless common practice in many RTM parameterizations. Table  provides an overview of the selected parameters for the vegetation RT simulations. The parameters were selected on the basis of their relevance to surface reflectance within the visible and near-infrared wavelength ranges. Their individual relevance was estimated in a sensitivity study by .

Table 2

Selected configuration of the SCOPE2.0 simulations.

Description	Symbol	Setting	Unit
Leaf chlorophyll concentration	$C_{\mathrm{ab}}$	40	$\mathrm{\mu }$g cm−2
Leaf carotenoid concentration	$C_{\mathrm{ca}}$	10	$\mathrm{\mu }$g cm−2
Leaf water equivalent layer	$C_{\mathrm{w}}$	0.009	cm
Leaf structure parameter	$N$	2.1	Unitless
BSM model parameter for soil brightness	$B$	0.5	Unitless
Vegetation height	$h_{\mathrm{c}}$	20	m
Output height	$h_{\mathrm{out}}$	40	m

3 Results and discussion

The one-dimensional simulations combining libRadtran and SCOPE2.0 are to be interpreted as synthetic measurements of $F^{\downarrow }\left(\mathit{\lambda }\right)$, $I^{\uparrow }\left(\mathit{\lambda }\right)$, and $\mathit{\rho }\left(\mathit{\lambda }\right)$ under varying atmospheric conditions, ranging from cloud-free to stratiform clouds, with values of $\mathit{\tau }$ between 0 and 40, and $\mathit{\theta }$ between 25 and 70$\mathit{°}$. By assuming various combinations of ${\mathit{\tau }}_{\mathrm{cal}}$, which prevailed during calibration, and ${\mathit{\tau }}_{\mathrm{mea}}$, present during actual measurements, the effects of changes in cloud conditions between RP measurements on ${\mathit{\rho }}_{\mathrm{mea}}\left(\mathit{\lambda }\right)$ and estimated VIs were determined. This was performed using simulated $I^{\uparrow }({\mathit{\tau }}_{\mathrm{mea}},\mathit{\lambda })$ and $F^{\downarrow }({\mathit{\tau }}_{\mathrm{cal}},\mathit{\lambda })$ in Eq. (), which is synonymous with the assumption of constant values of ${\mathit{\kappa }}_{\mathrm{cal}}$. The values of simulated ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$ range from 0 to 40 for the liquid water cloud and the ice water cloud. Two aspects influenced the chosen ranges. First, simulated ice clouds with $\mathit{\tau }$ up to 6 can be considered as a high-level cirrus. The natural variation in non-spherical ice particle shapes was taken into account by selecting aggregated ice particles (see Sect. ). Second, liquid water clouds were considered by simulating mid-level, continental stratus clouds, which are known to often reach values of up to $\mathit{\tau }=\mathrm{40}$ . The subsequent analysis focuses primarily on the spherical LAD, which is considered the general case. The effect of the selected LAD is discussed in Appendix .

3.1 Influence of clouds on the spectral reflectance

have shown the influence of molecules and aerosols on direct and diffuse $F^{\downarrow }\left(\mathit{\lambda }\right)$, the effects on $F^{\uparrow }\left(\mathit{\lambda }\right)$, and the resulting albedo effects over vegetated areas based on analytical equations. used coupled atmosphere–vegetation radiative transfer models to investigate these effects, explicitly considering clouds with spectrally dependent scattering and absorption in the atmosphere. An increase in $\mathit{\tau }$ leads to a decrease in $F_{\mathrm{dir}}^{\downarrow }\left(\mathit{\lambda }\right)$, while the response of $F_{\mathrm{dif}}^{\downarrow }\left(\mathit{\lambda }\right)$ depends on $\mathit{\tau }$. For values of $\mathit{\tau }$ less than 4 to 6, $F_{\mathrm{dif}}^{\downarrow }\left(\mathit{\lambda }\right)$ first increases and then decreases as $\mathit{\tau }$ increases further. The total $F^{\downarrow }\left(\mathit{\lambda }\right)$ and $f_{\mathrm{dir}}\left(\mathit{\lambda }\right)$ both continuously decrease as $\mathit{\tau }$ increases. In addition, $F^{\uparrow }\left(\mathit{\lambda }\right)$ became less sensitive to changes in $\mathit{\theta }$. Lastly, the presence of clouds modulates the incoming radiation spectrally by shifting the incoming radiation towards shorter wavelengths, as clouds primarily scatter radiation at shorter wavelengths and absorb radiation at longer wavelengths. also showed that radiative interactions between the canopy and the cloud base increase $F_{\mathrm{dif}}^{\downarrow }\left(\mathit{\lambda }\right)$ and albedo compared to cloud-free conditions. The present paper focuses on the related effects on $I^{\uparrow }\left(\mathit{\lambda }\right)$ and $\mathit{\rho }\left(\mathit{\lambda }\right)$.

3.1.1 Effect of diffuse radiation on spectral reflectance above vegetation

First, we consider situations, where the cloud conditions are similar during the RP calibration and the measurements (${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}$). Figure a shows $\mathit{\rho }\left(\mathit{\lambda }\right)$ for $\mathit{\theta }=\mathrm{25}\mathit{°}$ with the lowest $\mathit{\rho }\left(\mathit{\lambda }\right)$ under cloud-free conditions (${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}=\mathrm{0}$, black line). With increasing values of $\mathit{\tau }$, $\mathit{\rho }\left(\mathit{\lambda }\right)$ also increases, especially between 750 and 1300 nm wavelengths, where vegetation is generally characterized by high reflectivity and reflectance. The sensitivity of $\mathit{\rho }\left(\mathit{\lambda }\right)$ to $\mathit{\tau }$ is highest for small values of $\mathit{\tau }$ and quickly approaches an asymptotic value when $F_{\mathrm{dif}}^{\downarrow }\left(\mathit{\lambda }\right)$ dominates the radiation field. Figure b shows the ratio of $\mathit{\rho }\left(\mathit{\lambda }\right)$ between cloudy and cloud-free conditions. Independent of $\mathit{\tau }$, the ratio has a distinct spectral dependency with pronounced water absorption bands. Under cloudy conditions ($\mathit{\tau }=\mathrm{10}$) and for $\mathit{\theta }=\mathrm{25}\mathit{°}$, $\mathit{\rho }\left(\mathit{\lambda }\right)$ decreases by up to 30 % for wavelength below 750 nm, while for longer wavelengths $\mathit{\rho }\left(\mathit{\lambda }\right)$ increases by up to 8 % at about 780 nm. An exception is the dip in $\mathit{\rho }\left(\mathit{\lambda }\right)$ at about 1450 nm wavelength and wavelengths greater than 1800 nm. For $\mathit{\theta }>\mathrm{60}\mathit{°}$, an increase in $\mathit{\tau }$ leads to an increase in $\mathit{\rho }\left(\mathit{\lambda }\right)$ throughout the simulated wavelength range.

Figure 2

(a) Simulated spectral reflectance $\mathit{\rho }\left(\mathit{\lambda }\right)$ caused by an ice cloud for a solar zenith angle $\mathit{\theta }$ of 25$\mathit{°}$, and constant cloud optical thickness during calibration and the actual synthetic measurements (${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}$). The simulations base on a spherical LAD. (b) Ratios of spectral reflectance with ${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{cal}}>\mathrm{0}$ with respect to the spectral reflectance under cloud-free conditions. Panels (a) and (b) share the same legend. The gray marked areas highlight the Sentinel-2 bands B2, B4, B8, B8a, and B11.

[Figure omitted. See PDF]

The response of $\mathit{\rho }\left(\mathit{\lambda }\right)$ with increasing $\mathit{\tau }$ results from the transition from only direct to only diffuse radiation, which is controlled by the combination of $\mathit{\theta }$ and $\mathit{\tau }$. Two theoretical extremes are distinguished: (i) only direct radiation, also called “black-sky”, and (ii) only diffuse radiation, also called “white-sky”. Natural conditions typically lie between these two extremes; they are referred to as “blue-sky” . The amount of direct radiation controls two effects. The first effect acts on the canopy level, where diffuse radiation can penetrate deeper into the canopy and interacts with leaves that would be shaded under black-sky conditions. Consequently, for the same LAD and LAI, a greater total leaf area interacts with the incoming radiation . Model simulations by showed an increase in broadband solar canopy albedo with increasing $f_{\mathrm{dir}}$ . Furthermore, the extinction of radiation is sensitive to the incident angle, which is equal to $\mathit{\theta }$ for direct radiation and approaches an effective value of about 60$\mathit{°}$ under overcast conditions, due to the increasing contribution from diffuse radiation . The second factor acts on the individual leaf level, where the inherent directional reflection of radiation on surfaces is relevant. In the case of non-isotropic surfaces, direct radiation is scattered in all directions with a significant portion of the radiation being scattered in a specific solid angle, while diffuse radiation is scattered almost equally over the entire hemisphere. . Both effects are directly linked to the directional reflectance under blue-sky conditions that is described by the hemispherical–directional reflectance factor HDRF, , which includes the influence of the diffuse component in $F^{\downarrow }\left(\mathit{\lambda }\right)$. An example of a HDRF for a spherical LAD is given in Appendix .

3.1.2 Effect of cloud changes on spectral reflectance above vegetation

The effects presented above, do not include changes in cloud conditions between RP calibration and actual measurement over vegetated surfaces, where ${\mathit{\tau }}_{\mathrm{cal}}\ne {\mathit{\tau }}_{\mathrm{mea}}$ and $F_{\mathrm{cal}}^{\downarrow }\left(\mathit{\lambda }\right)\ne F_{\mathrm{mea}}^{\downarrow }\left(\mathit{\lambda }\right)$. The effects of cloud changes between RP measurements to $F_{\mathrm{mea}}^{\downarrow }\left(\mathit{\lambda }\right)$ are quantified and shown in Fig. , providing the illumination ratio $F_{\mathrm{mea}}^{\downarrow }\left(\mathit{\lambda }\right)/F_{\mathrm{cal}}^{\downarrow }\left(\mathit{\lambda }\right)$. The illumination ratio is calculated between the irradiance $F_{\mathrm{cal}}^{\downarrow }\left(\mathit{\lambda }\right)$ that prevailed during the RP calibration and the irradiance $F_{\mathrm{mea}}^{\downarrow }\left(\mathit{\lambda }\right)$ that is present during the above-vegetation measurements. The given illumination ratio represents the conditions below an ice cloud with ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{1}$ and three values of ${\mathit{\tau }}_{\mathrm{mea}}$ of 0, 4, and 10. In general, measurements where ${\mathit{\tau }}_{\mathrm{cal}}<{\mathit{\tau }}_{\mathrm{mea}}$ yield an illumination ratio less than one, because the reflectance at cloud top and the absorption inside the cloud become more intense, reducing the total $F^{\downarrow }\left(\mathit{\lambda }\right)$ below the cloud. For opposite conditions, where ${\mathit{\tau }}_{\mathrm{cal}}>{\mathit{\tau }}_{\mathrm{mea}}$, the illumination ratio is greater than one, since scattering and absorption are less during actual measurement compared to the RP calibration. Clearly visible are the ice absorption features, for example at 1500 nm wavelength, and a generally larger sensitivity of the illumination ratio towards longer wavelengths. Similar illumination ratios are calculated for liquid water clouds but with lower sensitivity for same values of $\mathit{\tau }$. This is due to the smaller cloud droplet size of 10 $\mathrm{\mu }$m compared to the ice particle size used in the simulations and the difference in the single-scattering phase function between liquid water droplets and ice crystals. In addition, differences in the imaginary part of the refractive indices lead to a slight spectral shift in the absorption features .

Figure 3

Ratio of downward irradiance $F_{\mathrm{mea}}^{\downarrow }\left(\mathit{\lambda }\right)$ present during the synthetic measurements and downward irradiance $F_{\mathrm{cal}}^{\downarrow }\left(\mathit{\lambda }\right)$ present during the calibration. The ratio is calculated for ice clouds, where ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{1}$ and ${\mathit{\tau }}_{\mathrm{mea}}$ is set to 0, 4, and 10. As an example, a solar zenith angle $\mathit{\theta }$ of 25$\mathit{°}$ is selected. The gray marked areas highlight the Sentinel-2 bands B2, B4, B8, B8a, and B11.

[Figure omitted. See PDF]

The spectral distortion in $F^{\downarrow }\left(\mathit{\lambda }\right)$ due to cloud changes has an immediate effect on $\mathit{\rho }\left(\mathit{\lambda }\right)$ that is calculated with Eq. (). For illustration, Fig.  shows $\mathit{\rho }\left(\mathit{\lambda }\right)$ of a vegetated surface for an intermediate value $\mathit{\theta }$ of 45$\mathit{°}$ and for four different combinations of ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$. The ground truth $\mathit{\rho }\left(\mathit{\lambda }\right)$, as it would be obtained from satellites, with ${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}=\mathrm{0}$, is given by the black line. Only in the trivial case, when the airborne observations are performed under the same cloud-free conditions, the same value of $\mathit{\rho }\left(\mathit{\lambda }\right)$ would be measured. For reference, $\mathit{\rho }\left(\mathit{\lambda }\right)$ under constant cloud conditions during RP calibration and measurement, where ${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}\ne \mathrm{0}$, is indicated by the gray line.

Figure 4

Panels (a, c): Synthetic spectral reflectance $\mathit{\rho }\left(\mathit{\lambda }\right)$ measurements, when cloud conditions change from ${\mathit{\tau }}_{\mathrm{cal}}$ = 0 to ${\mathit{\tau }}_{\mathrm{mea}}=$ 1 and 4, respectively. Panels (b) and (d): Same as (a) and (b) but for ${\mathit{\tau }}_{\mathrm{cal}}$ = 4 and values of ${\mathit{\tau }}_{\mathrm{mea}}$ of 1 and 4, respectively. In all panels: The black line represents $\mathit{\rho }\left(\mathit{\lambda }\right)$ under cloud-free conditions (${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}=\mathrm{0}$) and the gray line represents $\mathit{\rho }\left(\mathit{\lambda }\right)$ under constant cloudy conditions (${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}\ne \mathrm{0}$). Red lines represent ice water clouds and blue lines represent liquid water clouds. The gray marked areas highlight the Sentinel-2 bands B2, B4, B8, B8a, and B11.

[Figure omitted. See PDF]

The cases shown in Fig. a and c represent cloud-free conditions during the RP calibration (${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{0}$), while clouds are present during the actual measurements with ${\mathit{\tau }}_{\mathrm{mea}}=\mathrm{1}$ and ${\mathit{\tau }}_{\mathrm{mea}}=\mathrm{4}$, respectively. Due to the cloud change with ${\mathit{\tau }}_{\mathrm{cal}}$ $>$ ${\mathit{\tau }}_{\mathrm{mea}}$, the amused $\mathit{\rho }\left(\mathit{\lambda }\right)$ is lower compared to the reference. The difference between measured $\mathit{\rho }\left(\mathit{\lambda }\right)$ and the reference increases with increasing difference $\mathrm{\Delta }\mathit{\tau }={\mathit{\tau }}_{\mathrm{mea}}-{\mathit{\tau }}_{\mathrm{cal}}$ and is more pronounced for ice clouds than for liquid water clouds of the same $\mathit{\tau }$. For the ice cloud, $\mathit{\rho }\left(\mathit{\lambda }\right)$ at about 842 nm wavelength (Sentinel-2 B8) is reduced from about 0.38 (${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}=\mathrm{0}$) to about 0.33 (${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{0},{\mathit{\tau }}_{\mathrm{mea}}=\mathrm{1}$) and 0.26 (${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{0},{\mathit{\tau }}_{\mathrm{mea}}=\mathrm{4}$). Similarly, the opposite situation is possible, where $\mathit{\tau }$ decreases after calibration (${\mathit{\tau }}_{\mathrm{cal}}$ $>$ ${\mathit{\tau }}_{\mathrm{mea}}$), which is shown in the right column in Fig. . Since ${\mathit{\tau }}_{\mathrm{cal}}>{\mathit{\tau }}_{\mathrm{mea}}$, the measured $\mathit{\rho }\left(\mathit{\lambda }\right)$ overestimates the expected ground truth $\mathit{\rho }\left(\mathit{\lambda }\right)$, with the greatest bias in $\mathit{\rho }\left(\mathit{\lambda }\right)$ at about 1610 nm wavelength (Sentinel-2 B10) in the case of an ice cloud.

The examples given in Fig. b and d show that the estimated $\mathit{\rho }\left(\mathit{\lambda }\right)$ under constant cloud conditions (gray lines) are slightly enhanced compared to $\mathit{\rho }\left(\mathit{\lambda }\right)$ under cloud-free conditions but the change in $\mathit{\rho }\left(\mathit{\lambda }\right)$ due to differences between ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$ are much greater and dominating. Therefore, the subsequent analysis primarily focuses on the contribution of changing cloud conditions.

3.2 Effect of clouds and cloud changes on two-band vegetation indices

Figure  showed that the spectral distortion due to differences in ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$ affects certain wavelengths stronger than others. Subsequently, all wavelength combinations of $\mathit{\rho }\left({\mathit{\lambda }}_{\mathrm{1}}\right)$ and $\mathit{\rho }\left({\mathit{\lambda }}_{\mathrm{2}}\right)$, with ${\mathit{\lambda }}_{\mathrm{1}},{\mathit{\lambda }}_{\mathrm{2}}\in [\mathrm{400},\mathrm{2400}]$, that might be used in two-band VIs following Eq. () are examined with respect to their sensitivity to cloud changes. A similar approach was used, for example, by to determine the effect of cirrus clouds on upward radiances used for cloud remote sensing.

Figure 5

Absolute difference $\mathrm{\Delta }\mathit{\gamma }={\mathit{\gamma }}_{\mathrm{mea}}-\mathit{\gamma }$ for change in illumination conditions, when cloud optical thickness is different during the reflectance panel measurement (${\mathit{\tau }}_{\mathrm{cal}}$) and actual condition during the measurement (${\mathit{\tau }}_{\mathrm{mea}}$). Absolute differences in $\mathit{\gamma }$ are given for combinations of ${\mathit{\lambda }}_{\mathrm{1}}$ ($x$-axis) and ${\mathit{\lambda }}_{\mathrm{2}}$ ($y$-axis) simulated for an ice cloud, and for solar zenith angles $\mathit{\theta }$ of 25$\mathit{°}$ (top) and 50$\mathit{°}$ (bottom). Left column: Absolute differences, when ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{2}$ and ${\mathit{\tau }}_{\mathrm{mea}}=\mathrm{0}$. Right column: Absolute differences, when ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{2}$ and ${\mathit{\tau }}_{\mathrm{mea}}=\mathrm{4}$. Red areas show an overestimation of $\mathit{\gamma }$, while blue areas indicate an underestimation of $\mathit{\gamma }$. Wavelength combinations of $\mathit{\gamma }$ affected by strong water vapor absorption between 1350 and 1400 nm as well as 1830 and 1920 nm have been masked.

[Figure omitted. See PDF]

Figure  shows the effect of cloud changes on $\mathit{\gamma }$ expressed as $\mathrm{\Delta }\mathit{\gamma }={\mathit{\gamma }}_{\mathrm{mea}}-\mathit{\gamma }$, with ${\mathit{\gamma }}_{\mathrm{mea}}$ obtained over vegetation using the latest RP calibration that was performed under ${\mathit{\tau }}_{\mathrm{cal}}$. The expected ground truth $\mathit{\gamma }$ is the value that would be expected if an immediate estimate of $\mathit{\kappa }$ would be available or cloud conditions would not have changed (${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}\ne \mathrm{0}$). A negative $\mathrm{\Delta }\mathit{\gamma }$ therefore indicates an underestimation of the true $\mathit{\gamma }$ and vice versa. The combinations of ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$ were chosen to represent optically thin cirrus with changes in $\mathit{\tau }$ that could occur between RP calibration measurements, e.g., 10 min apart. Each panel in Fig.  is divided into four quadrants, starting with the first quadrant (Q1) in the lower right part and turning counter-clockwise. It follows from Eq. () that $\mathrm{\Delta }\mathit{\gamma }$ is point-symmetric with respect to the diagonal from the origin to the upper right corner. Thus, for any combination of ${\mathit{\lambda }}_{\mathrm{1}}$ and ${\mathit{\lambda }}_{\mathrm{2}}$ that causes an overestimation or underestimation, the inverse wavelength combination ${\mathit{\lambda }}_{\mathrm{2}}$ and ${\mathit{\lambda }}_{\mathrm{1}}$ will cause the same bias $\mathrm{\Delta }\mathit{\gamma }$, but with opposite sign. In addition to NDVI and NDWI, selected two-band VIs that follow Eq. () and use SWIR wavelengths were added as examples to the plot, namely: the modified normalized difference water index MNDVI; , the modified normalized difference water index MNDWI; , and the normalized difference infrared index NDII; . The list is not comprehensive and many more VIs exist, e.g., tabulated by , , and . Wavelengths affected by strong water vapor absorption between 1350 and 1400 nm as well as 1830 and 1920 nm wavelengths have been masked because $F^{\downarrow }\left(\mathit{\lambda }\right)$ and $F^{\uparrow }\left(\mathit{\lambda }\right)$ are close to zero causing numerical unstable values of $\mathit{\rho }\left(\mathit{\lambda }\right)$.

The largest absolute values of $\left|\mathrm{\Delta }\mathit{\gamma }\right|$, indicated by dark red or dark blue colors, occur for wavelength combinations with the greatest spectral distance $\left|\mathrm{\Delta }\mathit{\lambda }\right|=|{\mathit{\lambda }}_{\mathrm{1}}-{\mathit{\lambda }}_{\mathrm{2}}|$. Also wavelength combinations with one wavelength greater than 1400 nm, i.e., in Q1 and Q3, are affected by cloud transitions, since the spectral slope in $F^{\downarrow }\left(\mathit{\lambda }\right)$ is most pronounced towards longer wavelengths (see Fig. ).

In contrast, independent of the selected combination of ${\mathit{\tau }}_{\mathrm{cal}}$, ${\mathit{\tau }}_{\mathrm{mea}}$, and $\mathit{\theta }$ generally low values of $\left|\mathrm{\Delta }\mathit{\gamma }\right|$ are found for wavelength combinations with $\mathit{\lambda }<\mathrm{1400}$ nm (Q4). Additionally, small values of $\mathrm{\Delta }\mathit{\gamma }$ also occur along the diagonal from the origin to the upper right corner, i.e., the smaller $\left|\mathrm{\Delta }\mathit{\lambda }\right|$ becomes. While small $\left|\mathrm{\Delta }\mathit{\lambda }\right|$ would minimize the cloud effect, the proximity of ${\mathit{\lambda }}_{\mathrm{1}}$ and ${\mathit{\lambda }}_{\mathrm{2}}$ limits the information content that can be extracted from spectral ratios. A trade-off between information content and small cloud influence is required and Fig.  provides guidance to choose suitable wavelength combinations.

Figure a illustrates $\mathrm{\Delta }\mathit{\gamma }$ for $\mathit{\theta }=\mathrm{25}\mathit{°}$, where ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{2}$ decreases to ${\mathit{\tau }}_{\mathrm{mea}}=\mathrm{0}$. The decrease in $\mathit{\tau }$ after the RP calibration leads to predominantly negative $\mathrm{\Delta }\mathit{\gamma }$ (blue colors) in Q3 and positive $\mathrm{\Delta }\mathit{\gamma }$ (red colors) in Q2, indicating a respective underestimation and overestimation of the true $\mathit{\gamma }$. In the selected case, $\mathrm{\Delta }\mathit{\gamma }$ of NDVI (black circle) is small compared to the full potential range of NDVI (see Table ). Greater values of $\mathrm{\Delta }\mathit{\gamma }$ are calculated for ${\mathrm{NDWI}}_{\mathrm{1640}}$ with up to 0.184, resulting from the second wavelength being located in the SWIR region that is subject to the spectral slope in $F^{\downarrow }\left(\mathit{\lambda }\right)$.

Table 3

Difference $\mathrm{\Delta }\mathit{\gamma }={\mathit{\gamma }}_{\mathrm{mea}}-\mathit{\gamma }$ between estimated vegetation index ${\mathit{\gamma }}_{\mathrm{mea}}$ and ground truth vegetation index $\mathit{\gamma }$ for four vegetation indices for the four example cases given in Fig. a–d. The calculations assume a spherical LAD. Sentinel-2 band ratios were taken from .

Vegetation index	Sentinel-2 bands ratios	$\mathrm{\Delta }\mathit{\gamma }$
		Case (a)	Case (b)	Case (c)	Case (d)
NDVI	(B8 $-$ B4)$/$(B8 $+$ B4)	0.043	$-$0.011	0.009	$-$ 0.007
NDWI1240	(B3 $-$ B8)$/$(B3 $+$ B8)	$-$ 0.048	0.018	$-$ 0.022	0.015
NDWI1640	(B8a $-$ B11)$/$(B8a $+$ B11)	0.184	$-$0.212	0.216	$-$ 0.267
MNDVI	(B8 $-$ B12)$/$(B8 $+$ B12)	0.111	$-$0.127	0.112	$-$ 0.162
MNDWI	(B3 $-$ B11)$/$(B3 $+$ B11)	0.08	$-$0.113	0.156	$-$ 0.148
NDII	(B8 $-$ B11)$/$(B8 $+$ B11)	0.173	$-$0.205	0.208	$-$ 0.264

The second example in Fig. b shows the opposite transition from an optically thin cloud during the RP calibration with ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{2}$ to an optically thicker cloud during measurement ${\mathit{\tau }}_{\mathrm{mea}}=\mathrm{4}$. This leads to an inverted pattern of $\mathrm{\Delta }\mathit{\gamma }$. Although the change in cloudiness $\left|\mathrm{\Delta }\mathit{\tau }\right|={\mathit{\tau }}_{\mathrm{cal}}-{\mathit{\tau }}_{\mathrm{mea}}$ is similar to the example in Fig. a, the magnitude of $\left|\mathrm{\Delta }\mathit{\gamma }\right|$ for all wavelength combinations in Q2 and Q4 are slightly smaller, while the effect is greater for wavelength combinations in Q1 and Q3. This shows that the bias $\left|\mathrm{\Delta }\mathit{\gamma }\right|$ from changing cloud conditions is affected by $\mathrm{\Delta }\mathit{\tau }$ but also depends on the absolute value of ${\mathit{\tau }}_{\mathrm{cal}}$ during calibration. For this example, NDVI and ${\mathrm{NDWI}}_{\mathrm{1640}}$ are subject to biases $\mathrm{\Delta }\mathit{\gamma }$ of 0.011 and $-$0.212, respectively.

Figure c shows $\mathrm{\Delta }\mathit{\gamma }$ for the same combination of ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$ that is given in Fig. a but for a greater value of $\mathit{\theta }$ of 50$\mathit{°}$. This leads to greater absolute values of $\left|\mathrm{\Delta }\mathit{\gamma }\right|$ in Q1–Q3. An exception is Q4, which is characterized by $\mathrm{\Delta }\mathit{\gamma }$ around 0 for all wavelengths with less pronounced water vapor absorption features compared to $\mathit{\theta }$ of 25$\mathit{°}$. This indicates that with increasing $\mathit{\theta }$, changes in $\mathit{\tau }$ lead to reduced biases in $\mathit{\gamma }$ for wavelengths less than 1400 nm, while larger biases in $\mathit{\gamma }$ are expected for VIs using wavelengths beyond 1400 nm. In general, the effects of changing cloud conditions on $\mathit{\rho }\left(\mathit{\lambda }\right)$ and $\mathrm{\Delta }\mathit{\gamma }$ are less pronounced for the planophile LAD, while for the erectophile LAD greater effects were determined. This results primarily from the higher sensitivity of $\mathit{\rho }\left(\mathit{\lambda }\right)$ on $\mathit{\theta }$ in case of the erectophile LAD . An overview of $\mathrm{\Delta }\mathit{\gamma }$ for all example two-band VIs derived for the four cases that are marked in Fig.  are given in Table .

Subsequently, the effects of changing $\mathit{\tau }$ between RP calibration and measurement on three selected VIs are investigated.

3.2.1 Effect of cloud changes on the normalized differential vegetation index (NDVI)

The estimated $\mathrm{NDVI}$ for three values of $\mathit{\theta }$ depending on the combination of ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$ is shown in the top row of Fig. .

Figure 6

Top row: Absolute values of NDVI as obtained for solar zenith angles $\mathit{\theta }$ of 25, 50, and 70$\mathit{°}$. Subsequent rows same as top row but for ${\mathrm{NDWI}}_{\mathrm{1640}}$ and EVI. Color-coded is the cloud optical thickness ${\mathit{\tau }}_{\mathrm{cal}}$ that was present during the reflectance panel measurement. The cloud optical thickness ${\mathit{\tau }}_{\mathrm{mea}}$ during the measurement is given on the $x$-axis. Simulations for the liquid water cloud are given by solid lines and simulations for the ice cloud are given by the dotted lines. All simulations base on the assumption of a spherical LAD. The variability in VI due to the presence of clouds (${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}$) and the associated change in the surface reflectance, is highlighted in gray.

[Figure omitted. See PDF]

First, we consider the trivial case of cloud-free conditions with ${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}\approx \mathrm{0}$, where $\mathrm{NDVI}$ of about 0.87, 0.9, and 0.92 are calculated for $\mathit{\theta }$ of $\mathrm{25}$, $\mathrm{50}$, and $\mathrm{70}\mathit{°}$, respectively. These cases are marked by the dark-blue dots and represent the reference $\mathrm{NDVI}$. The increase of $\mathrm{NDVI}$ with increasing $\mathit{\theta }$ is related to scattering and absorption at gas molecules and aerosol particles. Next, cloudy conditions are considered, where ${\mathit{\tau }}_{\mathrm{cal}}\approx {\mathit{\tau }}_{\mathrm{mea}}\ne \mathrm{0}$ (colored dots), showing that the NDVI is enhanced by the presence of the cloud (gray highlighted area), which influences $f_{\mathrm{dir}}\left(\mathit{\lambda }\right)$ and $\mathit{\rho }\left(\mathit{\lambda }\right)$. The greatest variability is found for $\mathit{\theta }=\mathrm{25}\mathit{°}$, where $\mathrm{NDVI}$ increases from 0.87 to 0.91 with increasing $\mathit{\tau }$, which could be interpreted as an overestimation of vegetation health. Smaller effects are found for $\mathit{\theta }=\mathrm{50}\mathit{°}$ and for $\mathit{\theta }=\mathrm{70}\mathit{°}$ the NDVI even decreases. It should be noted that measuring NDVI under cloud-free conditions at different times of day, i.e., different $\mathit{\theta }$, causes the same variability as measuring NDVI at a fixed time with $\mathit{\theta }=\mathrm{25}\mathit{°}$ but under different cloud conditions, for example, on consecutive days.

For cases where cloud conditions differ between ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$, the change of ${\mathit{\tau }}_{\mathrm{mea}}$ for fixed ${\mathit{\tau }}_{\mathrm{cal}}$ is given by the colored lines. Moving to the right along lines of same ${\mathit{\tau }}_{\mathrm{cal}}$ can be understood as the advection of an optically thicker cloud during the measurement, while moving left represents the advection of an optically thinner cloud. The largest differences between measured and reference $\mathrm{NDVI}$ generally occur for $\mathit{\theta }=\mathrm{25}\mathit{°}$, and become successively smaller with increasing $\mathit{\theta }$. Thus, the NDVI is subject to the largest biases from cloud transitions when $\mathit{\theta }$ is small. The advection of an optically thicker cloud after the RP measurement (${\mathit{\tau }}_{\mathrm{cal}}<{\mathit{\tau }}_{\mathrm{mea}}$) results in a decrease in $\mathrm{NDVI}$ and an underestimation of the expected value. For ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{0}$ an extreme increase of ${\mathit{\tau }}_{\mathrm{mea}}$ from 0 to 40 results in a decrease in $\mathrm{NDVI}$ from 0.87 to 0.84, which could be interpreted as an underestimation of vegetation health. The advection of an optically thinner cloud after the RP measurement results in an overestimation of the actual $\mathrm{NDVI}$ value. For example, the combination of ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{40}$ and ${\mathit{\tau }}_{\mathrm{mea}}=\mathrm{0}$ increases the $\mathrm{NDVI}$ from 0.92 to 0.94. The aforementioned extreme changes in $\mathit{\tau }$ from 0 to 40, and vice versa, between RP overflights are unlikely in the case of stratiform clouds. These extremes of $\mathit{\tau }$ have been selected to estimate the maximum envelope for biases in $\mathrm{NDVI}$. Biases in $\mathrm{NDVI}$, retrieved under more homogeneous conditions, will be smaller. However, even seemingly stratiform clouds vary in $\mathit{\tau }$ and therefore bias the retrieved values of $\mathrm{NDVI}$. These variations have the greatest impact, when the RP overflight is performed under cloud-free conditions (${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{0}$) or under conditions with small values of ${\mathit{\tau }}_{\mathrm{cal}}$, since the slope of the curves is greatest in these situations (see, for example, the blue lines in Fig. a–c).

Comparing the responses of $\mathrm{NDVI}$ obtained below liquid and ice water clouds shows that both cloud types cause similar biases. However, the magnitude is generally greater for the ice cloud with the largest difference for small values of $\mathit{\theta }$ and for ${\mathit{\tau }}_{\mathrm{mea}}$ between 10 and 30. However, it is acknowledged that the calculated deviations between the measured and expected $\mathrm{NDVI}$ are small, considering the selected differences between ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$, and the typical range of $\mathrm{NDVI}$ between 0 and 1.

3.2.2 Effect of cloud changes on the normalized differential water index (NDWI)

Remote sensing of ${\mathrm{NDWI}}_{\mathrm{1640}}$ is based on the Sentinel-2 bands B8a and B11, a combination that is more sensitive to changes in $\mathit{\tau }$ compared to the NDVI, because one wavelength is located in the SWIR (Q3 in Fig. ). The middle row in Fig.  shows the change of NDWI induced by cloud changes.

Regardless of $\mathit{\theta }$, the bias in ${\mathrm{NDWI}}_{\mathrm{1640}}$ is under cloudy conditions with ${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}$ is of similar magnitude compared to the NDVI. A variation in ${\mathrm{NDWI}}_{\mathrm{1640}}$ with maximal values of $\pm$0.02 (gray highlighted area) is calculated, which is small compared to the influence of changes in cloudiness between RP calibrations. Here, the change in cloud conditions is discussed for $\mathit{\theta }=\mathrm{25}\mathit{°}$. When the RP measurements was performed for ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{0}$ the advection of an optically thicker cloud, represented by an increase in ${\mathit{\tau }}_{\mathrm{cal}}$ from 0 to 40, leads to an increase in ${\mathrm{NDWI}}_{\mathrm{1640}}$ from 0.28 to 0.89 ($\mathrm{\Delta }\mathit{\gamma }=\mathrm{0.61}$), which is relevant considering the typical range between $-$1 and 1. In such scenarios, the interpreted true health status would be well overestimated using ${\mathrm{NDWI}}_{\mathrm{1640}}$. For the same $\mathit{\theta }$, a calibration under cloudy sky with ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{40}$ and a subsequent decrease in ${\mathit{\tau }}_{\mathrm{mea}}$ from 40 to 0 results in a decrease in ${\mathrm{NDWI}}_{\mathrm{1640}}$ from 0.59 to $-$0.43 ($\mathrm{\Delta }\mathit{\gamma }=-\mathrm{1.02}$). Even though $\mathrm{\Delta }\mathit{\tau }$ is equal in both cases, the resulting $\mathrm{\Delta }\mathit{\gamma }$ differs, which emphasizes the dependence of the bias $\mathrm{\Delta }\mathit{\gamma }$ on the absolute ${\mathit{\tau }}_{\mathrm{cal}}$. As for the $\mathrm{NDVI}$, the transitions in $\mathit{\tau }$ from 0 to 40, or vice versa, are extreme cases. For practical applications, variations in ${\mathit{\tau }}_{\mathrm{mea}}$ for a given ${\mathit{\tau }}_{\mathrm{cal}}$ are more relevant. For the ${\mathrm{NDWI}}_{\mathrm{1640}}$, the slopes in Fig. d–f are shifted in absolute terms but are constant for the same values of $\mathit{\theta }$. Therefore, variations in ${\mathit{\tau }}_{\mathrm{mea}}$ between RP overflights will cause the same bias in ${\mathrm{NDWI}}_{\mathrm{1640}}$, independent of ${\mathit{\tau }}_{\mathrm{cal}}$. Similar to $\mathrm{NDVI}$, the effect of ice clouds on ${\mathrm{NDWI}}_{\mathrm{1640}}$ is qualitatively similar to the effect from liquid water clouds, however the magnitude is almost twice as large.

3.3 Effect of clouds and cloud changes on the three-band enhanced vegetation index (EVI)

As an example for a three-band VI, the effects of changing cloud conditions between RP calibration and measurement are shown for the $\mathrm{EVI}$ (see Fig. g–i). The effect of changes in $f_{\mathrm{dir}}\left(\mathit{\lambda }\right)$ under constant conditions (${\mathit{\tau }}_{\mathrm{cal}}={\mathit{\tau }}_{\mathrm{mea}}$) is generally small (gray area) compared to the variation in $\mathrm{EVI}$ that is associated with the mismatch between ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$. Irrespective of $\mathit{\theta }$, the advection of an optically thicker cloud (${\mathit{\tau }}_{\mathrm{cal}}<{\mathit{\tau }}_{\mathrm{mea}}$) causes an increase in estimated $\mathrm{EVI}$, which leads to an overestimation of the true $\mathrm{EVI}$. Conversely, the advection of an optically thinner cloud (${\mathit{\tau }}_{\mathrm{cal}}>{\mathit{\tau }}_{\mathrm{mea}}$) leads to and underestimation of the true $\mathrm{EVI}$. The bias gets more pronounced with increasing difference between ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$, and with increasing absolute value of ${\mathit{\tau }}_{\mathrm{cal}}$. The magnitude of these biases are more pronounced for the ice cloud than for the liquid water cloud. Furthermore, the biases become larger with increasing $\mathit{\theta }$, as well as with increasing difference between ${\mathit{\tau }}_{\mathrm{cal}}$ and ${\mathit{\tau }}_{\mathrm{mea}}$, which is particularly pronounced in combination with small ${\mathit{\tau }}_{\mathrm{cal}}$. Thus, estimated EVI are less susceptible to changes in cloud conditions at low values of $\mathit{\theta }$, e.g., around noon or at low latitudes, compared to cloud changes during measurements with the Sun close to the horizon or generally at higher latitudes. As an example, a calibration performed under $\mathit{\theta }=\mathrm{25}\mathit{°}$ and cloud-free conditions followed by the advection of an optically thicker ice cloud with ${\mathit{\tau }}_{\mathrm{mea}}=\mathrm{4}$ results in an increase in $\mathrm{EVI}$ from the expected value of 0.67 to 0.75. For the same transition in $\mathit{\tau }$ but for $\mathit{\theta }=\mathrm{70}\mathit{°}$ an increase from 0.67 to 1.04 is estimated. In both cases, the inferred ground-truth vegetation health would be overestimated.

The different response of $\mathrm{EVI}$ to cloud changes compared to the $\mathrm{NDVI}$ is related to two factors. First, the equation to calculate the EVI is fundamentally different from the two-band VIs. The use of an additional third spectral band at 490 nm wavelength makes the EVI generally more responsive to the spectral slope in $F^{\downarrow }\left(\mathit{\lambda }\right)$ caused by absorption. Second, the spectral slope effect is amplified by the pre-factors in Eq. (), since re-writing Eq. () results in: 10$$\mathrm{EVI}=G\cdot {\displaystyle \frac{a-b}{a+C_{\mathrm{1}}\cdot b-C_{\mathrm{2}}+L/{\mathit{\rho }}_{\mathrm{B}\mathrm{2}}}},$$ with the reflectance ratios $a={\mathit{\rho }}_{\mathrm{B}\mathrm{8}}/{\mathit{\rho }}_{\mathrm{B}\mathrm{2}}$ and $b={\mathit{\rho }}_{\mathrm{B}\mathrm{4}}/{\mathit{\rho }}_{\mathrm{B}\mathrm{2}}$ at Sentinel-2 bands B2, B4, and B8 (see Table ). The ratios $a$ and $b$ are affected by the spectral slope in the illumination ratio (see Fig. ). The spectral slope is greater for $a$ than for $b$, since $\mathrm{\Delta }\mathit{\lambda }$ is greater in $a$ than in $b$. The differences between $a$ and $b$ do not cancel out and are amplified by the pre-factor $C_{\mathrm{1}}$. The term $L/{\mathit{\rho }}_{\mathrm{B}\mathrm{2}}$ also contributes, since the constant $L$ is inversely scaled with ${\mathit{\rho }}_{\mathrm{B}\mathrm{2}}$ and band B2 is sensitive to atmospheric scattering.

3.4 Implication of biases in vegetation indices on estimated biophysical properties

In vegetation remote sensing, VIs are used to estimate biophysical properties such as LAI, gross primary production (GPP), fresh and dry biomass, and vegetation water content . For example, proposed an empirical linear regression given by: 11$$\mathrm{LAI}=\mathrm{3.618}\cdot \mathrm{EVI}-\mathrm{0.118},$$ where the LAI scales linearly with the EVI. Due to the direct relationship between $\mathrm{EVI}$ and LAI, biases in $\mathrm{EVI}$ are linearly scaled by the pre-factor 3.618 and lead to biased estimates of LAI. For example, when a calibration for $\mathit{\theta }=\mathrm{50}\mathit{°}$ was performed under cloud conditions with ${\mathit{\tau }}_{\mathrm{cal}}=\mathrm{0}$ but the measurement is influenced by an ice cloud with ${\mathit{\tau }}_{\mathrm{mea}}=\mathrm{10}$, the EVI is biased by 0.25, which corresponds to a bias in LAI of about $\mathrm{\Delta }\mathrm{LAI}=\mathrm{0.25}\cdot \mathrm{3.618}=\mathrm{0.9}$.

Similarly, the gross primary product is an important measure in the context of the global carbon cycle, since it indicates how quickly an ecosystem accumulates biomass . Several attempts have been made to estimate GPP based on NDVI and EVI; for example by , , or . Correlations have been estimated using in-situ carbon uptake flux measurements and remotely sensed NDVI or EVI. Although the correlations between GPP and EVI or NDVI vary by site and crop type, stable correlations could be identified. Using the correlations from or , a variation in EVI of $\mathrm{\Delta }\mathrm{EVI}=\mathrm{0.2}$ would yield variations in estimated GPP of about 2.6 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. also provided relationships for NDVI, where a variation of $\mathrm{\Delta }\mathrm{NDVI}=\mathrm{0.2}$ would result in a variation in GPP of about 0.51 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. Considering that the total range of GPP spans from 0 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, when the NDVI or EVI is zero, to 8 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, when NDVI or EVI is close to one, the bias in estimated GPP from biases in EVI is of relevance.

Several studies have attempted to estimate vegetation water content (VegWC) based on NDVI . The relationships derived between NDVI and VegWC vary greatly depending on the year, location, and crop type. Therefore, the identified cloud-induced biases in NDVI are expected to play only a minor role in the total uncertainty of estimated VegWC. aimed to estimate fresh and dry biomass, and LAI based on derived NDVI and EVI. These relationships are strongly dependent on vegetation type. Therefore, cloud-induced variations in NDVI are expected to contribute only a small amount to the total uncertainty in estimated fresh and dry biomass. However, since EVI is more sensitive to cloud-induced biases than NDVI, the uncertainty in EVI can be of a similar order of magnitude compared to the uncertainty in the relationships.

4 Summary and conclusions

This paper presented results of coupled atmosphere-vegetation radiative transfer (RT) simulations, using the coupled RT models libRadtran and SCOPE2.0, to systematically investigate biases in remotely sensed vegetation indices (VIs) due to changing cloud conditions. Simulations were performed for a stratiform liquid water cloud and an ice cloud representative of high-level cirrus. The optical thickness $\mathit{\tau }$ of the cloud varied between 0 and 40, and the solar zenith angles $\mathit{\theta }$ ranged from 25$\mathit{°}$ to 70$\mathit{°}$. The canopy optical properties were represented by spherical, erectophile, and planophile leaf angle distributions (LADs) and a leaf area index (LAI) of 3. The simulations were designed to resemble below-cloud measurements of downward irradiance $F^{\downarrow }\left(\mathit{\lambda }\right)$ and spectral surface reflectance $\mathit{\rho }\left(\mathit{\lambda }\right)$ used for remote sensing of vegetation from airborne observations. The synthetic measurements mimic the typical observation strategy during field measurements, where reflectance panels (RPs) are used to calibrate the surface reflectance measurements.

Field measurements can be performed below clouds, which reduce the direct fraction $f_{\mathrm{dir}}\left(\mathit{\lambda }\right)$ in $F^{\downarrow }\left(\mathit{\lambda }\right)$, which is a controlling factor of the surface reflectance. Clouds also change the downward irradiance $F^{\downarrow }\left(\mathit{\lambda }\right)$ spectrally and in absolute terms through spectral dependent scattering and absorption. Changes in cloud conditions, expressed as cloud optical thickness $\mathit{\tau }$, between periodic RP calibrations are expected to cause spectral distortions in the estimated $\mathit{\rho }\left(\mathit{\lambda }\right)$ and thus introduce biases in VIs estimated in the presence of clouds.

The synthetic observations allowed to separate the effect from changes in $f_{\mathrm{dir}}\left(\mathit{\lambda }\right)$ from the changes in $\mathit{\tau }$ between RP calibrations. For a solar zenith angle $\mathit{\theta }=\mathrm{25}\mathit{°}$, a reduction in $f_{\mathrm{dir}}\left(\mathit{\lambda }\right)$ led to an increase in $\mathit{\rho }\left(\mathit{\lambda }\right)$ for wavelengths greater than 650 nm compared to cloud-free conditions, while for wavelengths below 650 nm $\mathit{\rho }\left(\mathit{\lambda }\right)$ was reduced. For $\mathit{\theta }>\mathrm{55}\mathit{°}$, $\mathit{\rho }\left(\mathit{\lambda }\right)$ below 650 nm wavelengths was also increased with decreasing $f_{\mathrm{dir}}\left(\mathit{\lambda }\right)$. The change in $\mathit{\rho }\left(\mathit{\lambda }\right)$ under cloudy conditions therefore depends on the combination of $\mathit{\theta }$ and $\mathit{\tau }$, and the resulting $f_{\mathrm{dir}}\left(\mathit{\lambda }\right)$. The effect of changes in $f_{\mathrm{dir}}\left(\mathit{\lambda }\right)$ on estimated VIs was found to be small compared to the effect of changes in $\mathit{\tau }$ between RP calibrations.

The influence of cloud changes on two-band VIs was investigated for all possible wavelength combinations between 400 and 2400 nm, particularly for the normalized difference vegetation index (NDVI) and the normalized water index (NDWI). The influence of cloud changes on the enhanced vegetation index (EVI), which is representative for three-band VIs, was also investigated. The effect of cloud changes on narrow-band VIs was found to be small for VIs using wavelength combinations below 1400 nm and for decreasing spectral distance between wavelength combinations. However, an optimum must be found between the information content obtained from the wavelength ratio and the sensitivity of the respective wavelengths pair to cloud changes. Guidance in selecting wavelength-pairs was provided by presenting the cloud-induced bias for potential two-band wavelength combinations between 400 and 2400 nm. For wavelengths greater than 1400 nm, the sensitivity to cloud changes was found to increase and is particularly pronounced in the proximity and within the water absorption bands. For the NDVI, a generally low sensitivity was found. For an intermediate value of $\mathit{\theta }=\mathrm{50}\mathit{°}$, the transition in $\mathit{\tau }$ from 0 to 10 led to a bias of about $-$0.035. With increasing $\mathit{\theta }$ and for same transition in $\mathit{\tau }$, the bias in NDVI increased, while a decrease in $\mathit{\theta }$ resulted in lower biases. The NDWI is subject to greater biases of up to 0.26 for an intermediate $\mathit{\theta }=\mathrm{50}\mathit{°}$ and a transition in $\mathit{\tau }$ from 0 to 10. The biases in NDWI generally increase with increasing $\mathit{\theta }$. The EVI was also found to be sensitive to changes in $\mathit{\tau }$. A transition in $\mathit{\tau }$ from 0 to 10 in combination with $\mathit{\theta }=\mathrm{50}\mathit{°}$ led to an overestimation of the EVI by 0.25. For the same transitions in $\mathit{\tau }$, the bias in EVI decreases with decreasing $\mathit{\theta }$. The LAI estimated from EVI using empirical equations is directly affected by potential biases in EVI. For example, biases in EVI of 0.25 would cause a bias in LAI of 0.9. Similarly, biases in EVI and NDVI of about $\pm$0.2 can lead to biases in estimated gross primary product of up to $\pm$2 $\mathrm{gC}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. The impact for other LADs, such as the erectophile and planophile, has been investigated. It was found that the erectophile LAD is more susceptible to biases in VIs under constant cloud conditions than the spherical LAD, while the bias in VIs with underlying planophile LAD is almost unaffected. Although the absolute value of the VI biases are determined by the LAD, the relative change due to a transition in $\mathit{\tau }$ between RP measurements was found to be constant among the three LADs.

The presented analysis showed that the practice of using relative measurements with RPs is prone to uncertainties. With the improvement of drone technology and their ability to carry heavier payloads, and in combination with the advancements in sensor technology, it would be advantageous to measure $F^{\downarrow }\left(\mathit{\lambda }\right)$ directly instead of relying on relative measurements of $\mathit{\rho }\left(\mathit{\lambda }\right)$ using RPs. Uncertainties associated with changing cloud conditions could be minimized.

It is emphasized that the simulations here are limited in their representation because the full natural variability in vegetation and canopy types could not be covered. Furthermore, the assumed default values used in the atmosphere and vegetation RT simulations can influence the results. In particular, the assumed LAI, LAD, plant dry matter, and soil properties in the vegetation RT simulations, and the assumed atmospheric profile in the atmosphere RT simulations can influence the results . Since the natural variability cannot be covered in a single study, the presented work is to be interpreted as a conceptual study and to highlight the potential impacts of clouds during field observations. It is also emphasized that the presented simulations are based on one-dimensional RT only and lack a more detailed representation of the three-dimensional nature of RT below heterogeneous cloud fields. This is particularly problematic in the vicinity of clouds, where a nearby cloud casts a shadow while no cloud is in the zenith. Nevertheless, we argue that the presented study can be used as a first approximation for the transition between cloud-free and cloudy regions.

Appendix A Effect of cloud changes on the NDVI, NDWI1640, and EVI for erectophile and planophile LADs

This section provides an overview of how an erectophile or planophile leaf angle distribution (LAD) affects biases in normalized differential vegetation index (NDVI), normalized difference water index (NDWI1640), and the enhanced vegetation index (EVI), compared to the spherical LAD presented in the main text. First, we consider constant cloud optical thickness during the reflectance panel (RP) overflight with τcal and the actual measurement with τmea, such that τcal=τmea. Comparing Figs.  and reveals a greater variation in the NDVI bias for the erectophile LAD than for the spherical LAD (gray highlighted area). The biases are greatest for small solar zenith angle θ, when incoming radiation can penetrate deep into the canopy. In the case of θ=25°, the NDVI bias varies between 0.67 and 0.89 for the erectophile LAD, compared to a variation between 0.75 an 0.91 in case of the spherical LAD. In this particular case, the variation in NDVI exceeds the bias that results from a transition in cloud conditions, when τcal≠τmea. However, when θ is greater than 25° and the penetration depth of the radiation into the canopy is reduced, the bias from changing cloud conditions between RP overflights again dominates the NDVI bias. Comparing Figs.  and shows that for the planophile LAD and τcal=τmea, the variation in NDVI is negligible, with a nearly constant value of around 0.89. Similar effects are found for NDWI1640 and EVI. In general, the greatest biases occur for the erectophile LAD, independent of θ. Within on LAD the greatest biases occur for the lowest values of θ.

Taking into account changes in cloud conditions, i.e., when τcal≠τmea, the comparison of Figs. , , and shows that for the same values of θ, the relative variations in VI due to cloud changes are the same for all three LADs and VIs, which is indicated by a constant sensitivity of the VIs on τmea independent of the value of τcal. However, the absolute values in the VIs are shifted along the y-axis depending on τcal and the selected LAD.

Figure A1

Same as Fig.  but for the erectophile leaf angle distribution.

[Figure omitted. See PDF]

Figure A2

Same as Fig.  but for the planophile leaf angle distribution.

[Figure omitted. See PDF]

Appendix B Simulated reflectance factors and reflectance functions

The intrinsic reflectivity properties of a surface are given by its bidirectional reflectance distribution function (BRDF). The BRDF quantifies the reflection and scattering of incident radiation on the surface from one direction of the hemisphere to another. The spectral BRDF fBRDF, in units of sr⁻¹, gives the ratio between the reflected radiance Ir↑(θi,φi;θr,φr;λ) with respect to the incident irradiance Fi↓(θi,φi;λ). The resulting fBRDF depends on the zenith (θi) and the azimuth angle (φi) of the incoming radiation, and the zenith (θr) and azimuth angle (φr) of the reflected radiation. The spectral BRDF fBRDF also depends on the wavelength λ and is given by: B1fBRDF=dIr↑(θi,φi;θr,φr;λ)dFi↓(θi,φi;λ), given in the unit steradian (sr⁻¹).

Under atmospheric conditions, where F↓(λ) is composed of a direct and a diffuse fraction, the BRDF cannot be measured . What can be observed is the spectral reflectance ρ(λ) and the hemispherical-directional reflectance factor (HDRF), both of which include radiation from the entire hemisphere and not​​​​​​​​ ​​​​​​​​ just from a single direction. The hemispherical-directional reflectance factor is formalized by: B2RHDRF=dIr↑(θi,φi,2π;θr,φr)dF0↓(θi,φi,2π)=RHDRF(θi,φi;θr,φr)⋅fdir+R(2π;θr,φr)⋅(1-fdir). Thus, RHDRF represents the ratio of the reflected incident radiation from a surface in relation to an ideal, lambertian surface. Furthermore, in Eq. () it is assumed that the incident, diffuse radiation is isotropic .

Figure a and e show RHDRF at 550 nm wavelength for cloud-free conditions with τ=0 for θ of 25 and 70°, respectively. Scattering by gas molecules and aerosol particles lead to direct fraction fdir(λ)=0.94 just below 1. Under these conditions the first term in Eq. () dominates and the directional effect of the incoming radiation on the reflected radiation is most pronounced leading to the clear development of the hot spot that follows the solar zenith angle of the Sun (see white markers). With increasing values of τ and decreasing fdir(λ), the second term in Eq. () dominates since the radiation reaching the surface is predominantly diffuse and almost isotropic. The directional component of the surface reflectance vanishes, causing the hot spot to disappear and leading to a general smoothing of RHDRF, which is also confirmed in the cross-sections in panels (c), (d), (g), and (h).

Figure B1

All sub-panels show the normalized hemispherical–directional reflectance factor RHDRF at 550 nm wavelength. Panels (a), (b), (e), and (f) present polar plots of RHDRF for combinations of solar zenith angle θ and cloud optical thickness τ. Panels (c), (d), (g), and (h) present cross-sections of RHDRF along lines of 0–180 and 90–270° azimuth.

[Figure omitted. See PDF]

Data availability

The simulated spectra of radiance, irradiance, and vegetation albedo are made available via NetCDF files. The data are available on the Zenodo platform via 10.5281/zenodo.15275610 .

Author contributions

KW designed and implemented the model coupling, performed the simulations, and drafted the manuscript. EJ, AE, and MW, contributed to the preparation and revisions of the manuscript. MS, HF, and AH supported during the model set-up and provided suggestions for the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Acknowledgements

We thank the German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, which is a research center of the Deutsche Forschungsgemeinschaft (DFG). We would also like to acknowledge the two anonymous reviewers whose comments helped to improve the final manuscript.

Financial support

This research has been supported by the Sächsisches Staatsministerium für Wissenschaft und Kunst (grant no. 3-7304/44/4-2023/8846).

Review statement

This paper was edited by Cornelius Senf and reviewed by two anonymous referees.