Semi-analytical modeling and parameterization of particulates-in-water phase function for forward angles

Sanjay Kumar Sahu; Palanisamy Shanmugam

doi:10.1364/OE.23.022291

1. Introduction

Accurate modeling of the propagation of light in water, by solving the Radiative Transfer Equation, is important for many applications such as underwater environment monitoring, ecosystem and water quality assessment using remote sensing data, and underwater optical communication system development. In this approach, the medium is characterized by three parameters, namely, the absorption coefficient 'a', scattering coefficient 'b' and phase function $\tilde{β} (θ)$ . The phase function is the angular distribution of light intensity scattered by a particle at a given wavelength (here wavelength is ignored for brevity). It is defined as

\tilde{β} (θ) = \frac{β (θ)}{b} .

where

β (θ)

is the volume scattering function (VSF) which is the scattered intensity per unit incident irradiance per unit volume of water. Integrating β(θ) over all directions (solid angles) gives the total scattering coefficient,

b = 2 π \int_{0}^{π} β (θ) \sin θ d θ .

Equation (2) is often divided into two parts, forward scattering, 0 ≤ θ < π/2, and backward scattering, π/2 ≤ θ ≤ π. The angular distribution of scattered light intensity in the backward directions (π/2 ≤ θ ≤ π) is important for modeling the physical phenomena using satellite remote sensing, whereas the angular distribution of scattered light in the forward directions (0 ≤ θ < π/2) plays an important role in the problems involving multiple scattering and in implementing optical links for underwater communications and imaging systems. In particular, for under water optical communication, the channel capacity cannot be modeled without the exact information with regard to the forward volume scattering function.

There are two major factors that contribute to the values of the volume scattering function and hence phase function, i.e., (1) water itself and (2) particulate matter. The contribution due to pure seawater has been well studied theoretically by Morel in 1974 [1]. On the other hand, the theoretical study of particle contribution has not been made successfully because of the existence of a wide range of particles with different sizes, shapes and optical properties. From the observation point of view, the experimental measurement of VSF has been very difficult at very small and large angles because of the need of separating the original light beam and scattering from optical surfaces at these angles. Moreover, the magnitude of the scattered intensity typically increases by five or six orders of magnitude in going from θ = 90° to θ = 0.1° for a given natural water sample, and scattering at a given angle θ can vary by two orders of magnitude among the water samples. This implies that the dynamic range of an instrument measuring the VSF should be sufficiently high. The rapid change in VSF at small scattering angles requires very precise alignment of the optical elements, which is next to impossible on a rolling ship [2]. These are the reasons why the measurement of VSF and hence phase function have not been frequently made and very small number of measurement data have been reported in the literature. One among these data, which is most cited and used, was reported by Petzold in the early 70's [3] in San Diego Harbor. Later, Mobley et al. took average of three different types of waters from Petzold's data and named it as Petzold average particle phase function [2]. Many formulas have been suggested as the empirical fits to these measured data [4–6]. One among them is the Henyey-Greenstein phase function originally employed in studies related to astrophysics [4]. This phase function, which was rather popular because of its simplest form, was later used to estimate the scattering of light by oceanic particulates. However, it did not have scattering maxima at near forward and far backward angles, as of typical marine scatterers. To overcome this problem, Haltrin developed a phase function, following Kattawar, which consists of the weighted sum of two Henyey-Greenstein ones with maximum for forward and backscattering angles [6]. The advantage of these empirically fitted phase functions lies in their simplicity and analyticity, which allow some progress in treating the difficult multiple-scattering problems. However, the parameters used and the forms chosen for these empirical relations are not directly based on the physics of the problem. Therefore, information about physical parameters such as refractive index or particle size distribution cannot be extracted from these empirical relations. On the other hand, some work have been made on fitting experimental phase functions to the exact Mie scattering solutions numerically integrated over assumed or measured particle size distributions [7, 8]. These exact Mie solutions can be used to extract certain information about the physical parameters from the phase function, but they are purely numerical and not appropriate to use when dealing with light propagation problems involving multiple scattering. Another parameterized analytical phase function, developed to fit the Petzold average particle phase function, is the Fournier-Forand (FF) phase function [9]. This phase function, derived from anomalous diffraction theory with some set of approximations, was expressed as a function of two physically meaningful quantities, namely refractive index (n) and particle size distribution slope (ξ). Although the FF phase function is related to some physically meaningful quantities, it is based on some approximations and hence is away from the exact solution of scattering problem. The present work is devoted to model and parameterize the particulates-in-water phase function for forward angles using some physically meaningful quantities and the exact numerical solution obtained using Mie theory.

The significant particulate contribution to phase function can be modeled using the Mie theory if particles are assumed to be homogeneous spheres and the particle size distribution and refractive index values are known [10]. According to Holland et al., the spherical approximation of irregularly shaped oceanic particles does not give much error to the light scattering calculations for angles less than about 80° [11]. Thus, in this paper we present a parameterized phase function for forward angles (0.1° to 90°), derived by a semi analytical approach using Mie theory. An extensive simulation was performed on the basis of bulk refractive index 'n ' and particle size distribution slope 'ξ ' for the particle volume scattering function (VSF), particle scattering coefficient (b_p), particle backscattering coefficient (b_bp) and particle attenuation coefficient (c_p) (subscript 'p' on any symbol denotes the particle contribution to the total quantity). The values of 'n' varied from 1.04 to 1.3, with an increment of 0.02, and the values of 'ξ ' varied from 3 to 5, with an increment of 0.2. These ranges chosen for the values of 'n' and 'ξ ' are typical of oceanic environment [12, 13]. The wavelength was chosen to be 530 nm as the phase function we consider as a reference was measured at 530 nm by Petzold [3]. The particle distribution was taken to be Junge type (power law), typical of oceanic environment. The effect of variation in refractive index (relative to water), PSD slope, minimum and maximum particle sizes on backscattering ratio was also studied. In particular, the variation of backscattering ratio as a function of refractive index and PSD slope was modeled, which was later inverted to obtain the bulk refractive index from measured parameters (in situ). After modeling the backscattering ratio, the phase function was also modeled as a function of 'n' and 'ξ ' using the simulated data. For this work, the entire range of forward angles (0.1°-90°) is divided into two ranges, namely a small angle range 0.1°-5° and an intermediate angle range 5°-90°. This division boundary is decided by looking at the fact that all the phase function models cross each other near 5° as reported by Mobley et al. [14]. The phase function was then modeled for both the regions and the model equations were derived with different coefficients.

