Graziani F. (editor) Computational Methods in Transport

Подождите немного. Документ загружается.

Transport Theory for Optical Oceanography

N.J. McCormick

University of Washington, Mechanical Engineering Department, Seattle,

Washington 98195-2600

∗∗

Abstract. A general introduction to the ﬁeld of ocean optics is presented, includ-

ing features that make optical oceanography problems similar to and diﬀerent from

other transport problems. Methods for solving time-independent, one-dimensional

problems are discussed that are appropriate for passive illumination conditions.

1 Introduction

1.1 Background on Optical Oceanography

Non-military applications of ocean optics are a sub-discipline of biological

oceanography because photon absorption drives photosynthesis which drives

phytoplankton (“primary”) production at the bottom of the food chain. The

correlation of optical radiance measurements from satellites with other ocean

physics measurements (wave characteristics, temperature, salinity, etc.) is be-

coming more commonplace with recent enhanced eﬀorts in ecological moni-

toring. Such applications include both remote and in situ sensing of ocean wa-

ters, typically done for a passive time-independent surface illumination, which

is the emphasis of the discussion here. The solution of the time-dependent

transport equation is needed, on the other hand, for active illumination imag-

ing applications such as mine detection or communication by pulsed optical

signal propagation.

The types of in-water problems being solved with the classic radiative

transfer equation include a) forward problems for determining the angle-

dependent and angle-integrated light ﬁeld, and b) inverse problems for de-

termining, for example, the scattering and absorption properties in order to

monitor the biological primary production of water or the ecology of coral

reefs.

Good introductions to radiative transfer for oceanographic forward prob-

lem applications are available. The text by Kirk [Kir94] is for people more

interested in the biological aspects of radiative transfer. Mobley [Mob94] gives

details of how ocean optical transfer calculations can be done, as well as a

physicist’s approach to biological oceanography. The text by Thomas and

∗∗

mccor@u.washington.edu; http://faculty.washington.edu/mccor/

152 N.J. McCormick

Stamnes [TS99] and the monograph by Walker [Wal94] each contain an in-

troduction to atmospheric and oceanic optics, with the former giving more

details about atmospheric line-by-line radiative transfer calculations and the

latter containing an extensive coverage of sea-surface and refracted light sta-

tistics. The monograph by Bukata et al. [BJK95] pertains to water optical

properties. For oceanographic inverse problem applications, an excellent re-

view article is that by Gordon [Gor02].

The key features of optical transport in seawater are:

1. Light in the 400–700 nm range is of primary interest, where most biolog-

ical activity occurs and the optical transmission is the greatest.

2. Ocean waters are diﬃcult to characterize because

• they consist of water plus dissolved organic and inorganic matter not

well-characterized as to type and location,

• the “microscopic cross sections,” from which the eﬀective absorption

and scattering properties are determined, can only be idealized (e.g.,

with spherical or other regular shapes),

• the air-water interface condition is diﬃcult to simulate because of

refractive eﬀects and surface waves.

3. Ocean waters usually are modeled with the plane parallel approximation

for natural light illuminations because

• the water column is assumed to be layered (i.e., there are minimal

horizontal variations of the water constituents),

• the bottom is assumed ﬂat and of uniform composition or very deep,

• the incident radiation from the atmosphere is assumed to be uniform

over the sea surface.

The reasons ocean waters are diﬃcult to characterize are:

1. The wavelength-dependent scattering of phytoplankton depends on the

composition and particle size that can vary from approximately 0.7 µm

to ∼100 µm. An approximate model for the number size distribution n(x)

of biological particles with an equivalent diameter x, per unit volume in

open ocean waters, is [Mob94]

n(x) ∝ x

−s

for a slope of ln n(x)=−s with s ≈ 4.

2. The suspended sediments are predominately scatterers, rather than ab-

sorbers, although strong absorption features have been observed in iron-

rich sediment minerals.

3. Colored dissolved organic matter (CDOM) absorbs most strongly in ul-

traviolet wavelengths with various absorption peaks.

The similarities to other transport problems include:

1. The linear Boltzmann equation for neutral particles is the governing equa-

tion.

Transport Theory for Optical Oceanography 153

2. Active illumination problems (e.g., with a laser beam) usually are three-

dimensional.

3. The coeﬃcients of individual, optically-active constituents of seawater

are assumed additive so, for example, for the absorption coeﬃcient a (in

units of m

−1

)is

total

= a

seawater

+ a

phytoplankton

+ a

other particles

+ a

CDOM

+ ··· .

