Fast and accurate model of underwater scalar irradiance for stratified Case 2 waters

doi:10.1364/OE.14.001703

« Show journal navigation

Fast and accurate model of underwater scalar irradiance for stratified Case 2 waters

Cheng-Chien Liu »View Author Affiliations

Optics Express, Vol. 14, Issue 5, pp. 1703-1719 (2006)
http://dx.doi.org/10.1364/OE.14.001703

View Full Text Article

Acrobat PDF (437 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Radiative transfer in aquatic environments
3. Theory of approximation
4. Results
5. Summary
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (32)
Cited By
Figures (5)
Metrics

Abstract

This paper is devoted to the derivation of a fast and accurate model of scalar irradiance for stratified Case 2 waters. Five strategies are formulated and employed in the new model, including (1) reallocating the sky radiance, (2) approximating the influence of the air-water interface, (3) constructing a look-up table of average cosine based on the single-scattering albedo and the backscatter fraction, (4) calculating the phase function of surrogate particles in Case 2 waters, and (5) using the average cosine as an index to cope with stratified waters. A comprehensive model-to-model comparison shows that the new model runs more than 1,400 times faster than the commercially-available Hydrolight model, while it limits the percentage error to 2.03% and the maximum error to less than 6.81%. This new model can be used interactively in models of the oceanic system, such as biogeochemical models or the heat budget part of global circulation models.

1. Introduction

The quantitative study of the interactions of light with the aquatic environment on earth, such as oceans, estuaries, lakes, rivers and other water bodies, plays an essential role in hydrologic optics [1

C D. Mobley, Light and water: radiative transfer in natural waters (Academic Press, San Diego, CA, 1994), p. 592.

]. In particular, an accurate simulation of the vertical profile of underwater scalar irradiance E ₀(z) is a prerequisite for modeling hydrologic systems, such as the photosynthesis process in biogeochemical models or the heat budget part of ocean circulation models [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

]. After being scattered and absorbed by the atmosphere, the photons penetrate through the wavy boundary of air-water interface from all possible directions and interact with water bodies. A bipartite classification is generally used to categorize the water bodies as Case 1 and Case 2 waters [3

A. Morel and L. Prieur, “Analysis of variations in ocean color,” Limnol. Oceanogr. 22, 709–722 (1977). [CrossRef]

, 4

H. R. Gordon and A. Y. Morel, Remote assessment of ocean color for interpretation of satellite visible imagery: a review , Lecture Notes on Coastal and Estuarine Studies (Springer-Verlag, New York, 1983), Vol. 4, p. 114.

]. By definition, phytoplankton with their accompanying and covarying retinue of material of biological origin is the principal component responsible for variations in optical properties of Case 1 waters, while inorganic particles in suspension and yellow substances all influence the optical properties of Case 2 waters. There has been a lot of research interest in investigating Case 2 waters since over 60% of the human population, and 90% of the world’s fish catch are directly linked to the coastal zones of Case 2 waters [5

J. C. Pernetta and J. D. Milliman, “Land-ocean interactions in the coastal zone implementation plan,” IGBP Report 33 (Stockholm, 1995).

]. However, the study of radiative transfer in Case 2 waters is quite challenging because of the complex constitution of their optical components [6

S. Sathyendranath, ed., Remote sensing of ocean colour in coastal, and other optically-complex, waters , IOCCG Report (MacNab Print, Dartmouth, Canada, 2000), Vol. 3, p. 140.

, 7

A. Albert and C. D. Mobley, “An analytical model for subsurface irradiance and remote sensing reflectance in deep and shallow case-2 waters,” Opt. Express 11, 2873–2890 (2003). [CrossRef] [PubMed]

This paper aims at deriving a fast yet accurate approximation for calculating the underwater scalar irradiance for stratified Case 2 waters. In essence this is the extension of earlier work on the case of homogeneous Case 1 waters [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

]. A comprehensive model-to-model comparison shows that the new model runs more than 1,400 times faster than the full Hydrolight code, while it limits the percentage error to 2.03% and the maximum error to less than 6.81%. This new model can be used interactively in models of the oceanic system, such as biogeochemical models or the heat budget part of global circulation models.

2. Radiative transfer in aquatic environments

For a body of water with a constant index of refraction, the radiative transfer equation (RTE) can be derived from more fundamental equations [1

C D. Mobley, Light and water: radiative transfer in natural waters (Academic Press, San Diego, CA, 1994), p. 592.

], such as Maxwell’s electromagnetic equation [8

R. W. Preisendorfer, Hydrologic optics (U. S. Department of Commerce, Seattle, WA, 1976).

, 9

R. L Fante, “Relationship between radiative-transport theory and Maxwell’s equations in dielectric media,” J. Opt. Soc. Am. 71, 460–468 (1981).

], quantum mechanics [10

R. M. Measures, Laser remote sensing: fundamentals and applications (Kreiger, Malabar, Florida, 1992), p. 510.

] and the Boltzmann equation [11

K. Stamnes, “The theory of multiple scattering of radiation in plane parallel atmospheres,” Rev. Geophys. 24, 299–310 (1986). [CrossRef]

]. RTE states that the rate of change in radiance L through a distance r along the direction

\tilde{ξ}

at a spatial position

\vec{x}

is the linear combination of loss due to attenuation -cL and gain due to scattering L _* and internal emission

L_{*}^{S}

\frac{DL (\vec{x}; \hat{ξ}; λ)}{Dr} = - c (\vec{x}; \hat{ξ}; λ) L (\vec{x}; \hat{ξ}; λ) + L_{*} + L_{*}^{S},

(1)

where c is the beam attenuation coefficient. The radiative energy transfers within the water body at a constant speed v that equals the speed of light divided by the index of refraction. Therefore, the total rate of change along r is

\frac{D}{Dr} = \frac{1}{v} \frac{D}{Dt} = \frac{1}{v} (\frac{\partial}{\partial t} + \vec{v} ∙ \nabla) = \frac{1}{v} \frac{\partial}{\partial t} + \hat{ξ} ∙ \nabla = \hat{ξ} ∙ \nabla,

(2)

where

\hat{ξ}

≡

\vec{v}

/v is the unit vector along the path r. Note that the time scale for reaching the steady state in studies of hydrologic optics is generally much shorter than the sampling time, therefore, the time-independent assumption is valid and the term of partial derivative with t is dropped. L _* can be further divided into inelastic scattering

L_{*}^{I}

and elastic scattering

L_{*}^{E}

. For the elastic scattering, the volume scattering phase function β describes the probability of a photon scattered to the direction

\hat{ξ}

from all directions

\hat{ξ}

′ in a unit sphere Ξ. For the inelastic scattering,

L_{*}^{I}

is usually considered as an extra source that elevates the energy level at this wavelength. The common practice is to combine

L_{*}^{I}

and the emission term

L_{*}^{S}

into an effective source function S. Liu et al. [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

] tested the sensitivity of various factors in simulating the vertical profile of underwater scalar irradiance. They concluded that S plays a negligible role. Therefore, the general form of the monochromatic RTE is

\hat{ξ} ∙ \nabla L (\vec{x}; \hat{ξ}; λ) = - c (\vec{x}; λ) L (\vec{x}; \hat{ξ}; λ) + \int_{\hat{ξ′} \inΞ} L (\vec{x}; \hat{ξ′}; λ) β (\vec{x}; \hat{ξ'} \to \hat{ξ}; λ) d Ω (\hat{ξ'})

(3)

Plass and Kattawar [12

G. N. Plass and G. W. Kattawar, “Monte Carlo calculations of radiative transfer in the Earth’s atmosphere-ocean system. I. Flux in the atmosphere and ocean,” J. Phys. Oceanogr. 2, 139–145 (1972). [CrossRef]

] first applied the Monte-Carlo (MC) method to solve the general form of RTE in hydrologic optics. The basic principle of the MC method is to simulate a beam of light by using a very large number of photons. Following the path of each photon, a series of random numbers can be used to determine the photon’s life history according to different probabilities for different phenomena. The final light field is the cumulative contribution of all the photons. Based on the recent work of Liu and Woods [13

C.-C. Liu and J. Woods, “Prediction of ocean colour: Monte-Carlo simulation applied to a virtual ecosystem based on the Lagrangian Ensemble method,” Int. J. Remote Sens. 25, 921–936 (2004). [CrossRef]

], a 3D backward MC model was developed to simulate the light field in coastal waters with uneven bottom geometry [14

K. L. Carder, C.-C. Liu, Z. P. Lee, D. C. English, J. Patten, R. F. Chen, J. E. Ivey, and C. O. Davis, “Illumination and turbidity effects on observing faceted bottom elements with uniform Lambertian albedos,” Limnol. Oceanogr. 48, 355–363 (2003). [CrossRef]

]. Their model is able to provide solutions to the most general form of RTE. However, if the information required is not merely the surface reflectance but the detailed light field under the sea surface, the MC approach would require quite considerable computer resources that would therefore limit its practical applications [1

C D. Mobley, Light and water: radiative transfer in natural waters (Academic Press, San Diego, CA, 1994), p. 592.

Many efforts have been made to find a fast and accurate solution of the RTE without running the full MC simulation, and one way it can be done is by further simplifying the RTE. For the body of water with a flat bottom and where the horizontal variance of inherent optical properties (IOPs) is much less than their vertical variance, the assumption of plane-parallel waters can be made to reduce the spatial variable

\vec{x}

(vector) to one scalar variable, depth z. The general form of RTE, therefore, can be simplified as an integro-differential equation

\cos θ \frac{dL (z; θ; ϕ; λ)}{dz} = - c (z; λ) L (z; θ; ϕ; λ) + \int_{ϕ′ = 0}^{2 π} \int_{θ′ = 0}^{π} L (z; θ'; ϕ'; λ) β (z; θ' \to θ; ϕ' \to ϕ; λ) \sin θ' dθ' dϕ' .

(4)

By introducing the assumption of quad-averaging, the integro-differential equation can be further simplified to a set (typically hundreds) of ordinary differential equations. The invariant imbedding technique can be employed to imbed the air-sea and bottom boundary conditions to the set of ODEs, and the detailed distribution of L(z;θ;ϕ;λ) can then be solved simultaneously. The associated computer code was initially developed in 1979 and has since gone through a series of modifications [15–17

R. W. Preisendorfer and C. D. Mobley, “Unpolarized irradiance reflectances and glitter patterns of random capillary waves on lakes and seas, by Monte Carlo simulation,” NOAA Technical Memorandum, ERL PMEL-63 (NOAA Pacific Marine Environmental Laboratory, Seattle, WA, 1985).

]. Version 3.0 was funded by the US Office of Naval Research and renamed as Hydrolight in year 1995 [18

C. D. Mobley, “Hydrolight 3.0 Users’ Guide,” (SRI International, Menlo Park, CA, 1995).

]. The latest version of Hydrolight, 4.3, was released in March 2003 as commercial software. However, it is not always practical to use Hydrolight to calculate global underwater light fields, as the computational time required for running a full Hydrolight simulation interactively in biogeochemical models or global circulation models is often too long [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

]. Because of this, the general approach still relies on the empirical algorithm of K to calculate the underwater light field, resulting in considerable errors [19

C.-C. Liu, J. D. Woods, and C. D. Mobley, ”Optical model for use in oceanic ecosystem models,“ Appl. Opt. 38, 4475–4485 (1999). [CrossRef]

]. The progress of modeling the oceanic system is impeded somewhat in that while we know that error exists, and have the tool to correct the error, we simply cannot afford to use that tool. As a result, we all tolerate and live with error in simpler but faster models, despite that Hydrolight has been available for several years.

3. Theory of approximation

Liu et al. [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

] employed four strategies to accelerate the simulation of the scalar irradiance E ₀(z) without losing accuracy compared with Hydrolight. Their model runs more than fourteen thousand times faster than the full Hydrolight code, while limiting the percentage error to 2.20% and the maximum error to less than 4.78%. However, their model was limited to (1) vertically homogeneous water columns, (2) Case 1 waters or Case 2 waters that happen to be gelbstoff rich, (3) spectral ranges from 400 to 700 nm, and (4) chlorophyll concentration (Chl) from 0 to 10.0 mg∙m^-3. This paper modifies one strategy described in Liu et al. [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

] and proposes three new strategies to remove those limitations, resulting in a new model of E ₀(z) that can be applied to the stratified Case 2 waters. Details of these strategies are given as follows.

3.1 Reallocation of sky radiance

In the realm of hydrologic optics, the incident radiative energy originates from the Sun. After being scattered and absorbed by the atmosphere, the photons penetrate through the wavy boundary of the air-water interface from all possible directions (θ;ϕ) and interact with the bodies of water. Because the quantity of concern, E ₀(z), is computed from integrals of the radiance over the entire sphere Ξ,

E_{0} (z) = \int_{Ξ} L (z; \hat{ξ}) d Ω = \int_{ϕ = 0}^{2 π} \int_{θ = 0}^{π} L (z; θ; ϕ) \sin θd θdθ .

(5)

The incident photons which arrive at the sea surface (z = 0⁺) from all possible direction (θ;ϕ) can be reallocated to the plane of the Sun (θ;ϕ ₀) and the same integrated result of E ₀(z) still can be obtained. In other words, the planar irradiance incident on the surface can be expressed as

E_{d} (0^{+}) = \int_{ϕ = 0}^{2 π} \int_{θ = 0}^{\frac{π}{2}} L (0^{+}; θ; ϕ) \cos θ \sin θd θdϕ = \int_{θ = 0}^{\frac{π}{2}} \bar{L} (0^{+}; θ; ϕ_{0}) \cos θ \sin θd θ,

(6)

where the cumulative radiance in the azimuthal direction is

\overset{̅}{L} (0^{+}; θ; ϕ_{0}) = \int_{ϕ = 0}^{2 π} L (0^{+}; θ; ϕ) dϕ .

(7)

Applying the superposition principle at each depth z, the quantity of concern E ₀(z) can be quickly calculated by summing the contribution from each light source in the plane of the Sun (θ;ϕ ₀) by

E_{0} (z) = \int_{θ = 0}^{\frac{π}{2}} {\overset{̅}{E}}_{0} (z; θ; ϕ_{0}) \bar{L} (0^{+}; θ; ϕ_{0}) \cos θ \sin θd θ .

(8)

where E̅₀(z;θ;ϕ ₀) is the profile of E ₀ obtained by placing in a black sky a unit radiance in (θ;ϕ ₀). As long as the distribution of incident radiance L(0⁺;θ;ϕ) is given and E ₀(z;θ;ϕ ₀) is calculated and stored in advance, solving E ₀(z) is simply a procedure of multiplication and summation that can be done with very high efficiency [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

]. Note that no approximation needs to be made at this stage.

3.2 Approximation of the influence of the air-water interface

The air-water interface is usually roughened by wind, and the influence of this on simulating E ₀(z) can be regarded as redistributing the radiance field after the beam of light penetrates the interface. To consider this effect of redistribution, the profile of Ē₀(z;θ;ϕ ₀) can be scaled by multiplying a corrective factor for wind speed CW that is defined as

CW (θ; ϕ_{0}) ∣_{V_{wind} \to {V'}_{wind}} \equiv \frac{{\bar{E}}_{0} (0^{-}; θ; ϕ_{0}) ∣_{{V'}_{wind}}}{{\bar{E}}_{0} (0^{-}; θ; ϕ_{0}) ∣_{V_{wind}}},

(9)

where Ē₀(0^-;θ;ϕ ₀)|_{V_wind} and E ₀(0^-;θ;ϕ ₀)|_{V'_wind} are the subsurface values of E ₀ obtained by placing in a black sky a unit light source in (θ;ϕ ₀) under the condition of surface wind speed V_wind and V́_wind , respectively. Given CW(θ;ϕ ₀)|_{V_wind→V'wind}, the vertical profile of E ₀(z) calculated at wind speed V_wind (Eq. 8) can be easily extended to E ₀(z) at wind speed V′_wind by

E_{0} (z) ∣_{{V'}_{wind}} = \int_{θ = 0}^{\frac{π}{2}} {\bar{E}}_{0} (z; θ; ϕ_{0}) \bar{L} (0^{+}; θ; ϕ_{0}) CW (θ; ϕ_{0}) ∣_{V_{wind} \to {V'}_{wind}} \cos θ \sin θd θ .

(10)

To be more precise, CW is influenced by water constituents as well. Liu et al. [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

] demonstrated that CW is only a weak function of water constituents. It is therefore feasible to construct a look-up table (LUT) based on the situation of clear water for quick reference to the CW. Liu et al. [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

] calculated CW at 31 fixed wavelength from 400 nm to 700 nm at steps of 10 nm. To have a flexibility on the wavelength, the main variables of the new LUT are selected as ω ₀(0.01 – 0.99), V _wind (0.0 – 15.0) and θ (0.0 – 90.0). A standard volume scattering phase function for water molecules is obtained from the analytic Fournier-Forand phase function [20

G. R. Fournier and J. L. Forand, ”Analytic phase function for ocean water,“ presented at the SPIE: Ocean Optics XII, 1994.

] by giving the backscatter fraction (BF) a value of 0.5. A section of the LUT is denoted as LUT_CW and given in Table 1. Note that the table is for pure water with an assumed absorption thus that ω ₀ can be varied from 0.01 to 0.99.

Table 1. A section of LUT_CW for quick reference to the corrective factor CW based on wind speed V_wind and albedo ω ₀.

Quad		15°~25°			35°~45°
ω ₀	V_wind	1.0	3.0	10.0	1.0	3.0	10.0
0.01		1.000477	1.000185	1.000058	1.005423	1.002095	0.998207
0.03		1.000433	1.000163	1.000038	1.005430	1.002088	0.998149
0.05		1.000399	1.000142	1.000038	1.005444	1.002089	0.998093
0.07		1.000366	1.000122	1.000028	1.005447	1.002090	0.998039
0.09		1.000324	1.000102	1.000009	1.005457	1.002090	0.997987
0.11		1.000284	1.000073	0.999991	1.005455	1.002080	0.997920
0.13		1.000253	1.000063	0.999982	1.005460	1.002087	0.997872
0.15		1.000205	1.000036	0.999955	1.005478	1.002084	0.997794
0.17		1.000158	1.000018	0.999947	1.005470	1.002080	0.997743
0.19		1.000113	0.999991	0.999922	1.005467	1.002076	0.997678
0.21		1.000068	0.999966	0.999897	1.005462	1.002071	0.997609
0.23		1.000025	0.999941	0.999874	1.005454	1.002064	0.997543
0.25		0.999959	0.999917	0.999851	1.005458	1.002065	0.997480
0.27		0.999918	0.999894	0.999820	1.005437	1.002064	0.997420
0.29		0.999871	0.999871	0.999791	1.005421	1.002055	0.997348
0.31		0.999818	0.999850	0.999762	1.005409	1.002045	0.997281
0.33		0.999758	0.999821	0.999727	1.005387	1.002034	0.997210
0.35		0.999701	0.999801	0.999701	1.005356	1.002022	0.997136
0.37		0.999639	0.999767	0.999669	1.005328	1.002015	0.997074
0.39		0.999571	0.999741	0.999630	1.005297	1.002001	0.997002
0.41		0.999499	0.999709	0.999600	1.005257	1.001986	0.996935
0.43		0.999423	0.999679	0.999565	1.005207	1.001969	0.996867
0.45		0.999343	0.999650	0.999538	1.005154	1.001951	0.996797
0.47		0.999253	0.999609	0.999500	1.005099	1.001933	0.996727
0.49		0.999167	0.999577	0.999463	1.005034	1.001906	0.996662

3.3 Look-up table of average cosine based on ω₀ and BF

The key to achieving a fast and accurate simulation of E ₀(z) is the existence of a proper mathematical model that is able to approximate the vertical profile of the average cosine μ̄(z) with a set of only five parameters (B ₀, B ₁, P, B ₂, Q)

\frac{1}{\bar{μ} (ζ)} = B_{0} + B_{1} \exp (- Pζ) + B_{2} \exp (- Qζ) .

(11)

yet which still retains a very high accuracy [21

N. J. McCormick, ”Mathematical models for the mean cosine of irradiance and the diffuse attenuation coefficient,“ Limnol. Oceanogr. 40, 1013–1018 (1995). [CrossRef]

], where the optical depth ζ = cz and

\bar{μ}

(z) is

\bar{μ} (z) \equiv \frac{E (z)}{E_{0} (z)} .

(12)

Liu et al. [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

] demonstrated that a LUT can then be constructed for quick reference to set of five parameters (B ₀, B ₁, P, B ₂, Q) from the given IOPs. With this LUT, the net diffuse attenuation coefficient K_NET (z) can be calculated from μ̄(z) by use of Gershun’s equation [22]

K_{NET} (z) = \frac{a (z)}{\bar{μ} (z)} .

(13)

The definition of K_NET (z) is

K_{NET} (z) \equiv - \frac{d In E (z)}{dz} = - \frac{d In [E_{d} (z) - E_{u} (z)]}{dz},

(14)

where E_d and E_u are the downwelling and upwelling planar irradiance, respectively. Therefore, K_NET (z) gives us the profile of E(z) that can be converted to E ₀(z) using the definition of μ̄(z) (Eq. 12).

The major limitation of the LUT approach is the parameterization of IOPs. The earlier work [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

] was based on the biooptical model of Case 1 waters that uses the Chl, λ and the backscatter fraction (BF) as the main variables. As a result, the model of Liu et al. [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

] cannot be applied to Case 2 waters. The main variable Chl is limited to the range from 0 to 10.0 mg∙m^-3, and λ is limited to the range from 400 to 700 nm as well. Considering the fact that all IOPs can be derived from two principal IOPs: the absorption coefficient a and the phase scattering function β(ψ) [23

V. I. Haltrin, ”Chlorophyll-Based Model of Seawater Optical Properties,“ Appl. Opt. 38, 6826–6832 (1999). [CrossRef]

], where ψ is the scattering angle, the two principal IOPs should be selected as the main variables in the LUT. However, β(ψ) describes the probability of scattering in the ψ direction, which varies from a symmetric distribution (molecular scale) to a heavily peak distribution in the forward direction (large particle). To consider a large range of the distribution of β (ψ), it is necessary to find an accurate and efficient way to parameterize β(ψ). Fournier and Forand (FF) [20

G. R. Fournier and J. L. Forand, ”Analytic phase function for ocean water,“ presented at the SPIE: Ocean Optics XII, 1994.

] proposed an analytical phase function for ocean water, which parameterized β(ψ) as a function of the backscatter fraction (BF), defined as the backscattering coefficient b_b divided by the scattering coefficient b. A large range in β(ψ) can be obtained by varying the value of BF. For example, BF=0.0183 provides a very good fit to the Petzold’s average particle phase scattering function, while BF=0.5 yields the pure water phase function. Mobley el al. [24

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

] did a thorough sensitivity test of phase function effects on the oceanic light field. They suggested that the physically based FF phase function can be used to generate realistic phase functions having any desired BF. Therefore, BF is a good choice for parameterizing β(ψ).

Based on the work of Liu et al. [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

], we reconstructed a new LUT with two non-dimensional variables BF (0.0001 – 0.5) and ω ₀ (0.01 – 0.99) for quick reference to a set of parameters (B ₀, B ₁, P, B ₂, Q) used by the McCormick five-parameter model [21

N. J. McCormick, ”Mathematical models for the mean cosine of irradiance and the diffuse attenuation coefficient,“ Limnol. Oceanogr. 40, 1013–1018 (1995). [CrossRef]

] of μ̄(ζ). The single-scattering albedo ω ₀ is defined as b divided by c. A section of the LUT is denoted as LUT_μ and given in Table 2, where the Pearson correlation coefficient PCC shows how good the McCormick five-parameter model fits the vertical profiles of the average cosine. Throughout the entire LUT, the average value of PCC is 0.999577 and the lowest value is 0.980118.

Implementation of the new LUT is explained as follows. Given the IOPs of Case 2 waters, ω ₀ and BF_P can be calculated. The bulk property of BF is given by

BF (λ) \equiv \frac{b_{w} (λ)}{b (λ)} \times {BF}_{w} + \frac{b_{p} (λ)}{b (λ)} {BF}_{p},

(15)

where BF_w =0.5. Referencing to the new LUT, ω ₀ and BF_P give us a profile of μ̄_p(ζ), while ω ₀ and BF give us another profile of μ̄′(ζ) that can be regarded as the value for a pure medium with a backscatter fraction BF. However, neither μ̄_p(ζ) nor μ̄′(ζ) represents the actual profile of μ̄(ζ), because both particles and water molecules contribute to the effect of scattering. The earlier work [2

C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]

] suggested that a constant offset Δ exists between μ̄_p(ζ) and μ̄(ζ). Figure 1 illustrates that the assumption of a constant offset only held at depths where μ̄(ζ) has reached its asymptotic solution. To give a better approximation of μ̄(ζ) from μ̄_p(ζ) and μ̄′(ζ), this research proposes a new approach to describe the vertical profile of offset Δ(ζ). The offset near the surface is given by

Δ (0^{-}) = (1 - \frac{{BF}_{p}}{BF}) [{\bar{μ}}_{p} (0^{-}) - \bar{μ'} (0^{-})],

(16)

while the offset at infinite depth is

Δ (ζ_{\infty}) = \frac{1}{2} [{\bar{μ}}_{p} (ζ_{\infty}) - {\bar{μ'}}_{p} (ζ_{\infty})] .

(17)

μ̄(ζ) can be estimated from μ̄_p(ζ) by adding the offset Δ(ζ), i.e.,

\bar{μ} (ζ) = {\bar{μ}}_{p} (ζ) + Δ (ζ) .

(18)

Where the vertical profile of offset Δ(ζ) is given by

Δ (ζ) = [Δ (0^{-}) - Δ (ζ_{\infty})] [\frac{{\bar{μ}}_{p} (ζ) - {\bar{μ}}_{p} (ζ_{\infty})}{{\bar{μ}}_{p} (0^{-}) - {\bar{μ}}_{p} (ζ_{\infty})}] .

(19)

Table 2. A section of the LUT_μ for quick reference to a set of parameters (B ₀, B ₁, P, B ₂, Q) used by the McCormick five-parameter model [21

N. J. McCormick, ”Mathematical models for the mean cosine of irradiance and the diffuse attenuation coefficient,“ Limnol. Oceanogr. 40, 1013–1018 (1995). [CrossRef]

] of μ̄(ζ). This new LUT is based on two non-dimensional variables BF (0.0001 – 0.5) and ω ₀(0.01 – 0.99).

quad	BF	ω ₀	B ₀	B ₁	P	B ₂	Q	PCC
0°~5°	0.018	0.31	1.119133	0.086067	0.087585	0.017142	0.589276	0.999981
5°~15°	0.018	0.31	1.119114	0.073322	0.107362	0.016351	0.647022	0.999977
15°~25°	0.018	0.31	1.119486	0.062006	0.236405	-0.001879	0.236448	0.999166
25°~35°	0.018	0.31	1.114392	0.057128	0.324473	-0.040133	0.039546	0.998720
35°~45°	0.018	0.31	1.118637	-0.179313	0.089442	0.129340	0.269664	0.999598
45°~55°	0.018	0.31	1.119143	-0.261702	0.100093	0.128522	0.350085	0.999912
55°~65°	0.018	0.31	1.121492	-0.561937	0.134360	0.337626	0.281429	0.999969
65°~75°	0.018	0.31	1.121764	-0.582089	0.138694	0.221943	0.434598	0.999989
75°~85°	0.018	0.31	1.124118	-0.885466	0.170924	0.363883	0.454816	0.999959
85°~90°	0.018	0.31	1.125022	-1.013489	0.187065	0.304746	0.660470	0.999960

To illustrate the validity of Eq. (18), a total of six cases with various particle sizes and concentrations were selected by setting Chl (1.0, and 10.0 mg∙m^-3) and BF_P (0.00915, 0.0183 and 0.0366), respectively. Other conditions of computation were set as typical Case 1 waters and are listed in the caption of Fig. 1. Both Hydrolight and our approach are employed to simulate the profiles of μ̄(ζ). The corresponding profiles of μ̄′(ζ) and μ̄_p(ζ) are also plotted for comparison. The deviation between μ̄′(ζ) and μ̄_p(ζ) is apparent as the size of the particles and/or the concentration of particles is reduced, such as the case of Fig. 1(a). For all six cases, Eq. (18) gives a sound estimation of μ̄(ζ) that corresponds closely to the result from Hydrolight.

Fig. 1. Illustration of the validity of Eq. (18). The computational conditions are set as typical Case 1 waters with θ _s=30°, Cloudiness=0.0, V_wind =5.0 ms^-1, and (a) Chl=l.0 mg∙m^-3, BF_p =0.00915, (b) Chl=10.0 mg∙m^-3, BF_p =0.00915, (c) Chl=1.0 mg∙m^-3, BF_p =0.0183, (d) Chl=10.0 mg∙m^-3, BF_p =0.0183, (e) Chl=l.0 mg∙m^-3, BF_p =0.0366, (f) Chl=l0.0 mg∙m^-3, BF_p =0.0366. The values of the atmospheric parameters used are the default settings in Hydrolight. Both Hydrolight and the new model are employed to simulate μ̄(ζ).

3.4 Phase function of surrogate particles in Case 2 waters

Based on their optical characteristics, the Case 2 waters are generally categorized into four components with different parameterizations of IOPs, including pure water, large particles (assumed to absorb light like Chl and scatter light like large particles), CDOM (assumed to absorb light but not scatter it) and small particles. Note the component of small particles represents minerals in waters that can both absorb and scatter light. The amount of each component does not need to be controlled by Chl and the backscattering fraction for large particles BF_l and small particles BF_s do not need to be constant values either. To extend our approach to the general situation of Case 2 waters, it is equivalent to find a surrogate particle with backscattering fraction BF_p that is able to exhibit the same optical properties as those two types of particles presented. The normalized phase scattering function β̃[1

C D. Mobley, Light and water: radiative transfer in natural waters (Academic Press, San Diego, CA, 1994), p. 592.

] for the four-component Case 2 waters is

\tilde{β} (ψ; λ) \equiv \frac{b_{w} (λ)}{b (λ)} {\tilde{β}}_{w} (ψ) + \frac{b_{l} (λ)}{b (λ)} {\tilde{β}}_{l} (ψ; λ) + \frac{b_{s} (λ)}{b (λ)} {\tilde{β}}_{s} (ψ; λ) .

(20)

where β̃_l (ψ,λ) and β̃_s(ψ,λ) can be obtained from the analytic Fournier-Forand phase function [20

G. R. Fournier and J. L. Forand, ”Analytic phase function for ocean water,“ presented at the SPIE: Ocean Optics XII, 1994.

] by giving the values of BF_l and BF_s , respectively. Assuming a surrogate particle implies that

{\tilde{β}}_{p} (ψ; λ) = \frac{b_{l} (λ)}{b_{l} (λ) + b_{s} (λ)} {\tilde{β}}_{l} (ψ; λ) + \frac{b_{s} (λ)}{b_{l} (λ) + b_{s} (λ)} {\tilde{β}}_{s} (ψ; λ) .

(21)

Therefore, BF_p is the value that the resulting phase function gives the best fit to β̃_p(ψ,λ).

Fig. 2. (a) The analytical phase functions for ocean water proposed by Fournier and Forand [20

G. R. Fournier and J. L. Forand, ”Analytic phase function for ocean water,“ presented at the SPIE: Ocean Optics XII, 1994.

]. A large range in β(ψ) can be obtained by varying the value of BF. (b) An example of applying the technique of optimization to calculate BF_p for simulating the phase function of surrogate particles in Case 2 waters.

An example of calculating BF_p based on the technique of optimization is illustrated in Fig. 2. A set of β̃_p(ψ,λ) for ocean water is obtained by varying BF_p from 0.0001 to 0.5 (dashed lines). The thick solid line gives the corresponding β̃_p(ψ,λ) of a four-component Case 2 waters, with which BF_l =0.001, BF_s = 0.08 and ratio_l = 0.31. Note that

{ratio}_{l} \equiv \frac{b_{l} (λ)}{b_{l} (λ) + b_{s} (λ)} .

(22)

Based on the result of optimization, setting BF_p at a value of 0.053316 gives the best fit (thin solid line).

The process of searching for the optimized solution of BF_p consumes considerable computer time. To avoid repetitive calculation, the concept of constructing a LUT is again adopted. Three variables, ratio_l , BF_l and BF_s are selected, with ratio_l ranging from 0 to 1, BF_l ranging from 0.0001 to 0.02 to represent large particles, and BF_s ranging from 0.018 to 0.3 to represent small particles. A section of the LUT is denoted as LUT_BF and given in Table 3. Throughout the entire LUT_BF, the average value of PCC is 0.989577 and the lowest value is 0.930118. Note that the inelastic scattering effects were not included in the Hydrolight runs used to generate the three LUTs needed for the model of this paper.

Table 3. A section of the LUT_BF for quick reference to the value of backscattering fraction BF_p simulating the phase function of surrogate particles in Case 2 waters. The LUT_BF is based on three variables: ratio_l (0.0 – 1.0), the backscattering fraction for large particles BF_l (0.0001 – 0.0 the backscattering fraction for small particles BF_s (0.018 – 0.3).

ratio_i	BF_t	BF_S	BF_p	*PCC*
0.31	0.001	0.018	0.012388	0.999982
0.31	0.001	0.020	0.013687	0.999975
0.31	0.001	0.025	0.016998	0.999952
0.31	0.001	0.030	0.020336	0.999919
0.31	0.001	0.035	0.023545	0.999878
0.31	0.001	0.040	0.026882	0.999821
0.31	0.001	0.050	0.033383	0.999671
0.31	0.001	0.060	0.039866	0.999464
0.31	0.001	0.080	0.053316	0.998891
0.31	0.001	0.100	0.066865	0.998132
0.31	0.001	0.150	0.103043	0.994607
0.31	0.001	0.200	0.137565	0.990413
0.74	0.018	0.018	0.018000	1.000000
0.74	0.018	0.020	0.018520	1.000000
0.74	0.018	0.025	0.019801	1.000000
0.74	0.018	0.030	0.021077	1.000000
0.74	0.018	0.035	0.022331	0.999999
0.74	0.018	0.040	0.023562	0.999997
0.74	0.018	0.050	0.026020	0.999989
0.74	0.018	0.060	0.028401	0.999976
0.74	0.018	0.080	0.033082	0.999931
0.74	0.018	0.100	0.037621	0.999862
0.74	0.018	0.150	0.048219	0.999605
0.74	0.018	0.200	0.057609	0.999279

3.5 Use average cosine as an index to cope with stratified waters

The stratified water column can be divided into many thin layers, each with its own constitution of IOPs. Since a highly efficient solution method for homogeneous waters has been developed, it is advantageous overall to solve the RTE within each layer, and then to combine these “layer solutions” to obtain the solution of the RTE for an inhomogeneous body of water [1

C D. Mobley, Light and water: radiative transfer in natural waters (Academic Press, San Diego, CA, 1994), p. 592.

]. However, to determine the attenuation of radiative energy between each layer requires information about the radiance distribution at that particular depth. Recall that μ̄(ζ) is a crude but useful one-parameter measure of the directional structures of the light field in spatially uniform waters. It can be used as an index to provide rough information about the radiance distribution at that depth, and hence the attenuation of light within that layer. This idea is implemented and examined in Fig. 3.

Fig. 3. Illustration of using average cosine as an index to cope with stratified waters. The computational conditions are θ_S = 30°, Cloud = 0.0, λ = 440 nm and V_wind = 5m∙s^-1. The atmospheric parateters are those default settings in Hydrolight without consideration of inelastic scattering. (a) A profile of chlorophyll with one meter intervals. (b) A set of μ̄(z) that are calculated by considering eight cases of homogeneous bodies of water, each with one of those chlorophyll concentrations under the same conditions of incident light. (c) Combining these “layer solutions” of μ̄(z) to approximate the full solution of μ̄(z) for the corresponding case of an inhomogeneous body of water under the same conditions of incident light. (d) Comparison of the approximated E ₀ to the Hydrolight simulated E ₀.

Lewis et al. [25

M. R. Lewis, J. J. Cullen, and T. Platt, ”Phytoplankton and thermal structure in the upper ocean: consequences of nonuniformity in chlorophyll profile,“ J. Geoph. Res. 88, 2565–2570 (1983). [CrossRef]

] suggested that a typical profile of chlorophyll concentration for Case 1 waters can be approximated as a background value plus a Gaussian

C (z) = C_{o} + \frac{h}{s \sqrt{2 π}} \exp [- \frac{1}{2} {(\frac{z - z_{max}}{s})}^{2}] .

(23)

Similar parameterization of vertical profile was employed by Albert and Mobley [7

] as well. Figure 3(a) shows a profile of chlorophyll with one meter intervals, which was obtained by setting in Eq. (23) the parameters C_o , h, s, z _max values of 1, 20, 1.5, 7, respectively. This profile can also be treated as a combination of stratified layers with a total of eight concentrations of chlorophyll. Considering eight cases of homogeneous water bodies, each has one of those chlorophyll concentrations under the same conditions of incident light (as listed in caption). A set of μ̄_i(z) is calculated and illustrated in Fig. 3(b), where the subscript i is the index of the layer. These “layer solutions” of μ̄_i(z) can be combined to approximate the full solution of μ̄(z) for the corresponding case of inhomogeneous bodies of water under the same conditions of incident light (Fig. 3c). Note that the deviation of μ̄(z) between the approximation and the full Hydrolight run of the inhomogeneous case is small, and that the deviation of E ₀(z) is even lower (Fig. 3d). This verifies that the average cosine is a good index to extend the results of earlier work on homogeneous waters to the case of stratified waters.

Table 4. The randomly specified values of parameters for the 20 cases that are used in the model-to-model comparison to examine the accuracy, flexibility and applicability of the new model. The parameters include the solar zenith angle θ_S , cloudiness, surface wind speed V_wind , the backscattering fraction for large particles BF_l , the backscattering fraction for mineral particles BF_S , and the exponential coefficient of CDOM absorption γ.

No.	θ_S	Cloudiness	V_wind	BF_l	BF_S	γ
1	36.4	0.59	7.2	0.0056	0.0325	0.0100
2	7.5	0.99	6.7	0.0074	0.0291	0.0105
3	27.3	0.78	12.0	0.0076	0.0225	0.0127
4	38.0	1.00	15.0	0.0100	0.0754	0.0152
5	12.3	0.59	10.4	0.0039	0.0385	0.0101
6	27.4	0.57	5.4	0.0072	0.0721	0.0136
7	22.7	0.15	14.2	0.0146	0.0695	0.0130
8	8.0	0.82	7.1	0.0055	0.0323	0.0149
9	1.6	0.52	9.9	0.0074	0.0794	0.0107
10	13.4	0.74	8.5	0.0088	0.0903	0.0127
11	28.2	0.66	3.0	0.0103	0.0800	0.0132
12	33.6	0.25	2.2	0.0048	0.0622	0.0150
13	3.3	0.44	3.0	0.0049	0.0581	0.0145
14	15.1	0.18	14.9	0.0039	0.0322	0.0142
15	9.1	0.63	9.1	0.0045	0.0484	0.0114
16	40.7	0.24	2.8	0.0046	0.0589	0.0111
17	39.0	0.91	9.2	0.0115	0.0605	0.0148
18	14.7	0.08	1.1	0.0046	0.0296	0.0151
19	43.2	0.87	10.9	0.0081	0.0788	0.0144
20	22.7	0.23	4.4	0.0080	0.0260	0.0139

4. Results

A comprehensive model-to-model comparison was made to examine the accuracy, flexibility and applicability of the new model that includes all the aforementioned strategies. A total of 20 cases were compared by randomly specifying values to parameters, including the solar zenith angle θ_S , cloudiness, surface wind speed V_wind , the backscattering fraction for large particles BF_l , and the backscattering fraction for mineral particles BF_S (Table 4). The stratified layers were composed by randomly specifying values to four parameters of Eq. (23) to describe the vertical profiles of chlorophyll concentration Chl(z) and the mineral particle concentration M(z), respectively (Table 5). The CDOM absorption at a reference wavelength a _g,440(z) was set to be varied linearly from the surface value

a_{g, 440}^{Surface}

to the bottom value

a_{g, 440}^{Bottom}

, as listed in Table 5. Note that the spectral values of CDOM absorption are approximated as an exponential decay function

a_{g} (z, λ) = a_{g, 440} (z) \exp [- γ (λ - 440)],

(24)

where the value of γ is randomly specified as well. Figure 4(a), 4(b) and 4(c) show the resulting profiles of Chl(z), a _g,440(z) and M(z), respectively.

Table 5. The randomly specified values of parameters in Eq. (23) for describing the vertical profiles of chlorophyll concentration Chl(z) and the mineral particle concentration M(z), respectively. The CDOM absorption at a reference wavelength a _g,440(z) is set to be varied linearly from the surface value

a_{g, 440}^{Surface}

to the bottom value

a_{g, 440}^{Bottom}

, which are specified randomly as well.

No.	Chlorophyll				CDOM		Mineral particle
No.	C_o	h	s	z _max	$a_{g, 440}^{Surface}$	$a_{g, 440}^{Bottom}$	C_o	h	s	z _max
1	0.70	35.84	24.68	29.86	0.05	0.26	0.21	5.14	9.12	0.60
2	0.24	0.19	0.27	15.12	0.16	0.17	0.18	6.07	4.99	26.52
3	1.04	12.08	26.28	29.07	0.29	0.28	0.16	1.42	13.86	9.41
4	1.22	15.70	7.99	11.89	0.25	0.01	0.11	0.93	20.32	2.25
5	1.68	29.06	14.55	8.22	0.22	0.14	0.14	9.49	22.33	4.33
6	0.30	9.00	12.75	32.11	0.16	0.30	0.23	3.46	5.07	26.29
7	0.28	36.20	20.79	12.12	0.13	0.02	0.29	6.83	4.60	35.09
8	0.31	20.16	21.96	16.22	0.08	0.17	0.20	7.56	21.66	19.01
9	0.85	4.19	28.48	36.85	0.16	0.10	0.14	3.75	25.41	12.68
10	0.39	30.45	25.18	15.91	0.15	0.27	0.01	9.95	17.18	2.02
11	1.68	4.93	3.30	29.72	0.09	0.28	0.09	3.36	4.21	29.32
12	0.00	2.44	24.19	34.10	0.06	0.03	0.17	0.14	3.41	18.18
13	1.39	11.61	13.10	9.30	0.17	0.16	0.19	1.60	15.12	38.52
14	0.91	39.92	2.93	25.01	0.03	0.13	0.28	0.48	26.84	11.60
15	0.90	18.65	17.93	25.39	0.26	0.25	0.19	7.21	16.97	15.01
16	1.21	27.94	17.54	14.05	0.15	0.02	0.22	6.12	18.61	27.64
17	1.46	1.73	20.03	39.06	0.09	0.17	0.09	1.74	3.26	34.76
18	1.28	32.80	16.35	17.93	0.12	0.09	0.14	5.01	4.58	12.92
19	0.60	37.76	3.82	2.63	0.24	0.16	0.18	9.56	2.17	35.02
20	1.84	22.05	19.88	4.58	0.15	0.11	0.15	7.93	15.28	15.29

The historical chlorophyll-based bio-optical models were employed to calculate the absorption coefficient of the chlorophyll-bearing particles [26

L. Prieur and S. Sathyendranath, ”An optical classification of coastal and oceanic waters based on the specific spectral absorption curves of phytoplankton pigments, dissolved organic matter, and other particulate materials,“ Limnol. Oceanogr. 26, 671–689 (1981). [CrossRef]

, 27

A. Morel, ”Light and marine photosynthesis: a spectral model with geochemical and climatological implications,“ Prog. Oceanogr. 26, 263–306 (1991). [CrossRef]

]

a_{c} (z; λ) = 0.06 a_{c}^{*} (λ) Chl {(z)}^{0.65},

(25)

and the scattering coefficient [27

A. Morel, ”Light and marine photosynthesis: a spectral model with geochemical and climatological implications,“ Prog. Oceanogr. 26, 263–306 (1991). [CrossRef]

, 28

H. R. Gordon, D. K. Clark, J. W. Brown, O. B. Brown, R. H. Evans, and W. W. Broenkow, ”Phytoplankton pigment concentrations in the Middle Atlantic Bight: Comparison of ship determinations and CZCS estimates,“ Appl. Opt. 22, 20–36 (1983). [CrossRef] [PubMed]

]

b_{c} (z; λ) = 0.30 (\frac{550}{λ}) Chl {(z)}^{0.62},

(26)

where the normalized chlorophyll-specific absorption coefficient

a_{c}^{*}

(λ) is given in the paper of Morel [29

A. Morel, ”Optical modelling of the upper ocean in relation to its biogenous matter content (case 1 water),“ J. Geoph. Res. 93, 10749–10768 (1988). [CrossRef]

]. The spectral absorption a_w (λ) and scattering coefficient b_w (λ) for pure water are given in the papers of Pope and Fry [30

R. M. Pope and E. S. Fry, ”Absorption spectrum (380-700 nm) of pure water. II. Integrating cavity measurements,“ Appl. Opt. 36, 8710–8723 (1997). [CrossRef]

] and Smith and Baker [31

R. C. Smith and K. Baker, ”Optical properties of the clearest natural waters,” Appl. Opt. 20, 177–184 (1981). [CrossRef] [PubMed]

], respectively. The mineral absorption coefficient is given by

a_{m} (z; λ) = M (z) a_{m}^{*} (λ),

(27)

and the mineral scattering coefficient is

b_{m} (z; λ) = M (z) b_{m}^{*} (λ) .

(28)

Where the mass-specific mineral absorption coefficient

a_{m}^{*}

(λ) and scattering coefficient

b_{m}^{*}

(λ) are taken from published data [32

R. P. Bukata, J. H. Jerome, K. Y. Kondratyev, and D. V. Pozdnyakov, Optical properties and remote sensing of inland and coastal waters (CRC Press, 1995), p. 384.

]. With the set of bio-optical models described from Eq. (24) to Eq. (28), the vertical profiles of Chl(z), a _g,440(z) and M(z) can be converted to the vertical profiles of total absorption coefficient a(z,λ) and total scattering coefficient b(z,λ).

Fig. 4. Twenty cases of stratified Case 2 waters with vertical profiles of (a) Chl, (b) a _g,440 and (c) M. The corresponding profiles of IOPs at wavelength 440 nm are illustrated as (d) a, (e) b and (f) BF_p .

Given the values b_m (z,λ) and b_c (z,λ), the vertical profile of ratio_l (z,λ) can be calculated from Eq. (20). Following the procedure described earlier, the backscattering fraction BF_p (z,λ) of surrogate particles can be quickly referred to in the LUT from the given values of ratio_l (z,λ), BF_l and BF_S . For illustration, Fig. 4(d), 4(e) and 4(f) show the vertical profiles of a(z), b(z) and BF_p (z) respectively at wavelength 440 nm.

For each of the 20 cases described above, the vertical profiles of wavelength-integrated scalar irradiance in the Photosynthetical Active Radiation (PAR) range E _0,PAR(z) were simulated by use of the commercially available 4.3 version of the Hydrolight model and the new model developed in this work, respectively. A total of 35 wave bands were simulated, which ranged from 350 to 700 nm in steps of 10 nm. Figure 5 shows the comparison of E _0,PAR(z) from 0 to 50 m with one meter intervals. To quantitate the comparison, the correlation coefficient r, the maximum relative error ε _max, the average relative error ε _average and the percentage error ε% in the euphotic zone were calculated and given in Fig. 5. The definition of ε% is given by

ε_{%} = 10^{{RMSE}_{\log 10}} - 1,

(29)

where the root-mean-square-error in log 10 scale RMSE_log10 is defined by

{RMSE}_{\log 10} = \sqrt{\frac{\sum_{n = 1}^{N} {(\log_{10} E_{0, PAR}^{Model} - \log_{10} E_{0, PAR}^{Hydrolight})}^{2}}{N}} .

(30)

This kind of analysis puts the same emphasis on underestimates and overestimates. Comparing the results of our model to the results of the Hydrolight run, a very high correlation (r = 0.999924) as well as a large computational speed ratio (CSR) of 1402.8 is obtained for our model. The percentage error ε% is 2.03% and the maximum relative error ε _max is not more than 6.81%. The average relative error ε _average is as low as 1.62%. Figure 5 demonstrates that the new model can indeed provide an accurate and fast simulation of the vertical profile of E _0,PAR(z) for stratified Case 2 waters with a variety of compositions of the inherent optical properties.

5. Summary

This paper is devoted to the derivation of a fast and accurate model of scalar irradiance for stratified Case 2 waters. Over 60% of the human population and 90% of the world’s fish catch are directly linked to the coastal zone of Case 2 waters. The study of the radiative transfer in Case 2 waters is yet more challenging because of the complex constitution of their optical components. By definition, phytoplankton with their accompanying and covarying retinue of material of biological origin, as well as inorganic particles in suspension and yellow substances all influence the optical properties of Case 2 waters. Although the three-dimension MC optical model based on the backward ray-tracing technique is able to provide the solution to the most general form of radiative transfer equation, the computational resources required mean it is impractical to implement. A much faster calculation can be achieved by the commercially-available Hydrolight model that is based on the technique of invariant imbedding and the assumption of plane-parallel waters. However, it is not always practical to apply these numerical models to calculate global underwater light fields; the computational time required for running a full simulation interactively in biogeochemical models or global circulation models is often too long.

Extending earlier approaches for the case of homogeneous Case 1 waters, this paper describes a new model that runs more than 1,400 times faster than the full Hydrolight code, while it limits the percentage error to 2.03% and the maximum error to less than 6.81%. Having removed all limitations of earlier works, this new model is now fully comparable to the Hydrolight model in simulating the underwater scalar irradiance in Case 2 waters, yet still retains a very high computational speed. This new model can be used interactively in models of the oceanic system, such as biogeochemical models or the heat budget part of global circulation models.

Fig. 5. Comparison of accuracy and speed in simulating E ^0,PAR(z) (W m^-2) between the new model and Hydrolight. A very high correlation (r = 0.999924) as well as a large CSR of 1402.8 was obtained for our model. The percentage error ε% is 2.03% and the maximum relative error ε _max is not more than 6.81%.

Acknowledgments

This research was supported by National Science Council of Taiwan through grants NSC 94-2611-M-006-002. Thanks to Dr. Richard Miller for helps in Hydrolight simulation while the author held a National Research Council Research Associateship at NASA John C. Stennis Space Center.

References and links

1.	C D. Mobley, Light and water: radiative transfer in natural waters (Academic Press, San Diego, CA, 1994), p. 592.
2.	C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, “Fast and accurate model of underwater scalar irradiance,” Appl. Opt. 41, 4962–4974 (2002). [CrossRef] [PubMed]
3.	A. Morel and L. Prieur, “Analysis of variations in ocean color,” Limnol. Oceanogr. 22, 709–722 (1977). [CrossRef]
4.	H. R. Gordon and A. Y. Morel, Remote assessment of ocean color for interpretation of satellite visible imagery: a review , Lecture Notes on Coastal and Estuarine Studies (Springer-Verlag, New York, 1983), Vol. 4, p. 114.
5.	J. C. Pernetta and J. D. Milliman, “Land-ocean interactions in the coastal zone implementation plan,” IGBP Report 33 (Stockholm, 1995).
6.	S. Sathyendranath, ed., Remote sensing of ocean colour in coastal, and other optically-complex, waters , IOCCG Report (MacNab Print, Dartmouth, Canada, 2000), Vol. 3, p. 140.
7.	A. Albert and C. D. Mobley, “An analytical model for subsurface irradiance and remote sensing reflectance in deep and shallow case-2 waters,” Opt. Express 11, 2873–2890 (2003). [CrossRef] [PubMed]
8.	R. W. Preisendorfer, Hydrologic optics (U. S. Department of Commerce, Seattle, WA, 1976).
9.	R. L Fante, “Relationship between radiative-transport theory and Maxwell’s equations in dielectric media,” J. Opt. Soc. Am. 71, 460–468 (1981).
10.	R. M. Measures, Laser remote sensing: fundamentals and applications (Kreiger, Malabar, Florida, 1992), p. 510.
11.	K. Stamnes, “The theory of multiple scattering of radiation in plane parallel atmospheres,” Rev. Geophys. 24, 299–310 (1986). [CrossRef]
12.	G. N. Plass and G. W. Kattawar, “Monte Carlo calculations of radiative transfer in the Earth’s atmosphere-ocean system. I. Flux in the atmosphere and ocean,” J. Phys. Oceanogr. 2, 139–145 (1972). [CrossRef]
13.	C.-C. Liu and J. Woods, “Prediction of ocean colour: Monte-Carlo simulation applied to a virtual ecosystem based on the Lagrangian Ensemble method,” Int. J. Remote Sens. 25, 921–936 (2004). [CrossRef]
14.	K. L. Carder, C.-C. Liu, Z. P. Lee, D. C. English, J. Patten, R. F. Chen, J. E. Ivey, and C. O. Davis, “Illumination and turbidity effects on observing faceted bottom elements with uniform Lambertian albedos,” Limnol. Oceanogr. 48, 355–363 (2003). [CrossRef]
15.	R. W. Preisendorfer and C. D. Mobley, “Unpolarized irradiance reflectances and glitter patterns of random capillary waves on lakes and seas, by Monte Carlo simulation,” NOAA Technical Memorandum, ERL PMEL-63 (NOAA Pacific Marine Environmental Laboratory, Seattle, WA, 1985).
16.	C. D. Mobley and R. W. Preisendorfer, “A numerical model for the computation of radiance distributions in natural waters with wind-roughened surfaces,” NOAA Technical Memorandum, ERL PMEL-75 (NOAA Pacific Marine Environmental Laboratory, Seattle, WA, 1988).
17.	C. D. Mobley, “A numerical model for the computation of radiance distributions in natural waters with wind-roughened surfaces,” Limnol. Oceanogr. 34, 1473–1483 (1989). [CrossRef]
18.	C. D. Mobley, “Hydrolight 3.0 Users’ Guide,” (SRI International, Menlo Park, CA, 1995).
19.	C.-C. Liu, J. D. Woods, and C. D. Mobley, ”Optical model for use in oceanic ecosystem models,“ Appl. Opt. 38, 4475–4485 (1999). [CrossRef]
20.	G. R. Fournier and J. L. Forand, ”Analytic phase function for ocean water,“ presented at the SPIE: Ocean Optics XII, 1994.
21.	N. J. McCormick, ”Mathematical models for the mean cosine of irradiance and the diffuse attenuation coefficient,“ Limnol. Oceanogr. 40, 1013–1018 (1995). [CrossRef]
22.	A. Gershun, ”The light field,“ J. Math. Phys. 18, 51–151 (1939).
23.	V. I. Haltrin, ”Chlorophyll-Based Model of Seawater Optical Properties,“ Appl. Opt. 38, 6826–6832 (1999). [CrossRef]
24.	C. D. Mobley, L. K. Sundman, and E. Boss, ”Phase function effects on oceanic light fields,“ Appl. Opt. 41, 1035–1050 (2002). [CrossRef] [PubMed]
25.	M. R. Lewis, J. J. Cullen, and T. Platt, ”Phytoplankton and thermal structure in the upper ocean: consequences of nonuniformity in chlorophyll profile,“ J. Geoph. Res. 88, 2565–2570 (1983). [CrossRef]
26.	L. Prieur and S. Sathyendranath, ”An optical classification of coastal and oceanic waters based on the specific spectral absorption curves of phytoplankton pigments, dissolved organic matter, and other particulate materials,“ Limnol. Oceanogr. 26, 671–689 (1981). [CrossRef]
27.	A. Morel, ”Light and marine photosynthesis: a spectral model with geochemical and climatological implications,“ Prog. Oceanogr. 26, 263–306 (1991). [CrossRef]
28.	H. R. Gordon, D. K. Clark, J. W. Brown, O. B. Brown, R. H. Evans, and W. W. Broenkow, ”Phytoplankton pigment concentrations in the Middle Atlantic Bight: Comparison of ship determinations and CZCS estimates,“ Appl. Opt. 22, 20–36 (1983). [CrossRef] [PubMed]
29.	A. Morel, ”Optical modelling of the upper ocean in relation to its biogenous matter content (case 1 water),“ J. Geoph. Res. 93, 10749–10768 (1988). [CrossRef]
30.	R. M. Pope and E. S. Fry, ”Absorption spectrum (380-700 nm) of pure water. II. Integrating cavity measurements,“ Appl. Opt. 36, 8710–8723 (1997). [CrossRef]
31.	R. C. Smith and K. Baker, ”Optical properties of the clearest natural waters,” Appl. Opt. 20, 177–184 (1981). [CrossRef] [PubMed]
32.	R. P. Bukata, J. H. Jerome, K. Y. Kondratyev, and D. V. Pozdnyakov, Optical properties and remote sensing of inland and coastal waters (CRC Press, 1995), p. 384.

OCIS Codes
(010.4450) Atmospheric and oceanic optics : Oceanic optics
(010.7340) Atmospheric and oceanic optics : Water
(030.5620) Coherence and statistical optics : Radiative transfer

ToC Category:
Atmospheric and Oceanic Optics

History
Original Manuscript: December 15, 2005
Revised Manuscript: February 19, 2006
Manuscript Accepted: February 23, 2006
Published: March 6, 2006

Virtual Issues
Vol. 1, Iss. 4 Virtual Journal for Biomedical Optics

Citation
Cheng-Chien Liu, "Fast and accurate model of underwater scalar irradiance for stratified Case 2 waters," Opt. Express 14, 1703-1719 (2006)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-14-5-1703

Sort: Author | Year | Journal | Reset

References

C. D. Mobley, Light and water: radiative transfer in natural waters (Academic Press, San Diego, CA, 1994), p. 592.
C.-C. Liu, K. L. Carder, R. L. Miller, and J. E. Ivey, "Fast and accurate model of underwater scalar irradiance," Appl. Opt. 41, 4962-4974 (2002). [CrossRef] [PubMed]
A. Morel and L. Prieur, "Analysis of variations in ocean color," Limnol. Oceanogr. 22, 709-722 (1977). [CrossRef]
H. R. Gordon and A. Y. Morel, Remote assessment of ocean color for interpretation of satellite visible imagery: a review, Lecture Notes on Coastal and Estuarine Studies (Springer-Verlag, New York, 1983), Vol. 4, p. 114.
J. C. Pernetta and J. D. Milliman, "Land-ocean interactions in the coastal zone implementation plan," IGBP Report 33 (Stockholm, 1995).
S. Sathyendranath, ed., Remote sensing of ocean colour in coastal, and other optically-complex, waters, IOCCG Report (MacNab Print, Dartmouth, Canada, 2000), Vol. 3, p. 140.
A. Albert and C. D. Mobley, "An analytical model for subsurface irradiance and remote sensing reflectance in deep and shallow case-2 waters," Opt. Express 11, 2873-2890 (2003). [CrossRef] [PubMed]
R. W. Preisendorfer, Hydrologic optics (U. S. Department of Commerce, Seattle, WA, 1976).
R. L. Fante, "Relationship between radiative-transport theory and Maxwell's equations in dielectric media," J. Opt. Soc. Am. 71, 460-468 (1981).
R. M. Measures, Laser remote sensing: fundamentals and applications (Kreiger, Malabar, Florida, 1992), p. 510.
K. Stamnes, "The theory of multiple scattering of radiation in plane parallel atmospheres," Rev. Geophys. 24, 299-310 (1986). [CrossRef]
G. N. Plass and G. W. Kattawar, "Monte Carlo calculations of radiative transfer in the Earth's atmosphere-ocean system. I. Flux in the atmosphere and ocean," J. Phys. Oceanogr. 2, 139-145 (1972). [CrossRef]
C.-C. Liu and J. Woods, "Prediction of ocean colour: Monte-Carlo simulation applied to a virtual ecosystem based on the Lagrangian Ensemble method," Int. J. Remote Sens. 25, 921-936 (2004). [CrossRef]
K. L. Carder, C.-C. Liu, Z. P. Lee, D. C. English, J. Patten, R. F. Chen, J. E. Ivey, and C. O. Davis, "Illumination and turbidity effects on observing faceted bottom elements with uniform Lambertian albedos," Limnol. Oceanogr. 48, 355-363 (2003). [CrossRef]
R. W. Preisendorfer and C. D. Mobley, "Unpolarized irradiance reflectances and glitter patterns of random capillary waves on lakes and seas, by Monte Carlo simulation," NOAA Technical Memorandum, ERL PMEL-63 (NOAA Pacific Marine Environmental Laboratory, Seattle, WA, 1985).
C. D. Mobley and R. W. Preisendorfer, "A numerical model for the computation of radiance distributions in natural waters with wind-roughened surfaces," NOAA Technical Memorandum, ERL PMEL-75 (NOAA Pacific Marine Environmental Laboratory, Seattle, WA, 1988).
C. D. Mobley, "A numerical model for the computation of radiance distributions in natural waters with wind-roughened surfaces," Limnol. Oceanogr. 34, 1473-1483 (1989). [CrossRef]
C. D. Mobley, "Hydrolight 3.0 Users' Guide," (SRI International, Menlo Park, CA, 1995).
C.-C. Liu, J. D. Woods, and C. D. Mobley, "Optical model for use in oceanic ecosystem models," Appl. Opt. 38, 4475-4485 (1999). [CrossRef]
G. R. Fournier and J. L. Forand, "Analytic phase function for ocean water," presented at the SPIE: Ocean Optics XII, 1994.
N. J. McCormick, "Mathematical models for the mean cosine of irradiance and the diffuse attenuation coefficient," Limnol. Oceanogr. 40, 1013-1018 (1995). [CrossRef]
A. Gershun, "The light field," J. Math. Phys. 18, 51-151 (1939).
V. I. Haltrin, "Chlorophyll-Based Model of Seawater Optical Properties," Appl. Opt. 38, 6826-6832 (1999). [CrossRef]
C. D. Mobley, L. K. Sundman, and E. Boss, "Phase function effects on oceanic light fields," Appl. Opt. 41, 1035-1050 (2002). [CrossRef] [PubMed]
M. R. Lewis, J. J. Cullen, and T. Platt, "Phytoplankton and thermal structure in the upper ocean: consequences of nonuniformity in chlorophyll profile," J. Geoph. Res. 88, 2565-2570 (1983). [CrossRef]
L. Prieur and S. Sathyendranath, "An optical classification of coastal and oceanic waters based on the specific spectral absorption curves of phytoplankton pigments, dissolved organic matter, and other particulate materials," Limnol. Oceanogr. 26, 671-689 (1981). [CrossRef]
A. Morel, "Light and marine photosynthesis: a spectral model with geochemical and climatological implications," Prog. Oceanogr. 26, 263-306 (1991). [CrossRef]
H. R. Gordon, D. K. Clark, J. W. Brown, O. B. Brown, R. H. Evans, and W. W. Broenkow, "Phytoplankton pigment concentrations in the Middle Atlantic Bight: Comparison of ship determinations and CZCS estimates," Appl. Opt. 22, 20-36 (1983). [CrossRef] [PubMed]
A. Morel, "Optical modelling of the upper ocean in relation to its biogenous matter content (case 1 water)," J. Geoph. Res. 93, 10749-10768 (1988). [CrossRef]
R. M. Pope and E. S. Fry, "Absorption spectrum (380-700 nm) of pure water. II. Integrating cavity measurements," Appl. Opt. 36, 8710-8723 (1997). [CrossRef]
R. C. Smith and K. Baker, "Optical properties of the clearest natural waters," Appl. Opt. 20, 177-184 (1981). [CrossRef] [PubMed]
R. P. Bukata, J. H. Jerome, K. Y. Kondratyev, and D. V. Pozdnyakov, Optical properties and remote sensing of inland and coastal waters (CRC Press, 1995), p. 384.

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper