Canopy modeling of aquatic vegetation: A radiative transfer approach

RSE-09360; No of Pages 20 Remote Sensing of Environment xxx (2015) xxx–xxx

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Canopy modeling of aquatic vegetation: A radiative transfer approach Guanhua Zhou, Chunyue Niu ⁎, Wujian Xu, Wenna Yang, Jiwen Wang, Huijie Zhao School of Instrument Science and Opto-electronics Engineering, Beihang University, Beijing 100191, China

a r t i c l e

i n f o

Article history: Received 29 July 2014 Received in revised form 20 March 2015 Accepted 21 March 2015 Available online xxxx Keywords: Aquatic vegetation Radiative transfer Optical shallow water Monte Carlo BRF

a b s t r a c t Aquatic vegetation in shallow waters can absorb nitrogen and phosphorus from eutrophic waters, and effectively control water bloom. However, our knowledge of the interaction between solar radiation and aquatic vegetation is limited. This restricts the protection of aquatic vegetation. This paper proposes a radiative transfer approach for homogeneous aquatic vegetation canopy, based on the theories for terrestrial vegetation and Case II water. The effects of rough water surface and bottom reflection are also taken into consideration. Validations against an independent Monte Carlo model and field experiment data exhibit well consistencies on the bidirectional reflectance of homogeneous aquatic vegetation canopy. As expected, the results for sparse canopy degenerate. This model can be used for canopy modeling of emergent, floating or submerged aquatic vegetation. Potential applications include retrieval of biophysical or biochemical parameters based on hyperspectral data and physical process. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Much of polluted water is poured into coastal and inland waters due to human activities. As a consequence, rapid increase of total nitrogen and phosphorus creates eutrophic waters. This in turn causes a massive reproduction of algae. Water bloom may appear more frequently and unexpectedly (An et al., 2013). Fortunately, aquatic vegetation can absorb nitrogen and phosphorus from these eutrophic waters, and effectively control water bloom (Abe, Komada, & Ookuma, 2008). Aquatic vegetation usually grows in the transition regions between land and water, typically a wetland ecosystem. However, these regions are generally difficult for humans to set foot due to the presence of marshes and mires. As a non-contact detection method, remote sensing has unique advantages in the monitoring and information extraction of aquatic vegetation. Various applications based on multispectral or hyperspectral data have been reported, e.g. classification (Schmidt & Skidmore, 2003), mapping (Gilmore et al., 2008) or monitoring (Barducci, Guzzi, Marcoionni, & Pippi, 2009) of aquatic vegetation, and retrieval of biomass (Proisy, Couteron, & Fromard, 2007) or leaf area index (LAI) (He, Quan, & Xing, 2013). Detailed reviews can be found in Adam, Mutanga, and Rugege (2010) and Klemas (2011, 2013). But the majority of the present methods belong to the empirical or semiempirical category and greatly rely on field survey and measurement. ⁎ Corresponding author at: School of Instrument Science and Opto-electronics Engineering, Beihang University, No. 37 Xueyuan Road, Haidian District, Beijing 100191, China. E-mail addresses: [email protected] (G. Zhou), [email protected] (C. Niu), [email protected] (H. Zhao).

The relationship between reflectance spectra and parameters (canopy structure parameters, etc.) is not described in detail. Further research on the radiative process in aquatic vegetation is also needed (Adam et al., 2010). The habitats of aquatic vegetation usually have a shallow water depth and the incident radiation can propagate to the bottom. This is the so-called optical shallow water (Ackleson, 2003). Based on the relative position of canopy and water, aquatic vegetation can be divided into three categories, namely emergent aquatic vegetation (EAV), floating aquatic vegetation (FAV) and submerged aquatic vegetation (SAV). Aquatic vegetation together with its habitat has the following one or more characteristics compared with Case I ocean water and terrestrial vegetation: (1) the simultaneous appearance of canopy and water results in a more complicated radiative transfer process, and the mixed spectra differ from those when either exists separately; (2) the water is relatively turbid, and the influences of water components (suspended particle, etc.) and bottom reflection on the spectra measured above water surface cannot be neglected (especially SAV vs. Case I water); (3) either part of or the entire canopy is immersed in water periodically or permanently, and the water surface has significant directional reflectance characteristics (especially EAV vs. terrestrial vegetation) (Vanderbilt et al., 2002); and (4) aquatic vegetation communities have high spatial complexity and temporal variability (Klemas, 2013), and their reflectance signals are easily affected by adjacent land, especially for airborne or spaceborne sensors. Thus the canopy modeling of aquatic vegetation needs to adhere to relevant theories, and grasp its own characteristics as well. The canopy modeling of aquatic vegetation cannot be completely beyond the scope of the modeling of vegetation canopy. Various

http://dx.doi.org/10.1016/j.rse.2015.03.015 0034-4257/© 2015 Elsevier Inc. All rights reserved.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

2

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

approaches have been summarized in Qin and Liang (2000), Disney, Lewis, and North (2000) and Chen, Li, Nilson, and Strahler (2000). Actually some early researches are based on the traditional vegetation canopy models. An example is the submerged vegetation canopy model by Suits (1984) (noted as Suits84) which is based on the Suits canopy model (Suits, 1972). The media below water surface was divided into two parts: turbid water (the upper layer) and mixed vegetation and water. Calculation of attenuation and scattering coefficients for mixed medium was given (Suits, 1984). Similarly, Ackleson and Klemas (1986) also built a radiative transfer model for submerged vegetation canopy by extending the Suits canopy model, namely Extended Suits Model (ESM). ESM uses matrix formulation and can be applied to multilayer media. However, the canopy morphology of Suits model is simplified to horizontal and vertical leaf area projections, and this leads to deviations of bidirectional reflectance profiles (Verhoef, 1984). Furthermore, only submerged canopy was discussed in Suits84 and ESM. Despite such deficiencies, these models lay a foundation for further research. Lee, Carder, Mobley, Steward, and Patch (1998) built a semi-analytical (SA) model for optical shallow water. The contribution of water column to remote sensing reflectance Rrs is separated from that of the bottom. This SA model was validated by the simulation results of Hydrolight (Mobley, 1994). However, bottom albedo is only a boundary condition in the SA model. No descriptions about the bidirectional reflectance distribution function (BRDF) of the bottom can be found. Mobley, Zhang, and Voss (2003) analyzed the effects of BRDF of non-Lambertian bottoms on Rrs, and errors were quantified if these bottoms were replaced by Lambertian ones. The Rahman model (Rahman, Pinty, & Verstraete, 1993) was used as an ersatz aquatic vegetation BRDF. However, as noted by the authors, BRDFs of aquatic vegetation, e.g. seagrass, may differ from those described by the Rahman model (Mobley et al., 2003). Zimmerman (2003) built a parameterized twoflow model of plane irradiance distribution for seagrass canopy. Impacts of canopy architecture on irradiance distribution were explored, and field measurements were used for validation. But this model pays more attention to irradiance distribution and photosynthesis rather than the BRDF of submerged canopy. Hedley (2008) proposed a threedimensional (3D) radiative transfer model for shallow water environments based on spatial and directional discretization. This model can solve radiative transfer problems in anisotropic scattering media, and achieve good performances in a plane parallel configuration. In a later paper (Hedley & Enríquez, 2010), this 3D model was applied to seagrass canopies. Their BRDFs were demonstrated to be non-Lambertian and exhibited different features with those used in Mobley et al. (2003). However, computer simulation models are usually difficult to retrieve despite their high accuracy and capability of handling complicated 3D problems. It is notable that most of the models mentioned above have a spectrum mainly within the visible bands (400–700 nm). Although this is a general configuration in water color remote sensing, it is insufficient for exploring aquatic vegetation. Based on field experiments of several emergent vegetation species, Kearney, Stutzer, Turpie, and Stevenson (2009) demonstrated significant reflectance reductions in near infrared bands (900–1100 nm) with progressive inundation. Turpie (2012) added a shallow water background (Lee et al., 1998) to a two-layer canopy reflectance model (ACRM) (Kuusk, 2001) and developed the Wetland Canopy Reflectance Model (WCRM). WCRM can calculate the reflectance spectra of aquatic vegetation in visible and near infrared (VNIR) bands (400–1000 nm). Simulation results indicated a shift in the red-edge position, but a simple relationship to the relative inundation level could not be confirmed (Turpie, 2013). Moreover, WCRM can simulate the bidirectional reflectance characteristics of emergent vegetation, the hot spot of canopy and the specular reflection of water surface, simultaneously (Turpie, 2012). However, water surface together with the media beneath is regarded as background in WCRM, and increasing inundation levels are simulated by decreasing the LAI values above water. Interactions between radiation and submerged canopy are not

