1. Introduction

The average cosine $\overline{\mu }$ of the underwater light field is a single parameter describing the angular distribution of light at a particular depth, defined as:

The average cosine of the oceanic light field in the asymptotic regime, ${\overline{\mu }}_{\infty }$, is an inherent optical property (IOP), i.e., it is not dependent on the shape of the incident light field but on the absorption and scattering properties of the water column [2–6]. Asymptotic theory is based on the principle that the shape of the light field with depth gradually transforms from being dependent on the incident surface light field to being constant and azimuthally symmetric (so L is only a function of θ). In the asymptotic regime, the diffuse attenuation coefficient for irradiance reaches an asymptotic value $K_{\infty }$, where ${\overline{\mu }}_{\infty }=\frac{\text{a}}{{\text{K}}_{\infty }}$, and a is the absorption coefficient. This is the asymptotic representation of Gershun’s equation, obtained by integrating the radiative transfer equation in all directions. Attenuation coefficients for all aspects of the light field, i.e., for all radiances and therefore all irradiances, are equivalent in the asymptotic regime and are also IOPs [7,8].

Absorption and scattering control ${\overline{\mu }}_{\infty }$ in a complex manner. Absorption elongates the radiance distribution with depth about the vertical direction because absorption increases with longer pathlength; the vertical direction is the shortest path and thus experiences the least effect from absorption. Scattering conversely disperses light in all directions, broadening the radiance distribution and decreasing${\overline{\mu }}_{\infty }$. If scattering were isotropic, i.e., equal in all directions, $\overline{\mu }$would approach 0 with no absorption. However, the angular distribution of scattering in natural waters is strongly forward-peaked, which reduces such broadening.

Light fields in the asymptotic regime have received considerable attention, particularly in ocean optics literature from 1970’s to 1990’s (see reviews by [6] and [9]), because of the direct link ${\overline{\mu }}_{\infty }$ and $K_{\infty }$ provide between AOPs and IOPs. In a series of papers, Zaneveld [2,3,10], developed radiative transfer approximations (RTAs) and an inversion to derive IOPs from radiometric measurements in the asymptotic regime that was adapted in [2] to an RTA for remote sensing reflectance R_RS ( = L_u/E_d) in surface waters in terms of the IOPs. In this latter case, surface water IOPs define a conceptual ${\overline{\mu }}_{\infty }$ extrapolated to a hypothetically infinite, homogeneous water column. Even though the light field in surface waters is far from asymptotic, a critical parameter for remote sensing, the attenuation coefficient for upwelling nadir radiance K_Lu in surface waters, was hypothesized by Zaneveld [2] to be approximately equivalent to $K_{\infty }$ because of the decoupling of upwelling radiances to downwelling radiance distributions. Note this would not apply in actual light fields in deep waters where influence from Raman inelastic scattering becomes apparent in the green and red [11].

Previous attempts to represent ${\overline{\mu }}_{\infty }$ in terms of IOPs have focused on the single scattering albedo ω ( = b/c), where b is total scattering and c is attenuation [2,3,5,6,12–14]. While convenient in that the ω range is bounded between 0 and 1, a drawback of this approach in the context of RTAs for ocean color remote sensing applications is b and c are not parameters that are readily derived from R_RS (e.g., [15]). For algorithm testing and validation, these parameters are also difficult to directly measure without large errors because typical collimated collectors for c measurements have significant acceptance of scattered light at near-forward angles [16,17]. Total scattering b is normally derived from c – a, with the current convention of measuring in situ a and c with WET Labs ac devices [18]. The c measurement error is typically >25% for a device with ~1° acceptance [16,17] and there is no easily applied correction without high quality angular scattering measurements in the near-forward direction.