2. Basic Mie theory and simulation

The light scattering problem in ocean has some analytical solutions, according to Mie theory, with certain assumptions. Particulates in water should be (1) spherical, (2) homogenous, and (3) sufficiently diluted. It is traditional in light scattering theory to assume the shape of the marine particles as spherical. Particles with other shapes such as ellipsoid can also be handled, but in that case the scattering will also depend on the orientation of the particle and hence the scattering function must be averaged over all orientations. These complexities demand that the uniform spherical model must be employed to study the scattering of light by marine particles. On the other hand, composition of a particle relates to the light scattering properties through its refractive index, which may be a function of position within the particle itself. But it is very difficult to handle this intra-particulate variation of refractive index theoretically. Therefore, the homogeneity assumption is necessary for theoretical study of light scattering in underwater.

The interaction of the electromagnetic field with a sphere can be described in terms of a vector spherical harmonic expansion satisfying the Maxwell’s equations and has been well described in literature [10, 15]. With above mentioned assumptions and by applying the appropriate boundary conditions of the electric and magnetic fields on the surface of a sphere, the expressions for scattering efficiency (Q_bp), backscattering efficiency (Q_bbp) and extinction efficiency (Q_ext) can be derived. These are defined as follows,

Q {}_{b p}= \frac{2}{α^{2}} \sum_{n = 1}^{\infty} (2 n + 1) ({| a_{n} |}^{2} + {| b_{n} |}^{2}) .

Q {}_{e x t}= \frac{2}{α^{2}} \sum_{n = ​ ​ ​ 1}^{\infty} (2 n + 1) Re (a_{n} + b_{n}) .

Q {}_{b b p}= \frac{1}{α^{2}} {| \sum_{n = 1}^{\infty} (2 n + 1) {(- 1)}^{n} (a_{n} - b_{n}) |}^{2} .

where ' a_n ' and ' b_n' are Mie-scattering coefficients and ' α ' is the size parameter given by,

α = \frac{π D}{λ_{m}} .

'λ_m' is the wavelength of light in the host medium and 'D' is the diameter of the particle.

According to energy conservation theorem, the absorption efficiency 'Q_abs' can be calculated from 'Q_bp' and 'Q_ext' using the relation below [16]

Q_{e x t} = Q_{b p} + Q_{a b s} .

All infinite series can be truncated after n_max terms. For this number Bohren and Huffman [17] proposed the value

n_{\max} = α + 4 α^{\frac{1}{3}} + 1.

To calculate the detailed shape of the angular scattering pattern in the far field zone, one must calculate the scattering amplitude functions S₁ and S₂ given by,

S_{1} (\cos θ) = \sum_{n = 1}^{\infty} \frac{2 n + 1}{n (n + 1)} (a_{n} π_{n} + b_{n} τ_{n}) .

S_{2} (\cos θ) = \sum_{n = 1}^{\infty} \frac{2 n + 1}{n (n + 1)} (a_{n} τ_{n} + b_{n} π_{n}) .

where 'n' is an integer and 'θ ' is the angle between the incident electric field and scattering plane. The functions π_n and τ_n describe the angular scattering patterns of the spherical harmonics and are related to Legendre Polynomials. For an extensive and detailed description about Mie theory one can refer to Van de Hulst [10] or Bohren and Huffman [17].

With S₁ and S₂ values, the phase function can be found which is proportional to S₁₁(cosθ) and is given by

S_{11} (\cos θ) = \frac{1}{2} ({| S_{1} (\cos θ) |}^{2} + {| S_{2} (\cos θ) |}^{2}) .

However, the Eq. (11) is not the normalized phase function and the normalized phase function is given by

S_{11} (\cos θ) = \frac{({| S_{1} (\cos θ) |}^{2} + {| S_{2} (\cos θ) |}^{2})}{\sum_{n = 1}^{\infty} (2 n + 1) ({| a_{n} |}^{2} + {| b_{n} |}^{2})} .

For a poly-dispersion the optical properties can be calculated using the following relations,

b_{p} = \frac{λ^{3}}{8 π^{2}} \int_{α_{\min}}^{α_{\max}} Q_{b p} (α) f (α) α^{2} d α .

b_{b p} = \frac{λ^{3}}{8 π^{2}} \int_{α_{\min}}^{α_{\max}} Q_{b b p} (α) f (α) α^{2} d α .

B_{p} = \frac{\int_{α_{\min}}^{α_{\max}} Q_{b b p} (α) f (α) α^{2} d α}{\int_{α_{\min}}^{α_{\max}} Q_{b p} (α) f (α) α^{2} d α} .

\tilde{β} (θ) = \int_{α_{\min}}^{α_{\max}} S_{11} (\cos θ) f (α) α^{2} d α .

