1. Introduction

Any ground-truth study in optical remote sensing of the ocean or a study comparing in-water radiative transfer calculations with measured water leaving radiance needs a formula linking inherent optical properties (IOPs) of sea water such as absorption coefficient a and backscattering coefficient b _b with sea surface reflectance R. Such a formula, first introduced in this form by Morel and Prieur [1] but based on earlier functions (among them Refs. 2 & 3 which discuss even earlier formulations), has the following form:

where f is a coefficient which to a first order of approximation may be treated as dependent only on the shape of the scattering phase function and the angular distribution of downwelling radiance just below the sea surface [4]. The formula, universally used in ocean optics, is actually a surprisingly good approximation if one notes that it is based on a form of single scattering approximation (identification of upward scattering with backscattering becomes more defective with every order of scattering). For remote sensing reflectance (R _rs), the formula is analogous to (1) but with the coefficient f/Q in place of f, where Q is the ratio of upwelling irradiance to upwelling radiance (also to a first approximation a function of phase function shape but generally also of scattering to absorption ratio).

Of course, this approach works only in a vertically homogeneous sea. Any stratification calls for use of a depth weighting function. One such function was proposed in 1980 by Gordon and Clark [5]:

where z ₉₀ is the depth at which 90% of the surface downwelling irradiance had been attenuated and

Analogous formulas are used for depth averaging other remote sensing retrieved products like absorption or chlorophyll concentration values. The symbol “2K” is usually given as doubled diffuse attenuation of downwelling radiation (2 K _d) but more correctly it should represent K _u+K _d (according to Ref. 5 the resulting error “is not serious in the present application”). This makes the weighting function proportional to round-trip irradiance transmittance to and from a given depth (or “round-trip attenuation” [6]). In the case of R _rs, “2K” is usually replaced by K _d+c (where c is beam attenuation) [6] while for lidar remote sensing of the ocean, the practice is to use 2 c (which is a direct application of the the well known lidar equation [7]).

Recently Zaneveld et al. [6] proposed a modified form of the depth weighting, using the first derivative of “round-trip attenuation” as the weighting function:

where

Formula (4) is “self normalizing” as long as it is integrated over the whole “infinitely” deep ocean which explains the change of the integral upper limit and the lack of the denominator integral necessary in formula (2).

 

Fig. 1. A simple calculation illustration of the Gordon and Clark (3) and Zaneveld et al. (5) weighting functions for an optically active submerged layer (represented by increase of K _d between 5 and 14 m depths)

Download Full Size | PDF

The motivation for looking for a better weighting function was errors in reflectance values calculated with formula (2) for stratified oceans [8]. This made Zaneveld et al. look at the problem again and mathematically derive an improved formula. Their depth weighting formula (5), unlike formula (3) is not a monotonously decreasing function. The qualitative difference of the functions is shown in Fig. 1 showing both functions for an optically active submerged layer (represented simply by increase of K _d). Function (5) has a maximum at the top of the active layer while function (3) overestimates the optically inactive surface layer making the resulting product overdependent on the optical properties of the surface layer.

Zaneveld et al. did not present any calculation examples showing the impact of the weighting approach on reflectance values which may explain why their paper has largely gone unnoticed in the community. Therefore we decided to calculate, using a radiance transfer model, the reflectances (both R and R _rs) for two cases of submerged layers: bubble and chlorophyll rich. We also decided to test how different “K” choices in the weighting functions and the method of their estimation (K _d and K _u are never directly measured) influence the computed reflectance values.

2. Methods

The study used a Monte Carlo radiative transfer code previously used for studying ocean optics phenomena [9,10] and instrumentation [11,12]. The spectral IOPs for the homogeneous ocean case and for the background water in the stratified runs were calculated for chlorophyll concentration Chl = 0.05 mg m^-3 according to formulas from Ref. 8 (default Hydrolight formulas). However, unlike Stramska and Stramski, we used stepwise layers to decrease the complication and run times. A chlorophyll rich (Chl = 3 mg m^-3) layer between 20 and 40 m depths represented a case of a high absorption layer, while a bubble layer was simulated between 10 m and 20 m with b _bub = 1.5 m^-1 of additional scattering and with the “bub43” phase function from Ref. 12. We do not claim such fully submerged bubble clouds are a frequent phenomenon. However it certainly is a good “extreme” test for the weighting functions.

In all Monte Carlo runs, incoming irradiation was modeled with the sun at 30 deg zenith angle and with 30% of uniform diffused irradiation (“sky glow”). Sea surface was modeled with the Cox-Munk parameterization [13] for wind speed 5 m s^-1. We did all runs for 24 wavelengths placed quasi-uniformly (to cover all remote sensing relevant values) between 300 and 800 nm. For each wavelength, 10⁸ virtual “photons” were traced. No inelastic effects (Raman scattering or fluorescence) were included.

Because the f and f/Q coefficients are IOP independent only to a first approximation [4], we calculated them for each wavelength from the Monte Carlo results for the homogeneous ocean (we could not assume a priori any form of depth weighting function while testing them).

