Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548

Contents lists available at SciVerse ScienceDirect

Journal of Quantitative Spectroscopy &

Radiative Transfer
journal homepage: www.elsevier.com/locate/jqsrt

Radiative transfer solutions for coupled atmosphere ocean systems

using the matrix operator technique
Andre Hollstein n, Jurgen
Fischer
Freie Universit
at, Institute for Space Sciences, Berlin, Germany

a r t i c l e in f o abstract

Article history: Accurate radiative transfer models are the key tools for the understanding of radiative
Received 7 November 2011 transfer processes in the atmosphere and ocean, and for the development of remote
Received in revised form sensing algorithms. The widely used scalar approximation of radiative transfer can lead
16 January 2012
to errors in calculated top of atmosphere radiances. We show results with errors in the
Accepted 17 January 2012
Available online 26 January 2012
order of 7 8% for atmosphere ocean systems with case one waters. Variations in sea
water salinity and temperature can lead to variations in the signal of similar magnitude.
Keywords: Therefore, we enhanced our scalar radiative transfer model MOMO, which is in use at
Radiative transfer Freie Universit
at Berlin, to treat these effects as accurately as possible. We describe our
Bio-optical ocean model
one-dimensional vector radiative transfer model for an atmosphere ocean system with
Matrix operator method
a rough interface. We describe the matrix operator scheme and the bio-optical model
Ocean optics
Polarization for case one waters. We discuss some effects of neglecting polarization in radiative
MERIS transfer calculations and effects of salinity changes for top of atmosphere radiances.
OLCI Results are shown for the channels of the satellite instruments MERIS and OLCI from
412.5 nm to 900 nm.
& 2012 Elsevier Ltd. All rights reserved.

1. Introduction and Zhai [8] to name a few. The work described in this
paper is based on the radiative transfer model MOMO
An accurate and exible remote sensing scheme has a which is itself based on the work of Fischer and Grassl [9],
broad range of possible applications in the eld of atmo- Fell and Fischer [5] and Bennartz and Fischer [10]. It has a
spheric and oceanic research. Virtually, all analyses of long tradition of successfully developed remote sensing
measurements made by radiance sensors need radiative applications, including the sensing of lakes [11], analysis
transfer (RT) calculation results to derive meaningful of hyper spectral data to derive surface uorescence
physical quantities. In this paper we describe a radiative signals [12], the analysis of ocean color data from MERIS
transfer scheme which is able to calculate the vector measurements [13], and the retrieval of land surface
radiance eld in an atmosphere ocean system (AOS) with pressure from MERIS data [14]. We decided to upgrade
a wind blown interface. We assume that the system has the MOMO FORTRAN code to account for polarization in
no horizontal, but arbitrary vertical structure. Hence, the order to base the development of future remote sensing
scheme is a one-dimensional vector radiative transfer algorithms on accurate RT calculations.
solver. Similar systems have been described in the past In Sections 2 and 3 we introduce the radiative transfer
and recent literature, such as the works from Kattawar equation and the matrix operator method. In Sections 4
and Adams [1], Nakajima and Tanaka [2], Takashima [3], and 5 we describe the models for pure ocean water and
Chami [4], Fell and Fischer [5], Chowdhary et al. [6], He [7] the bio-optical model for in water constituents. Section 6
is devoted to the validation of the code and in Section 7
n
Corresponding author. Tel.: 49 30 838 56656. we describe rst applications as mentioned in the
E-mail address: andre.hollstein@fu-berlin.de (A. Hollstein). abstract.

0022-4073/$ - see front matter & 2012 Elsevier Ltd. All rights reserved.
doi:10.1016/j.jqsrt.2012.01.010
 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548 537

2. Radiative transfer equation multiplying with the Gauss Lobatto weights ci:
Z X
k
The differential radiative transfer equation (RTE) given dm f m  f mi ci : 8
in Eq. (1) states that the change of the diffuse light eld i1

@t Lt with respect to the optical thickness t is propor- We dene matrices that contain the Gaussian points,
tional to both the light eld itself, and the diffuse sources weights, phase matrix values, and source term values:
Jt at this optical depth:
c diagc1 , . . . ,ck , 9
m@t Lt Lt Jt: 1
M diagm1 , . . . , mk , 10
The light eld is described by a real four-dimensional
Stokes vector [15,16] (and references therein). To nd 7, 7
Pm i,j P m 7 mi , 7 m0j , i,j 2 1, . . . ,k, 11
unique solutions, it is necessary to dene boundary
conditions that dene the top and the bottom of the
Jm
i,j o0 S0 P m i,j et=mi , 12
atmosphere. Eq. (2) states that there is no diffuse down-
ward directed radiation at the top of the atmosphere, and J 
m i,j o0 S0 P m i,j e
t=mi , 13
Eq. (3) states that the upward directed radiation at the
bottom of the AOS is given by the reection of the where d0m is the Kronecker delta. Dening the matrices
downward directed radiation. The surface reection is Gm == = and S = as abbreviations:
0
modeled using a real 4  4 reection matrix Rm, f, m0 , f Gm M1 1o0 p1 d0m Pm c, 14
0
which depends on the direction of incidence (m , f ) and
0