where, f(α) is the probability density function for the PSD. Since particle size distributions (PSD) are typically continuous distributions and Mie-scattering properties can only be calculated for a monodisperse system, the PSD is fractionated and Mie-scattering properties are calculated for each fraction. These Mie scattering properties are then combined with the weight for each fraction to derive the bulk optical properties. As the number of fractions is unknown and needs to be optimized for each calculation, the entire range of PSD is kept fractionated until the desired properties (b, b_b and β) converge to stable values as suggested by Aernouts et al. [18]. Since these PSD functions are often defined for sizes going from zero to infinity, the minimum and maximum size should be defined such that all the excluded sizes have no significant effect and the included size window is minimal to reduce the calculation time.

The simulation was carried out by choosing the wavelength λ = 530 nm and the refractive index variation from n = 1.04 to n = 1.3 with an increment of 0.2. The range of values of refractive index is typical for oceanic particles and chosen according to the previous work. The particle size distribution was taken to be Junge type, since the size distribution of the total particle suspension and of the pelagic organisms in aquatic ecosystem has been shown to be well represented by such a distribution [10, 12, 19]. The probability density function for Junge type distribution is given as

f (D) = K D^{- ξ} .

where, ' D ' is the particle diameter and

K = \frac{(ξ - 1)}{(D_{\min}^{1 - ξ} - D_{\max}^{1 - ξ})} .

D_min and D_max are the lower and upper limits of the size range under consideration and 'ξ ' is the PSD slope.

3. Model description

3.1 Modeling of the backscattering ratio

A model for calculating the backscattering ratio is derived from simulation results using the least square technique. It considers the two major factors responsible for light scattering, namely refractive index and PSD slope, as the inputs. The variation of backscattering ratio values with refractive index for some fixed values of PSDs is shown in Fig. 1(a). From this figure it is intuitive that the backscattering coefficient can be modeled as a function of 1/(n-1)² by a power law equation, for a particular value of PSD slope. The model equation satisfying this relationship between 'B_p' and 'n' for all the values of PSD slope in the range 3 to 5 is expressed as

B_{p} = P_{1} \cdot {(x)}^{P_{2}} .

where,

x = \frac{1}{{(n - 1)}^{2}} .

Clearly, the intermediate parameters P₁ and P₂ are functions of the PSD slope (Fig. 1 (b)). This variation is modeled using a least square technique, as expressed below.

P_{m} = a_{m} {(ξ - 3)}^{2} + b_{m} (ξ - 3) + c_{m} .

where, m = 1, 2 and 'ξ ' is the PSD slope. In fitting all these equations care was taken so that the value of R² (correlation coefficient) is nearly equal to one and the value of RMSE is nearly zero. The corresponding coefficients 'a_m', 'b_m' and 'c_m' are given in Table 1.

Fig. 1 (a) Variation of the Backscattering Ratio (B_P) with 1/(n-1)². (b) Variation of the intermediate parameters P₁ and P₂ with PSD slope (ξ).

Download Full Size | PDF

Table 1. Model parameters for calculating the backscattering ratio.

View Table | View all tables in this article

Equation (19) can be inverted to give the values of bulk refractive index using easily measurable parameters 'ξ ' and 'B_P '. The equation for predicting “n” is expressed as,

n = 1 + {(\frac{P_{1}}{B_{P}})}^{\frac{1}{2 P_{2}}} .

3.2 Modeling of the phase function

The main difference between the phase functions calculated using Mie theory and real phase functions of random natural grains of sea water lies in the fact that, oscillations are present in the former and disappear for the later. The main scattering peaks of these two are similar [20]. Also, oscillations present in the Mie scattering solution for a monodisperse particulate system start disappearing when a polydisperse particulate system (typical of oceanic environment) is considered [9]. For lower values of PSD the oscillation is more and it starts disappearing towards higher values of PSDs. Therefore, to include the entire range of PSD slopes in a single phase function model only forward angles are considered in the present work. Moreover, the spherical approximation of irregularly shaped oceanic particles in Mie theory gives good results for forward angles [11].

To model the phase function from 0.1° to 90°, the entire range is divided into two parts, one from 0.1° to 5° and the other from 5° to 90°. Such a division was necessary as the phase function is highly nonlinear with a change in slope around 5° [21] and it is very difficult to model the entire range by a single power equation. Also, all the existing phase function models cross each other near this angle [14]. Then the least square method is applied to fit the simulated data of phase function to derive model equations separately for both the regions.

For the range 0.1° to 5° the model equation as a function of PSD slope and refractive index is given below,

\log (\tilde{β} (θ)) = P_{1} {(\ln (θ))}^{2} + P_{2} (\ln (θ)) + P_{3} .

where,

P_{m} = a_{m} \exp (- x) + b_{m} (x) + c_{m}

and

\begin{array}{l} \begin{array}{l} a_{m} = \frac{d_{m}}{y^{2}} + e_{m} \sin (5 y) + f_{m} \\ b_{m} = \frac{h_{m}}{y^{2}} + i_{m} \sin (5 y) + j_{m} \\ c_{m} = \frac{k_{m}}{y^{2}} + l_{m} \sin (5 y) + o_{m} \end{array}] . \\ x = ξ - 3. \\ y = n - 1 a n d m = 1, 2, 3 \end{array}

For the range 5° to 90°, the model equation is given below,

\log (\tilde{β} (θ)) = P_{1} {(\ln (θ))}^{3} + P_{2} {(\ln (θ))}^{2} + P_{3} (\ln (θ)) + P_{4} .

where,

P_{m} = a_{m} \exp (- x) + b_{m} (x) + c_{m} .

and

\begin{array}{l} \begin{array}{l} a_{m} = \frac{d_{m}}{y^{2}} + e_{m} \sin (5 y) + f_{m} \\ b_{m} = \frac{h_{m}}{y^{2}} + i_{m} \sin (5 y) + j_{m} \\ c_{m} = \frac{k_{m}}{y^{2}} + l_{m} \sin (5 y) + o_{m} \end{array}] . \\ x = ξ - 3. \\ y = n - 1 a n d m = 1, 2, 3, 4 \end{array}