Gordon et al. [19] have shown, however, that R_RS can be closely related to b_b / (a + b_b), where b_b is the backscattering coefficient representing integrated scattering in the backward hemisphere. Furthermore, in cases where a is much greater than b_b, the relationship can be simplified to b_b/a. It can be easily shown that b/c = b_b / (aB + b_b) where B is the backscattering ratio, or b_b/b, describing the general shape of the volume scattering function (VSF) [20,21]. Therefore b_b / (a + b_b) and b_b/a may have good potential to also describe${\overline{\mu }}_{\infty }$, with the advantage that these relationships 1) are directly applicable to remote sensing RTAs and 2) may be directly measured with accuracies within a few percent with commercially in situ available instrumentation [22,23]. Such relationships may have potential for developing novel ocean color inversion algorithms for next generation imagers, such as the Plankton, Aerosol, Cloud, and ocean Ecosystem (PACE) instrument with scheduled launch in 2022, which may also include a polarimeter with angular scattering information (http://pace.gsfc.nasa.gov/). With this application in mind, the relationship between ${\overline{\mu }}_{\infty }$ and b_b/a is investigated across a representative range of VSF shapes.

2. RT simulations

Hydrolight (HL; www.sequoiascientific.com; see [1]) version 5.2.2 was used to solve the radiative transfer equation. HL obtains ${\overline{\mu }}_{\infty }$ and $K_{\infty }$ by solving the following [12,24]:

All computations here were performed in the mid-visible at 532 nm. The interest here is in application of asymptotic light field theory in RTAs to use in remote sensing RTAs and inversion algorithms in the surface ocean. As such, we do not include Raman inelastic effects in the following. The effects of Raman on mesopelagic light fields have been rigorously addressed in other works [11,25]. Results are applicable to other wavelengths once water contributions are taken into account.

To compile the input matrix for IOPs, b_b was first set to 0.001, 0.0013, 0.0020, 0.0040, 0.01, 0.1 m⁻¹. Since backscattering of pure seawater b_bw(532) at 20°C and 35 PSU is 0.00098 m⁻¹ [23], these values correspond to η_bb = b_bw / b_b = 0.98, 0.75, 0.49, 0.25, 0.098, 0.0098, where η_bb represents the fractional seawater content of backscattering (cf. Morel and Gentili [26]). While in oceanic waters the contribution of pure seawater to total scattering is negligible, the seawater contribution is usually significant in the backward direction [27,28]. Values for b_b/a were then set from 0.001 to 1 in twenty, logarithmically equal increments. Total absorption was calculated from b_b/(b_b/a) for each permutation. Pure water absorption was set to 0.0441 m⁻¹ at 532 nm after [29].

Fournier-Forand (FF) analytical phase functions for the particulate fraction (β_p(θ)/b_p [30,31], latest version in [32]) were computed from bulk refractive index n_p, ranging from 1.04 to 1.24 in 0.02 steps, and particle size distribution (PSD) slope γ, ranging from 3.1 to 4 in 0.1 steps. A subset of aggregate phase functions (particles and seawater) used in HL are plotted in Fig. 1. See [33] for description of discretization of the phase functions and input in HL.

 

Fig. 1 Selected total phase functions used in radiative transfer computations, (A) varying fractional water content for backscattering η_bb, keeping n_p and γ fixed at 1.1 and 3.4, respectively, for the FF phase function representing the particulate fraction, and (B) varying FF phase function input parameters (n_p, γ in legend) while keeping η_bb fixed at 0.0098 (i.e., particle dominated).

Download Full Size | PDF

For computations of vertical profiles of K, simple boundary conditions were chosen in the model runs, i.e. 1 m/s wind speed, a semi-empirical sky model based on RADTRAN-X, with 0% cloud coverage. The ocean was assumed to be infinitely deep and optically homogeneous. No other boundary conditions were considered because the asymptotic light field depends only on the IOPs.

3. Results

Figure 2 shows ${\overline{\mu }}_{\infty }$ as a function of inputs (n_p and γ) for the FF phase function and b_b/a for the extreme values of η_bb considered, 0.98 and 0.0098. The largest b_b/a assessed for η_bb = 0.98 was 0.0227, as b_b was fixed at 0.001 m⁻¹ and a could not be less than a_w = 0.0441 m⁻¹. At low b_b/a values, the impact of absorption was dominant for all cases and ${\overline{\mu }}_{\infty }$ tended toward 1, with negligible effect from phase function shape. As b_b/a increased, ${\overline{\mu }}_{\infty }$decreased as the effects of scattering became more significant. This effect became more significant for phase functions with lower n_p and γ, i.e., for “softer” particle fields with higher relative fractions of large particles. However, in general, phase function shape had minimal impact on${\overline{\mu }}_{\infty }$. This has also been noted by Gordon et al. [6] who used Henyey-Greenstein phase function shapes.

 

Fig. 2 Asymptotic average cosine ${\overline{\mu }}_{\infty }$ as a function of refractive index n_p and PSD slope γ inputs to Fournier-Forand phase function model. Horizontal grayscale surfaces represent fixed b_b/a values (smallest and largest labeled with arrows at right side of each graph). Note 0.0227 was the largest b_b/a possible at η_bb = 0.98 for the input IOP matrix (see text). Subset of results are shown for fractional seawater backscattering η_bb = 0.98 (left panels; b_b fixed at 0.001 m⁻¹) and 0.0098 (right panels; b_b fixed at 0.1 m⁻¹). Colored vertical surfaces in upper panels correspond to fixed n_p panels below.

Download Full Size | PDF

Results from all IOP matrix inputs are plotted in Fig. 3(a), with phase function shape variability manifested as the spread of points at any specific b_b/a for a specific η_bb (see Fig. 3(a) legend). A 2nd order polynomial provided a satisfactory fit $\left[{\overline{\mu }}_{\infty }={\displaystyle \sum }_{n=0}^2{p}_n{\left(\text{Log}\frac{b_b}{a}\right)}^n\right]$ to individual η_bb relationships (Table 1). With decreasing η_bb, absolute error of the fit increased slightly due to the increasing impact of phase function variability.

 

Fig. 3 (a) All results for0${\overline{\mu }}_{\infty }$ plotted with respect to b_b/a, including full range of phase function shapes. A curve derived from the relationship obtained by Berwald et al. [5] is also shown (see text). (b) 4th order polynomial fits to all data and for η_bb ≤ 0.49. Vertical lines show maximum error for each b_b/a for each fit.

Download Full Size | PDF

Table 1. Fitting statistics for the average cosine for each η_bb as a function of b_b/a after 2nd order polynomial fits, including full range of phase function shapes.

View Table | View all tables in this article

Good agreement was observed for results with η_bb ≤ 0.49, representative of coastal waters and the majority of surface waters in the open ocean [27,28]. Combining all results and results for η_bb ≤ 0.49, 4th order polynomial fits $\left[{\overline{\mu }}_{\infty }={\displaystyle \sum }_{n=0}^4{p}_n{\left(\text{Log}\frac{b_b}{a}\right)}^n\right]$were satisfactory (Fig. 3(b)), with fitting parameters provided in Table 2. Percent absolute error %δ for the fit to all data was a reasonable 3.4%. Residuals of ${\overline{\mu }}_{\infty }$ to fitted values were approximately normally distributed with standard deviation of 0.04 (not shown), indicating insignificant contribution of outliers and a highly representative relationship. A previous relationship between ${\overline{\mu }}_{\infty }$ and ω found by Berwald et al. [5] could also be overlaid in Fig. 3 by knowing B = 0.0183 for the Petzold phase function used in their HL simulations, where b_b/a = 0.0183 / (ω⁻¹ – 1).

Table 2. Fitting statistics for the average cosine for combined η_bb relationships as a function of b_b/a after 4th order polynomial fits, including full range of phase function shapes.

View Table | View all tables in this article

A further relationship was also developed that explicitly included the η_bb parameter for all data, which could be useful if η_bb could be known or estimated a priori. For example, estimation of η_bb from chlorophyll concentration derived from ocean color R_RS would be comparable to the current convention for choosing a bidirectional reflectance distribution function (BRDF) correction to retrieve normalized water-leaving radiance from satellite imagery, following [34]. The resulting fit for ${\overline{\mu }}_{\infty }$ for all results was as follows:

4. Discussion

A key motivation of this work was to assess the dependence of ${\overline{\mu }}_{\infty }$ on b_b/a for use in RTAs for ocean color R_RS algorithms such as those proposed by Zaneveld [2]. Results in Table 2 show ${\overline{\mu }}_{\infty }$ is a strong function of b_b/a with minimal dependence on variability in the shape of phase functions for representative particle fields. Varying contribution of molecular scattering from pure seawater for different turbidity levels complicates the relationship, but %δ was still only 3.4% when considering the full range of turbidities, i.e., the full range of η_bb. If η_bb can be estimated a priori [34], then errors may be further reduced, particularly for very clear waters (Table 1). If remotely sensed multi-angle polarimetry becomes available, as is proposed for the future PACE mission, then it may be possible to glean enough phase function shape information with that sensor to directly estimate η_bb, as the angular shapes of the water and particle contributions in the backward trend in opposite directions (e.g [35].). For the majority of waters, however, the strong agreement in results for η_bb ≤ 0.49 will be representative.

HL computation results here also allow assessment of a central assumption in the Zaneveld approximations, namely that K_Lu near the surface is approximately equivalent to$K_{\infty }$. Zaneveld [2] argued these parameters should be in close agreement based on measurements made in the 1960’s and 1970’s of upward radiance fields, showing a near constant shape and attenuation rate with depth [9,36–38]. Figure 4(a) shows profiles of the % difference between K_Lu just below the surface (0^-) and $K_{\infty }$ for a relatively turbid case (η_bb = 0.098), for a range of b_b/a and particle phase function shapes. For these phase functions, B varied between 0.5% and 2%. Over the full range of variables tested, maximal percent difference was within 3% for the majority of conditions, with the sun at zenith. Another recent study by [11], also employing HL computations, found similar results in the asymptotic regime for water columns with vertical structure in IOPs, but only in the blue spectral region, as their work included the effects of inelastic Raman scattering. Figure 4(b) shows the effects of varying solar zenith angle θ_z for a range of η_bb, b_b/a and particle phase function shapes. As θ_z increased, % differences also increased, but in a predictable way for all cases (4th order polynomial fit provided in Fig. 4(b)). Good predictability implies K_Lu can be directly tied to $K_{\infty }$ and thus to b_b and a with low uncertainties.

 

Fig. 4 Percent difference between the attenuation of upwelling nadir radiance K_Lu and asymptotic attenuation $K_{\infty }$ (a) for η_bb = 0.098 and γ = 3.4, with a range of n_p and b_b/a given in legend; (b) as a function of θ_z for a range in η_bb, just below the air-sea interface (0^-), with symbol colors matching legend in (a). Fit provided for η_bb = 0.098 (solid black curve). The physical approximation from Eq. (4) is also plotted (dashed red curve).

Download Full Size | PDF

The general relationship found in Fig. 4(b) has direct physical meaning. We may assume K_d and K_Lu may be similar, as nadir radiance is derived from backscattered downward irradiance. For a flat water surface, we may assume K_d(θ_w) ≈K_d(0°) / cos(θ_w), where θ_w is the solar zenith angle in the water, so that K_Lu(θ_w) ≈K_Lu(0°) / cos(θ_w). Setting θ_w = 0° (as in Fig. 4(a)), results in K_Lu(0°) ≈$K_{\infty }$, so that:

Results thus provide justification for the targeted application of $K_{\infty }$ and ${\overline{\mu }}_{\infty }$ in ocean color RTAs for the surface ocean. Assessing the performance of ocean color RTAs using the above representations of ${\overline{\mu }}_{\infty }$ in terms of b_b/a with associated uncertainties is the subject of other ongoing work.