4. Phytoplankton and CDOM have the ability to ﬂuoresce by re-emitting

light within distinct wavebands that are somewhat longer than the ab-

sorbed light, which means a coupled wavelength analysis may be neces-

sary.

5. Raman (inelastic) scattering events cause a wavelength increase (i.e.,

downscattering in energy).

6. Bioluminescence

is analogous to an external source from neutron-gamma

reactions in gamma transport or spontaneous ﬁssion in neutron transport

analyses.

The diﬀerences from most other linear transport problems include:

1. “Inverse problems” (e.g., to characterize properties of medium) are very

important, more so than for analyses of designed systems (e.g., nuclear

engineering applications).

2. A plane-parallel approximation usually is suﬃcient for waters illuminated

by natural light so the radiance L (i.e., the angular dependence of the

radiation ﬁeld) depends only on the water depth z and the polar angle

θ =cos

−1

µ and azimuthal angle ϕ.

3. The optical sensors most often used enable a simpliﬁed plane-parallel

approximation so the azimuthal angular dependence need not be deter-

mined.

4. The scattering is very strongly forward peaked (e.g., more like the scat-

tering of atmospheric aerosols than neutrons), with peak forward-to-

backward scattering ratios of O(10

–10

) [Pet72,Mob94].

5. The refractive index mis-match at the air-water interface causes angle-

dependent internal reﬂection for some surface emerging directions.

6. Air-water interface waves are hard to characterize.

7. Problems with polarization eﬀects have four radiance components for

each wavelength.

It is worth emphasizing that ocean waters are dynamically changing sys-

tems that are hard to characterize, and hence the challenge of obtaining good

input data for computations means the old adage “garbage in, garbage out”

is a big concern. Furthermore, the environment for performing optical exper-

iments in the ﬁeld is often diﬃcult, which leads to appreciable measurement

Bioluminescence results from an exothermic internal chemical reaction releas-

ing energy as light, either spontaneously or in response to mechanical stimulation,

and occurs with some species of bacteria, dinoﬂagellates, zooplankton, and ﬁsh.

154 N.J. McCormick

uncertainties. For these reasons computational results are not needed with

more precision than measurement uncertainties.

1.2 The Transfer Equation Used

Because the water optical properties typically vary little in the horizontal

direction, and because the surface illumination often is nearly uniform over

the sea surface, ocean optics problems often can be treated spatially with

only the depth variable provided the eﬀects of surface waves are averaged

out. In this case, in the customary notation of the ocean optics literature,

the governing Boltzmann integrodiﬀerential for the radiance (i.e., “angular

ﬂux”) L(z,µ,ϕ, λ) (in Wm

−2

−1

)is

∂L(z, µ, ϕ,λ)

∂z

+ c(z, λ)L(z, µ, ϕ, λ)

= b(z, λ)



2π



dϕ





−1

dµ



β(z, µ



,µ,ϕ



,ϕ,λ)L(z, µ, ϕ, λ)+S(z,λ) , (1)

where µ is the cosine of the polar angle between the downward distance z

from the sea surface and the direction of interest and ϕ is the azimuthal

angle in the horizontal plane. The sum of the absorption and scattering coef-

ﬁcients is the so-called beam attenuation coeﬃcient c(z,λ)=a(z, λ)+b(z, λ)

(in m

−1

). These coeﬃcients are “inherent optical properties” (IOP). The

scattering phase function is

β(z, µ



,µ,ϕ



,ϕ,λ) (in sr

−1

), also is an IOP. The

source from Raman scattering or ﬂuorescence is eﬀectively isotropic,

S(z,λ)=



∞

dλ





2π



dϕ





−1

dµ



k(z,µ



,ϕ



,λ



,λ)L(z,µ



,ϕ



,λ



)(2)

for a conversion kernel k(z,µ



,ϕ



,λ



,λ)(inm

−1

). Often it is suﬃ-

cient to use only a simpliﬁed transport equation that does not depend on the

source S(z, λ) because bioluminescent sources are negligibly small during the

daytime hours for most ocean waters and Raman scattering and ﬂuorescence

eﬀects also are small in comparison to the sea surface illumination except for

deeper depths (τ>∼30) and longer visible wavelengths (e.g., λ>∼520 nm).

In such cases the radiative transfer equation can be solved independently for

each wavelength λ.

Another simpliﬁcation that often can be made when attempting to predict

subsurface measurements is to use a transport equation that depends only

on µ and not ϕ, even though the sea surface illumination depends on both

angle variables. The reason the single-angle equation can be used is that the