discussed. Beget, Bettachini, Di Bella, and Baret (2013) built a radiative transfer model for flooded vegetation named SAILHFlood. Two vegetation layers, the emerged vegetation layer and the submerged vegetation layer, are included. However, a flat water surface is used in SAILHFlood, and the vegetation is assumed to be submerged into clear water without particles. Until now, seldom general purpose model has been reported on the directional reflectance spectra of aquatic vegetation, much less appropriate descriptions of the interaction between incident radiation and coupled vegetation–water system. In this paper, we propose a radiative transfer approach for homogeneous aquatic vegetation, namely the Aquatic Vegetation Radiative Transfer model (AVRT). The AVRT model can calculate the bidirectional reflectance spectra of emergent, floating and submerged aquatic vegetation in visible and near infrared bands (400–1000 nm). The submerged part of canopy is also taken into consideration besides canopy architecture, optical active components of water, rough water surface, etc. Basic principles and relative models used in the AVRT are given in Section 2. Implementations are given in Section 3, including a schematic diagram. Then a brief description of the Aquatic Vegetation Monte Carlo model (AVMC) and the corresponding validation tests, a model-to-model comparison and a primary field experiment, are given in Section 4. Finally, the conclusions are given in Section 5. 2. Methodology According to the above review, radiative transfer modeling of aquatic vegetation is recommended to provide a comprehensive description of the interaction between incident radiation and the coupled vegetation–water system at least in VNIR bands. The effects of rough water surface and bottom reflection cannot be overlooked. Model versatility and a good balance between accuracy and efficiency are also important for future applications. Based on these principles, the AVRT takes the following approximations: (1) the coupled vegetation–water system is regarded as homogeneous continuous media, and variations of its optical properties are only allowed in vertical direction; thus (2) the system can be divided into medium layers, and every elemental layer has homogeneous optical properties (only three kinds of medium layers are considered, i.e. canopy layer, water layer and mixed layer of canopy and water); (3) the average height of water surface is taken as zero, and the thickness of water surface is assumed to be zero; and (4) a Lambertian bottom is used for now. The basic framework of the AVRT is based on the methodology of the SAIL model (Verhoef, 1984, 1985). Reflectance and transmittance of a single leaf are calculated by the PROSPECT model (Jacquemoud & Baret, 1990). The PROSPECT and SAIL families are undergoing continuous improvement, and have gone through extensive validations and various applications (Feret et al., 2008; Jacquemoud, Bacour, Poilve, & Frangi, 2000; Jacquemoud et al., 2009; Kallel, Verhoef, Le Hégarat-Mascle, Ottlé, & Hubert-Moy, 2008; Pinty et al., 2004; Verhoef & Bach, 2007). In addition, the hot spot effect (Verhoef, 1998; Kuusk, 2001) of the emergent part of aquatic vegetation canopy, and reflection and transmission across a rough water surface (Cox & Munk, 1956) are also considered. 2.1. PROSPECT At present, the emergent and submerged leaves are all treated as biLambertian. The PROSPECT-5 (Feret et al., 2008) model is used to calculate the reflectance and transmittance of a single leaf in the 400– 1000 nm range (available on the website http://teledetection.ipgp. jussieu.fr/prosail/). The needed leaf biochemical parameters include leaf mesophyll structure (N), pigment contents (Cab for chlorophyll a and b, Ccx for total carotenoid, other pigments such as anthocyanin are not considered), water depth (Cw or EWT for equivalent water thickness) and dry matter content (Cm or LMA for leaf mass per area)

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

(Feret et al., 2008; Jacquemoud & Baret, 1990). For lack of measured data, parameters for common reed (Phragmites australis) are taken from the LOPEX'93 database (Hosgood et al., 1995) (available on the website http://opticleaf.ipgp.fr/index.php?page=database). Default values and their dimensions are listed in Table 1.

(d), upward flux (u), top of the layer (t) and bottom of the layer (b) (Verhoef, 1985). Nomenclaturea Radiation flux Es E− E+ Eo

2.2. Radiative transfer in vegetation canopy 2.2.1. Radiative transfer in turbid medium The AVRT model mainly refers to the SAIL model (Verhoef, 1984, 1985) to solve the radiative transfer problem in a single medium layer. The two-flow approximation is employed to handle the diffuse radiation. The radiative transfer equations in a medium layer can be written as (Verhoef, 1985, 1998) 3 2 Es k 0 6 6 d 6 E− 7 7 ¼ 6 −s dx 4 Eþ 5 4 s Eo w 2

0 α β v

0 −β −α 0 v

32 3 Es 0 7 6 0 76 E− 7 7: 0 54 Eþ 5 Eo −K

3 2 Es ð−1Þ τ ss 6 E− ð−1Þ 7 6 τsd 6 7¼6 4 Eþ ð0Þ 5 4 ρsd Eo ð0Þ ρso

0 τdd ρdd ρdo

0 ρdd τdd τ do

ð1Þ

3 32 Es ð0Þ 0 6 E− ð0Þ 7 0 7 7: 76 0 54 Eþ ð−1Þ 5 Eo ð−1Þ τ oo

Reflectance and transmittance of medium layer or interface ρ(θs,φs;θo,φo) Bidirectional reflectance distribution function (BRDF) τ(θs,φs;θo,φo) Bidirectional transmittance distribution function (BTDF) τss Direct transmittance for solar flux Direct transmittance in observation direction τoo ρsd Directional–hemispherical reflectance for solar flux τsd Directional–hemispherical transmittance for solar flux ρdo Hemispherical–directional reflectance in observation direction τdo Hemispherical–directional transmittance in observation direction ρddc Bi-hemispherical reflectance τddc Bi-hemispherical transmittance



  d E ðbÞ ¼ Td u Rt E ðtÞ

Rb Tu



d

E ðtÞ u E ðbÞ

Absorption and scattering coefficients of waterd a Absorption coefficient b Scattering coefficient bb Back scattering coefficient Vegetation canopy L Leaf area index (LAI) h Height of canopy or a medium layer s‘ Hot spot size parameter

ð2Þ

Eq. (2) can be noted as Eout = ZEin, where Ein and Eout are the matrix of incident and outward fluxes, respectively. Matrix Z is the so-called layer scattering matrix (Verhoef, 1985). Annotations of radiation fluxes include 0 and − 1 for the top and bottom of the layer, respectively. Coefficients in matrix Z are divided into reflectance ρ and transmittance τ. Double subscripts are taken to indicate types of the incident and the outward radiation flux, i.e. s for Es, d for E− or E+ and o for Eo (Verhoef, 1985) (described in the Nomenclature table). Vector E can be substituted by two 2 × 1 subvectors, and matrix Z by four 2 × 2 submatrices. The formulism of block matrices will facilitate the following calculation (Section 3.3). Thus Eq. (2) can be written as (Verhoef, 1985)  ð3Þ

where R and T are matrices of reflectance factors and transmittance factors, respectively. Annotations refer to downward flux

Table 1 Leaf parameters for common reed.

Geometry and others θ Zenith angle φ Azimuth angle Ω Projected solid angle s Unit vector in the direction of solar incidence o Unit vector in the direction of observation n Unit normal vector of wave facet or leaf a

Only part of the parameters are listed, mainly those used in more than one sections. May have subscripts c or w to indicate contributions from canopy or water. c For water surface, two more subscripts are added to indicate the direction of downward (aw, from air to water) and upward (wa, from water to air) diffuse flux. d Have subscripts for water (w) and its components pure water (pw), phytoplankton (ph), nonalgal particles (NAP), colored dissolved organic matter (CDOM) and particles (p, for scattering coefficients only). b

2.2.2. Hot spot of vegetation canopy The original SAIL model (Verhoef, 1984) regards vegetation canopy as turbid medium, and the sizes of scattering elements such as leaf are not considered. Thus, it cannot simulate the hot spot effect of canopy. Correction is added in a later version named SAILH (Verhoef, 1998). For a single layer vegetation canopy, the single scattering contribution to the bidirectional reflectance can be written as (Verhoef, 1998) s

r so ¼ wc

Parameters

N

Cab

Ccx

EWT

LMA

Dimensions

–

μg·cm−2

μg·cm−2

cm

g·cm−2

Value

1.58

98.80

28.35

0.0076

0.0077

Direct solar irradiance Diffuse downward irradiance Diffuse upward irradiance Flux-equivalent radiance Eo = πLo, where Lo is the radiance in the direction of observation

Extinction, attenuation and scattering coefficients (EASCs) of medium layerb k Extinction coefficient for direct flux in the direction of solar incidence K Extinction coefficient for direct flux in the direction of observation α Attenuation coefficient for diffuse flux β Back scattering coefficient for diffuse flux s Scattering coefficient from Es to E+ s′ Scattering coefficient from Es to E− v Scattering coefficient from E− to Eo v′ Scattering coefficient from E+ to Eo w Bidirectional scattering coefficient from Es to Eo

In order to describe the radiative transfer process in water as well, absolute height is used as the vertical coordinate. Positive direction is taken upward. The radiation fluxes and the coefficients in Eq. (1) are described in the Nomenclature table. We will call all these coefficients the EASCs for simplicity. For a mixed layer, subscript c is added to distinguish the contributions of vegetation canopy from those of water (represented by subscript w). Calculation of EASCs of vegetation canopy can be found in Verhoef (1984). Those of water and mixed layer will be given in Section 3.1. The EASCs of canopy and water have the same dimension (m− 1) in the AVRT model. The solution of Eq. (1), or the relationship between incident and outward fluxes, can be expressed as (Verhoef, 1985) 2

3

Z

0

−h

P so ðxÞdx þ r s P so ð−hÞ

ð4Þ

where rs is the reflectance of bottom and h is the canopy height. Pso(x) is the joint probability of the presence of sunlight and the free line of sight from outside the canopy (Verhoef, 1998), also named the bidirectional

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

4

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

gap probability in the canopy at level x (Kuusk, 2001). The following Pso(x) is used (Verhoef & Bach, 2007) h i pffiffiffiffiffiffi qx  P so ðxÞ ¼ exp ðK þ kÞx þ Kk 1−e =q

reference wavelength 500 nm is taken, the scattering coefficient of pure water can be expressed as (Mobley, 1994; Morel & Prieur, 1977; Smith & Baker, 1981)

ð5Þ

where q = Δ / (s‘h) · 2D / (K + k) is a factor controlling the width of the hot spot effect. Parameter Δ is the angle distance measure (Verhoef & Bach, 2007), s‘ is the hot spot size parameter (Verhoef, 1998) and D is the leaf area density defined by D = L / h (Verhoef, 1984). Correction factor 2D / (K + k) is applied to compensate the effect of shadow lengthening in the horizontal plane that occurs for inclined leaves at larger zenith angles (Verhoef & Bach, 2007). Because absolute height and D are used, the form of q given here is not identical with that in Verhoef and Bach (2007). 2.3. Bio-optical model of Case II water The optical active components of Case II water include pure water, phytoplankton, nonalgal particles (NAP) and colored dissolved organic matter (CDOM) (Babin et al., 2003). For lack of measured data and the consideration of model versatility, the representative models from literature are temporarily used. The absorption coefficient, scattering coefficient and back scattering coefficient of water are noted as aw, bw and bbw, respectively. They can be expressed as (Babin et al., 2003; Lee et al., 1998)

 bpw ðλÞ ¼ 0:0022

λ 500

−4:32

:

ð12Þ

The scattering of pure water is simply assumed to be Rayleigh scattering. So the back scattering ratio of pure water can be regarded as 0.5. The scattering coefficient of particles can be expressed as (Lee et al., 1998) bp ðλÞ ¼ EC

0:62

ð13Þ

550=λ

where C is still the concentration of chlorophyll a and E is an empirical value. Lee et al. (1998) used the values 0.3, 1.0 and 5.0 to simulate a range from normal to highly turbid waters. A fixed value E = 1.0 is used in the validation tests, and a back scattering ratio 1.9% is taken due to the scattering phase function of particles used in the AVMC (described in Section 4.1). The scattering phase function is not considered in the present AVRT. The absorption, scattering and back scattering coefficients of water are described in this section. Calculations of the EASCs for water and mixed layer will be given in Section 3.1.

aw ðλÞ ¼ apw ðλÞ þ aph ðλÞ þ aNAP ðλÞ þ aCDOM ðλÞ

ð6Þ

2.4. Reflectance and transmittance of rough water surface

bw ðλÞ ¼ bpw ðλÞ þ bp ðλÞ

ð7Þ

bbw ðλÞ ¼ bbpw ðλÞ þ bbp ðλÞ:

ð8Þ

Reflection and transmission happen when solar radiation goes through water surface. For a flat surface, reflectance and transmittance can be calculated by using the Fresnel Law. However, natural water surface cannot be totally flat under general conditions. When the plants penetrate through the water surface, surface tension will curve the surface around the plants. The shape of the curved surface between two vertical plates is described in Vella and Mahadevan (2005). But the actual situation is much more complicated. The transection of a plant could have a specific shape. Large numbers of plants are located in the scene randomly or following certain distribution. Besides, the plants could have different inclination angles. In addition, there probably are wind-generated waves in the same time. All these seem to make the surface shape unintelligible. However, a noteworthy fact is that the height and the stretch of the curved surface are very short, like about 10 mm (Vella & Mahadevan, 2005). These subtle features could probably be buried in the capillary waves or shadowed by the canopy, for a non-imaging spectrometer or a nadir-view imager with a resolution of more than a few meters. Moreover, the reflection of water surface is usually something that we try to avoid in most of the applications. The situation seems really complicated and unfortunately, we cannot figure out a realistic description of such a water surface now. So the interaction between aquatic vegetation and waves is temporarily not considered, and the widely used slope statistics model is taken. It is our obligation to remind the readers that in certain applications, especially close observation of aquatic vegetation, the following approximation of water surface may be not appropriate and need improvement. The practical method is to regard the rough wave surface as a collection of facets. The probability density function (PDF) of facet slopes depends on wind speed and direction. Details of the slope statistics model can be found in Cox and Munk (1956). The BRDF ρ(θs,φs;θo,φo) and BTDF τ(θs,φs;θo,φo) of wave surface can finally be expressed as (He, Bai, Zhu, & Gong, 2010)

The wavelength λ will be omitted for simplicity. Subscripts pw, ph, NAP, CDOM and p indicate the contributions of pure water, phytoplankton, NAP, CDOM and particles, respectively (see the Nomenclature). The absorption coefficient of pure water refers to both Buiteveld, Hakvoort, and Donze (1994) (400–800 nm) and Palmer and Williams (1974) (800–1000 nm) and is interpolated to make a smooth transition. The absorption coefficient of NAP can be expressed as (Babin et al., 2003) aNAP ðλÞ ¼ aNAP ð443Þ  0:75e

−0:0123ðλ−443Þ

ð9Þ

where the reference wavelength is 443 nm. aNAP(443) is related to the mass of concentration of suspended particulate matter (SPM, in g·m−3) by aNAP(443) = 0.041 m2·g−1·SPM (Babin et al., 2003). A fixed value SPM = 5 g·m−3 is used in the validation tests. The absorption coefficient of CDOM can be expressed as (Volpe, Silvestri, & Marani, 2011) −0:0192ðλ−375Þ

aCDOM ðλÞ ¼ aCDOM ð375Þe

:

ð10Þ

The value 1.25 m−1 is used for aCDOM(375), the absorption coefficient at reference wavelength 375 nm. The specific absorption coefficient of phytoplankton a⁎ph (m2·mg−1) in 400–700 nm can be expressed as (Staehr & Markager, 2004) 

ln aph ðλÞ ¼ −BðλÞ ln C− ln AðλÞ

ð11Þ

where C is the concentration of chlorophyll a (mg·m−3). A(λ) and B(λ) are wavelength dependent parameters determined by fitting the measurement data. Absorption spectra are collected in estuarine, coastal and oceanic waters (covering a chlorophyll concentration range of 0.03–88.1 mg·m−3). According to the specific absorption coefficient of chlorophyll provided by PROSPECT-5 (Feret et al., 2008), the absorption of chlorophyll above approximately 750 nm need not be considered. So we simply assume that the absorption coefficient of phytoplankton has a linear decrease in 700–750 nm, and is zero in 750–1000 nm. If a

ρðθs ; φs ; θo ; φo Þ ¼





r ðδÞ 2 p zx ; zy S θs ; θo ; σ 4 4 cos θs cos θo cos θn

ð14Þ

τðθs ; φs ; θo ; φo Þ ¼





nw t ðδÞ cosδ cosδt 2 p zx ; zy S θs ; θo ; σ 4 c cos θs cos θo cos θn

ð15Þ

2

2

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

where r(δ) and t(δ) are the Fresnel reflectance and transmittance with incident angle δ, respectively; δt is the refraction angle relative to wave facet; nw is the refractive index of water and c is a parameter (Mobley, 1994). θ and φ are the zenith angle and azimuth angle, respectively. The subscripts s, o and n indicate the sun, the observer and the normal vector of the wave facet. p is the probability of a specific facet slope (Cox & Munk, 1956). σ is the root mean square slope. The values for different surface conditions are given in Cox and Munk (1956) and several other values are summarized in Kay, Hedley, and Lavender (2009). The shadowing factor S is used to correct the abnormal output under large zenith angles. Description of this factor can be found in Ottaviani et al. (2008) and the references therein. The roles of these equations will be described in Section 3.2. 3. Implementation

3.1. The EASCs of water and mixed layer As explained at the beginning of Section 2, only three kinds of medium layers are considered in the AVRT, i.e. canopy, water and mixed layer. Three kinds of aquatic vegetation are simulated by the arrangement of different medium layers and water surface. Some aquatic vegetation scenes generated by the AVMC (see Section 4.1 and Appendix A) are given in Fig. 1. If variable X can represent any of the EASCs, then the EASCs of canopy and water can be noted as Xc and Xw, respectively. The EASCs of mixed layer can be divided into contributions from canopy and water (Ackleson & Klemas, 1986; Suits, 1984) X ¼ Xc þ F wXw

ð16Þ

where the weighting parameter Fw can be interpreted as the volume fraction occupied by water in the mixed medium and ranges between 0 and 1 (Ackleson & Klemas, 1986). Calculations of X c can be found in Verhoef (1984) and X w can be calculated by (Suits, 1984) 8 α w ¼ bbw þ 1:4aw > > > > βw ¼ bbw > > > > > sw ¼ bbw sec θs > 0 > > < sw ¼ ðbw −bbw Þ sec θs kw ¼ ðaw þ bw Þ sec θs > > K w ¼ ðaw þ bw Þ sec θo > > > > > vw ¼ bbw secθo > > 0 > > > : vw ¼ ðbw −bbw Þ secθo ww ¼ bbw secθs secθo

effect of ww in the calculated or measured reflectance. So we still follow the approximation in Suits (1984) at present. Once the EASCs are determined, the layer scattering matrix of three kinds of medium layers can be solved by following the procedure in the SAIL model (Verhoef, 1985). 3.2. Scattering matrix for air–water interface Till now the interaction between radiation and medium layer has been expressed by the layer scattering matrix. If the interaction between radiation and air–water interface can be described with a similar matrix, the whole calculation will be better facilitated by using the adding method. A scattering matrix for air–water interface is defined by imitating the layer scattering matrix in Eq. (2). 2

The implementation of the AVRT model will be explained in this section. First, the EASCs of water and mixed layer are discussed in Section 3.1. Then a scattering matrix for rough water surface is described in Section 3.2. Calculation of the bidirectional reflectance for a multilayer media based on adding method is given in Section 3.3. Finally, a schematic diagram of the AVRT is provided in Section 3.4.

ð17Þ

where aw, bw and bbw are the absorption coefficient, scattering coefficient and back scattering coefficient of water as described in Section 2.3. The mean optical path length of diffuse radiation flux in water is greater than that of direct radiation flux, so an empirical coefficient of 1.4 is introduced in α w to account for this (Suits, 1984). Only scattering of pure water is considered in Suits84, so the equations in Eq. (17) are not identical with those in Suits84. Most of the variables are also used in the radiative transfer modeling of water body and their meanings are distinct, except for the w. The variable w in Eq. (1) is originally used in the modeling of vegetation to account for the single scattering contribution. The idea of the single scattering contribution of water body may seem awkward. For the EAV, the whole submerged part only has a minor contribution to the reflectance. We suppose that it is difficult to distinguish the

5

3 2 Es ð0−Þ τ ss 6 E− ð0−Þ 7 6 τsd 6 7¼6 4 Eþ ð0þÞ 5 4 ρsd Eo ð0þÞ ρso

0

0

τddaw ρddaw ρdo

ρddwa τ ddwa τdo

3 32 Es ð0þÞ 0 6 7 0 7 76 E− ð0þÞ 7 0 54 Eþ ð0−Þ 5 Eo ð0−Þ τ oo

ð18Þ

where annotations 0 − and 0 + are for radiation just below and just above the water surface, respectively. Annotations for reflectance and transmittance are omitted and will be explained in the following calculations. Because of total internal reflection, the reflectance and transmittance for upward diffuse flux differ from those for downward diffuse flux. So two more subscripts are added to indicate the diffuse flux, i.e. aw for downward (from air to water) and wa for upward (from water to air). The scattering matrix for interface is named the interface scattering matrix. However, some problems remain to be solved if we want to take advantage of the matrix formulism. The scattering matrix is originally introduced for the vegetation. The incident fluxes are usually divided into the direct radiation and the diffuse radiation. The assumption of isotropic diffuse radiation is reasonable for vegetation in the NIR. Discussions can be found in Verhoef (1985). If the medium has a strong absorption, like the water body in the NIR, the diffuse radiation is very small and has minor effect on the output. So the assumption of isotropic diffuse radiation is also taken in the AVRT model. But water surface is not identical with the vegetation. The incident radiation will be refracted and change its propagating direction when crossing the water surface. The parallel incident radiation will be dispersed by the rough water surface. This could be illustrated by Fig. 7 in D'Alimonte, Zibordi, Kajiyama, and Cunha (2010). The spatial distribution of downward radiation is related with the surface shape and the distance to water surface. But a regression value is reasonable, and the regression value just below water surface is very close to the value of the incident radiation just above water surface. Besides, Haltrin (2006) gives the ratio of the direct irradiance to the diffuse irradiance just below the water surface and the transmittance of direct solar radiation through wind-generated rough water surface. If the wind speed or the incident zenith angle is not too large, the transmittance of direct radiation is very close to the Fresnel transmittance. So we assume that, especially for smooth water surface under small wind speed, the direct radiation will not be evidently dispersed just below the water surface, and the Fresnel transmittance is taken as a reference value. Besides, the definitions of reflectance factors and transmittance factors and the energy conservation law are taken as the principles. The approaches are listed as follows. For a conservative scattering case, the energy conservation law can be expressed as 8 < τ ss ðθs ; φs Þ þ τ sd ðθs ; φs Þ þ ρsd ðθs ; φs Þ ¼ 1 ρ ðθ ; φ Þ þ τdo ðθo ; φo Þ þ τ oo ðθo ; φo Þ ¼ 1 : : do o o ρdd þ τdd ¼ 1

ð19Þ

Given the unit incident radiation from (θs,φs), the sum of reflected and transmitted energy should be equal to the incident energy. In an

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

6

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

above-water incidence case, the energy conservation law is expressed as Z

Z Πþ

ρðθs ; φs ; θo ; φo ÞdΩo þ

Π−

τðθs ; φs ; θo ; φo ÞdΩo

¼ Rðθs ; φs Þ þ T ðθs ; φs Þ ¼ Sðθs ; φs Þ ¼ 1

ð20Þ

where Ω is the projected solid angle and is related to solid angle ω by dΩ = cosθdω = cosθsinθdθdφ (Nicodemus, Richmond, Hsia, Ginsberg, & Limperis, 1977). The upper and lower hemispheres are indicated by Π + and Π− , respectively. Integration results of total reflected energy and total transmitted energy are noted by R(θs,φs) and T(θs,φs), respectively. The sum of R(θs,φs) and T(θs,φs) is noted by S(θs,φs). The directional–hemispherical reflectance factor (DHRF) can be written as Z ρsd ðθs ; φs Þ ¼ Rðθs ; φs Þ ¼

Πþ

ρðθs ; φs ; θo ; φo ÞdΩo :

ð21Þ

The following approximation is taken to calculate the direct and diffuse transmittance factors for direct incidence. T(θs,φs) is determined according to the energy conservation law (Eq. 20) and compared with t(θs). If T(θs,φs) is larger than t(θs), then we get 

τss ðθs ; φs Þ ¼ t ðθs Þ : τsd ðθs ; φs Þ ¼ T ðθs ; φs Þ−τ ss ðθs ; φs Þ

ð22Þ

τss ðθs ; φs Þ ¼ T ðθs ; φs Þ : τsd ðθs ; φs Þ ¼ 0

ð23Þ

The field of view (FOV) of sensor or observer is not considered in present simulations. So the bidirectional reflectance factor (BRF) of rough water surface is related to BRDF by ρso ðθs ; φs ; θo ; φo Þ ¼ πρðθs ; φs ; θo ; φo Þ:

ð24Þ

The diffuse radiation is assumed to be isotropic in the AVRT model. For observers above water, the hemispherical–directional reflectance factor (HDRF) of rough water surface can be expressed as (Nicodemus et al., 1977)

1 π

Z

Z Π Π∓

ρðθs ; φs ; θo ; φo ÞdΩo dΩs :

ð29Þ

According to the energy conservation law, the bi-hemispherical transmittance factor (BHTF) can be determined by ρdd þ τ dd ¼ 1:

ð30Þ

However, Monte Carlo simulation shows that there is little difference between the diffuse transmittance for rough surface and that for flat surface, and a numerical value of 0.93 is given (Turpie, 2012). Thus a flat surface is also used in the AVRT model when calculating the bi-hemispherical reflectance and transmittance. Then analytical formulae can be used. Given a refractive index of water nw = 1.34, the BHRF for diffuse incidence from air will be ρddaw = 0.0675 according to Preisendorfer (1976). So we can get τddaw = 0.9325, which is consistent with the result of Turpie (2012). The effect of total internal reflection needs to be considered for the diffuse incidence from water, and a value of ρddwa = 0.4807 can be obtained according to the formula in Lyzenga (1977). Then we can get τddwa = 0.5193. Again these parameters are based on assumptions that have divergences with the real water surface. Potential errors caused by this need further evaluation in certain applications.

ρdo ðθo ; φo Þ ¼ Rðθo ; φo Þ ¼

The optical properties of infinitesimal medium layer will change in vertical direction for the coupled vegetation–water system. So the whole system must be divided into medium layers to assure homogeneous optical properties in each layer. According to Sections 2.2 and 3.2, the relationships between radiation and medium layer or interface have been expressed by the scattering matrix. Thus the bidirectional reflectance of the whole system can be easily calculated based on the adding method (Cooper, Smith, & Pitts, 1982; Verhoef, 1985). The whole system is divided into N + 1 layers (including the air– water interface and bottom) as shown in Fig. 2. The reflectance properties of the bottom can be described by the surface reflectance matrix (Verhoef, 1985) 

Z Πþ

ρðθs ; φs ; θo ; φo ÞdΩs :

ð25Þ

Similarly, the hemispherical–directional transmittance factor (HDTF) is defined by Z τdo ðθo ; φo Þ ¼ T ðθo ; φo Þ ¼

ρdd ¼

3.3. Adding method

Otherwise we get 

For diffuse incidence, the bi-hemispherical reflectance factor (BHRF) can be calculated by

Π−

τðθs ; φs ; θo ; φo ÞdΩs :

ð26Þ

Rs ð1Þ ¼

r sd r so

 r dd : r do

ð31Þ

So matrix Rs(i) is for the bottom surface of Layer i, i.e. the reflectance properties of the ensemble from Layer 0 to Layer i − 1. The upward and downward fluxes at the bottom of Layer i can be related by Rs(i) u

d

E ði; bÞ ¼ R s ðiÞE ði; bÞ:

ð32Þ

Thus Eq. (3) can be rewritten with the index of layer According to the reciprocity principle, we can get ρsd ðθ; φÞ ¼ ρdo ðθ; φÞ:

 ð27Þ

However, the reciprocity equation is not valid for transmission because of the total internal reflection. Numerical results also demonstrate that the sum of R(θo,φo) and T(θo,φo) is generally smaller than unit for observers above water. According to the energy conservation law, we assume that the rest of the energy received by the observer comes from the direct transmission in observation direction. Then we have τoo ðθo ; φo Þ ¼ 1−ρdo ðθo ; φo Þ−τdo ðθo ; φo Þ:

ð28Þ

  d E ði; bÞ ¼ Td ðiÞ u Rt ðiÞ E ði; tÞ

Rb ðiÞ Tu ðiÞ



 d E ði; tÞ : u E ði; bÞ

ð33Þ

Then we can get the following equation by solving Eqs. (32) and (33) (Verhoef, 1985)     u  d r ðiÞ r dd ðiÞ d E ði; tÞ E ði; tÞ ¼ Rt ðiÞE ði; tÞ ¼ sd   r so ðiÞ rdo ðiÞ h i −1 d ¼ Rt ðiÞ þ Tu ðiÞðI−Rs ðiÞR b ðiÞÞ Rs ðiÞTd ðiÞ E ði; tÞ

ð34Þ

where I is the unit matrix. The reflectance properties of the ensemble from Layer 0 to Layer i are described by R∗t (i): directional–hemispherical

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

reflectance r∗sd(i), bi-hemispherical reflectance r∗dd(i), bidirectional reflectance r∗so(i) and hemispherical–directional reflectance r∗do(i). Obviously we arrive at R∗t (i) = Rs(i + 1). The adding begins from the bottom. After a layer is added, the present position is calculated and compared with the position of water surface. If it reaches the water surface, the scattering matrix of the interface is used in the adding. Then the adding continues to the next layer. Therefore, the reflectance matrix of the whole system R∗t (N) can be easily calculated by iteration of Eq. (34) (Verhoef, 1985). Finally we obtain the bidirectional reflectance of aquatic vegetation under both direct and diffuse incident radiations     r so ¼ r so ðNÞEs ðN; t Þ þ rdo ðNÞE− ðN; t Þ =ðEs ðN; t Þ þ E− ðN; t ÞÞ:

ð35Þ

3.4. Overall schematic of the AVRT model The details of the AVRT model have been given and an overall schematic is shown in Fig. 3. All the input parameters are represented by the document symbol. Models or methods are given by blue rectangles. White rectangles stand for temporary data or output results. In our calculations, the whole system is first divided into medium layers according to the need of a specific problem. The reflectance and transmittance of a single leaf are calculated by the PROSPECT model. The optical properties of water are calculated by the bio-optical model of water. Thus the scattering matrix of each medium layer can be calculated based on the algorithms in the SAILH model. The scattering matrix of air–water interface can be obtained based on the Cox–Munk model. With all these scattering matrices and bottom reflectance matrix, the reflectance matrix R∗t of aquatic vegetation can be calculated based on the adding method. Then the ratio of diffuse flux is considered and we can get the bidirectional reflectance of aquatic vegetation. 4. Validation Validation of the AVRT model will be given in this section. In Section 4.1, a Monte Carlo model for aquatic vegetation is briefly described. Then the model-to-model comparison validation is given in Section 4.2. In Section 4.3, a primary field experiment is described and simulation results of the AVRT model are compared with measured data.

7

4.1. Monte Carlo simulation Structure elements of real vegetation canopies have finite sizes and specific spatial distributions. The AVRT model is mainly based on the assumption of turbid medium and the sizes of scattering elements are neglected. However, the turbid medium is still a good approximation, especially for homogeneous dense canopies. So the Aquatic Vegetation Monte Carlo (AVMC) model with randomly positioned equilateral triangle leaves is used as a transition between the turbid medium and the real canopy scene. The results of the AVMC are used for comparison and validation. In the AVMC, the inclination angles of leaves fulfill certain probability density functions (PDF) and cumulative distribution functions (CDF). Leaves are regarded as bi-Lambertian. Optical properties of water and water surface are the same as those described in Sections 2.3 and 2.4, respectively. The Petzold's San Diego scattering phase function is taken from Leathers, Downes, Davis, and Mobley (2004) and the water bottom is assumed to be Lambertian. When measured dataset or better information is available, the relative parts in the AVMC model can be modified. Some relative details of the AVMC are given in Appendix A. In this paper, we focus on the bidirectional reflectance properties of EAV. Results under various conditions are given. The effects of the side length of leaf (SLL), incident zenith angle, leaf inclination distribution function (LIDF) and LAI are also analyzed. Three kinds of SLLs are used in the validation tests. They are 0.10 m, 0.05 m and 0.01 m, and abbreviated as SL10, SL5 and SL1, respectively. The emergent vegetation scene is divided into two layers. The lower layer is a mixed layer with a thickness of 0.4 m and a LAI value of 1. The upper layer is a canopy layer with a thickness of 0.6 m and a default LAI value of 2. The dimensions of the scene are 2 m × 2 m for SL10 and SL5, and 1 m × 1 m for SL1 to reduce data size. The number of leaves in each layer is determined according to scene area, LAI and SLL. The barycenter of a triangle leaf is randomly located in each layer. Then the finite scene is repeated outwards in eight neighborhoods to simulate an infinite vegetation scene. Look-up tables (LUTs) for the CDF of leaf inclination angles are built according to the LIDF of the SAIL model (Verhoef, 1998). The LIDF is determined by two parameters LIDF.a and LIDF.b. Three kinds of LIDFs are used in the validation tests: planophile (PLA, with LIDF.a = 1), uniform

Fig. 1. Examples of aquatic vegetation scene: EAV (a), SAV (b) and FAV (c, d and e). FAV may only have the emergent part (d) or submerged part (e). Green triangles represent for leaves. Spaces bounded by blue boxes are filled with water. For simplicity, a plane z = 0 is used to indicate the average position of rough water surface. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

8

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

then counted according to their directions. The hemisphere above the scene is divided into patches of equal zenith angle width dθ = 2° and equal azimuth angle width dφ = 5°. The number of patches is N = 2880. The weight of photons is used in the AVMC (described in Appendix A). Then the BRF in the direction of a specific patch can be calculated by (Govaerts & Verstraete, 1998) r MC ¼ πW o =ðdΩo W s Þ

ð36Þ

where dΩo is the corresponding projected solid angle of the patch in observation direction. The sum of the weights of total incident photons is Ws, and that of exit photons in the observation solid angle is Wo. The result of the AVRT in Eq. (35) is rewritten as rRT, and rMC is regarded as the reference. Then the relative error is defined by ε r ¼ ðr RT −r MC Þ=r MC : Fig. 2. Medium layers of emergent aquatic vegetation (a) and submerged aquatic vegetation (b) after division. They are layers 0 to N from bottom to top. The air–water interface is regarded as a special layer and its thickness is exaggerated. The vertical green lines indicate vegetation and areas in blue indicate water. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

(UNI, with LIDF.a = 0) and erectophile (ERE, with LIDF.a = −1). LIDF.b is always assigned to be 0. Only parallel incidence is considered in the tests. All the incident photons have the same incident direction but each has a random position on top of the canopy (TOC). The number of incident photons used in each test is 2 × 108. Three incident zenith angles 150°, 135° and 120° are used and the default incident azimuth angle is 180°. Incident directions are noted as Ts30, Ts45 and Ts60 for practice. Since a simulated infinite scene is used, the exit photons from TOC are identified and

ð37Þ

The absolute value of the relative error is written as |εr|. Two wavelengths are chosen in the tests according to Adam et al. (2010) and Turpie (2012). They are 675 nm in RED band and 920 nm in NIR band. Reflectance rl and transmittance tl of a single leaf are listed in Table 2. Results under different conditions are given in the following sections, mainly BRF and εr or |εr| in principle plane (PP), orthogonal plane (OP) and the upper hemisphere. 4.2. Model comparison and validation 4.2.1. Side length of leaf The AVRT can be regarded as a hybrid model of turbid medium and canopy hot spot correction. Turbid medium is the ultimate state of

Fig. 3. Overall schematic of the AVRT model. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx Table 2 Reflectance and transmittance of a single leaf. Band

rl

tl

RED NIR

0.0397 0.4653

0.0001 0.4597

infinitesimal scattering elements. However, the hot spot effect of vegetation canopy is caused by finite discrete scattering elements such as leaves. The triangle leaves in the AVMC have a specific size and therefore infinitesimal leaves would be impractical. Thus, the influence of SLL is first evaluated. The default incident zenith angle is Ts45. Three SLLs are tested, namely SL1 (0.01 m), SL5 (0.05 m) and SL10 (0.10 m). The corresponding hot spot size parameters s‘ used in the AVRT are 0.005, 0.02 and 0.03, respectively. This setting gives good fit around the canopy hot spot direction as shown in Fig. 4. There is another hot spot in forward direction due to the specular reflection of water surface. This is named the forward hot spot or mirror hot spot in this paper. The backward direction is distinguished from the forward direction by negative view zenith angles. As shown in Fig. 4, the width of canopy hot spot becomes larger with an increasing SLL. Meanwhile, an evident rise of rMC around the mirror hot spot direction is observed. However, the situation is the opposite except for the hot spot directions. The decrease of rMC is obvious in OP. All these indicate a larger gap probability in the canopy with the same LAI but larger leaves. Relative errors are generally within ± 2% for SL1, and get larger for SL5 and SL10. Errors larger than ±10% are not plotted. The distribution of |εr| in the upper hemisphere for SL1 is given in Fig. 12.b. For SL1, about 96.7% of the 2880 points have |εr| values smaller than 2%. Large errors mainly distribute around the hot spot directions. Unfortunately, the results for SL10 are poor. Considering the small noises in rMC, SL1 gives fine results. Smaller leaves will consume more computer resources. So SL1 is used as a default setting in the rest of the validation tests.

9

4.2.2. Leaf inclination distribution function Triangle leaves in the AVMC are generated according to specific LIDF. The distribution of leaf inclination angle has significant effects on the bidirectional reflectance properties of vegetation canopy. So the influence of LIDF is evaluated is this section. The default incident zenith angle is still Ts45. Three kinds of LIDFs are tested, i.e. PLA, UNI and ERE. The corresponding values of the hot spot size parameter s‘ used in the AVRT are 0.005, 0.01 and 0.2, respectively. As shown in Fig. 5, the results for ERE around canopy hot spot direction are not that satisfying. With the LIDF changing from PLA to ERE, more leaves get larger inclination angles. The BRF values increase at large view zenith angles in PP, while a sharp decrease is observed around the nadir in both PP and OP. The summit of canopy hot spot also decreases. The values of rMC become larger than that of rRT at large forward zenith angles in PP. Distinct divergences between rMC and rRT are observed at large zenith angles in OP. The AVMC gives a higher mirror hot spot for UNI and ERE. This also indicates a larger gap probability. Although the results for PLA are fine, the AVRT underestimates the BRF values at large view zenith angles for UNI and ERE. Distributions of rMC and rRT under different kinds of LIDFs are given in Fig. 6. This underestimation is evident for ERE (Fig. 6.c and f). Distributions of |εr| in the upper hemisphere for UNI and ERE are given in Fig. 12.g and h. Results for UNI are still not too bad at small zenith angles and backward directions. About 38.8% of the values are within 2%, and 60.7% within 5%. But the results for ERE are rather unacceptable except for those at some backward directions. The corresponding numbers are 6.9% and 13.8%, respectively. These underestimations imply more contributions from multiple scattering. The two-flow approximation used in the AVRT cannot handle the multiple scattering well when dealing with erect-dominated leaves. So the BRFs are divided into contributions from single scattering and multiple scattering to verify this assumption. Each interaction between a photon and the canopy, water or water surface is defined as scattering. So a photon transmitted into water

Fig. 4. BRFs and relative errors in PP and OP under various SLLs. The dot lines and solid lines are for rMC and rRT, respectively. Their labels are on the left y axis. The crosses are for εr and corresponding labels are on the right y axis. View zenith angles range from −79° to 79° with an interval of 2°, and negative values indicate backward directions relative to the solar incidence. The first row shows results in PP and the second shows results in OP. Columns 1, 2 and 3 are for SL1, SL5 and SL10, respectively.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

10

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

Fig. 5. BRFs and relative errors in PP and OP under different kinds of LIDFs. Lines, symbols and settings are the same with Fig. 4. The first row shows results in PP and the second shows results in OP. Columns 1, 2 and 3 are for PLA, UNI and ERE, respectively.

Fig. 6. BRFs in the upper hemisphere under different kinds of LIDFs. Rows 1 and 2 are results of AVMC and AVRT, respectively. Columns 1, 2 and 3 are for PLA, UNI and ERE, respectively. Figures in the same column share the same color bar. View zenith angles range from 1° to 79°. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

must be scattered more than once before it reaches a sensor. Then Eq. (4) can be used to calculate the single scattering contribution in the AVRT, with water surface at the bottom of the canopy layer. Scattered times of a photon can be counted in AVMC. The other part is regarded as the multiple scattering contribution after the single scattering contribution is determined. The single scattering and multiple scattering contribution under different kinds of LIDFs are shown in Figs. 7 and 8, respectively. For PLA, the AVRT slightly underestimates the single scattering contribution, as shown in Fig. 7.a and d. However, the situation is the opposite for ERE, as shown in Fig. 7.c and f. The differences are generally not significant except backward directions for UNI (Fig. 7.b and e) and ERE (Fig. 7.c and f). This is due to the large s‘ used in the AVRT to get a better fit of total BRF. In other words, the single scattering contribution for UNI and ERE is overestimated when total BRF fits well around the canopy hot spot direction. This also indicates an underestimation of the multiple scattering contribution, as shown in Fig. 8. The multiple scattering contributions calculated by the AVRT for UNI and ERE (Fig. 8.e and f) are evidently smaller than the results of the AVMC (Fig. 8.b and c) especially at large zenith angles. However, the differences for PLA (Fig. 8.a and d) are not significant. This demonstrates the AVRT's shortage in proper handling of multiple scattering for erect-dominated leaves. Thus a better approximation of multiple scattering should be one of the emphases in the future AVRT. 4.2.3. Leaf area index Compared with sparse canopy, dense canopy is a better approximation of turbid medium. Therefore canopy density, or LAI for a fixed area, is also an important factor affecting the bidirectional reflectance of aquatic vegetation. The default LAI of the canopy layer is 2 in the above tests. It is changed to 1 and 3, and the performances of the

11

AVRT are evaluated. The incident zenith angle is still Ts45. The s‘ value used in the AVRT is 0.005 for SL1. The BRFs and relative errors are shown in Fig. 9. Most BRFs get larger values in both PP and OP with an increasing canopy LAI. However, the mirror hot spot becomes much weaker, and nearly vanishes for LAI = 3. No significant reduction of εr is observed when LAI increases from 2 to 3. But values of εr become a little larger especially around the mirror hot spot when LAI decreases from 2 to 1. The biggest εr is about −10% in mirror direction. This is evidently shown in Fig. 12.e and f, the distributions of |εr| for LAI = 1 and LAI = 3. For LAI = 3, about 97.7% of the values are within 2%. However, the number for LAI = 1 is just 16.2%. It indicates that a LAI of 1 is insufficient for the AVMC to simulate a turbid medium properly. However, the mirror hot spot is not significant for a LAI of 3. So the default LAI of 2 is still used in the rest of the tests. 4.2.4. Incident zenith angle The BRFs of aquatic vegetation have been calculated under different kinds of SLL, LIDF and LAI. A fine fit between models can be obtained with the settings of SL1, PLA and LAI = 2. In this scene, the reflectance characteristics of aquatic vegetation are significant. So these settings are used as the default in the rest of the tests. Influences from other factors can then be evaluated. The s‘ value used in the AVRT is 0.005 for SL1. The bidirectional reflectance of vegetation canopy or water is affected by incident zenith angle. For water surface, a sharp increase of reflectance can be observed with an increasing incident zenith angle. So the reflectance properties of aquatic vegetation also need to be investigated under different incident zenith angles. Three kinds of incident zenith angles are tested: Ts30, Ts45 and Ts60. The results are shown in Fig. 10. The BRFs and relative errors exhibit no significant variations in OP under these incident zenith angles. However, BRFs around mirror direction and at large zenith angles in PP become

Fig. 7. Single scattering contribution of BRF under different kinds of LIDFs. Settings are the same with Fig. 6.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

12

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

Fig. 8. Multiple scattering contribution of BRF under different kinds of LIDFs. Settings are the same with Fig. 6.

Fig. 9. BRFs and relative errors in PP and OP under different LAIs of the canopy layer. Lines, symbols and settings are the same with Fig. 4. The first row shows results in PP and the second shows results in OP. Columns 1, 2 and 3 are for canopy LAI of 1, 2 and 3, respectively.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

13

Fig. 10. BRFs and relative errors in PP and OP under different incident zenith angles. Lines, symbols and settings are the same with Fig. 4. The first row shows results in PP and the second shows results in OP. Columns 1, 2 and 3 are for Ts30, Ts45 and Ts60, respectively.

larger when the incident zenith angle changes from Ts30 to Ts60. The sharp increase of BRF due to reflection of water surface is obvious for Ts60. Because the sample view zenith angles have an interval of 2°, those around the canopy hot spot for Ts30 and Ts60 are not the exact hot spot direction. So a summit like that for Ts45 does not appear. Relative errors are generally within ±2% except around mirror direction for Ts45 and Ts60. Distributions of |εr| for Ts30 and Ts60 are shown in Fig. 12.a and c. The errors are generally very small except around mirror direction. For Ts30, about 96.0% of the values are within 2%. The number for Ts60 is 94.8%. The divergences around mirror direction become

larger with increasing incident zenith angles. This also indicates a larger gap probability. 4.2.5. Wavelength Healthy leaves usually have high reflectance and transmittance values in NIR band. However, their reflectance and transmittance are generally low in VIS band. This implies different bidirectional reflectance properties of the canopy. Vegetation indexes are commonly applied to extract information of vegetation canopies from remote sensing images. Usually one reflectance in RED band and another in NIR

Fig. 11. BRFs and relative errors in PP (left) and OP (right) at 675 nm. Lines, symbols and settings are the same with Fig. 4.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

14 G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

Fig. 12. Distributions of |εr| for SL1 under different conditions: (a) Ts30_PLA; (b) Ts45_PLA; (c) Ts60_PLA; (d) Ts45_PLA_675; (e) Ts45_PLA_LAI = 1; (f) Ts45_PLA_LAI = 3; (g) Ts45_UNI; (h) Ts45_ERE. If wavelength or LAI is not noted, that indicates the default values of 920 nm and LAI = 2. View zenith angles range from 1° to 79°. Gray areas indicate values larger than 10.0% in (a)–(f), and 15.0% in (g) and (h). Black areas are actually not negative values but due to smoothing.

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

band are used. So the performances of the AVRT in RED band are evaluated in this section. BRFs and relative errors are shown in Fig. 11. BRF values are significantly smaller than those of 920 nm (column SL1 in Fig. 4). The bowl shape is no longer seen. The mirror hot spot is higher than the canopy hot spot. The AVRT again underestimates the BRFs around the mirror direction. Values of |εr| for 675 nm (Fig. 12.d) are larger than those for 920 nm (Fig. 12.b). Only 22.0% of the values are within 2%, and 96.4% within 5%. However, the numbers for 920 nm are 96.7% and 99.8%, respectively. 4.2.6. Summary and discussion Due to the difficulties encountered in field experiments, the results from the AVMC are used temporarily as the comparison and validation. Distributions of BRF and relative error under different conditions have been described and analyzed. Distributions of |εr| are summarized in Fig. 12. Statistics are given in this section to get an overall performance evaluation of the AVRT. Only conditions with SL1 are selected. Correlation coefficient (R) and root mean square error (RMSE) defined by Eqs. (38) and (39) are used as the statistics. N X ðrRT ðiÞ−rRT Þðr MC ðiÞ−r MC Þ i¼1 R ¼ vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N N X uX 2 2 t ðr RT ðiÞ−r RT Þ  ðr MC ðiÞ−rMC Þ i¼1

ð38Þ

i¼1

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N u1 X 2 ðr ðiÞ−r MC ðiÞÞ RMSE ¼ t N i¼1 RT

ð39Þ

15

where the number of samples is N = 2880. The average value of rRT is noted as r RT , and that of rMC is r MC . Different conditions are labeled by the incident zenith angle and LIDF. If the wavelength or LAI is not the default value, i.e. 920 nm and LAI = 2, the value used will also be listed. Results of the model to model comparison are shown in Fig. 13. Results for PLA generally have small RMSE (b1%) except when LAI = 1 (Fig. 13.e). But the correlation coefficients are not that satisfying (about 90% or even lower) in NIR band (Fig. 13.a, b, c and f). This is mainly due to the large divergence around hot spot and the noises. Because of the lower canopy hot spot and higher mirror hot spot, values of R get larger in RED band (Fig. 13.d) and with LAI = 1 (Fig. 13.e). This implies that R is affected by the abnormal values of rRT around canopy hot spot (circles above the text in Fig. 13.b and f). The correlation coefficients are fine for UNI (Fig. 13.g) and ERE (Fig. 13.h). However, values of RMSE are larger, and the AVRT underestimates the majority of all the BRF values. Besides the noises in the AVMC, other factors may also influence the results. (1) Firstly it is the hot spot size parameter s‘ used in the AVRT. Although we compared the results under several different s‘ values before making a selection, the values used at present are not the optimal. Moreover, the AVRT and the AVMC cannot give close match in the canopy hot spot direction due to the differences of their principles. (2) Secondly, the LUT for rough water surface is expedient to avoid a complicated implementation such as triangular network and keep the same with the AVRT in the characteristics of water surface. Partition of the sphere with equal angle width is insufficient around the glitters and creates redundant data. An adaptive partition or a realistic implementation of rough water surface may give better results around mirror direction and improve the efficiency. (3) Thirdly, sampling of positions, i.e. the barycenter of triangle leaves and initial position of

Fig. 13. Comparisons of rMC and rRT for SL1 under different conditions. Conditions are the same with Fig. 12. All the x axes have the same label rMC, and y axes have the label rRT. Solid lines are the 1:1 line. Correlation coefficient (R) and root mean square error (RMSE) are given in the figure. Each figure has its own scale.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

16

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

Fig. 14. Location of the experiment site. The red dots indicate the location on a map of Beijing (a) and an image from Google Earth™ (b). Locations of the two aquatic vegetation communities are marked by red rectangles on a zoom image (c). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

photons on TOC, also makes a difference. Basic random sampling is used temporarily in the AVMC. That is to say, intersections of leaves are possible and the initial positions of some photons might be too close. Some space-filling sampling strategies may be helpful to improve the efficiency (Janssen, 2013). 4.3. Field experiment In order to further evaluate the performance of the AVRT model, field experiments were carried out in the Beijing Olympic Park. Two kinds of aquatic vegetation, Scirpus fluviatilis and Scirpus triqueter, were chosen as the main objects because relatively homogeneous communities could be found. The experiment site is shown in Fig. 14 and photos of the two vegetation communities are shown in Fig. 15. A

simple multiple-angle observing platform was set up and a hand-held ASD radiometer was used to collect the multiple-angle spectra of the vegetation. Some other parameters such as LAI were also measured. However, the optical properties of water and the biochemical data of leaves are not available. The basic information of the two vegetation communities is listed in Table 3. A micro camera was fixed on the radiometer to capture quasisynchronous images of the object when measuring the spectra. Example multiple-angle images of S. fluviatilis are shown in Fig. 16. Due to lack of some auxiliary data, some parameters must be fitted according to the measured spectra. Only BRFs in the NIR band are analyzed here because the reflected spectra are less affected by water body and pigments. A comparison of measured data and model simulation results at 920 nm is shown in Fig. 17.

Fig. 15. Photos of the two aquatic vegetation communities: (a) S. fluviatilis and (b) S. triqueter.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

5. Conclusions

Table 3 Basic information of the two vegetation communities. Vegetation type

LAIa

Mean tilt anglea (°)

Canopy heightb(m)

Water depthb (m)

S. fluviatilis S. triqueter

2.42 1.26

51 61

1.0 0.7

0.2 0.2

a b

17

Measured by LAI-2000. Approximation.

The measured data and simulation results exhibit fine overall consistency as shown in Fig. 17, but divergences caused by various factors are also evident. In Fig. 17.a, distinct divergences can be found around 15° where the measured data has a hollow curve and their values are lower than the simulated ones. This can be attributed to the inhomogeneity of the canopy, as shown in Fig. 16.b. Theoretically speaking, two hot spots are expected on the BRF pattern of the EAV: a backward one due to the scattering of vegetation canopy and a forward one due to the mirror reflection of water surface (Fig. 16.a). Unfortunately, the measured data in Fig. 17.a do not seem to provide satisfying evidences. No evident peaks around ±40° can be identified in the measured data. Because the sensor has a specific FOV and the output is an averaged value, the reflectance peaks may be smoothed or even buries due to the inhomogeneity of the scene (the dark region in Fig. 16.a) or the shadows of the holder and the equipments (Fig. 16.c). Under large view zenith angles, the AVRT model underestimates the BRF values. This is a drawback of the present model because of the usage of the two-flow approximation and has been discussed in Section 4.2.2. In addition, the height of the present holder is within 2.0 m which may be insufficient because the objects in the sensor FOV will change distinctly under various view zenith angles. This may also contribute to the divergences under large view zenith angles. The S. triqueter canopy is more homogeneous and has smaller LAI value than the S. fluviatilis canopy, so better consistency between measured data and simulation results can be found around nadir and the peak due to reflection of water surface can also be identified around 30°, as shown in Fig. 17.b. However, the AVRT model overestimates the BRFs around the mirror direction. Under large view zenith angles, the situation is similar with that in Fig. 17.a. Due to the limited availability of relative equipments and auxiliary data, the present validation stops here but more work should be done in the future. A more complete dataset which includes the optical properties of water and the biochemical data of leaves is highly recommended and the measured and simulated spectra of various aquatic vegetation should be compared. Besides, it is worthwhile to evaluate the influences of observation distance and sensor FOV on the BRFs, especially around mirror direction. In certain applications, the bidirectional reflectance of water surface may need a better description. Moreover, other challenges still exist in the field experiments of aquatic vegetation, such as the lack of an established protocol (Turpie, 2012).

A radiative transfer model of aquatic vegetation (the AVRT model) is proposed in this paper. Aquatic vegetation and its habitat, i.e. water and bottom, are regarded as a whole system. Theoretically, the AVRT model can calculate the bidirectional reflectance spectra of emergent, floating and submerged aquatic vegetation. The bidirectional reflectance properties of emergent aquatic vegetation are validated here by a Monte Carlo model (the AVMC model) and compared with the measured data of a primary field experiment. Randomly located triangle leaves are used in the AVMC to simulate a homogeneous canopy. The two hot spots caused by canopy and water surface are clearly shown in both models. Comparison and validation generally exhibit close consistencies and small errors, with settings of planophile leaf inclination distribution and small leaf size. However, results are not so satisfying under certain conditions such as erect-dominated leaves. Besides, the simulation results have divergences with measured data because of drawbacks of the present model and uncertainties in the field experiment. Although it seems to be complicated, canopy modeling of aquatic vegetation based on a general approach such as radiative transfer is indeed feasible. The AVRT model can be applied to spectra simulation and sensitivity analysis of aquatic vegetation. The radiative transfer process may be analyzed to promote our understanding of the interactions between radiation and aquatic vegetation. Possible applications may include classification and information extraction, e.g. retrieval of chlorophyll or LAI in a shallow lake with aquatic vegetation. It might be a little earlier to draw a conclusion on the performance of the AVRT model now. Further validations, especially complete field experiments, are still needed in the future. A better handling of multiple scattering is highly recommended, because erect-dominated leaves are commonly seen in aquatic vegetation communities. Precise description of the water surface in aquatic vegetation region remains an unsolved problem which might cause potential problems in certain applications. Evaluation of the potential problems and a better description of water surface are included in future objectives, as least for us. Considering the high spatial complexity of aquatic vegetation, a geometrical optical (GO) model or hybrid GO-RT model may be useful in future applications. Available hyperspectral data with high spatial resolution are helpful for model validation and applications. Moreover, a more realistic 3D computer model may be better for comparison and validation than randomly positioned leaves. Acknowledgment This work was supported by the National Natural Science Foundation of China (NSFC) under Grant No. 40901168. The authors thank Dr. Feng Zhao, Dr. Yongming Du and Xu Dai for helpful suggestions; Mr. Jackson Turner for language improvement; and Peng Zhang and Yiqing Guo for discussions. The authors are also grateful to the anonymous reviewers for their valuable comments.

Fig. 16. Images of S. fluviatilis under various view zenith angles: (a) 40°, (b) 10° and (c) −40°.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

18

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

Fig. 17. Comparison of measured data and model simulation results at 920 nm for S. fluviatilis (a) and S. triqueter (b). Measured data are marked by triangles and the corresponding simulation results are marked by circles. Negative angle values indicate backward directions relative to the solar incidence. R is correlation coefficient and N is the number of measured data.

Appendix A. Aquatic Vegetation Monte Carlo model The principles of Monte Carlo (MC) simulation and applications to water and vegetation can be found in Leathers et al. (2004), Gimond (2004), Govaerts and Verstraete (1998), and Disney et al. (2000). The settings of the AVMC model have been described in Section 4.1. Details of some characteristics are given below. Construction of the scene A brief description of the emergent vegetation scene is given in Section 4.1. Generation of a triangle leaf and ray tracing accelerating strategy used in the AVMC will be described in this section. Generation of the equilateral triangle leaf is shown in Fig. A1. Firstly an equilateral triangle is generated in the xOy plane. The origin is its barycenter and a vertex A is on +x axis. Its unit normal vector nl is on +z axis. Then the triangle and nl are rotated around z axis by a certain angle φ1 (①). The vector nl will not change but vertex A gets a

random azimuth angle in the xOy plane. A certain CDF of leaf inclination angle is calculated according to the SAIL model (Verhoef, 1998). nl is rotated around y axis by a specific inclination angle θl determined by the CDF and a random number (②). Then nl is rotated around z axis again by another azimuth angle φ2 (③). In these two steps, the triangle is also rotated in the same way with vector nl. Azimuth angles φ1 and φ2 are between 0° and 360°, and the inclination angle θl is between 0° and 90°. After these steps, vector nl may point to any direction in the upper hemisphere and θl fulfills a specific LIDF. Finally a translation vector Vl for the barycenter is randomly generated according to the dimensions of the scene and all the vertexes are translated (④). The normal vector nl will not change in this translation. In a ray tracing model, calculation of ray-object intersection is time consuming. There are a large number of triangle leaves in the AVMC, and a proper accelerating strategy can evidently improve the efficiency. Because the leaves are uniformly distributed in the scene, acceleration can be achieved based on space subdivision, typically the Three Dimensional Digital Differential Analyzer (3DDDA) method (Fujimoto, Tanaka, & Iwata, 1986). The basic idea is to divide the space into 3D grids, namely the voxels. An auxiliary data structure is built to set up mappings between voxels and the leaves therein. A ray propagates in the voxels along its path. Intersection test is only needed between the ray and the leaves in the present voxel. Once intersection is detected, there is no need to check the rest of the voxels on the path. Thus only a few leaves need to be examined with the intersection test. With this the efficiency is improved. Relative details, e.g. traversal of the grids, intersection calculation between a line and a triangle, can be found in monographs for ray tracing or rendering techniques (Shirley & Morley, 2003). The AVMC model is written in C++ language. Our test platform is a Windows PC with 2.7 GHz CPU and 1 G RAM. When there are about 7 × 104 triangle leaves in the scene, it takes about 12 h to trace 5 × 108 photons for EAV in the NIR band. The consumed time slightly varies with leaf number, incident zenith angle, LIDF, canopy height, etc. LUT for rough water surface

Fig. A1. Generation of the equilateral triangle leaf.

When the photon is intercepted by water surface, its exit direction needs to be determined. Reflected and transmitted energies are concentrated around the specular reflection and transmission direction, respectively. The BRDF and BTDF patterns depend on incident direction, wind speed and direction. So LUTs are built based on Sections 2.4 and

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

3.2 to determine the exit direction. The energy conservation law (Eq. 20) can be rewritten as Z 0Z 2π π

0

f ðθs ; φs ; θo ; φo Þ cos θo sin θo dθo d φo ¼ 1

ðA1Þ

where f(θs,φs;θo,φo) is the BRDF or BTDF of rough water surface. Then we can get the CDF of energies in different exit directions under different incident directions. The sphere is divided into grids with equal angle width 1° × 1°, and exit energy in each grid is calculated. Incident zenith angles range from 0° to 180° with an interval of 1°. Numerical integration has a certain error especially under grazing incidence. So the CDF values are normalized after this primary calculation. Then small increments are removed according to a given threshold (like 10−6) to reduce the number of data in the LUT. An example is shown in Fig. A2 with an incident zenith angle of 150°. Arrow A in Zoom A and Arrow B in Zoom B indicate the refracted glitter and reflected glitter, respectively. Larger threshold value will reduce the number of data in the LUT but may smooth away the subtle features, e.g. Arrow B will become Arrow C. However, smaller threshold value indicates more data in the LUT and a longer searching time. For any given incident zenith angle, LUTs for the nearby two integral zenith angles are used to determine the exit direction. Each LUT gives an exit direction, and the average is regarded as the final result. Photon weight and energy conservation The single scattering albedo of water is generally very small especially in NIR band. If only a single photon is traced at a time, most of the photons may be absorbed in water and cannot reach the sensor. In other words, these photons are wasted. So the photon packet is introduced in MC simulation of water to reduce the variance (Mobley, 1994). The method is to replace the above mentioned photon by a packet of many identical photons. The weight of a packet is usually set to be unit when it is generated. The weight of the packet is multiplied by the single scattering albedo of water every time it is scattered in water. Similarly, the weight of the packet is multiplied by the single scattering albedo of a leaf when it is intercepted by the leaf. This process continues until the packet reaches a sensor or its weight is smaller than a given threshold. When the latter happens, a roulette is taken to guarantee energy conservation.

Fig. A2. An example of the LUT for rough water surface with an incident zenith angle of 150°. The labels of x and y axes are the same in main figure and the subfigures. Arrow A in Zoom (a) indicates a rapid accumulation of exit energies, corresponding to the refracted glitter. Arrow B in Zoom (b) indicates a locally rapid accumulation, corresponding to the reflected glitter.

19

When a photon is intercepted by leaves or water surface, it may be reflected or transmitted if it is not absorbed by leaves. A photon packet can also interact with leaves and water surface for more than once. The algorithm will become much more complicated if the packet is split into two parts and they are traced separately. So the following methods are taken in the AVMC. (1) If the packet is intercepted by leaves, the behavior (reflected or transmitted) of the unabsorbed part is determined by the ratio rl/(rl + tl). The leaves are regarded as bi-Lambertian. Then the exit direction can be determined according to the incident direction and the normal of the leaf. (2) If the packet is intercepted by water surface, the exit direction is determined by the LUT described in Fig. A2. Absorption of water surface is neglected. So the weight of a packet will not change when it interacts with water surface. The LUT has guaranteed that the probability of the reflected or transmitted packet is consistent with the ratio of reflected or transmitted energy. Only water, leaves and bottom can absorb energy from packets in the AVMC. If they are set to be nonabsorbing, validation tests demonstrate that total weights of the packets exiting from the scene equal to that of the incident packets (1 million) under different incident zenith angles (180°, 150°, 120° and 100°). So energy conservation is fulfilled in the AVMC. References Abe, K., Komada, M., & Ookuma, A. (2008). Efficiency of removal of nitrogen, phosphorus, and zinc from domestic wastewater by a constructed wetland system in rural areas: A case study. Water Science and Technology, 58, 2427–2433. http://dx.doi.org/10. 2166/wst.2008.570. Ackleson, S. G. (2003). Light in shallow waters: A brief research review. Limnology and Oceanography, 48, 323–328. http://dx.doi.org/10.4319/lo.2003.48.1_part_2.0323. Ackleson, S. G., & Klemas, V. (1986). Two-flow simulation of the natural light field within a canopy of submerged aquatic plants. Applied Optics, 25, 1129–1136. http://dx.doi. org/10.1364/AO.25.001129. Adam, E., Mutanga, O., & Rugege, D. (2010). Multispectral and hyperspectral remote sensing for identification and mapping of wetland vegetation: A review. Wetlands Ecology and Management, 18, 281–296. http://dx.doi.org/10.1007/ s11273-009-9169-z. An, S., Tian, Z., Cai, Y., Wen, T., Xu, D., Jiang, H., et al. (2013). Wetlands of Northeast Asia and High Asia: An overview. Aquatic Sciences, 75, 63–71. http://dx.doi.org/10.1007/ s00027-012-0281-4. Babin, M., Stramski, D., Ferrari, G. M., Claustre, H., Bricaud, A., Obolensky, G., et al. (2003). Variations in the light absorption coefficients of phytoplankton, nonalgal particles, and dissolved organic matter in coastal waters around Europe. Journal of Geophysical Research, 108, 3211. http://dx.doi.org/10.1029/2001JC000882. Barducci, A., Guzzi, D., Marcoionni, P., & Pippi, I. (2009). Aerospace wetland monitoring by hyperspectral imaging sensors: A case study in the coastal zone of San Rossore Natural Park. Journal of Environmental Management, 90, 2278–2286. http://dx.doi.org/ 10.1016/j.jenvman.2007.06.033. Beget, M. E., Bettachini, V. A., Di Bella, C. M., & Baret, F. (2013). SAILHFlood: A radiative transfer model for flooded vegetation. Ecological Modelling, 257, 25–35. http://dx. doi.org/10.1016/j.ecolmodel.2013.02.025. Buiteveld, H., Hakvoort, J., & Donze, M. (1994). Optical properties of pure water. Ocean Optics XII (pp. 174–183). International Society for Optics and Photonics. Chen, J. M., Li, X., Nilson, T., & Strahler, A. (2000). Recent advances in geometrical optical modelling and its applications. Remote Sensing Reviews, 18, 227–262. http://dx.doi. org/10.1080/02757250009532391. Cooper, K., Smith, J. A., & Pitts, D. (1982). Reflectance of a vegetation canopy using the adding method. Applied Optics, 21, 4112–4118. http://dx.doi.org/10.1364/AO.21. 004112. Cox, C., & Munk, W. (1956). Slopes of the sea surface deduced from photographs of sun glitter. Bulletin of the Scripps Institution of Oceanography, 6, 401–488 (https:// escholarship.org/uc/item/1p202179). D'Alimonte, D., Zibordi, G., Kajiyama, T., & Cunha, J. C. (2010). Monte Carlo code for high spatial resolution ocean color simulations. Applied Optics, 49, 4936–4950. http://dx. doi.org/10.1364/AO.49.004936. Disney, M. I., Lewis, P., & North, P. (2000). Monte Carlo ray tracing in optical canopy reflectance modelling. Remote Sensing Reviews, 18, 163–196. http://dx.doi.org/10. 1080/02757250009532389. Feret, J., François, C., Asner, G. P., Gitelson, A. A., Martin, R. E., Bidel, L. P., et al. (2008). PROSPECT-4 and 5: Advances in the leaf optical properties model separating photosynthetic pigments. Remote Sensing of Environment, 112, 3030–3043. http://dx.doi. org/10.1016/j.rse.2008.02.012. Fujimoto, A., Tanaka, T., & Iwata, K. (1986). Arts: Accelerated ray-tracing system. Computer Graphics and Applications, IEEE, 6, 16–26. http://dx.doi.org/10.1109/MCG. 1986.276715. Gilmore, M. S., Wilson, E. H., Barrett, N., Civco, D. L., Prisloe, S., Hurd, J. D., et al. (2008). Integrating multi-temporal spectral and structural information to map wetland vegetation in a lower Connecticut River tidal marsh. Remote Sensing of Environment, 112, 4048–4060. http://dx.doi.org/10.1016/j.rse.2008.05.020.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015

20

G. Zhou et al. / Remote Sensing of Environment xxx (2015) xxx–xxx

Gimond, M. (2004). Description and verification of an aquatic optics Monte Carlo model. Environmental Modelling & Software, 19, 1065–1076. http://dx.doi.org/10.1016/j. envsoft.2003.11.010. Govaerts, Y. M., & Verstraete, M. M. (1998). Raytran: A Monte Carlo ray-tracing model to compute light scattering in three-dimensional heterogeneous media. IEEE Transactions on Geoscience and Remote Sensing, 36, 493–505. http://dx.doi.org/10. 1109/36.662732. Haltrin, V. I. (2006). Absorption and scattering of light in natural waters. In A. A. Kokhanovsky (Ed.), Light scattering reviews: Single and multiple light scattering (pp. 445–486). Springer Berlin Heidelberg. http://dx.doi.org/10.1007/3-54037672-0_10. He, X., Bai, Y., Zhu, Q., & Gong, F. (2010). A vector radiative transfer model of coupled ocean–atmosphere system using matrix-operator method for rough sea-surface. Journal of Quantitative Spectroscopy and Radiative Transfer, 111, 1426–1448. http:// dx.doi.org/10.1016/j.jqsrt.2010.02.014. He, B., Quan, X., & Xing, M. (2013). Retrieval of leaf area index in alpine wetlands using a two-layer canopy reflectance model. International Journal of Applied Earth Observation and Geoinformation, 21, 78–91. http://dx.doi.org/10.1016/j.jag.2012.08.014. Hedley, J. (2008). A three-dimensional radiative transfer model for shallow water environments. Optics Express, 16, 21887–21902. http://dx.doi.org/10.1364/OE.16.021887. Hedley, J., & Enríquez, S. (2010). Optical properties of canopies of the tropical seagrass Thalassia testudinum estimated by a three-dimensional radiative transfer model. Limnology and Oceanography, 55, 1537–1550. http://dx.doi.org/10.4319/lo.2010.55. 4.1537. Hosgood, B., Jacquemoud, S., Andreoli, G., Verdebout, J., Pedrini, G., & Schmuck, G. (1995). Leaf optical properties experiment 93 (LOPEX93) report EUR-16095-EN. Jacquemoud, S., Bacour, C., Poilve, H., & Frangi, J. (2000). Comparison of four radiative transfer models to simulate plant canopies reflectance: Direct and inverse mode. Remote Sensing of Environment, 74, 471–481. http://dx.doi.org/10.1016/S00344257(00)00139-5. Jacquemoud, S., & Baret, F. (1990). PROSPECT: A model of leaf optical properties spectra. Remote Sensing of Environment, 34, 75–91. http://dx.doi.org/10.1016/00344257(90)90100-Z. Jacquemoud, S., Verhoef, W., Baret, F., Bacour, C., Zarco-Tejada, P. J., Asner, G. P., et al. (2009). PROSPECT+ SAIL models: A review of use for vegetation characterization. Remote Sensing of Environment, 113, S56–S66. http://dx.doi.org/10.1016/j.rse.2008. 01.026. Janssen, H. (2013). Monte-Carlo based uncertainty analysis: Sampling efficiency and sampling convergence. Reliability Engineering & System Safety, 109, 123–132. http://dx. doi.org/10.1016/j.ress.2012.08.003. Kallel, A., Verhoef, W., Le Hégarat-Mascle, S., Ottlé, C., & Hubert-Moy, L. (2008). Canopy bidirectional reflectance calculation based on Adding method and SAIL formalism: AddingS/AddingSD. Remote Sensing of Environment, 112, 3639–3655. http://dx.doi. org/10.1016/j.rse.2008.05.014. Kay, S., Hedley, J. D., & Lavender, S. (2009). Sun glint correction of high and low spatial resolution images of aquatic scenes: A review of methods for visible and near-infrared wavelengths. Remote Sensing, 1, 697–730. http://dx.doi.org/10.3390/rs1040697. Kearney, M. S., Stutzer, D., Turpie, K., & Stevenson, J. C. (2009). The effects of tidal inundation on the reflectance characteristics of coastal marsh vegetation. Journal of Coastal Research, 25, 1177–1186. http://dx.doi.org/10.2112/08-1080.1. Klemas, V. (2011). Remote sensing of wetlands: Case studies comparing practical techniques. Journal of Coastal Research, 27, 418–427. http://dx.doi.org/10.2112/ JCOASTRES-D-10-00174.1. Klemas, V. (2013). Remote sensing of emergent and submerged wetlands: An overview. International Journal of Remote Sensing, 34, 6286–6320. http://dx.doi.org/10.1080/ 01431161.2013.800656. Kuusk, A. (2001). A two-layer canopy reflectance model. Journal of Quantitative Spectroscopy and Radiative Transfer, 71, 1–9. http://dx.doi.org/10.1016/S00224073(01)00007-3. Leathers, R. A., Downes, T. V., Davis, C. O., & Mobley, C. D. (2004). Monte Carlo radiative transfer simulations for ocean optics: A practical guide. http://oai.dtic.mil/oai/oai? verb=getRecord&metadataPrefix=html&identifier=ADA426624 Lee, Z., Carder, K. L., Mobley, C. D., Steward, R. G., & Patch, J. S. (1998). Hyperspectral remote sensing for shallow waters. I. A semianalytical model. Applied Optics, 37, 6329–6338. http://dx.doi.org/10.1364/AO.37.006329. Lyzenga, D. R. (1977). Reflectance of a flat ocean in the limit of zero water depth. Applied Optics, 16, 282–283. http://dx.doi.org/10.1364/AO.16.000282. Mobley, C. D. (1994). Light and water: Radiative transfer in natural waters. San Diego: Academic press. Mobley, C. D., Zhang, H., & Voss, K. J. (2003). Effects of optically shallow bottoms on upwelling radiances: Bidirectional reflectance distribution function effects. Limnology and Oceanography, 48, 337–345. http://dx.doi.org/10.4319/lo.2003. 48.1_part_2.0337. Morel, A., & Prieur, L. (1977). Analysis of variations in ocean color. Limnology and Oceanography, 22, 709–722. http://dx.doi.org/10.4319/lo.1977.22.4.0709.

Nicodemus, F. E., Richmond, J. C., Hsia, J. J., Ginsberg, I. W., & Limperis, T. (1977). Geometrical considerations and nomenclature for reflectance. http://physics.nist.gov/ Divisions/Div844/facilities/specphoto/pdf/geoConsid.pdf Ottaviani, M., Spurr, R., Stamnes, K., Li, W., Su, W., & Wiscombe, W. (2008). Improving the description of sunglint for accurate prediction of remotely sensed radiances. Journal of Quantitative Spectroscopy and Radiative Transfer, 109, 2364–2375. http://dx.doi. org/10.1016/j.jqsrt.2008.05.012. Palmer, K. F., & Williams, D. (1974). Optical properties of water in the near infrared. Journal of the Optical Society of America, 64, 1107–1110. http://dx.doi.org/10.1364/ JOSA.64.001107. Pinty, B., Widlowski, J., Taberner, M., Gobron, N., Verstraete, M. M., Disney, M., et al. (2004). Radiation Transfer Model Intercomparison (RAMI) exercise: Results from the second phase. Journal of Geophysical Research: Atmospheres (1984–2012), 109. http://dx.doi.org/10.1029/2003JD004252. Preisendorfer, R. W. (1976). Hydrologic Optics. Vol VI. Surfaces. Honolulu: US Dept. of Commerce, National Oceanic and Atmospheric Administration, Environmental Research Laboratories, Pacific Marine Environmental Laboratory. Proisy, C., Couteron, P., & Fromard, F. (2007). Predicting and mapping mangrove biomass from canopy grain analysis using Fourier-based textural ordination of IKONOS images. Remote Sensing of Environment, 109, 379–392. http://dx.doi.org/10.1016/j.rse. 2007.01.009. Qin, W., & Liang, S. (2000). Plane-parallel canopy radiation transfer modeling: Recent advances and future directions. Remote Sensing Reviews, 18, 281–305. http://dx.doi.org/ 10.1080/02757250009532393. Rahman, H., Pinty, B., & Verstraete, M. M. (1993). Coupled surface–atmosphere reflectance (CSAR) model: 2. Semiempirical surface model usable with NOAA advanced very high resolution radiometer data. Journal of Geophysical Research, [Atmospheres], 98, 20791–20801. http://dx.doi.org/10.1029/93JD02072. Schmidt, K. S., & Skidmore, A. K. (2003). Spectral discrimination of vegetation types in a coastal wetland. Remote Sensing of Environment, 85, 92–108. http://dx.doi.org/10. 1016/S0034-4257(02)00196-7. Shirley, P., & Morley, R. K. (2003). Realistic ray tracing (2nd ed.). Natick: AK Peters, Ltd. Smith, R. C., & Baker, K. S. (1981). Optical properties of the clearest natural waters (200– 800 nm). Applied Optics, 20, 177–184. http://dx.doi.org/10.1364/AO.20.000177. Staehr, P. A., & Markager, S. (2004). Parameterization of the chlorophyll a-specific in vivo light absorption coefficient covering estuarine, coastal and oceanic waters. International Journal of Remote Sensing, 25, 5117–5130. http://dx.doi.org/10.1080/ 01431160410001716932. Suits, G. H. (1972). The calculation of the directional reflectance of a vegetative canopy. Remote Sensing of Environment, 2, 117–125. http://dx.doi.org/10.1016/00344257(71)90085-X. Suits, G. H. (1984). A versatile directional reflectance model for natural water bodies, submerged objects, and moist beach sands. Remote Sensing of Environment, 16, 143–156. http://dx.doi.org/10.1016/0034-4257(84)90058-0. Turpie, K. R. (2012). Enhancement of a canopy reflectance model for understanding the specular and spectral effects of an aquatic background in an inundated tidal marsh. (PhD Dissertation) Maryland: University of Maryland, College Park. Turpie, K. R. (2013). Explaining the spectral red-edge features of inundated marsh vegetation. Journal of Coastal Research, 29, 1111–1117. http://dx.doi.org/10.2112/ JCOASTRES-D-12-00209.1. Vanderbilt, V. C., Perry, G. L., Livingston, G. P., Ustin, S. L., Barrios, D., BrEOn, F., et al. (2002). Inundation discriminated using sun glint. IEEE Transactions on Geoscience and Remote Sensing, 40, 1279–1287. http://dx.doi.org/10.1109/TGRS.2002.800233. Vella, D., & Mahadevan, L. (2005). The “Cheerios effect”. American Journal of Physics, 73, 817–825. http://dx.doi.org/10.1119/1.1898523. Verhoef, W. (1984). Light scattering by leaf layers with application to canopy reflectance modeling: The SAIL model. Remote Sensing of Environment, 16, 125–141. http://dx.doi. org/10.1016/0034-4257(84)90057-9. Verhoef, W. (1985). Earth observation modeling based on layer scattering matrices. Remote Sensing of Environment, 17, 165–178. http://dx.doi.org/10.1016/00344257(85)90072-0. Verhoef, W. (1998). Theory of radiative transfer models applied in optical remote sensing of vegetation canopies. (PhD Dissertation) Wageningen: Agricultural University. Verhoef, W., & Bach, H. (2007). Coupled soil–leaf–canopy and atmosphere radiative transfer modeling to simulate hyperspectral multi-angular surface reflectance and TOA radiance data. Remote Sensing of Environment, 109, 166–182. http://dx.doi.org/10.1016/ j.rse.2006.12.013. Volpe, V., Silvestri, S., & Marani, M. (2011). Remote sensing retrieval of suspended sediment concentration in shallow waters. Remote Sensing of Environment, 115, 44–54. http://dx.doi.org/10.1016/j.rse.2010.07.013. Zimmerman, R. C. (2003). A biooptical model of irradiance distribution and photosynthesis in seagrass canopies. Limnology and Oceanography, 48, 568–585. http://dx.doi.org/ 10.4319/lo.2003.48.1_part_2.0568.

Please cite this article as: Zhou, G., et al., Canopy modeling of aquatic vegetation: A radiative transfer approach, Remote Sensing of Environment (2015), http://dx.doi.org/10.1016/j.rse.2015.03.015