reection (m, f Gm  M1 o0 p1 d0m P m  c, 15

Lt 0, m o0 0, 2
G
m M
1
o0 p1 d0m Pm c, 16
Z 1 Z
Lt t0 , m 40, f dm0
0 0 0
df Rm, f, m0 , f Lt0 , m0 , f : G
m M
1
1o0 p1 d0m P 
m c, 17
0
3 Sm7 M1 Jm7 18
The complexity of the RTE comes from the coupling of we can insert them into Eq. (5) and write the result as a
the eld by the scattering source term J, which is shown compact matrix equation:
in Eq. (4). It consists of a scattering term for the direct !
! ! !
d L Gm Gm 
Lm Sm
solar radiation and a scattering term for the diffuse eld    : 19
  dt L Gm Gm Lm Sm
d
m 1 Lt, m, f o0 Pt, m, f, ms , fs et=ms S0
dt |{z}
single scattering term
Z 3. Matrix operator method
0 0 0
o0 dm0 df Pt, m, f, m0 , f Lt, m0 , f , 4
|{z} The method is based on the interaction principle
which has been described by Twomey et al. [17] and later
diffuse scattering term
by Grand [18]. It includes any order of scattering and is
where S0 is the solar constant, the solar position is set to applicable to systems with any optical thickness.
ms , fs , and o0 is the single scattering albedo. We assume The interaction principle states that the upward directed
0
that the scattering matrix P only depends on ff and light eld at a given optical thickness depends linearly on
expand L and P in a Fourier series with the expansion the transmitted light eld from a layer at higher optical
coefcient m. The equation then decouples into a series of thickness, and the downward directed intensity at the same
equations in Fourier space that are now independent of level. The interaction coefcients are called reection rij and
the viewing azimuth angle: transmission tij and a schematic is shown in Fig. 1. This
  holds analogously for the downward directed light eld at
d
m 1 Lm t, m o0 Pm m, ms et=ms the lower level:
dt
Z L t2 t 21 L t1 r 12 L t2 J 21 , 20
o0 p1 d0m dm0 P m m, m0 Lm t, m0 : 5
L t1 r 21 L t1 t 12 L t2 J 12 : 21
We discretize Eq. (5) for numerical treatment on a Stating the interaction principle for two consecutive
computer system, and split the light eld into parts for atmospheric layers with three boundaries, one can elim-
the upper and lower hemisphere: inate the transmission and reection operators of the

Lm t, m Lm t, m 4 0, 6 intermediate layer. By writing the resulting equations in
the same form as the interaction principle, the transmis-
sion and reection operators of the combined layers can
m t, m Lm t, m o0:
L 7
be expressed as [17,19,20,1,5]:
Integrations are then replaced by summing over the
integrand at Gaussian quadrature points mi and t31 t32 1r12 r32 1 t21 , 22
538 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548


r01 Cm dt, 29

r10 C
m dt: 30

4. Sea water optical model

Scattering of radiation by ocean water can be modeled

using a Rayleigh-like phase matrix with a depolarization
factor of d 0:039 [22]. Sources of scattering in the
homogeneous water bulk are thermodynamic uctuations
of sea water density and salt ion concentration. A detailed
discussion in terms of thermodynamics has been given by
Fig. 1. Interaction principle. Zhang and Hu [23]. As an option for the user we keep the
previously used model from Morel [24] which neglects
temperature and salinity effects. Morel described a model
with exponential spectral dependency with two constants
that have been derived from measurements in open
oceans:
 4:32
1 l
bMorel
i 0:00288 n : 31
m 500 nm
The absorption coefcient aw of pure sea water is taken
from the algorithm theoretical basis document (ATBD) of
the ESA project WATERRADIANCE [22] which includes
data from many sources [2531]. The absorption is
Fig. 2. Doubling and adding scheme. modeled as linear expansion with coefcients CS for
salinity, CT for temperature, and absorption measure-
ments at T 0 201 and S0 0 PSU:
t13 t12 1r32 r12 1 t23 , 23 aw T,S, l aw T 0 ,S0 , l TT 0 CT l SS0 CS l: 32
The refractive index of air relative to sea water can be
r31 r21 t12 r32 1r12 r32 1 t21 , 24
calculated using:
2
r13 r23 t32 r12 1r32 r12 1 t23 : 25 nair
sw T,S, l 2 300 nm,800 nm n0 n1 n2 T n3 T S
n5 n6 S n7 T n8 n9
Applying the algorithm for two layers with the same n4 T 2 2 3 33
l l l
optical properties is known as doubling; applying it to
layers with different optical properties as adding. with values of the ki from Ref. [22]. For longer wavelength
Fig. 2 shows how this concept is applied to the atmo- up to 4000 nm the temperature and salinity dependency
sphere ocean system. The initial optical thickness t0 after at 800 nm can be used.
n doublings grows exponentially like tn to 2n and should The volume scattering coefcient of sea water bsw is
be chosen such that nal result is independent of to . The the sum of contributions from density uctuations (bdf )
layers can be combined using the adding algorithm and and concentration uctuations (bcf ) and have been dis-
the interaction principle can be used to calculate the cussed by Zhang and Hu [23]:
radiances at inter layer boundaries. bsw l,T,S bdf l,T,S bcf l,T,S: 34
The mathematical foundation of this procedure was
published in a series of papers by Grant and Hunt In Fig. 3 we show relative differences of the model by
[18,21] and a result was the differential matrix form of Morel and by Zhang and Hu. The variation with salinity
the interaction principle, which has the same form as and temperature is shown as a set of gray curves and at
Eq. (19): 400 nm the relative differences are in the order of 10%.
! !
! The variation of the relative refractive index with
d L d L d J 01 respect to wavelength, salinity and temperature is shown
t 01 1 r 10 r 01 t 10  1 :
dt L dt L dt J 
10 in Fig. 4. On the scales of temperature and salinity
26 relevant to the Earths oceans, the changes with salinity
are more pronounced than those with temperature.
A comparison of the two equations yields the de-
nition of the elementary transmission and reection
5. Bio-optical model for case one waters
operators:

t10 1Cm dt, 27 The bio-optical model relates the chlorophyll concen-
tration to its spectrum of inherent optical properties
t01 1C
m dt, 28 (IOP). The latest measurements of IOPs available to us
 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548 539

Fig. 3. Volume scattering coefcient with log scale according to the Morel model (in black) and varying with temperature and salinity (in gray).
The dashed line represents (right scale) relative differences between the model for a temperature of 201 and salinity of 20 PSU.

refractive index 1.0

1.355
T 15C, varying S chlorophyll
0.9
1.350 S 20PSU, varying T concentration

1.345 0.8
0.1
,T,S

1.340 0.7 0.5

0

T 10C
n

1.335 T 20C 1.
0.6
5.
1.330
S 35PSU
S 25PSU 0.5 10.
1.325 S 15PSU
S 5PSU Bricau d et al. 50.
0.4 from Mie calc.
400 600 800 1000 1200 1400
in nm 500 600 700 800 900
in n m
Fig. 4. Relative refractive index with respect to wavelength, salinity and
temperature. Variations due to salinity changes at constant temperature Fig. 5. Single scattering albedo spectra from Bricaud 1998 (dashed
are shown with black lines, those with temperature at constant salinity lines), and best t results with penalty terms using a log normal size
are shown with gray dashed lines. distribution. The gray vertical columns indicate the position and the
width of the OLCI channels in that spectral band.

are from Bricaud et al. (from 2010 via personal commu-

nication, see also [32]) and include measurements of the
single scattering albedo o0 and the chlorophyll absorp- perfect t to the data, but the phase matrices show
tion coefcient achl; but not the phase matrix. unphysical high oscillating phase functions. To overcome
Information about the shape of chlorophyll phase this we add a penalty value to the norm used in the
functions has been given by Morel [33], Chowdhary [6], optimization scheme if the backward direction of the
Chami [4] and others. Spectral constant size distribution resulting phase function is either too small, or shows too
parameters and the complex refractive index spectra are strong oscillations. This describes a tradeoff between
input parameters for Mie calculations. We derive the rst achieving physical phase functions and a better t of the
guess of the imaginary part of the refractive index by single scattering albedo spectrum. Fig. 5 shows the resulting
using the following equation [32,34]: single scattering albedo spectra and the measurements for
chlorophyll concentrations from 0:1 mg=l to 50 mg=l.
achl ll
ni l , 35 Fig. 6 shows the spectral dependence of all absorption
4p
coefcients in the bio-optical model. Data for the chlor-
where achl is the chlorophyll absorption coefcient. This ophyll absorption is present where the symbols are
simple model assumes that the electric eld in a chlor- drawn. For a given chlorophyll concentration in the
ophyll particle is an evanescent wave. To compensate the Bricaud 1989 data set, the absorption of CDOM can be
shortcomings of this simple model, we introduce an offset calculated using the proposed exponential model:
to the complex refractive index. We use the Levenberg
Marquardt method to optimize the size distribution 1
aCDOM l aCDOM 440 nmeSl440 nm , S 0:011
parameters and the offset to nd a t to the measured nm
single scattering albedo spectra. The result is an almost 36
540 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548

100 chlorophyll absorbtion pure water absorbtion

101
absorbtion in m1

102

103
C 0.1 C 0.5
C 1 C 5
yellow substance absorbtion
104
400 450 500 550 600 650 700
in nm

Fig. 6. Absorption coefcients of the constituents of the bio-optical

ocean model. The black solid line shows the absorption of pure water,
the dashed lines show the absorption of CDOM and the straight lines
show the absorption of chlorophyll. Different gray tones represent
different chlorophylls and associated CDOM concentrations.

6. Model validation

The standard way of validating new implementa-

tions of RTE solvers is to compare results with indepen-
dent implementations for common test cases. Tables Fig. 7. Difference of MOMO calculations and Natrajs tables [40]. All
containing numeric values of the polarized radiation available cases are plotted in gray. The zenith resolved mean is shown in
eld for special cases have been published by several black, and the overall mean as number on the right scale. (a) Relative
authors [3541]. We used data from Natraj et al. [40] differences of the published tables and MOMO results and (b) signicant
differences of the published tables and MOMO results. Signicant differ-
and Kokhanovsky et al. [41] to validate MOMO. First,
ences indicate the number of equal digits of the compared numbers.
we compare Stokes vector tables for cases of pure
Rayleigh scattering. Second, we compare light elds
calculated with the SCIATRAN model [42] for scattering surface. As scatterers they have chosen the Rayleigh, an
by three scattering matrices (including Rayleigh) over a aerosol and a cloud-type phase matrix. The components
black surface. At last we look at physical aspects of the of the phase matrices under consideration are shown in
light eld. Fig. 8. Both phase functions (see m1 and m3 component)
are strongly peaked in the forward direction. The other
6.1. Rayleigh scattering phase matrix elements m5 and m6 are rich of features and
hence useful to test the accuracy of the models for given
The tables for Rayleigh scattering by Natraj et al. [40] resolutions.
have been calculated using a numerically more stable The mean relative deviations of MOMO and SCIATRAN
approach to the solution of the X and Y functions of results for all cases are shown in Fig. 9. Results of MOMO
Chandrasekhar and Coulsen [43,35]. The range of viewing calculations for the two cases together with the SCIATRAN
angles in the tables is from 01 to 901 (16 values), solar results are shown in Appendix in Fig. B2. Our results agree
positions range from 01 to 901 (7 values), and azimuth very well with the results from SCIATRAN. The dashed
positions range from 01 to 1801 (7 values). Surface albedos lines in the gure represent cases in which the phase
range from 0.0 to 1.0 (3 values) and Rayleigh optical depth functions have been truncated and straight lines repre-
range from 0.02 to 1.0 (5 values). sent results for the original phase functions. This approx-
The main result of our comparison is shown in Fig. 7. imation cannot reproduce the intensity in the forward
Panel (a) shows the relative differences with respect to scattering direction, but reproduces the components Q, U,
viewing angle for every point in the table. Panel (b) is V in this region quite well and results can be reproduced
the same type of graph, but for signicant differences. The with much less computational effort (see also Appendix
signicant difference of two numbers a and b is the B.1 and Fig. B2).
difference of the signicant in language independent
arithmetic form of a and b divided by 10 to the power 6.3. Conservation of ux
of the exponent of a. The value shows the number of
(nonzero) equal digits of the two numbers. Elastic scattering conserves the radiance ux. To test
only the effects of scattering, we calculated a number of
6.2. Mie scattering test cases where the direct solar source is only present in
the top of atmosphere layer. This layer is the only source
Kokhanovsky et al. [41] have published results from of diffuse radiation which propagates through the system.
the SCIATRAN model where three scattering cases involve With no internal sources and no absorption present, the
a homogeneous purely scattering layer over a black outward directed uxes from a single layer must balance
 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548 541

resolved radiances at a layer boundary:

10 Z
m m
1000 Aerosol
Cloud
F m,k t dm df mLm,k t, m, f: 39
100
10
A nonzero value of d has two possible origins. First,
due to errors in the implementation, and second due to
1
the numeric approximation of the integrals. Since it can
0.1
be almost impossible to distinguish the two effects d must
0 50 100 150 be small. For layers with Mie scattering d is smaller than
 in 1% when choosing 30 atmospheric Gaussian points, 36
Fourier terms, and azimuthal output at every 301. Using
0.15 m m
this procedure we can also test the atmosphere ocean
0.10
Aerosol
Cloud
interface and at the same resolution we nd deciencies
0.05
smaller than 3%.
0.00
0.05 6.4. Ocean surface reection
0.10
0.15 Tables of the radiation emerging from the ocean, to our
0 50 100 150
best knowledge, have not been published yet.1 In a paper by
 in
Mobley et al. [44], the authors showed a comparison of
Fig. 8. Benchmark phase matrices from Kokhanovsky et al. [41]. (a) M1 several models for oceanic waters. The scalar version of our
and M3 components and (b) M5 and M6 components. model agreed well with the published values [5]. However,
they do not include polarization or detailed radiance elds,
and hence are not be sufcient for a comparison. For this
1.02
benchmark: rayleigh aerosol cloud comparison, we consider an atmosphere with Rayleigh
1.015 Truncated:
Original: scattering with an optical thickness of 0.1, a rough sea
1.01 surface as lower boundary, no ocean body, and the wind
I

1.005 speed is set to 7 m/s. Via personal communication, the

1.0 author of the NASA GISS radiative transfer model Jacek
I

Chowdhary shared with us some results of his model [6,45].

0.995
0 20 40 60 80 The compared model output was the upward directed
zenith angle  in vector radiance at the top of the atmosphere. In Fig. 10 we
show results for the upward directed radiance and degree of
1.01 polarization for both models. The dashed lines represent
Chowdharys result, and the plot markers show the model
1.005
resolution of 0.1 in terms of zenith cosines. Since for a given
1.0 zenith resolution both of the models use different sampling
points we used MOMO with a high resolution of 60 zenith
S

0.995
angles to be able to compare with Chowdharys results. In
S

0.99 both panels we show results for different solar angles

0 20 40 60 80
indicated by different gray shades, and show solar positions
 in
with vertical dashed lines. The two models agree very well
Fig. 9. Mean of relative differences for MOMO and SCIATRAN results. with each other which holds for all other viewing geome-
The mean is taken over the azimuth values and up and downward tries. A comparison in which we used the same zenith
direction. Scatter types are indicated with gray tones. Straight lines
resolution produced similar results. The largest deviations
show results for the original phase matrices and dashed lines for the
truncated phase matrices. Absolute values are shown in Fig. B2. are seen in the direction of the horizon near sun glint for
(a) Intensity and (b) Q, U and V parameter. low solar angles.

6.5. Known properties

the inward directed uxes from the neighbor layers:
E F m tn F k tn F m tn1 F k tn 1 , 37 To test if the radiative transfer model is physically
|{z} |{z}
F out F in correct we look at special cases and verify if the model
reproduces our expectations. The sun glints directed into
F out F in the atmosphere and into the ocean are good proxies for
d 1
: 38
2F out F in the correct implementation of the reection and trans-
mission matrices of the interface. Both are dominant
The deciency value E should be small for all layers,
features, and show strong dependencies with zenith and
and can serve as a proxy for the accuracy of the model for
azimuth angle, and they are computed with the same
a given resolution. Here we use d to describe the ratio of
the deciency to the mean of the two uxes. To test the
Fourier expansion, we calculate the uxes from azimuthal 1
A statement also given by Zhai et al. in 2010 [8].
542 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548

78.5 66.4 53.1 36.9 0. opposite solar position

0.50 with sea surface
10 Rayleigh only case

0.20

I TOA
0.10 w 0.4m s
Intensity I

w 1.0m s
10
0.05 w 3.0m s
w 5.0m s

0.02 w 7.0m s

10 0.01
50 0 50
50 0 50
zenith angle  in deg
zenith angle  in deg

Fig. 11. Upward directed radiance at the top of the atmosphere in the
solar position 0. 36.9 53.1 66.4 78.5
1.0 principle plane. Solar position is indicated with a vertical line. The black
line with plot markers represents the result with a black surface and also
shows the zenith resolution. Result for different surface wind speeds is
0.8 shown as gray shaded lines. The solar constant was set to unity, so the
unit represents reection.
degree of polarization dop

0.6
Critical Angle
1.0
opposite solar position diffraction
0.4 for solar
0.8 position
upward directed radiance
0.2 just above sea surface
0.6

downward directed
an d I

0.4
0.0
radiance just
50 0 50 below sea surface
zenith angle  in deg 0.2
I

Fig. 10. Comparison of radiances and degree of polarization computed 0.0

with our and Chowdharys model. The viewing geometry is the principal
plane, and curves show the results for different solar positions. 0.2
(a) Comparison of radiances in the principal plane for different solar
positions and (b) comparison of degree of polarization in the principal 0.0 0.2 0.4 0.6 0.8 1.0
plane for different solar positions. zenith Angle 

matrices that are used for the diffuse reection and Fig. 12. Upward (black) and downward (gray) directed radiances just
above and below the ocean surface. The solar position and refracted
transmission. In Fig. 11 we show upward directed TOA solar position are shown with vertical dashed lines. The critical angle is
radiances for a Rayleigh atmosphere with a black surface, shown by the gray vertical line.
and cases with an ocean surface with varying wind speed.
The Rayleigh optical thickness for these cases was set to of 0.1 m/s we can clearly see, that the transmitted diffuse
0.1 (  545 nm) using a US standard atmosphere. As the radiance is almost entirely refracted into the Fresnel cone.
wind speed becomes smaller, the sun glint becomes more At the boundary of the cone we can see the brightening of
narrow and stronger as one would expect if the surface the underwater horizon. If the surface wind speed
roughness is changing from rippled to at. increases the boundary becomes more smooth and radi-
When the light eld propagates into the ocean it is ance is distributed to the outside of the Fresnel cone.
refracted according to Snells law. In Fig. 12 (t 0:1,w For a at ocean surface the Fresnel cone has a sharp
1 m=s) we show the upward and downward directed boundary and all transmitted radiance from the atmo-
radiance just above and below the ocean interface. The sphere is refracted into the cone. Upward directed radi-
directions are indicated by the signs of the radiances. We ance just below the ocean interface becomes total
show the position of the sun, the critical angle, and the internally reected if outside the Fresnel cone. In Fig. 14
position of the refracted solar position with vertical lines. As we show the downward directed radiance just below the
expected, the glint into the ocean is refracted, and the sun ocean interface; the sun is in zenith, the sky is black, and
glint has a steeper angle. the ocean is purely scattering. Different wind speeds are
An other way to test the refraction in the model is to shown with different shades of gray and one can see how
look at the downward directed radiance just below the the sharp boundary of the Fresnel cone becomes smother
ocean surface. We set the solar position to the zenith and for rougher ocean surfaces.
chose a high water absorption so that there is no upward
directed radiance that can contribute to the signal from 7. First applications
reection at the ocean surface. In Fig. 13 we show the
results of this simulation for different surface wind MOMO now accounts for polarization, sea water salinity,
speeds. When looking at the lowest chosen wind speed and sea water temperature and all three effects can have a
 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548 543

transmitted intensity much faster to compute and less difcult to implement,

0.08
but resulting relative deviations from vector calculations
can reach up to 20% for atmosphere ocean systems [8] and
0.06 will in general depend on wavelength and scene. For a
brief discussion we computed the effects of neglecting
polarization for the 412 nm OLCI channel where the sea
0.04 water is most transparent. Fig. 15 shows hemispheric
I

plots with radiances calculated with the vector model

w 0.1m s
w 1.0m s (see top of each panel) and the relative difference 1~ I=I
w 3.0m s
0.02
w 5.0m s
with respect to the scalar mode (bottom of each panel).
w 7.0m s We call 1~ I=I the relative polarization error and if it is
smaller than zero, the scalar calculations are
0.00
0.0 0.2 0.4 0.6 0.8 1.0
zenith Angle 

Fig. 13. Downward directed radiance just below the ocean surface. The
solar position is set to zenith and different surface wind speed cases are
shown with different gray shades. The ocean is highly absorbing and
reection from upward directed radiance from the ocean is not con-
tributing to the signal.

Fresnel cone

0.35

0.30 w 0.1m s
w 0.2m s
0.25 w 0.4m s
w 1.0m s
0.20 w 3.0m s
w 5.0m s
I

0.15
w 7.0m s
0.10

0.05

0.00
0.0 0.2 0.4 0.6 0.8 1.0
zenith Angle 

Fig. 14. Downward directed radiance just below the ocean surface
without contribution from transmitted radiance. Different surface wind
speeds are indicated with different shades of gray.

non-negligible impact on top of atmosphere and water

leaving radiance. In this section we discuss two applications
of the model for top of atmosphere and water leaving
radiances, for channels from 412.5 nm up to 900 nm of
the upcoming ocean color instrument OLCI (Ocean Land
Colour Instrument) onboard ESAs Sentinel 3 satellite.

7.1. Effects of neglecting polarization in case one waters

MOMO can be operated in a scalar mode which

neglects polarization, and therefore consumes less mem-
ory and is much faster than the vector mode. When
switching from polarization to scalar mode, the difference
in computing time and memory is slightly smaller than
Fig. 15. Top of atmosphere vector radiance eld in the 412.5 nm channel
1/(3  3) for linear, or 1/(4  4) for complete, Stokes vector
for an atmosphere ocean system with clear atmosphere and case one
calculations. The main reason for this effect is that higher waters. The chlorophyll concentration is 0:1 mg=l, the wind speed is 7 m/s,
order Fourier terms of the scattering matrices tend to and the salinity is 35 PSU. In panel (a) the sun is at 25.91 and in panel (b)
contain more zero-valued elements than those of lower the sun is at 50.31. The upper part of each panel shows the vector radiance
orders (depending on the actual physical constituent) and where the radiance direction on the unit sphere has been projected to the
equatorial plane. The gray colors indicate the value of the radiance in
that we use a sparse matrix multiplication approach
W/m2/nm/sr. The lower part of the panel uses the same projection
which is much faster for matrices containing zero-valued technique but shows the value of 1~ I =I to show the relative difference
elements. Using only scalar radiative transfer is usually of the two radiances in %. (a) Sun at 25.91 and (b) sun at 50.31.
544 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548

underestimating the real value. Results are shown for the pattern, therefore we see no easy way of correcting scalar
solar positions 25.91 (see panel (a)) and 50.31 (see panel radiative results without running a full vector model.
(b)). The salinity is 35 PSU, the wind speed is 7 m/s, and
the sea water temperature is 15 1C. For both cases
combined, the relative polarization error is in the order 7.2. Salinity
of 78%. A change of sign occurs and the pattern of this
change depends strongly on solar position. The highest The salinity of the ocean can vary from values as little
relative polarization error can be found in the principle as 5 PSU, as in the northern part of the Baltic sea, to a
plane opposite to the solar position (Fig. 16). maximum of 40 PSU. A standard value of 35 PSU for the
Fig. 15 shows results for water leaving radiances for open oceans can be assumed [46]. Figs. 17 and 18 show
the same case. The range of the relative polarization error the effects of a salinity change from 35 PSU to 5 PSU on
increases from 2.5% to  5% to values from 1.5% to  6% top of atmosphere and water leaving zenith radiance:
for the higher solar angle. With increasing solar angle a 1I5 PSU =I35 PSU . The water leaving radiance is dened as
sign change of the relative polarization error occurs. the upward directed radiance just above the ocean sur-
We can conclude that the effects of neglecting polariza- face, but without radiance contributions from diffuse and
tion for a realistic atmosphere ocean system with case one direct reection from the ocean surface. It is the radiance
waters depend strongly on viewing geometry and solar directly emerging from the ocean and is an useful starting
angle, and in the shown cases reached values up to 78%. point for ocean colour retrievals using radiative transfer
The relative polarization error can show a rather complex results for the ocean alone.
Both gures show values for the sun at zenith (top value
per cell), at 25.91 (bottom value per cell), and a chlorophyll
concentration ranging from zero (pure sea water) to
50 mg=l. The atmosphere was modeled aerosol free, and
gaseous absorption was taken into account. The root causes
of the salinity effect are changes in the sea water absorption,
the sea water bulk scattering coefcient, and the real part of
the sea water refractive index. The results for the pure sea
water case can be seen as upper limit to the salinity effect.
Increasing chlorophyll concentration leads to additional

Fig. 17. Salinity effect (1I35 PSU =I5 PSU ) for top of atmosphere zenith
radiance in %. Top value per cell represents results with the sun in zenith
position, and bottom values for the sun at 25.91. Results are shown for
the OLCI channel subset from 412.5 nm to 900 nm, and a chlorophyll
concentration range from pure sea water to 50 mg=l.

Fig. 18. Salinity effect (1I35 PSU =I5 PSU ) for zenith water leaving radi-
ance in %. Top value per cell represents results with the sun in zenith
position, and bottom values for the sun at 25.91. Results are shown for
Fig. 16. Same case as in Fig. 15, but for the water leaving radiance. the OLCI channel subset from 412.5 nm to 900 nm, and a chlorophyll
(a) Sun at 25.91 and (b) sun at 50.31. concentration range from pure sea water to 50 mg=l.
 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548 545

absorption and scattering by ocean components which are Appendix A. Rough atmosphere ocean interface
independent from the salinity, therefore decreasing the
salinity effect. A chlorophyll concentration of 0:1 mg=l The model of the surface reection matrix R and trans-
represents global mean value and for this concentration mission matrix T is based on papers by Nakajima and
the salinity effect is in the order of 3.114.85%. For top of Tanaka [2], Kattawar and Adams [1], and Sancer [47]. This
atmosphere radiances and channels up to 753.8 nm, the is modeled as a Gaussian distribution of surface facet
salinity effect is generally decreasing with increasing chlor- normals which is caused by surface roughening due to the
ophyll concentration, and remains almost constant for the wind blowing over the surface [48,49]. Shadowing effects of
channels with higher wavelength. This is caused by the the surfaces facets are modeled using a shadowing function
increase in sea water absorption (see Fig. 6). High sea water Z 1 Z 2p
absorption renders the water body almost black and the 1
Em 1 dm0 dfRm, m0 , f Tm, m0 , f: A:1
remaining cause for effects are changes of the refractive m 0 0

index, and therefore the reectivity of the atmosphere ocean

The deviation of the deciency term Em from unity
interface. When this effect is dominating the signal the
shows how the models conserve radiance ux. In Fig. A1
salinity effect becomes independent from the chlorophyll
we show the hemispheric reection and transmission for
concentration.
wind speeds of 5, 10, and 15 m/s for relative refractive
Fig. 18 shows the salinity effect for the water leaving
indexes 1.33 and 1.33  1. The top panel shows the values
radiance. For a chlorophyll concentration of 0:1 mg=l the
for light incident from the atmosphere and the bottom
salinity effect is ranging from 4.85% to 16.84%. The effect
panel for incident light from the ocean. The window of
is highest when the sea water absorption is lowest. For
total internal reection (Snells cone) can be clearly seen
pure sea water the effect is almost constant with variation
in panel (b) of Fig. A1. The deciency is larger for larger
from 17.59% to 19.26%, and with increasing chlorophyll
wind speeds which may be caused from multiple scatter-
concentration the salinity effect becomes smaller. For the
ing and shadowing effects which are modeled by a
channels in the NIR, due to the high sea water absorption,
shadowing function.
the chlorophyll signal dominates the water leaving radi-
ance, therefore the salinity effect becomes small.
We showed that the sea water salinity can have a 1.00
signicant effect on both, the top of atmosphere radiance,
and the water leaving radiance. The effect is caused by 0.50
changes of the sea surface reection, sea water absorption
and scattering. The changes for water leaving radiance
can reach up to 19.26%, and could have a signicant effect 0.20

for retrievals of ocean constituents if the actual value of

0.10
the salinity is neglected.

0.05

8. Conclusion and outlook

0.02
MOMO, now accounts for the polarization of radiation
0 20 40 60 80
in the atmosphere ocean system, and we are condent
angle of incidence in
that our implementation is free of major errors. During
the development of the vector version of the program, the
latest stable scalar version was used to implement the 1.00

treatment of inelastic Raman scattering. The result of

merging the two versions would be an almost feature 0.50

complete, one-dimensional, radiative transfer system for

atmosphere ocean systems. 0.20

0.10

Acknowledgments
0.05

We thank the DFG (Deutsche Forschungsgemeinschaft)

with the SPP 1294 (HALO) and the ESA (European Space 0.02

Agency) WATERRADIANCE project for providing funding for 0 20 40 60 80

this research. We want to thank Jacek Chowdhary and Peng- angle of incidence in
Wang Zhai for sharing the results of their radiative transfer
programs with us, and also thank for the fruitful discussions Fig. A1. Hemispheric reectivity and transmissivity for wind speeds 5,
10, and 15 m/s for a relative refractive index of 1.34 (from atmosphere)
with Rene Preusker during the implementation of this and 1.34  1 (from ocean). The deciency E is shown gray dashed lines for
model. We also want to thank the two anonymous reviewers all three wind speeds. (a) Radiation incident from atmosphere and
for their very much appreciated comments. (b) radiation incident from ocean.
546 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548

Appendix B. Numerical techniques B.2. Fourier series modication for the conservation
of radiation
B.1. Phase function truncation
The conservation of intensity for scattering can become
Phase function truncation is a widely used method to an issue if the zenith resolution of the RT computations is too
decrease the number of Fourier terms necessary for low. Incident radiation ~ Sm0 is scattered0 to all other direc-
m
azimuthal radiative transfer calculations. The general tions m described by the phase matrix: ~ S m Mm0 , m~
Sm0 .
procedure is shown in Fig. B1. The strong forward peak ~
In Fourier space, the zeroth matrix M 0 describes the mean of
of the original phase function (shown in black) is con- the angular distribution and we therefore describe the
tinuously replaced by a second order polynomial. To energy conservation in Fourier space as:
accurately represent the modied phase function far Z Z
fewer Fourier terms are necessary, but forward scattered m 0
I m0 dm ~
S 0 m ~ 0 m0 , m~
dmM Sm0 0 : B:1
radiation is effectively treated as unscattered. We show
this effect in Fig. B2 using the phase functions used in
the comparison with SCIATRAN. We show results from In the case of unpolarized incident solar radiation the
the SCIATRAN model, and results from MOMO using the matrix product simplies, and the conservation of inten-
original and the truncated phase function. In panel (b) we sity can be written as:
show the downward directed radiation. The strong diffuse Z
forward scattering could not be reproduced, but away 8m0 -1 ~ 11 m0 , m 0-Em0 :
dm M B:2
0
from this feature the results with truncated phase func-
tion agree well.
Due to the limited number of zenith angles, the left
To keep the polarization properties of the original
hand side of Eq. (B.2) may not vanish and is set to E.
phase matrix, we rescale all other elements of the matrix
Increasing the number of zenith angles would diminish
with the ratio of the truncated and the original phase
this residual but this may not be possible due to con-
function [4].
strains in available computation resources. For this reason
we modify the phase matrix:
10
1000
~ 11 m0 , m0 M 11 m0 , m0  Em0
M B:3
0 0
g m0
Phase Function

100 original phase function

10
As described in Appendix B.1, the phase matrix is then
1
0.1
truncated one modied to keep the polarization state of the scattered

0 50 100 150
Scatering Angle
10.00 up down
Fig. B1. Phase function and truncated phase function. Truncation is 5.00 Original
Modified
performed from 151 on. 1.00 Highres
0.50
I

0.10
0.05
3.00
azimuth angle: 0 90 180 0.01
2.00 Sciatran result: 0 20 40 60 80
1.50 Original phase matrix:
Truncated phase matrix: zenith angle  in deg
1.00
I

0.70
0.50
up down
0.30 0.20 Original
Modified
0 20 40 60 80 0.15 Highres
d op

zenith angle  in 0.10

0.05

0.00
5.0 0 20 40 60 80

2.0 zenith angle  in deg

I

1.0 Fig. B3. Effect of diagonal balancing of residuals to better conserve the
0.5 total intensity of the system. Comparing intensity and degree of
polarization for the original and balanced case using a cloud model
0.2
with optical thickness 5 and a black underlying surface. Azimuthal
0 20 40 60 80 averaged calculations where done using 180 (high res) and 60 zenith
zenith angle  in angles. For the cases with lower resolution the original and the modied
cases where calculated and the high resolution is seen as truth.
Fig. B2. Detailed comparisons of MOMO and SCIATRAN results. (a) Upward (a) Comparison of upward and downward azimuthal averaged intensity
directed radiance for the aerosol case and (b) downward directed radiance and (b) comparison of upward and downward azimuthal averaged
for the aerosol case. degree of polarization.
 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548 547

radiation unaltered, and the rescalation factor d becomes: retrieval from ocean color in case i waters. J Geophys Res 2003;108.
doi:10.1029/2002JC001638.
Em [14] Lindstrot R, Preusker R, Fischer J. The retrieval of land surface
dm 1 : B:4 pressure from MERIS measurements in the oxygen a band.
g m M 11 m, m0
J Atmos Oceanic Technol 2009;26(7):136777, doi:10.1175/2009
JTECHA1212.1.
The effect of these phase function modications for a [15] Tyo JS. Design of optimal polarimeters: maximization of signal-to-
noise ratio and minimization of systematic error. Appl Opt
cloud case is shown in Fig. B3. The highres case was 2002;41(4):61930, doi:10.1364/AO.41.000619.
calculated using 180 zenith Gaussian points and the [16] Aiello JPWA. Linear algebra for Mueller calculus, arXiv:math-ph/
unaltered phase matrix. The original and modied were 0412061v3.
[17] Twomey S, Jacobowitz H, Howell HB. Matrix methods for multiple-
calculated using 60 zenith Gaussian points but the origi- scattering problems. J Atmos Sci 1966;23:28998.
nal and the modied phase matrix. The surface is black [18] Grant IP, Hunt GE. Discrete space theory of radiative transfer. I.
and the optical thickness of the cloud is 5. The highres case fundamentals. Proc R Soc London, Ser A: Math Phys Sci
1969;313(1513):18397.
represents the truth. The modied version agrees very
[19] Grant IP, Hunt GE. Solution of radiative transfer problems using the
well with the highres results but needed signicantly less invariant sn method. Mon Not R Astron Soc 1968;141:2741.
computation effort. The usage of the original phase matrix [20] de Hulst V. A new look at multiple scattering, NASA Institute for
leads to errors in the intensity and degree of polarization. Space Studies; 1963.
[21] Grant IP, Hunt GE. Discrete space theory of radiative transfer. II.
The case was chosen so that the effect of the discussed stability and non-negativity. Proc R Soc London, Ser A: Math Phys
procedure becomes visible; for standard cases the effect Sci 1969;313:199216.
would not be as pronounced since real droplet size [22]
Rottgers
R, Doerfer R, McKee D, Schonfeld W. Pure water spectral
absorbtion, scattering, and real part of refractive index model. ESA
distributions tend to dilute the pronounced features of algorithm technical basis document; 2010.
such special cases. [23] Zhang X, Hu L. Estimating scattering of pure water from density
uctuation of the refractive index. Opt Express 2009;17(3):16718.
[24] Morel A, Prieur L. Analysis of variations in ocean color. Limnol
References Oceanogr 1977;22(4):70922.
[25] Pope RM, Fry ES. Absorption spectrum (380700 nm) of pure water. II.
[1] Kattawar GW, Adams CN. Stokes vector calculations of the submarine Integrating cavity measurements. Appl Opt 1997;36(33):871023.
light eld in an atmosphereocean with scattering according to a [26] Lu Z. Optical absorbtion of pure water in the blue and ultraviolet.
Rayleigh phase matrix: effect of interface refractive index on radiance PhD thesis, Texas A&M University; 2006.
and polarization. Limnol Oceanogr 1989;34(8):145372. [27] Wang L. Measuring optical absorption coefcient of pure water in
[2] Nakajima T, Tanaka M. Effect of wind-generated waves on the uv using the integrating cavity absorption meter. PhD thesis, Texan
transfer of solar radiation in the atmosphereocean system. J Quant A&M University; 2008.
Spectrosc Radiat Transfer 1983;29(6):52137, doi:10.1016/0022- [28] Segelstein DJ. The complex refractive index of water. PhD thesis,
4073(83)90129-2. Department of Physics, University of Missouri Kansas City; 1981.
[3] Takashima T. Polarization effect on radiative transfer in planetary [29] Wieliczka DM, Weng S, Querry MR. Wedge shaped cell for highly
composite atmospheres with interacting interface. Earth Moon absorbent liquids: infrared optical constants of water. Appl Opt
Planets 1985;33:5997. 1989;28(9):17149.
[4] Chami M, Santer R, Dilligeard E. Radiative transfer model for the [30] Max J-J, Chapados C. Isotope effects in liquid water by infrared
computation of radiance and polarization in an oceanatmosphere spectroscopy. III. h2 o and d2 o spectra from 6000 to 0 cm  1. J Chem
system: polarization properties of suspended matter for remote Phys 2009;131(18):184505, doi:10.1063/1.3258646.
sensing. Appl Opt 2001;40(15):2398416, doi:10.1364/AO.40. [31]
Rottgers R. Technical note, unpublished data; 2010.
002398. [32] Bricaud A, Morel A, Prieur L. Optical efciency factors of some
[5] Fell F, Fischer J. Numerical simulation of the light eld in the phytoplankters. Limnol Oceanogr 1983;28(5):81632.
atmosphereocean system using the matrix-operator method. [33] Morel A, Antoine D, Gentili B. Bidirectional reectance of oceanic
J Quant Spectrosc Radiat Transfer 2001;69(3):35188, doi:10.1016/ waters: accounting for raman emission and varying particle scat-
S0022-4073(00)00089-3. tering phase function. Appl Opt 2002;41(30):6289306, doi:10.1364/
[6] Chowdhary J, Cairns B, Travis LD. Contribution of water-leaving AO.41.006289.
radiances to multiangle, multispectral polarimetric observations [34] Morel A, Bricaud A. Theoretical results concerning light absorption
over the open ocean: bio-optical model results for case 1 waters. in a discrete medium, and application to specic absorption
Appl Opt 2006;45(22):554267, doi:10.1364/AO.45.005542. of phytoplankton. DeepSea Res Part A. Oceanogr Res Pap
[7] He X, Bai Y, Zhu Q, Gong F. A vector radiative transfer model of 1981;28(11):137593, doi:10.1016/0198-0149(81)90039-X.
coupled oceanatmosphere system using matrix-operator method [35] Coulson K, Dave J, Sckera Z. Tables related to radiation emerging
for rough sea-surface. J Quant Spectrosc Radiat Transfer 2010;111(10): from a planetary atmosphere with Rayleigh scattering. Berkeley:
142648, doi:10.1016/j.jqsrt.2010.02.014. University of California Press; 1960.
[8] Zhai P-W, Hu Y, Chowdhary J, Trepte CR, Lucker PL, Josset DB. [36] Garcia RDM, Siewert CE. A generalized spherical harmonics solu-
A vector radiative transfer model for coupled atmosphere and tion for radiative transfer models that include polarization effects.
ocean systems with a rough interface. J Quant Spectrosc Radiat J Quant Spectrosc Radiat Transfer 1986;36(5):40123, doi:10.1016/
Transfer 2010;111(78):102540. 0022-4073(86)90097-X.
[9] Fischer J, Grassl H. Radiative transfer in an atmosphereocean [37] Garcia RDM, Siewert CE. The FN method for radiative transfer
system: an azimuthally dependent matrix-operator approach. models that include polarization effects. J Quant Spectrosc Radiat
Appl Opt 1984;23(7):10329, doi:10.1364/AO.23.001032. Transfer 1989;41(2):11745, doi:10.1016/0022-4073(89)90133-7.
[10] Bennartz R, Fischer J. A modied k-distribution approach applied to [38] Evans KF, Stephens GL. A new polarized atmospheric radiative
narrow band water vapour and oxygen absorption estimates in the transfer model. J Quant Spectrosc Radiat Transfer 1991;46(5):
near infrared. J Quant Spectrosc Radiat Transfer 2000;66(6):53953. 41323, doi:10.1016/0022-4073(91)90043-P.
[11] Heege T, Fischer J. Mapping of water constituents in lake constance [39] Wauben W, Hovenier J. Polarized radiation of an atmosphere
using multispectral airborne scanner data and a physically based containing randomly-oriented spheroids. J Quant Spectrosc
processing scheme. Can J Remote Sensing 2004;30(1):7786, doi:10. Radiat Transfer 1992;47(6):491504, doi:10.1016/0022-4073(92)
5589/m03-056. 90108-G.
[12] Guanter L, Alonso L, Gomez-Chova L, Meroni M, Preusker R, Fischer [40] Natraj V, Li K-F, Yung YL. Rayleigh scattering in planetary atmo-
J et al. Developments for vegetation uorescence retrieval from spheres: corrected tables through accurate computation of x and y
spaceborne high-resolution spectrometry in the o2-a and o2-b functions. Astrophys J 2009;691(2):1909.
absorption bands, J Geo Res 2010;115, doi:10.1029/2009JD01371. [41] Kokhanovsky AA, Budak VP, Cornet C, Duan M, Emde C, Katsev IL et al.
[13] Zhang T, Fell F, Liu Z-S, Preusker R, Fischer J, He M-X. Evaluating the Benchmark results in vector atmospheric radiative transfer. J Quant
performance of articial neural network techniques for pigment Spectrosc Radiat Transfer 2010. doi:10.1016/j.jqsrt.2010.03.005.
548 A. Hollstein, J. Fischer / Journal of Quantitative Spectroscopy & Radiative Transfer 113 (2012) 536548

[42] Rozanov A, Rozanov V, Buchwitz M, Kokhanovsky A, Burrows J. [46] Antonov JI, Seidov D, Boyer TP, Locarnini RA, Mishonov AV, Garcia
Sciatran 2.0a new radiative transfer model for geophysical HE et al. World Ocean Atlas 2009, vol. 2: Salinity, NOAA Atlas
applications in the 1752400 nm spectral region. Adv Space Res NESDIS 68; 2009.
2005;36(5):10159. [47] Sancer M. Shadow-corrected electromagnetic scattering from a
[43] Chandrasekhar S. Radiative transfer. Courier Dover Publications; randomly rough surface. IEEE Trans Antennas Propag 1969;17:
1960. 57785.
[44] Mobley CD, Gentili B, Gordon HR, Jin Z, Kattawar GW, Morel A, et al. [48] Cox C, Munk W. Measurement of the roughness of the sea surface
Comparison of numerical models for computing underwater from photographs of the suns glitter. J Opt Soc Am 1954;44(11):
light elds. Appl Opt 1993;32(36):7484504, doi:10.1364/AO.32. 83850.
007484. [49] Ebuchi N, Kizu S. Probability distribution of surface wave slope
[45] Chowdhary J. Multiple scattering of polarized light in atmosphere derived using sun glitter images from geostationary meteorological
ocean systems: application to sensitivity analyses of aerosol satellite and surface vector winds from scatterometers. J Oceanogr
polarimetry. PhD thesis, Columbia University; 2000. 2002;58:47786, doi:10.1023/A:1021213331788.