To model the phase function in a simple and useful way, the dependency of the phase function on scattering angle is first modeled by a power law equation for both the ranges 0.1° to 5° and 5° to 90° as given in Eq. (22) and Eq. (25) respectively. The intermediate parameters P_m are modeled as a function of the PSD slope and refractive index which are the major parameters influencing the scattering. The variation of these intermediate parameters with PSD slope is shown in Fig. 2 for the case of phase function in the 5° to 90° range. The similar variation of all the intermediate parameters with PSD resulted in a single equation describing the dependency. This dependency was fitted by using the least square method and is given by Eq. (26). The coefficients a_m, b_m and c_m of Eq. (26) are then modeled as a function of refractive index 'n' and is given in Eq. (27) (these figures are not shown for brevity). The same procedure is employed to model the phase function in the range 0.1° to 5°.

Fig. 2 Variation of the intermediate parameters P_m (P₁- P₄) with PSD slope values.

Download Full Size | PDF

The corresponding model coefficients for the ranges 0.1° to 5° and 5° to 90° are given in Tables 2 and 3 respectively. Since the scattering phase function is highly nonlinear and has significantly larger variations in the forward angles, a large number of coefficients are required to model this phase function. However, it should be noted that the coefficients d_m, h_m and k_m are very small values and the model is straightforward and less complex compared to the other phase function models. The number of coefficients with our phase function model can be decreased but at the expense of compromising the exact shape of the function.

Table 2. Model parameters for computing the phase function in the range 0.1ᵒ-5ᵒ.

View Table | View all tables in this article

Table 3. Model parameters for computing the phase function in the range 5°- 90°.

View Table | View all tables in this article

3.3 Model comparison

In the past, many phase function models have been used for the numerical calculation of light fields in underwater. Many of these models are basically derived by fitting empirically to the experimentally measured Petzold average particle phase function. Other models are numerically derived based on Mie theory calculations. Here, performance evaluation of the new phase function model is accomplished by comparing its results with the Petzold average particle phase function derived from experimental measurements and other phase function models – such as the Henyey-Greenstein Phase Function (OTHG), Haltrin's Two-Term Henyey-Greenstein Phase Function (TTHG) and Fournier-Forand (FF) Phase Function. Besides, the present model is also compared independently with the measured data from Sokolov et al [21]. These models are briefly described in what follows.

3.3.1 Petzold Average Particle Phase Function

This phase function was reported by Petzold in 1972 [3] using three measurements of VSF in San Diego Harbor. The values of the phase function are given in Mobley et al. [22] and Mobley [2]. Numerical integration of this phase function over 90° ≤ θ ≤ 180° gives the particle backscattering ratio of B_p = 0.0183.

3.3.2 Fournier-Forand Phase Function

Fournier-Forand derived an approximate analytical phase function for an ensemble of particles that have a hyperbolic (Junge type) distribution, with each particle scattering according to the anomalous diffraction approximation to the exact Mie theory [9]. This phase function in its latest form is given by

\begin{array}{l} {\tilde{β}}_{F F} (θ) = \frac{1}{4 π {(1 - δ)}^{2} δ^{v}} {[ν (1 - δ) - (1 - δ^{v})] + [δ (1 - δ^{v}) - ν (1 - δ)] \sin^{- 2} (\frac{θ}{2})} \\ + \frac{1 - δ^{v}_{180}}{(16 π (δ_{180} - 1) δ^{v}_{180})} (3 \cos^{2} θ - 1) . \end{array}

where,

ν = \frac{3 - ξ}{2} . and δ = \frac{4}{3 {(n - 1)}^{2}} \sin^{2} (\frac{θ}{2}) .

Here 'n' is the real index of refraction of particles, 'ξ ' is the PSD slope and δ₁₈₀ is δ evaluated at θ = 180ᵒ. Eq. (28) can be integrated to obtain the particle backscattering ratio

B_{P} = 1 - \frac{1 - δ_{90}^{ν + 1} - 0.5 (1 - δ_{90}^{ν})}{(1 - δ_{90}) δ_{90}^{ν}} .

Here δ₉₀ is δ evaluated at θ = 90ᵒ. It is clear from Eq. (30) that one can have several sets of PSD and refractive index values for a specific value of backscattering ratio.

3.3.3 One -Term Henyey-Greenstein phase function (OTHG)

It is a single parameter analytical phase function [4], which is frequently used in radiative transfer calculations in astrophysics, atmospheric science and ocean optics, because of its mathematical simplicity. The expressions for the phase function and the corresponding backscattering ratio are given below.

{\tilde{β}}_{O T H G} (θ) = \frac{1}{4 π} \frac{1 - g^{2}}{{(1 + g^{2} - 2 g \cdot \cos θ)}_{}^{3 / 2}} .

B_{P} = \frac{1 - g}{2 g} (\frac{1 - g}{\sqrt{1 + g^{2}}} - 1) .

The parameter 'g ' is the mean cosine of scattering angles and a value of g = 0.9185 gives B_p = 0.0183.

3.3.4 Haltrin's Two-Term Henyey-Greenstein Phase Function (TTHG)

Haltrin used the weighted sum of OTHG phase function to give better values in the near forward angles [6]. The expression for this phase function is as given below.

{\tilde{β}}_{T T H G} (θ) = α {\tilde{β}}_{O T H G} (θ, g_{1}) + (1 - α) {\tilde{β}}_{O T H G} (θ, g_{2}) .

The value of α, 0 ≤ α ≤ 1, gives the relative contributions of the two OTHG phase functions to the TTHG phase function. Haltrin developed a version of the TTHG in which the 'α' and 'g₂' parameters are given as functions of 'g₁' as follows,