sensors for passive (i.e., ambient) light typically measure only the “apparent

optical properties” of the downward and upward “planar irradiances” (i.e.,

“partial currents”)

Transport Theory for Optical Oceanography 155

(z,λ)=



2π

dϕ



dµ µL(z, µ, ϕ, λ) , (3)

(z,λ)=



2π

dϕ



dµ µL(z, −µ, ϕ, λ) , (4)

and the vertically upward radiance

L(z,−1,λ)=



2π

L(z,−1,ϕ,λ)dϕ . (5)

Such sensors depend on sunlight and they are restricted to measurements

down to a depth z ∼150 m, depending on the water properties.

Polarization eﬀects for the radiation can be important near the air-water

interface, in which case the source-free form of (1) must be replaced by

∂I(z, µ, ϕ, λ)

∂z

+ c(z, λ)I(z, µ, ϕ, λ)

= b(z, λ)



2π

dϕ





−1

dµ



Σ(z,µ



,µ,ϕ



,ϕ,λ)I(z, µ



,ϕ



,λ)+S(z,λ)(6)

for the Stokes vector I =[I,Q,U,V ]

, a diagonal interaction matrix c(z,λ),

and the Stokes scattering phase matrix

Σ(z,µ



,µ,ϕ



,ϕ,λ). The eﬀects of

polarization for ocean optics applications has been investigated by Kattawar

and collaborators [AK93, RK98].

1.3 Fourier Uncoupling of the Transfer Equation

Many ocean optics computational programs take advantage of the relatively

simple geometry to convert (1) that depends on µ and ϕ into a set of uncou-

pled equations that depend only on µ, especially when only the downward

and upward planar irradiances and the vertically upward radiance are needed.

Consider the equation

∂L(τ,µ,ϕ)

∂τ

+ L(τ, µ, ϕ)=ω



2π

dϕ





−1

dµ



β(µ, µ



,ϕ,ϕ



)L(τ,µ



,ϕ



)(7)

in terms of the optical distance variable τ =



c(z)dz, the single-scattering

albedo ω = b/c, and an implicit single wavelength λ. An expansion of the

phase function as

β(µ, µ



,ϕ,ϕ



)=(2π)

−1



m=0

(2 − δ

m,0

)

(µ, µ



)cosm(ϕ −ϕ



)(8)

gives the set of uncoupled equations [Cha60,Mob94]

∂L

(τ,µ)

∂τ

+ L

(τ,µ)=ω



−1

(µ



,µ)L

(τ,µ



)dµ



(9)

156 N.J. McCormick

for m =0toM, with

L(τ,µ,ϕ)=



m=0