Both K _d and K _u were calculated from three approximate equations: from the definition but using finite differences K _d1 = -ΔE _d/Δz with irradiances calculated every Δz = 2 m; from IOPs: K _d2 = a/$\overline{\mu }$ _d and K _d3 = (a+b _b)/$\overline{\mu }$ _d where $\overline{\mu }$ d is the average cosine of downwelling irradiance (derived from the Monte Carlo calculations). Analogous approximate values K _u1, K _u2 and K _u3 (using average cosine of upwelling irradiance) were also calculated for K _u.

3. Results and discussion

Figure 2 presents the Monte Carlo calculated spectral R compared to reflectance values calculated using formulas (3) and (5) for the chlorophyll rich case. It can be seen that the Gordon and Clark weighting is almost blind to the chlorophyll layer (the calculated reflectances are much closer to the “no layer” homogeneous ocean case - practically indistinguishable from it - than to the correct “layer” value). The Zaneveld et al. functions are much closer to the Monte Carlo calculated “layer” value, with the K _d2+K _u2 version being an outlier. The K _d1+K _u1 version, which is based on the diffuse attenuation definition, has problems towards longer wavelengths, probably because the spatial step (2 m) of irradiance sampling is too large. The 2K _d variant (not shown) had in all tested cases values similar to K _d+K _u confirming the above quoted hypothesis of Gordon and Clark.

 

Fig. 2. Monte Carlo calculated spectral R for the “no layer” case and for a chlorophyll rich layer compared to calculated directly using Gordon & Clark (3) and Zaneveld et al. (5) functions. See text for details.

Download Full Size | PDF

Figure 3 shows remote sensing reflectance results for the same case. This figure corresponds to the second row of Fig. 2 in Ref. 8 (our “layer” case is close to the middle of the case range shown there with dotted lines). The reflectance overestimation within the blue region and underestimation for shorter wavelengths, comparing to Ref 8, is caused by not including inelastic processes in our study. However, this does not impact on the weighting function discussion. Actually our slightly unrealistic case makes it easier to compare formulas which were not derived for inelastic processes. The Fig. 3 results explain why Stramska and Stramski had problems using the Gordon-Clark weighting function. For R _rs all its variants are practically blind for the 30 m centered chlorophyll layer, as in the case of reflectance. The Zaneveld et al. functions give values closer to the calculated “layer” R _rs. It is interesting that all K _u+K _d versions fare better than the K _d+c ones (the best one being K _u2+K _d2). This is not so strange if one notices that c for radiance is not the counterpart of K _d for irradiance. A component of beam attenuation is forward scattering which does not usually strongly affect upwelling radiance as the upwelling light field is normally highly diffused.; its result is rather more diffusing than attenuation. Therefore K _u may be indeed a better proxy for an “upwelling radiance attenuation” than c.

 

Fig. 3. Monte Carlo calculated spectral R _sr for the “no layer” case and for a chlorophyll rich layer compared to calculated directly using Gordon & Clark (3) and Zaneveld et al. (5) functions. See text for details.

Download Full Size | PDF

One can expect that the massive multiple scattering inside a bubble layer will affect the possibility of retrieving depth averaged IOPs even more than the case of absorbing chlorophyll layer. Figure 4 confirms the difficulty. While the formula (5) variants seem on average better, there is large variability between them. One of the Gordon and Clark variants (K _d2+c) did not notice the bubble cloud at all (which is not surprising because attenuation K _d2 is not scattering dependent). Most of the other variants of both the weighting functions have similar values but the Zaneveld et al. ones seem to better reproduce the spectral shape of the “layer” Monte Carlo calculated spectral R _rs. One of the formula (5) curves even overshoots it. Because Monte Carlo scales well with the scattering order, we believe this must be the result of error in using formula (1) in such a massively scattering case. The spectral values of f/Q coefficient we used were calculated for much less scattering homogeneous ocean and are probably unfit for this case. In this case using c for upwelling radiance attenuation is even more dangerous because of the highly forward scattering phase functions (in bubble clouds up to 80% scattering happens at angles smaller than 1 deg [12]). Similarly in the case of a remote sensing lidar (not tested), a “2 c” based weighting function will also overestimate the influence of the bubble cloud because of the high amount of forward scattering (as shown in Fig. 4 of Gordon 1982 [14] for the highest values of single scattering albedo).

 

Fig. 4. Monte Carlo calculated spectral R _sr for the “no layer” case and for a highly scattering bubble layer compared to calculated directly using Gordon & Clark (3) and Zaneveld et al. (5) functions. See text for details.

Download Full Size | PDF

4. Conclusions

Monte Carlo modeling results confirmed that the Gordon and Clark depth weighting function is inadequate for calculating reflectance for the stratified ocean, giving too much weight to the surface layer. It should be replaced with the Zaneveld et al. weighting in all applications relating the IOPs or water component concentrations to reflectance or remote sensing reflectance values. This is true both in case of chlorophyll rich absorbing layers and highly scattering bubble clouds.

Our model results show that using 2K _d instead of K _d+Ku in the weighting functions leads to only minimal errors (<1% for the chlorophyll rich layer) as predicted by Gordon and Clark [5]. Calculating the diffuse attenuation values from irradiance profiles gives better results than from IOPs as long as the spatial sampling density is adequate. In the case of remote sensing reflectance of a stratified ocean, use of K _d+c in the depth weighting function leads to worse results than K _d+K _u, most probably because of the forward scattering component of beam attenuation.