g_{2} = - 0.30614 + 1.0006 g_{1} - 0.01826 g_{1}^{2} + 0.03644 g_{1}^{3} .

α = \frac{g_{2} (1 + g_{2})}{(g_{1} + g_{2}) (1 + g_{2} - g_{1})} .

Choosing g₁ = 0.9809 gives B_p = 0.0183 for the Haltrin phase function. Further details on these phase functions can be found in Mobley et al. [14].

3.3.5 Measured phase function by Sokolov et al

This phase function was derived from the measurements of VSF profiles in the black sea costal area during the summer periods of 2002, 2003 and 2004. More than 1000 VSF profiles were obtained using a VSF meter in four spectral channels 443 nm, 490 nm, 555 nm and 620 nm with bandwidth of 10 nm. These measurements were performed in between 0.6° to 178°. The mean VSF values calculated from these data at four wavelengths, with extrapolated values of VSF for angles less than 0.6ᵒ and more than 178ᵒ, are provided in Table 1 of Sokolov et al. [21]. The corresponding mean values of scattering and backscattering coefficients are also reported in their study.

4. Results and discussion

4.1 Influence of the optical properties on backscattering ratio

The backscattering ratio is estimated using Eq. (15) with varying parameters. Results are shown in Fig. 3, wherein the first subset shows the variation of backscattering ratio (in percentage) with maximum size of particles in the distribution (D_max) for fixed values of refractive index, D_min, PSD slope and wavelength (λ = 530) (Fig. 3(a)). It is found that particles greater than 100 μm contribute ~1% to the backscattering ratio. Figure 3(b) shows the effect of small particles on backscattering ratio for fixed values of other parameters. This shows that the B_p value increases significantly with decrease in D_min, particularly when D_min < 1 μm, which clearly illustrates the significant role of submicron particles in scattering. Further, it is found that the backscattering ratio is independent of D_min below 0.01μm. From these results, the minimum particle size (D_min) was chosen as 0.01 μm and maximum particle size (D_max) was chosen as 200 μm. Figure 3(c) shows that for a poly-dispersed distribution obeying power law (Junge type) the backscattering ratio is independent of wavelength over the visible range. This spectral invariance of particulate backscattering ratio was verified experimentally by Whitmire et al. [23]. Therefore, the wavelength of simulation was chosen to be 530 nm which is of primary interest for computations involving RTE in underwater. Figure 3(d) shows the variation of B_p as a function of the particle size distribution slope 'ξ ' for different refractive index values. The range of values of 'ξ ' was chosen based on the previous studies [12, 13]. From this plot it is clear that B_p is highly dependent on 'ξ ' and this dependency increases with increase in the refractive index 'n '. Waters with higher values of 'ξ ' (a more negative slope in a log-log plot) will have a higher backscattering ratio. This strong dependency of 'B_p' on 'ξ ' suggests that the underwater optical properties depend on the shape of the particle size distribution.

Fig. 3 Variation of the particle backscattering coefficient with (a) maximum particle diameter 'D_max' with fixed D_min = 0.01 μm, λ = 530 nm, m = 1.05-0.0001i and PSD slope (ξ) = 3, 3.5, 4, 5, (b) minimum particle diameter ' D_min' with fixed D_max = 200 μm, λ = 530 nm, m = 1.05-0.0001i and PSD slope (ξ) = 4, (c) wavelength ' λ ' with fixed D_min = 0.01 μm, D_max = 200 μm, and PSD slope (ξ) = 4, and (d) PSD slope with fixed D_min = 0.01 μm, D_max = 200 μm, λ = 530 nm and real refractive indices values n = 1.04- 0.02- 1.2.

Download Full Size | PDF

4.2 Evaluation of the inversion model

The performance of the model is evaluated based on the in situ data of B_P and c_p (λ) that were measured from turbid coastal waters off Point Calimere (on the southeast part of India) and relatively clear waters off Chennai during the periods 13 - 21 May 2012, 15 - 18 August 2012 and 24 August - 1 September 2013. The in situ backscattering ratio values (B_P) were derived from two independent measurements carried out in the field. For example, the particulate backscattering coefficients (b_bp(λ)) at nine discrete wavelengths (412, 440, 488, 510, 532, 595, 650, 676, 715 nm) were measured with a WET Labs BB9 sensor. The particulate scattering coefficients (b_p(λ)) throughout the visible and near-infrared spectrum (400-750 nm; in the interval of approximately 3-4 nm) were obtained from the measurements of beam attenuation (c_p(λ)) and absorption (a_p(λ)) coefficients made with a WET Labs ac-s meter (i.e., b_p(λ) = c_p(λ) - a_p(λ)). Prior to this analysis, data from the ac-s instrument were corrected for temperature, salinity and scattering effects (more details in Nasiha et al. [24]). The value of 'ξ ' was estimated from the beam attenuation spectrum c_p(λ), according to Diehl and Haardt [25], which follows the exact relation γ = ξ - 3, with 'γ' as the slope of the beam attenuation spectrum. According to Twardowski et al. [26], this relation is valid only for ξ ≥ 3.5 and for ξ values less than 3.5, nonlinearity is observed. Therefore, for harbor waters where the attenuation slope is less than 3.5, Eq. (29) of Nasiha et al. [24] is used to calculate the PSD slope. The estimated values of refractive index were examined for different water types, and are given in Table 4 along with the corresponding values of ' γ ', 'ξ ' and 'B_p'. For turbid coastal waters off Point Calimere (May 2012), the values of bulk refractive index range from 1.11 to 1.16. The higher values were observed at morning times and as the day passed the refractive index values decreased. This trend is because of a decrease in sediment concentration caused by diurnal effects of tides and currents in this region. A similar trend was observed for the time series data from the Jan. 2014 cruise. For harbor and lagoon waters, the bulk refractive index values ranged from 1.09 to 1.11 because of the significant contribution of detrital particles and phytoplankton relative to suspended sediments. For clear waters off Chennai, the refractive index values lie in the range 1.067-1.075, which was expected as the clear waters are dominated mainly by phytoplankton and its co-varying constituents. These ranges of bulk refractive index values for different water types are consistent with those reported by Twardowski et al. [26].

Table 4. Particulate optical properties obtained for different types of waters. Results also show the refractive index values derived for these waters.

View Table | View all tables in this article

4.3 Validation of the phase function model

From Eqs. (22) to (27) it is clear that the key inputs to the phase function model are 'n ' and 'ξ '. For validation of the new phase function, the Petzold average particle phase function is chosen as done in the previous studies [14, 21]. Since this phase function is based on experimental measurements and more frequently used by researchers, it is taken as the bench mark and is used for validating the new phase function model. The comparison results are shown in Fig. 4, wherein the Petzold average particle phase function is represented by the black line with a backscattering ratio value of B_p = 0.0183.

Fig. 4 Comparison of the phase functions obtained through different models for the scattering angles (a) from 0.1° to 90° and (b) 0.1° to 10° (extended 10° instead of 5° for better clarity) to emphasize the large variation and deviation of different phase functions at small angles.

Download Full Size | PDF

Figure 5 shows the different contours of B_p calculated from the backscattering model. These contour lines do not show the asymptotic behavior as the Fournier-Forand phase function does. This is due to the reason that in the analytical derivation of FF phase function, anomalous diffraction approximation of Mie scattering was taken and the infinite series of backscattering efficiency was approximated by a finite series. Further, this finite series was approximated, for simplicity, by an expression having asymptotic behavior at both small and large values of particle sizes [9]. On the other hand, in case of Mie theory, the infinite series is terminated according to Eq. (8). The same kind of non-asymptotic behavior can be found for the inversion model of Twardowski et al. [25]. From Fig. 5 it is clear that there is no unique set of values of 'n' and 'ξ ' for a particular B_p value. Several sets of values of 'n' and 'ξ ' are obtained for the backscattering ratio of 0.0183 using the backscattering model. However, it is found that for a particular set of 'n' and 'ξ ' (n = 1.16 and ξ = 3.4586), the new model provides a phase function closely matching with the Petzold average particle phase function, as shown in Fig. 4. It can be seen that the OTHG phase function levels out of the Petzold average particle phase function for small angles less then few degrees. Though the TTHG model gives improved values of $\tilde{β}$ at near forward angles, in comparison to the OTHG function model, it shows a considerable amount of deviation from the Petzold average particle phase function in the far forward angles. Looking at the FF phase function calculated with the same values of 'n' and 'ξ', it gives a good fit to the Petzold phase function in almost all angles while overestimating the values at small angles.

Fig. 5 Contours of the backscattering ratio B_p from Eq. (21).

Download Full Size | PDF

The quantitative measure of the differences of all these phase functions from a standard phase function (i.e., Petzold phase function) is obtained based on some standard statistical parameters: for example, Root Mean Square Error (RMSE), Mean Relative Error (MRE), slope, intercept and correlation coefficient (R²) which are commonly used for evaluating the model results [27]. The present model yielded MRE = −0.13, RMSE = 0.11, bias = −0.017, slope = 1, intercept = −0.042, and R² = 0.99, which are significantly better than those of the OTHG and TTHG models. The above statistical values of our model are nearly identical to those of the FF model. The lower values of error and higher values of slope and R² associated with the new model indicates its potential in predicting the phase function of arbitrary water type with bulk refractive index and particle size distribution slope as the inputs.

The best fit value n = 1.16 is consistent with the idea that the San-Diego Harbor water may be dominated by suspended sediments with an average refractive index of around 1.16 [28]. Moreover, according to Satyendranath et al. [19], oceanic particles having B_p values more than a few percent are mineral particles with indices of refraction in the order of about 1.15.

The sensitivity of our model with the bulk refractive index value of n = 1.16, which was determined based on the measurement data from San Diego harbor waters giving the best fit to the Petzold average phase function, is further verified by another independent measured data reported by Sokolov et al. [21]. As the present simulation was carried out for a wavelength of 530 nm, this independent comparison was performed for a wavelength of 555 nm. The corresponding value of backscattering is given as B_p = 0.0188. Several sets of values of n and ξ were obtained for this particular backscattering value using the backscattering model and the corresponding phase function values were calculated. It is found that for a particular set of values of 'n' and 'ξ ' (i.e., n = 1.16; ξ = 3.5188), the present model gives a phase function closely matching with the data reported by Sokolov et al. [21] (Fig. 6). A change in slope at 5° can be observed for both the phase functions. This best fit value of n = 1.16 is consistent with the value of average bulk refractive index reported by Chami et al. [29] for turbid coastal waters of the Black Sea. The reported bulk refractive index values varied from 1.12 to 1.2, with an average of 1.16 which is consistent with the mineral composition of particles commonly found in the coastal ocean (such as quartz with n = 1.15 or mica with n = 1.2).

Fig. 6 Comparison of the present study with the data presented by Sokolov et al. [21]

Download Full Size | PDF

5. Conclusion