(2 − δ

m,0

(τ,µ)cosmϕ + L

(τ,µ) . (10)

Here L

(τ,µ) accounts for that portion of the radiance not included in the

expansion (e.g., that may arise from a Dirac-delta type surface illumina-

tion). The reduced phase function can be expanded in terms of the associated

Legendre functions P

(µ)as

(µ



,µ)=



n=m

(µ



(µ) , (11)

with f

=1and

(2n +1)

(n − m)!

(n + m)!

,n≥ m.

As mentioned previously, often just the simplest form of the transport

equation for m = 0 is needed,

∂L(τ,µ)

∂τ

+ L(τ,µ)=(ω/2)



n=0

(2n +1)f

(µ)



−1

(µ



)L(z,µ



)dµ



(12)

as when only E

(τ), E

(τ), L(τ,−1), and the “scalar irradiance” (i.e., “total

ﬂux”)

(τ,µ)=



2π

dϕ



−1

dµL(τ,µ,ϕ)

are desired.

2 Aspects Requiring Special Computational Attention

2.1 Optical Properties

“Case 1” ocean waters are those for which the inherent optical properties are

determined primarily by phytoplankton and co-varying CDOM and detritus.

(It should be noted that the label is not a synonym for open ocean waters,

although many such waters are of the case 1 type). Examples for correlations

of the absorption and scattering coeﬃcients with the non-water “particle”

concentration C in mg m

−3

are [PS81]

(λ, C)=0.06A

(λ)C

0.602

, (13)

(λ, C) ≈ 0.5B

(550)



550



0.62

,C<10 . (14)

Transport Theory for Optical Oceanography 157

(Here A

(λ) is a tabulated coeﬃcient normalized to the particle absorption

at 440 nm and B

(550) ranges from 0.12 to 0.45 m

−1

.) “Case 2” water optical

properties are signiﬁcantly inﬂuenced by CDOM, detritus, mineral particles,

bubbles, and other substances whose concentrations do not co-vary with the

phytoplankton concentration. These waters, not necessarily coastal, are a

signiﬁcant complication for ocean optics computations compared to the many

other ﬁelds of transport where the input material parameters follow regular

laws of physics.

A diﬃculty often also arises when attempting to describe the strongly

anisotropic scattering phase function. The classic experimental scattering

data of Petzold [Pet72,Mob94] traditionally has been used to model case

1 waters. Another approach is the “Fournier-Forand model” [FF94] that is

for an ensemble of particles that have a hyperbolic particle distribution given

by F (r) ∝ r

−k

for radius r and slope parameter k. For a phase function

scattering angle ψ, the phase function P (ψ) of that model is [MSB02]

P (ψ)=

4π(1 − η)

[ν(1 − η) − (1 −η

)



[η(1 − η

) − ν(1 − η)] sin

−2





1 − η

180

16π(η

180

− 1)η

180

(3 cos

ψ − 1) , (15)

where

ν =

3 − k

,η=

3(n − 1)

sin





. (16)

Here n is the real index of refraction of the particles and η

180

is η evaluated

at ψ = 180

◦

. Knowledge of the shape of the phase function for a scattering

angle of ≤ 15

◦

is not necessary in many applications, however [Gor93].

An examination of the eﬀects of diﬀerent scattering models on the com-

puted radiance has shown that [MSB02]:

• Use of the correct phase function is just as necessary as use of the correct

absorption and scattering coeﬃcients.

• An accurate shape of the phase function is not so important if the overall

shape does not diﬀer greatly from the correct shape, but the backscatter-

ing fraction must be nearly correct.

• The Fourier-Forand phase function works well.

2.2 Index of Refraction Eﬀects

Light rays change direction when passing through the air-water interface

because of Snell’s law,

sin θ

air

= n

sin θ

water

, (17)

158 N.J. McCormick

for a relative index of refraction of water with respect to air, n

≈ 4/3. Such

an interface phenomenon does not occur for most other transport problems

(e.g., those involving neutrons). To analyze those light rays passing from

air to water through a ﬂat interface, only a simple change of direction is

required to obtain the direction beneath the surface. But for those within

the water that reach the interface, only those with a polar angle of θ

water

≤

∼48.6

◦

with respect to the downward normal can pass through the interface

while those with ∼48.6

◦

≤ θ

water

≤ 90

◦

are totally reﬂected back into the

water. This means that numerical techniques that do not discretize the polar

angle variable, such as the spherical harmonics method, are cumbersome to

implement.

The mis-match of the diﬀerent indices of refraction for air and water tends

to aﬀect the radiance mostly near the surface. At a few optical depths beneath

the surface, where almost all photons have undergone one or more scatter-

ing interactions, the directional dependence of the radiance is dominated by

the optical properties more than the directional dependence of the surface

illumination. For deep, spatially uniform waters containing no sources (e.g.,

inelastic scattering) the radiance approaches an asymptotic regime where the

directional dependence does not change with increasing depth, although the

magnitude of the radiance decreases exponentially [McC92].

2.3 Sea Surface Waves and Wave Focusing Eﬀects

If the radiance is to be accurately computed near the surface, then the ef-

fects of waves also must be considered. This requires incorporating the index

of refraction direction changes at the point on the wave surface where the

ray strikes. The Monte Carlo technique, the logical method of choice, has

been implemented for this purpose [Mob94]. But a major problem is how to

describe the shape of the waves. Traditionally the classic Cox-Munk model

[CM54a,CM54b,Mob94,Wal94] has been used. It is a wave-slope wind-speed

correlation for capillary (i.e., small surface) waves. The dimensionless horizon-

tal wave slopes in the alongwind and crosswind directions vary independently

and are normally distributed with zero means, with variances proportional

to the wind speed U in ms

−1

= a

U, where a

=3.16 × 10

−3

−1

= a

U, where a

=1.92 × 10

−3

−1

s .

The wave elevation also is distributed normally.

Wind produces a frictional drag at the interface that transmits energy

to the sea surface, so when a steady wind starts then small capillary waves

develop ﬁrst. Energy then is transferred to the longer-wavelength gravity

waves and the spectrum energy grows until equilibrium is reached such that

the input wind energy balances the energy dissipation. For a radial spectrum

S(k) and an angular speading function Φ(k, ϕ), the “Elfouhaily” statistical

Transport Theory for Optical Oceanography 159

model gives a sea surface roughness spectrum Ψ(k, ϕ)thatistheFourier

transform of the autocorrelation function of the surface height [ECK96],

Ψ(k, ϕ)=k

−1

S(k)Φ(k, ϕ) ≈ (2πk)

−1

S(k)[1 + a

cos(2ϕ)] . (18)

Another consequence of a wind-roughened surface is that wave focusing

occurs near the surface and the plane geometry approximation is no longer

valid [ZBB01]. There may be signiﬁcant implications of this for remote sensing

applications, but for in-water measurements the eﬀects are damped out with

increasing optical depth, like those due to the air-water index of refraction

mis-match.

3 Computational Programs

An excellent comparison of diﬀerent computational methods used for oceano-

graphic analyses has been given by Mobley et al. [MGG94]. A recent but more

limited comparison also is available [CCG03].

3.1 HYDROLIGHT

R

This is a commercially available program [MS00] in Fortran for solving the

sub-surface light ﬁeld given an incident illumination on the sea surface. For

the unit sphere of directions oriented with a downward vertical axis, it sepa-

rates out the azimuthal angle dependence with the Fourier decomposition of

Sect. 1.3. The polar angle range 0 ≤ θ ≤ 180

◦

is divided into equally-spaced

segments of ∆θ, and equally-spaced segments also are used for ∆ϕ. The result

is a set of surface areas or “quads”.

The program uses the invariant imbedding procedure to solve (9) for each

azimuthal component of the light ﬁeld [Mob94]. This procedure has an advan-

tage for calculations involving many layers of distinct optical properties. A

key feature of the program is that it contains a Monte Carlo ray-tracing capa-

bility for simulating waves at the surface to generate the in-water boundary

condition just beneath the sea surface. The program has an extensive collec-

tion of built-in data sets to model a wide variety of water optical properties

that are mixed according to user-speciﬁed constituent concentrations in the

water. The PC-computer version has a graphical user interface with many

input options.

The program also provides a means of coupling atmospheric radiative

transport with the in-water light ﬁeld because the inherent optical properties

for both water and atmosphere are available to cover the wavelength range

300 nm to 1000 nm. A nice feature is that the Fortran source code of the

program is available for personal modiﬁcation (e.g., a spherically-symmetric

problem has even been solved with it [MK00]). Because the program has been

developed commercially, it is not inexpensive; the current price is $10K that

includes updated versions from the authors.

160 N.J. McCormick

3.2 DISORT

This program is based on the classic discrete ordinates method [STW88,JS94,

TS99] for solving the set of azimuthally-independent transfer equations (9),

with extra ordinates to allow for the change in direction across the air-water

interface. The program uses an “iteration of the source function method”

(i.e., “post processing”) to compute the radiance in directions other than

the ordinates speciﬁed in the program input. The newest version for ocean

applications is the Coupled Atmosphere-Ocean (CAO-DISORT) program

[GSH03,YK03]. A vector version, VDISORT, also is available [SS00] for in-

vestigating polarized radiation problems.

3.3 Analytic Discrete Ordinates Method

This is a discretized version of the analytic eigenmode expansion of the ho-

mogeneous radiative transfer equation [Sie00]. The advantage of preserving

that functional form of the solution is that highly accurate benchmark-

quality results can be obtained for the radiance. With this method the

variables τ and µ of (9) for the phase function of (11) are separated with

(τ,µ)=φ

(ν

,µ)exp(−τ/ν

). With a suppressed superscript m for φ

and ν, the eigenfunctions φ(ν, µ

) satisfy



1 −



φ(ν, µ



n=m

(µ

)



k=1

(µ

)[φ(ν, µ

)

+(−1)

n−m

φ(ν, −µ

)] , (19)





φ(ν, −µ



n=m

(µ

)



k=1

(µ

)[(−1)

n−m

φ(ν, µ

)

+φ(ν, −µ

)] (20)

for quadrature nodes µ

on [0, 1], i =1toK, and integration weights w

From these equations the eigenvalues ±ν

, j =1toJ, can be computed so

that the solution in vector form for a region 0 ≤ τ ≤ τ

can be expanded

with

(τ)=



j=1

(ν

−τ/ν

+ B

∓

(ν

−(τ

−τ)/ν

]+L

p,±

(τ) (21)

for

(τ)=[L(τ, ±µ

),L(τ,±µ

),...L(τ,±µ

]

(ν

)=[φ(ν, ±µ

),φ(ν, ±µ

),...,φ(ν, ±µ

max

)]

, (22)

where the A

and B

are expansion coeﬃcients that depend on the surface

illumination and L

p,±

(τ) is the particular solution [BGS00] incorporating