In the present work, an extensive numerical computation for various inherent optical properties of aquatic particles was performed. The variation of 'B_p' with different optical properties of particles was investigated. From this study, it is concluded that particles with size less than 0.01 μm and greater than 200 μm have very little effects on light scattering in the visible range. However, the particle size distribution plays a significant role in deciding the intensity distribution of scattered light. In particular, particles with submicron size play a major role on the light scattering properties. The model of backscattering ratio (B_p) was inverted to form a relationship which correlates the refractive index to some physically meaningful and measurable optical parameters. This model was further employed to calculate the refractive indices for some selected ocean and coastal waters and expected values were obtained. In the second part, a parameterized relation for the scattering phase function of marine particles (for forward angles) was established purely from the numerically computed data. The two dependent variables of this phase function, namely the refractive index and PSD slope, are easily obtainable parameters from in situ measurements using the commercially available instrumentation. The present phase function model when compared with the Petzold average particle phase function and results from the existing models yielded encouraging results. Further, the model results were verified independently with another recently measured phase function data from Sokolov et al. [21]. It was found that the value of refractive index predicted by the present model is in agreement with that reported for turbid coastal waters. Some discrepancies were evident between the measured and modeled values of phase function due to approximations in Mie computation as discussed in section 2. However, results from this study clearly indicate that the present model predicts the phase function with the desired accuracy and hence will have important implications in marine optics and underwater optical studies, particularly in underwater imaging and channel modeling for underwater optical wireless communication systems.

Acknowledgments

This research work was supported by the Indian Institute of Technology (IIT) Madras, Chennai – 600 036. This research was supported in part by the Indian National Centre for Ocean Information Services (INCOIS) through the SATCORE program (OEC1314117INCOPSHA). We sincerely thank Dr. Gerard Wysocki, Associate Editor, and the two anonymous reviewers for their constructive comments and suggestions, which helped to improve the quality of this manuscript.

References and links

1. N. G. Jerlov and E. S. Nielsen, Optical Aspects of Oceanography (Academic Press, 1974), Chap. 1.

2. C. D. Mobley, Light and Water: Radiative Transfer in Natural Waters (Academic Press, 1994).

3. T. J. Petzold, “Volume Scattering Functions for selected ocean waters,” SIO Ref. 72–78, (Scripps Institute of Oceanography, (1972).

4. L. C. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]

5. L. O. Reynolds and N. J. McCormik, “Approximate two-parameter phase function for light scattering,” J.O.S.A. 70(10), 1206–1212 (1980). [CrossRef]

6. V. I. Haltrin, “One-parameter two-term Henyey-Greenstein phase function for light scattering in seawater,” Appl. Opt. 41(6), 1022–1028 (2002). [CrossRef] [PubMed]

7. O. B. Brown and H. R. Gordon, “Two component mie scattering models of sargasso sea particles,” Appl. Opt. 12(10), 2461–2465 (1973). [CrossRef] [PubMed]

8. M. Jonasz and H. Prandke, “Comparison of measured and computed light scattering in the Baltic,” Tellus 38B(2), 144–157 (1986). [CrossRef]

9. G. Fournier and J. L. Forand, “Analytic phase function for ocean water,” in Ocean Optics XII, J.S Jaffe, eds., Proc. SPIE 2258, 194–201 (1994).

10. H. C. Van De Hulst, Light Scattering by Small Particles (Dover Publications Inc., 1981).

11. A. C. Holland and G. Gagne, “The scattering of polarized light by polydisperse systems of irregular particles,” Appl. Opt. 9(5), 1113–1121 (1970). [CrossRef] [PubMed]

12. H. Bader, “The hyperbolic distribution of particle sizes,” J. Geophys. Res. 75(15), 2822–2830 (1970). [CrossRef]

13. J. C. Kitchen, “Particle size distributions and the vertical distribution of suspended matter in the upwelling region off Oregon,” 77–10 (Oregon State Univ., Corvallis, 1977), pp. 118.

14. C. D. Mobley, L. K. Sundman, and E. Boss, “Phase function effects on oceanic light fields,” Appl. Opt. 41(6), 1035–1050 (2002). [CrossRef] [PubMed]

15. M. Kerker, The Scattering of Light and Other Electromagnetic Radiation (Academic Press, 1969).

16. C. Mätzler, “MATLAB functions for Mie scattering and absorption,” 2002–08 (University of Bern, 2002), pp. 18.

17. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley, 1983).

18. B. Aernouts, R. Watté, J. Lammertyn, and W. Saeys, “A flexible tool for simulating the bulk optical properties of polydisperse suspensions of spherical particles in an absorbing host medium,” in Optical Modeling and Design II, Frank Wyrowski, John T. Sheridan, Jani Tervo, Youri Meuret, eds., Proc. of SPIE 8429,1–11 (2012).

19. O. Ulloa, S. Sathyendranath, and T. Platt, “Effect of the particle-size distribution on the backscattering ratio in seawater,” Appl. Opt. 33(30), 7070–7077 (1994). [CrossRef] [PubMed]

20. Y. C. Agrawal and O. A. Mikkelsen, “Empirical forward scattering phase functions from 0.08 to 16 deg. for randomly shaped terrigenous 1-21 µm sediment grains,” Opt. Express 17(11), 8805–8814 (2009). [CrossRef] [PubMed]

21. A. Sokolov, M. Chami, E. Dmitriev, and G. Khomenko, “Parameterization of volume scattering function of coastal waters based on the statistical approach,” Opt. Express 18(5), 4615–4636 (2010). [CrossRef] [PubMed]

22. C. D. Mobley, B. Gentili, H. R. Gordon, Z. Jin, G. W. Kattawar, A. Morel, P. Reinersman, K. Stamnes, and R. H. Stavn, “Comparison of numerical models for computing underwater light fields,” Appl. Opt. 32(36), 7484–7504 (1993). [CrossRef] [PubMed]

23. A. L. Whitmire, E. Boss, T. J. Cowles, and W. S. Pegau, “Spectral variability of the particulate backscattering ratio,” Opt. Express 15(11), 7019–7031 (2007). [CrossRef] [PubMed]

24. H. J. Nasiha, P. Shanmugam, and V. G. Hariharasudhan, “A new inversion model to estimate bulk refractive index of particles in costal oceanic waters: Implications for remote sensing,” IEEE J. Sel. Topics in App, Earth. Obs. Remote Sensing 7, 3069–3083 (2014).

25. P. Diehl and H. Haardt, “Measurement of the spectral attenuation to support biological research in a ‘plankton tube’ experiment,” Oceanol. Acta 3, 89–96 (1980).

26. M. S. Twardowski, E. Boss, J. B. Macdonald, W. S. Pegau, A. H. Barnard, and J. R. V. Zaneveld, “A model for estimating bulk refractive index from the optical backscattering ratio and the implications for understanding particle composition in case I and case II waters,” J. Geophys. Res. 106(C7), 14,129–14,142 (2001). [CrossRef]

27. Z. P. Lee, “Remote Sensing of Inherent Optical Properties: Fundamentals, Tests of Algorithms, and Applications,” 5 (IOCCG, Dartmouth, Canada, 2005), pp. 126.

28. D. Stramski, A. Bricaud, and A. Morel, “Modeling the inherent optical properties of the ocean based on the detailed composition of the planktonic community,” Appl. Opt. 40(18), 2929–2945 (2001). [CrossRef] [PubMed]

29. M. Chami, E. B. Shybanov, T. Y. Churilova, G. A. Khomenko, M. E. G. Lee, O. V. Martynov, G. A. Berseneva, and G. K. Korotaev, “Optical properties of the particles in the Crimea coastal waters (Black Sea),” J. Geophys. Res. 110(C11), C11020 (2005). [CrossRef]

m	*a_m*	*b_m*	*c_m*
1	0.03182	0.00416	0.1514
2	0.101	0.00372	−0.6116

m	*d_m*	*e_m*	*f_m*	*h_m*	*i_m*	*j_m*	*k_m*	*l_m*	*o_m*
1	0.000124	-0.1437	0.01982	0.000035	-0.03221	-0.06899	-0.00015	0.05137	0.06542
2	0.000218	-0.1622	-0.1345	0.000017	0.06596	0.2118	-0.00026	-0.04575	-0.58010
3	-0.000054	0.8617	-0.7334	0.000047	-0.08305	-0.64645	0.00032	-1.00100	2.79600

m	*d_m*	*e_m*	*f_m*	*h_m*	*i_m*	*j_m*	*k_m*	*l_m*	*o_m*
1	0.0012	0.4594	−0.5601	0.00029	0.214	−0.2145	−0.00088	−0.6234	0.6062
2	−0.0098	−3.1480	4.1990	−0.00250	−1.754	1.7921	0.00768	4.9780	−5.1440
3	0.0246	6.5966	−9.6520	0.00723	4.430	−4.3770	−0.02148	−11.640	11.930
4	−0.0198	−4.3123	6.2890	−0.00605	−3.486	2.9660	0.01853	7.7180	−6.6490

Measurement period / Water type	Station / No. of samples	Attenuation slope (γ)	PSD Slope	*Back Scattering ratio (B_p)*	Bulk Refractive Index(n)
May-2012 (Turbid water)	9 am	0.733	3.733	0.021	1.153
	10:30 am	0.816	3.816	0.017	1.143
	12:30 pm	0.928	3.928	0.018	1.128
	2 pm	1	3.820	0.018	1.118
	3:30 pm	1.01	3.911	0.015	1.116
Jan-2014(Turbid water)	8:30 am	0.671	3.671	0.025	1.190
	10 am	0.661	3.661	0.024	1.182
	11:30 am	0.744	3.744	0.023	1.165
	1 pm	0.734	3.734	0.021	1.152
	2:30 pm	0.785	3.785	0.021	1.145
	4 pm	0.65	3.650	0.019	1.148
2013 (Clear water)	8 am	0.836	3.836	0.011	1.077
	10 am	0.836	3.836	0.011	1.077
	12	0.933	3.933	0.011	1.068
	1:30 pm	0.86	3.860	0.011	1.075
	3 pm	0.856	3.856	0.011	1.075
2013 (Harbour water)	1	0.321	3.217	0.009	1.100
	2	0.19	3.190	0.011	1.103
	3	0.258	3.258	0.014	1.102
	4	0.391	3.391	0.015	1.095
2013 (Productive lagoon water)	1	0.924	3.924	0.017	1.105
	2	0.912	3.912	0.017	1.109
	3	0.793	3.793	0.017	1.123
	4	0.756	3.756	0.017	1.127
2014 (Productive lagoon water)	1	0.887	3.887	0.017	1.111
	2	0.703	3.703	0.017	1.132
	3	0.597	3.597	0.017	1.142
	4	0.618	3.618	0.017	1.141
	5	0.597	3.597	0.017	1.143

m	*a_m*	*b_m*	*c_m*
1	0.03182	0.00416	0.1514
2	0.101	0.00372	−0.6116

Semi-analytical modeling and parameterization of particulates-in-water phase function for forward angles

Abstract

1. Introduction

2. Basic Mie theory and simulation

3. Model description

3.1 Modeling of the backscattering ratio

3.2 Modeling of the phase function

3.3 Model comparison

3.3.1 Petzold Average Particle Phase Function

3.3.2 Fournier-Forand Phase Function

3.3.3 One -Term Henyey-Greenstein phase function (OTHG)

3.3.4 Haltrin's Two-Term Henyey-Greenstein Phase Function (TTHG)

3.3.5 Measured phase function by Sokolov et al

4. Results and discussion

4.1 Influence of the optical properties on backscattering ratio

4.2 Evaluation of the inversion model

4.3 Validation of the phase function model

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (6)

Tables (4)

Equations (36)

Optics Express