Abstract

Generation of random sea surfaces using wave variance spectra and Fourier transforms is formulated in a way that guarantees conservation of wave energy and fully resolves wave height and slope variances. Monte Carlo polarized ray tracing, which accounts for multiple scattering between light rays and wave facets, is used to compute effective Mueller matrices for reflection and transmission of air- or water-incident polarized radiance. Irradiance reflectances computed using a Rayleigh sky radiance distribution, sea surfaces generated with Cox–Munk statistics, and unpolarized ray tracing differ by 10%–18% compared with values computed using elevation- and slope-resolving surfaces and polarized ray tracing. Radiance reflectance factors, as used to estimate water-leaving radiance from measured upwelling and sky radiances, are shown to depend on sky polarization, and improved values are given.

© 2015 Optical Society of America

1. Introduction

Optical reflection and transmission by water surfaces are fundamental to radiative transfer calculations for oceanic and atmospheric applications. Reflection of Sun and sky radiance by the sea surface creates glint. In ocean color remote sensing of water-column conditions, glint is an undesired contribution to the measured upwelling radiance that, along with atmospheric path radiance, must be estimated and removed from the measured radiance in order to obtain the water-leaving radiance, which carries the information about the water column itself. In remote sensing of the surface wave state or wind speed and direction, glint is the signal because it depends on surface roughness. The water-leaving radiance itself comes from upwelling underwater radiance, which is transmitted through the sea surface from water to air. The underwater light field arises from light that is transmitted through the surface from air to water and is affected by underwater radiance that is reflected back downward by the underside of the sea surface.

Accurate calculation of surface optical properties has three requirements. First, the sea surface itself must be modeled in a way that accounts for wave elevation and slope. To first order, reflection and transmission of a light ray incident onto the sea surface are governed by the slope of the surface at the point of intersection. However, wave elevation can lead to wave shadowing for incident angles near the horizon, which refers to rays predominately striking the sides of the waves tilted toward the light source. In the ocean, wave elevation is dominated by the highest amplitude waves, which are generally gravity waves with wavelengths of tens to hundreds of meters. The waves with the largest slopes are often gravity waves with wavelengths of a meter or less and capillary waves with wavelengths of a few millimeters. Thus, the simulated sea surface must include the elevation and slope effects for spatial scales that vary by a factor of 105.

Second, the simulated light rays must account for polarization. Surface reflection and transmission induce polarization, even if the incident source (e.g., the Sun’s direct beam) is unpolarized. Incident sky radiance is usually partially linearly polarized. Thus, both reflected and transmitted light fields are, in general, partially polarized. Although many sensors are designed to be insensitive to polarization, the state of polarization of the incident light affects the total energy reflected or transmitted by the sea surface. Polarization also carries information that potentially can be used to improve algorithms for removal of surface glint, underwater visibility, or separation of particle types, if polarization-sensitive sensors are employed.

Third, the ray tracing must account for multiple scattering between light rays and waves. Although single scattering dominates ray-surface interactions, multiple scattering can occur for incident angles that are large relative to the local normal to the sea surface. This occurs most often when the incident ray is nearly horizontal, e.g., when the Sun is near the horizon.

Simulations of sea surfaces for optical calculations often employ Cox–Munk wind speed-wave slope statistics for generation of random surface realizations. Their equations [1]:

σa2=0.003+0.00316W12.5±0.002,
σc2=0.00192W12.5±0.004,
σ2=0.003+0.00512W12.5±0.004,
relate the along wind (σa2), cross wind (σc2), and total (σ2=σa2+σc2) variances of the nondimensional sea surface slopes to the wind speed W12.5 in meters per second at 12.5 m above sea level (the coefficients of W12.5 have units of sm1). Surfaces generated using these equations capture slope effects for the full range of wave spatial scales at the times of observation. The data upon which these equations rest were obtained for wind speeds of 1–14ms1. The sea states were not described but were probably mature but not fully developed. Preisendorfer and Mobley [2] used Cox–Munk surfaces and a ray tracing algorithm that included all orders of multiple scattering to study sea surface optical properties, but only for unpolarized light. Although the Cox–Munk model provides a reasonable simulation of surface slopes, it cannot reproduce wave elevations.

Both elevation and slope statistics can be captured if wave variance spectra and Fourier transforms are used for surface generation. These surfaces are referred to here as FFT surfaces because the fast Fourier transform (FFT) is used to evaluate the discrete Fourier transforms. Several previous studies [37] have used FFT surfaces but have not included polarization or multiple scattering. Tulldahl and Steinvall [8] used a two-step, two-scale model to simulate the effects of both gravity and capillary waves on wave slopes. In the first step of their calculations, they used a coarse spatial grid to resolve the gravity waves. In the second step, they added capillary-scale waves to the coarse grid. Their capillary-scale waves were generated using slope variances equal to the difference of the Cox–Munk slope variance and the slope variance of the gravity waves generated in the first step. Their surfaces, thus, accounted for slope effects due to the full range of wavelengths, and their ray tracing used the same multiple scattering algorithm as [2]. They did not, however, include polarization. McLean [5] accounted for unresolved slope variance in a study of lidar transmission through the surface by adding extra divergence to the lidar beam when it passed through the sea surface. A later similar study [6] used a grid resolution of 0.5 cm to capture the slopes of capillary-scale waves, but the longest gravity waves were not fully resolved. The recent study by Kay et al. [7] simulated the full range of spatial scales from 200 m to 3 mm needed to resolve the amplitude and slope variances, but their brute-force computations required 65536×65536 grid points in the FFTs, with consequently very long run times of 6 h to generate a single surface realization. They also did not include polarization.

A number of studies have modeled the polarization properties of sea surfaces but only for level [913] or Cox–Munk surfaces [1417]. The latter studies sometimes included an analytic function to model wave shadowing, but none accounted for multiple scattering between wave facets.

The previous studies just cited have, thus, included at best two of the three abovementioned requirements for accurate computation of sea surface optical properties. The first goal of the present paper is to introduce techniques for computing optical properties for sea surfaces that incorporate elevation and slope effects for the full range of gravity and capillary spatial scales, with all orders of multiple scattering between light rays and surface waves, and including polarization via the full Stokes vector formalism. The second goal is to compare predictions made with the complete formalism with those made using the simplifications of Cox–Munk surfaces, unpolarized ray tracing, and/or single scattering. Fortran 95 code was written to implement the surface-generation and ray-tracing algorithms described in the next two sections. That code can generate either Cox–Munk or FFT surfaces, and the ray tracing can be either polarized or unpolarized. The ray tracing keeps a separate tally for singly and multiply scattered rays. It is, thus, easy to compare results obtained for different sea surface models, for polarized versus unpolarized treatments of the light, and for single versus multiple scattering.

The next section describes the surface wave generation, including an algorithm to account for elevation and slope variances without the need for large numbers of grid points. The polarized ray tracing algorithm is outlined in the following section. Section 4 then shows examples of reflectance and transmittance for both water- and air-incident light for the simple case of an unpolarized, collimated incident beam. Reflection and transmission of a partially polarized Rayleigh sky radiance distribution are then illustrated in Section 5. Section 6 then compares exact calculations (i.e., using FFT surfaces and polarized ray tracing with multiple scattering) of sea surface irradiance reflectance with results obtained using various approximations. The use of radiance reflectance factors for surface glint removal is well established [18]. Finally, Section 7 compares these factors as computed using FFT surfaces and polarized ray tracing with the factors obtained with Cox–Munk surfaces and unpolarized ray tracing.

2. Sea Surface Generation

The initial part of this section reviews the relations between sea surfaces and the wave elevation and slope variance spectra needed for sea surface generation. The next part then shows how a wave variance spectrum is used to define random sea surfaces that are consistent with the chosen wave spectrum. Particular attention is paid to the normalization factors that are required for conservation of wave elevation variance (proportional to wave energy) in a round-trip calculation from a wave variance spectrum to a sea surface realization and back to a wave spectrum. The final section shows how a wave elevation variance spectrum can be rescaled so that a finite range of sampling of the spectrum in frequency space generates sea surfaces that reproduce the elevation and slope statistics of the real sea surface corresponding to the original variance spectrum.

A. Wave Variance Spectra

The Monte Carlo ray tracing begins with the generation of random realizations of sea surfaces. Let z(x)=z(x,y) be the sea surface elevation in meters at point x=(x,y) at a particular time. The spatial extent of the sea surface is 0x<Lx and 0y<Ly, with Lx and Ly in meters. A wind-centered coordinate system is used, with the +x direction chosen to be downwind; x is then upwind, and ±y are the cross-wind directions. This surface is sampled on a rectangular grid of Nx by Ny points, where both Nx and Ny are powers of 2 to enable use of FFTs. The spatial sampling points are then x(r)=rΔx,r=0,,Nx1 with Δx=Lx/Nx. A similar equation holds for the y(s)=sΔy points. A discrete sample of z(x,y) is denoted z(xrs) or z(r,s).

The two-dimensional (2D) forward discrete Fourier transform of z(r,s), which is evaluated using the FFT algorithm, converts the surface elevations to a 2D grid of complex Fourier amplitudes z^(u,v)=FFT2D{z(r,s)}. These amplitudes are defined at positive and negative discrete angular spatial frequencies kx(u)=uΔkx,u=(Nx/2+1),,Nx/2, with Δkx=2π/Lx. The fundamental frequency of magnitude Δkx corresponds to the longest resolvable wave. The positive frequency kx(Nx/2)=2π/(2Δx) is the Nyquist frequency, which corresponds to the shortest resolvable wavelength of 2Δx. Corresponding results hold for the ky(v) spatial frequencies. These spatial frequencies have units of radians per meter. A discrete frequency pair kuv=[kx(u),ky(v)] is labeled by (u,v). The sea surface z(xrs) is real valued, so the amplitudes are Hermitian: z^*(kuv)=z^(+kuv), where z^* denotes the complex conjugate of z^.

The absolute value squared of the discrete amplitudes is the discrete elevation variance spectrum corresponding to the surface sample z(r,s). This spectrum is denoted Ψ(u,v)=|z^(u,v)|2. Dividing the discrete spectrum by the sampling bandwidth gives a discrete estimate (a periodogram) of the continuous elevation variance spectral density: Ψ(kx,ky)=Ψ(u,v)/(ΔkxΔky). Averaging the periodograms from many independent samples of the sea surface averages out the statistical noise inherent in each individual periodogram and leads to the elevation spectral density of the sea surface. It is important to note that the discrete z2(r,s), |z^(u,v)|2, and Ψ(u,v) are all point functions with units of m2, whereas the continuous function Ψ(kx,ky) is a spectral density with units of m2/(rad/m)2. The integral of Ψ(kx,ky) over all spatial frequencies gives the variance of the zero-mean sea surface elevation:

var{z}=z2=Ψ(kx,ky)dkxdky,
where indicates expectation or ensemble average over many measurements of the sea surface. The energy contained in a surface wave is proportional to the wave amplitude squared, so variance spectra are often loosely called energy spectra, and conservation of wave energy corresponds to conservation of elevation variance.

The FFT of z(r,s) leads to a discrete two-sided (2S) spectrum, denoted here by Ψ2S(kuv). “Two-sided” means that positive and negative frequencies, which correspond to waves propagating in opposite directions, are present. With the choice of +x being downwind, frequencies with a positive kx value correspond to waves propagating generally downwind (|φ|<π/2, where φ=tan(ky/kx)); negative kx corresponds to waves propagating upwind (|φ|>π/2).

The resulting Ψ2S(k)=Ψ2S(kx,ky) is equal for directions +k and k. This symmetry arises because the FFT cannot determine which direction a wave is propagating because there is no time dependence in the surface sample z(r,s). In many applications in physical oceanography, the quantity of interest is just the total variance (i.e., total energy) per unit frequency interval contained in waves of a given frequency magnitude k=|k|=(kx2+ky2)1/2, without regard for the direction of wave propagation. In this case, a one-sided (1S) spectral density is used, which has only positive k values but is twice the magnitude of the two-sided spectrum: Ψ1S(k)=2Ψ2S(k). One-sided spectra are also known as folded spectra [19]. Wave spectral densities as presented in the oceanographic literature are almost always one-sided spectra. However, when generating sea surfaces as described below, two-sided spectra are needed, in which case the one-sided spectrum must be divided by two to convert it back to the proper magnitude of a two-sided spectrum. In reality, Ψ2S(k)Ψ2S(k) because much more energy propagates downwind than upwind at a given frequency, but symmetric spectra are consistent with the generation of time-independent surfaces as needed for the present study.

This study uses the one-sided elevation variance spectral density of Elfouhaily et al. [20], which they write in polar coordinates (k,φ) as

Ψ1S(k,φ)=1kS(k)Φ(k,φ).
Here, S(k) is the 1D omnidirectional spectral density with units of m2/(rad/m), and Φ(k,φ) is a nondimensional spreading function. Ψ1S(k,φ), thus, has the same units as Ψ1S(kx,ky). Their spreading function has the form:
Φ(k,φ)=12π[1+Δ(k)cos(2φ)].
This function is symmetric about φ=±π/2, i.e., the function gives symmetric spreading for downwind and upwind propagation. The division of Ψ1S by 2 when generating time-independent surfaces is, in essence, a division of the one-sided or folded spreading function to give equal energy propagation in the φ and φ+π directions.

The Ψ spectra give the variance of the sea surface elevations. The corresponding variance spectrum of the sea surface slopes is k2Ψ. Thus, the variance of the wave slopes (the mean square slope) in the along-wind (±x) direction is given by either of

σa2=kx2Ψ2S(kx,ky)dkxdky=0ππk2cos2φΨ2S(k,φ)kdkdφ.
A similar equation gives the cross-wind slope variance. The total slope variance is then
σ2=σa2+σc2=0k2S(k)dk.
These relations are derived in [21] and [22, Appendix A], both of which contain a general discussion of wave spectra.

B. Surface Generation Via FFTs

The next task is to generate a set of random discrete complex Fourier amplitudes that are physically consistent with the chosen wave elevation variance spectrum. These amplitudes must be defined for positive and negative frequencies, and they must be Hermitian, so that the inverse Fourier transform generates a real sea surface. In essence, the elevation variance spectrum is used to define the wave amplitudes as a function of spatial frequency, and random variables are used to define a random phase for each wave.

Following the formulation of Tessendorf [23], first define

z^o(kuv)12[ρ(kuv)+iσ(kuv)]×[Ψ1S(k=kuv)2ΔkxΔky]1/2,
=12[ρ(kuv)+iσ(kuv)]Ψ2S(kuv).
Here, ρ(kuv) and σ(kuv) are independent normal random variables with zero mean and unit variance, denoted N(0,1). A different random variable is drawn for each kuv value. Ψ1S(k=kuv) is a one-sided continuous spectral density evaluated at the discrete frequencies kuv. Note in Eq. (4) that the one-sided spectrum is divided by 2 in the process of generating the symmetric discrete two-sided spectrum of Eq. (5). This division by two reflects the assumption that the one-sided spectrum corresponds to a symmetric two-sided spectrum. Multiplication of the continuous spectral density by the sampling frequency bandwidths converts the spectral density to a discrete point function showing how much variance is contained in a finite frequency interval ΔkxΔky at each frequency kuv. Ψ2S(kuv) denotes the resulting two-sided, discrete variance spectrum.

The amplitudes z^o(kuv)=z^o(u,v) are random variables. The statistical expectation of |z^o(kuv)|2 is

z^o(u,v)z^o*(u,v)={12[ρ(u,v)+iσ(u,v)]Ψ2S(u,v)}×{12[ρ(u,v)iσ(u,v)]Ψ2S(u,v)}=Ψ2S(u,v)2[ρ2(u,v)+σ2(u,v)]=Ψ2S(u,v),
where it has been noted that ρ2(u,v)=σ2(u,v)=1 for independent N(0,1) random variables. Thus, z^o(u,v) is consistent with the chosen variance spectrum. However, z^o(u,v) is not Hermitian, so the inverse FFT would not give a real sea surface. However, the spectral amplitudes defined by
z^(kuv,t)12[z^o(+kuv)exp(iωuvt)+z^o*(kuv)exp(+iωuvt)]
are clearly Hermitian, so the inverse FFT applied to z^(kuv,t) will give a real-valued z(xrs,t). Time dependence is included in this definition for complete generality. Just as was done for |z^o(kuv)|2 above, it is easy to show that |z^(kuv,t)|2=12[Ψ2S(kuv)+Ψ2S(kuv)]. Thus, the amplitudes of Eq. (6) give a real sea surface that is consistent with the chosen energy spectrum. This result holds even if the two-sided spectrum is not symmetric for ±k. The time independence of the expectation shows that even though the shape of the wave surface depends on time, the total energy of the wave field does not.

If a sequence of independent surfaces is to be generated, as needed for the Monte Carlo simulations below, then the time can be set to 0 in Eq. (6). Each surface is then generated by drawing a new set of random numbers in Eq. (5). If the time-dependent propagation of a given surface is to be simulated, e.g., for creating a visual rendering of moving waves, then a single set of random numbers is drawn, and the time dependence is obtained from the exponentials of Eq. (6) and an appropriate dispersion relation to obtain the temporal angular frequency ωuv from the spatial angular frequency kuv. For deep-water gravity waves, the dispersion relation is ωuv2=gkuv, where g is the acceleration due to gravity. Note, however, that, when generating a time-dependent sequence of surfaces for a propagating wave field, it is necessary to use an asymmetric spreading function so that most of the energy propagates downwind. A commonly used spreading function for this purpose is the cosine-2s family [24] of the form Φ(k,φ)=C(s)cos2s(φ/2), where s is a spreading parameter that depends on k, and C is a normalization coefficient. These functions are asymmetric about φ=π/2 and give much stronger propagation downwind than upwind. In this case, essentially all of the variance is contained in the downwind half of Ψ2S(k), and no division by two is needed in Eq. (4). In any case, the random realization of the sea surface is obtained from the inverse FFT of the amplitudes: z(r,s,t)=FFT2D1{z^(kuv,t)}.

Equations similar to Eqs. (4) and (6) are widely used in the movie industry for computer-generated imagery of sea surfaces. However, the literature on those applications (e.g., [23]) always ignores the division of the one-sided variance spectrum by 2, as seen in Eq. (4), and the leading factor of 1/2, as seen in Eq. (6), is likewise missing. Multiplication by the frequency interval, as seen in Eq. (4), is also often ignored. The justification seems to be that scale factors do not matter because the generated surfaces subsequently will be distorted for artistic purposes, e.g., to make the waves appear larger. However, for scientific applications, it is crucial that the generated sea surfaces be physically correct. When generating time-independent surfaces with a symmetric spreading function such as that of Eq. (3), if Ψ1S is not divided by 2 in Eq. (4) and the 1/2 factor is omitted in Eq. (6), then the amplitudes will be a factor of 2 too large. This corresponds to a factor of four violation of conservation of variance (i.e., energy) in a round trip calculation from variance spectrum to sea surface and back to variance spectrum.

C. Accounting for Unsampled Variance

Generation of a sea surface requires specific values for Lx, Ly, Nx, and Ny in the above equations. This raises the question of how large a spatial region and how many points must be used in generating sea surface realizations. Not surprisingly, the answer depends on the application. The visual impression of a sea surface is, to first order, determined by the height of the waves, which, in turn, is governed by the largest gravity waves for the given wind speed. When creating sea surfaces for rendering into a visually appealing movie scene, Lx and Ly can be chosen large enough to resolve the gravity waves with the greatest elevation variance (estimated from the peak of the S(k) spectrum). The number of grid points can be small; 29=512 or 210=1024 grid points in each direction are usually adequate.

Accurate scientific calculations of sea surface reflectance require many more grid points in order to resolve the surface down to the millimeter scale where capillary waves make a significant contribution to the surface slope variance. It is now computationally possible to create 2D FFT surfaces with sufficiently large N values in the x and y directions, but run times are prohibitive. Kay et al. [7] created 200m×200m surfaces with 216=65536 points in each dimension. This allowed sampling of the variance spectrum for wavelengths from 200 m gravity waves to 3 mm capillary waves. However, it took 6 h to create just one surface on a 3 GHz computer. Many ray tracing applications require tens to hundreds of thousands of independent surface realizations to obtain satisfactory statistical estimates. Thus, it is necessary to finesse certain calculations so that large-N FFTs can be avoided.

The upper-left panel of Fig. 1 shows the omnidirectional elevation spectrum S(k) of Elfouhaily et al. [20] for a wind speed of 10ms1. The upper-right panel shows the corresponding slope spectrum k2S(k). Note that the slope spectrum falls off much more slowly with increasing k than does the elevation spectrum. The elevation and slope variances contained in these spectra can be computed by numerically integrating the spectra from some very low frequency k0 to a very high frequency kh, which cover the frequency range where the spectra are non-negligible in magnitude:

z2=k0khS(k)dk,
and
σ2=k0khk2S(k)dk.
Using k0=0.01 and kh=104 gives z2=0.4296m2 and σ2=0.06011rad2 for the spectrum of the present example.

 figure: Fig. 1.

Fig. 1. Example sampling of elevation and slope spectra. The red dots on the lines in the upper panels show the bounds of the elevation and slope spectra sampled using 1024 points with L=200m. The light red lines are the true spectra S, and the truncated heavy blue lines are the adjusted spectra S˜ obtained from Eqs. (9 )–(11). The dashed lines at 370rad/m indicate the conventional boundary between gravity and capillary waves. The lower panels show random realizations of short segments along the x direction of the sea surface elevation and slope for the true (light red line) and adjusted (heavy blue line) spectra.

Download Full Size | PPT Slide | PDF

Now suppose this spectrum is used to generate a surface with Lx=200m using Nx=1024 sampling points. The fundamental spatial frequency corresponding to the longest resolvable wavelength is then kf=2π/Lx=0.0314radm1 and the Nyquist frequency is kNy=16.085radm1; these frequencies are shown by the two red dots in the upper panels of Fig. 1. These two frequencies and the Nx2 evenly spaced frequencies in between are the frequencies at which the elevation variance spectrum is sampled. Using kf and kNy as the lower and upper limits in Eqs. (7) and (8) gives z2(N=1024)=0.4219m2 and σ2(N=1024)=0.02584rad2. Thus, the finite range of the sampled frequencies includes the fraction,

fEz2(N)z2=0.42190.4296=0.982,
of the total variance of the sea surface elevation. However, the corresponding fraction of the sampled slope variance is just fS=0.02584/0.06011=0.430. Thus, N=1024 is sampling 98% of the elevation variance but only 43% of the slope variance. This sampling is acceptable for creating a sea surface that looks realistic to the eye. However, the unsampled frequencies greater than the Nyquist frequency account for a large part of the optically important slope variance.

One way to account for the unsampled slope variance is to alias the variance of the waves with frequencies greater than the Nyquist frequency into the higher-frequency waves with frequencies less than the Nyquist frequency. The higher-frequency waves that are sampled will then contain too much variance, i.e., they will have amplitudes that are too large for their wavelengths, which will increase their slopes. This is done by adjusting or rescaling the elevation spectrum,

S˜(k)[1+δ(k)]S(k),
such that the integral of k2S˜(k) over the sampled region kf to kNy equals the integral of the true k2S(k) over the entire spectral range. Sampling the rescaled spectrum k2S˜(k) over the kf to kNy frequency range then gives the same slope variance, as would be obtained by sampling the true spectrum k2S(k) over the entire range of frequencies.

L is chosen sufficiently large for the low frequencies to be well sampled starting at the fundamental frequency kf. There is, thus, no need to modify the low-frequency part of the variance spectrum, which, if done, would adversely affect the total elevation variance. Only the high-frequency part of the spectrum needs modification. A simple approach is to take δ(k) to be a linear function of k between the spectrum peak kp and the highest sampled frequency kNy, and zero elsewhere:

δ(k){0ifkkpδNy(kkpkNykp)ifk>kp,
δNy is a parameter that depends on the spectrum (i.e., the wind speed), the size L of the spatial domain, and the number of samples N. δNy must be determined so that δ(k) “adds back in” the unresolved slope variance. The added variance will be zero at the peak frequency kp and largest at the Nyquist frequency. That is, the δ(k) function will make only a small change to the variance spectrum at the low frequencies, and the change will be largest near the highest sampled frequencies, which is consistent with the idea that the high-frequency waves have the largest slopes.

The parameter δNy is determined by noting that

σ2k0khk2S(k)dk=k0kNyk2S(k)dk+kNykhk2S(k)dkk0kNyk2S˜(k)dk=k0kNyk2S(k)dk+δNykpkNyk2(kkpkNykp)S(k)dk.

The right-most terms of these equations give (after recalling that δ(k)=0 for kkp)

δNy=kNykhk2S(k)dkkpkNyk2(kkpkNykp)S(k)dk.
These integrals are numerically evaluated for the given spectrum.

The heavy blue lines in the upper panels of Fig. 1 show the S˜(k) and k2S˜(k) spectra for the present example. It is clear that the δ(k) function has added progressively more variance to the higher frequencies. The rescaled variance spectrum does, of course, contain somewhat more variance over the sampled region than does the true spectrum. As the upper inset fE value shows, this rescaling has increased the fraction of sampled/true variance from 0.982 to 1.020. However, the upper fS number in the upper-right panel shows that 99.5% of the slope variance is now being sampled, as opposed to just 43% for the true spectrum. Having slightly too much total elevation variance is a reasonable trade-off for being able to model the optically important slope statistics without the need for very large numbers of grid points.

The lower-left panel of Fig. 1 shows a 10 m slice along x of a random realization of the surfaces generated from these two spectra (with the same set of random numbers). The surface elevations z(xr) from the true spectrum (shown as the light red line) and the rescaled spectrum (thick blue line) are almost indistinguishable at the scale of this figure. The lower-right panel shows the surface slopes computed from finite differences of the discrete surface heights. The Cox–Munk along-wind mean square slope of 0.0346 given by Eq. (1a) compares well with the value of 0.032 obtained with the rescaled spectrum. However, the value obtained from sampling the truncated true spectrum is only 0.022, or 64% of the Cox–Munk value. Thus, in this example, the δ(k) correction to the elevation variance spectrum reproduces the slope variance at the expense of a 2% error in the elevation variance, compared to the theoretical values that would be obtained if the spectra were sampled with enough points to fully resolve the elevation and slope spectra.

The rescaled spectrum, thus, reproduces the slope variance at the expense of some increase to the elevation variance. The question still remains as to how many FFT grid points should be used in practice. Figure 2 illustrates the approach to the full slope variance as the number of grid points increases when no adjustment is made to the wave spectra, versus the slope variances obtained with rescaled spectra. The filled red diamonds show the mean square slopes (mss) when an L×L=100×100m area of sea surface is simulated using an N×N FFT, for a wind speed of W=5ms1. The open red diamonds show the corresponding mean square slopes for the rescaled spectra. For N=27=128 points, the mss is only 0.0119, versus a value of mss=0.0328 for the rescaled spectrum. Not until N=215=32,768 does the mss for an unadjusted spectrum equal the value for the adjusted spectrum. For L=100m and N=215, the surface is being resolved down to a grid spacing of L/N=0.003m (or a Nyquist spatial frequency of kNy=1029radm1), which fully resolves the capillary waves. At this value, the rescaled and unscaled spectra are the same, i.e., δ=0 in Eq. (9). Larger values of N then make no difference in the mss values. The red dots show the elevation variances z2 for unadjusted (solid dots) and adjusted (open dots) spectra. For L=100m and W=5ms1, there is almost no difference in the elevation variances. The blue curves show the corresponding results for L=300m and W=10ms1. In that case, the unadjusted spectrum does not fully resolve the mss until N=217=131,072 points are used, which corresponds to a grid point spacing of 0.0023 m or a Nyquist frequency of 1373radm1. Now, however, the small N values give a substantial increase to the elevation variance, as seen by the open blue circles. For N=128, z2 for the adjusted spectrum is 1.54 times the value for the unadjusted spectrum. This difference decreases to an increase of 1.057 at N=1024 and 1.029 at N=2048. Numerical experimentation shows that N=1024 provides a reasonable balance between having enough grid points to minimize the increase in elevation variance and keeping the run times acceptably small. The dependence of surface-reflected radiance on the choice of N is discussed below in Fig. 10.

 figure: Fig. 2.

Fig. 2. Approach to fully resolved mean square slopes (diamonds, left ordinate) and elevations (circles, right ordinate) as a function of the number of grid points N, for winds speeds of W=5 (red) and 10ms1 (blue).

Download Full Size | PPT Slide | PDF

The surface generation techniques of this section give sea surfaces that correctly reproduce the elevation and slope statistics of a real sea surface corresponding to the chosen variance spectrum. Figure 3 shows an example of a 2D sea surface generated using these techniques. The along-wind, cross-wind, and total slope variances for this realization (computed by finite differences) are in good agreement with the respective Cox–Munk values of 0.035, 0.019, and 0.054. The significant wave height of H1/3=2.18m for this realization with Lx=100m is somewhat less than the theoretical value of H1/34z2=2.622m obtained by numerical integration of the elevation variance spectrum, as in Eq. (2). However, increasing Lx to 400 m accounts for the very longest waves and gives H1/3=2.613m. Table 1 compares the Cox–Munk slope variances of Eq. (1) with those obtained by random realizations of sea surfaces, as shown in Fig. 3, for various wind speeds. It should be noted that the along-wind and cross-wind mean square slopes (i.e., the ratio σa2/σc2) depend on choice of spreading function Φ(k,φ), but the total variance depends only on the omnidirectional spectrum S(k).

Tables Icon

Table 1. Mean Square Slopes as Given by Cox–Munk Formulas in Eq. (1) and by Random Realizations of FFT Surfacesa

 figure: Fig. 3.

Fig. 3. 2D example sea surface random realization for Lx=Ly=100m and Nx=Ny=512 and a wind speed of 10ms1 in the +x direction. Light color is high surface elevation; dark color is low elevation. Inset numbers give the significant wave height H1/3, mean square slopes (total, mss; along wind, mssx; and cross wind, mssy), and mean surface tilt angles from the vertical in the along-wind (θx) and cross-wind (θy) directions.

Download Full Size | PPT Slide | PDF

3. Polarized Ray Tracing

Ray tracing is most conveniently carried out if the surface is modeled as a grid of triangular wave facets. Resampling of rectangular FFT grids is described first. The process of ray tracing for polarized light is then outlined.

A. Resampling Rectangular Grids

The wave spectrum-FFT process gives sea surface elevations defined on a rectangular grid of points. The Monte Carlo computation of surface optical properties requires repeated determinations of the points of intersection of rays with surfaces and of the angles of incidence of the rays with the surface normals at the points of intersection. Those calculations are most easily performed if the sea surface is modeled as a grid of triangular, rather than rectangular, wave facets. The mathematics then amount to finding the intersection point of a line with a plane defined by the three vertices of a wave facet. The normal to the planar wave facet is easily computed from facet vertices.

There is no unique or one best way to map a rectangular FFT grid to a grid of triangles. However, the ray tracing algorithm that determines when a ray intersects the surface does depend on the particular triangulation scheme chosen because it is necessary to determine when a ray crosses the boundaries of the facet triangles. If the FFT surface is generated using Ny=Nx/2, then, essentially, every other x grid point can be omitted to generate a hexagonal grid of triangular wave facets as illustrated in Fig. 4 for the case of Nx=16 and Ny=8 (and Lx=Ly=8m). This resampling generates a hexagonal region with 3(Ny/2)(Ny/2+1)+1 facet vertices and 3Ny2/2 wave facets. Although a triangulation based on right triangles half the size of the isosceles triangles shown in Fig. 4 would use every FFT grid point, the present triangulation scheme was chosen because a well-proven computer code (used in HydroLight) for the ray tracing was already available. The Cox–Munk Eq. (1) can be used to generate surfaces of triangular wave facets with the chosen triangulation scheme, as described in [25, Section 4.3].

 figure: Fig. 4.

Fig. 4. Example mapping of a rectangular FFT grid (light blue lines) to a hexagonal grid of triangles as used in ray tracing. Red dots show the triangle vertices. FFT grid points at the right and top, shown by the dotted lines, are obtained by the inherent spatial periodicity of the surface as determined by FFT techniques. Thus, the surface elevation of point 35 is the same as that of point 27, point 57 is the same as point 1, etc. One of the triangles generated from the FFT grid is shaded in red.

Download Full Size | PPT Slide | PDF

B. Ray Tracing

Tracing the path of a light ray through 3D space, as it is reflected and refracted by the sea surface, is in principle no more complicated than determining where a straight line intersects a plane. Those calculations are purely geometric and do not depend on the polarization state of the light. However, the application of the basic geometric concepts to the irregular geometry of the wave facets comprising a random realization of the sea surface is rather tedious. The details of those geometric calculations are given in [26, Appendix B].

Polarized light is described by the four-component Stokes vector S̲=[I,Q,U,V]T, whose components describe the state of polarization in a particular coordinate system. The convention used here [27] is that positive Q represents vertical polarization (the electric field in the vertical, or meridian plane, normal to the mean sea surface) and negative Q represents horizontal polarization (electric field parallel to the sea surface). Positive U is +45deg polarization, and negative U is 45deg. Positive V is left circular polarization, and negative V is right circular. Incident rays approaching the sea surface before any interaction are described in a coordinate system based on the incident meridian plane, which is the plane containing the z axis (outward normal to the mean sea surface) and the direction ξi of the incident ray. The associated Stokes vector S̲i is specified relative to horizontal and vertical directions, hi and vi, which are, respectively, perpendicular and parallel to the incident meridian plane. Unit vectors hi and vi are determined so that “perpendicular” cross “parallel” is in the direction of light propagation. For an incident ray, this rule gives

hi=z×ξi|z×ξi|,vi=ξi×hi,ξi=hi×vi.

An incident ray always intersects a wave facet. The incident direction ξi and the wave facet normal n define the scattering plane containing the incident, reflected (ξr), and transmitted (ξt) rays. The incident Stokes vector must be “rotated” (transformed) from the incident meridian plane into a coordinate system whose axes are perpendicular (s) and parallel (p) to the scattering plane. Using the choices of incident direction cross facet normal equals the perpendicular vector, and perpendicular cross parallel equals direction of propagation, gives

s=ξi×n|ξi×n|,p=ξi×s,ξi=s×p.
Vector s resides in the plane of the wave facet and is normal to each of ξi, n, ξr, and ξt. Vector p resides in the scattering plane and is normal only to ξi.

The rotation angle 0αi<2π that rotates the incident Stokes vector into the scattering plane is the angle that rotates the initial perpendicular direction hi into the direction perpendicular to the scattering plane. (This same angle rotates vi into p.) The rotation angle can, therefore, be obtained from the dot product of hi and p. However, there is some ambiguity in the angle obtained this way. This is resolved by choosing a positive rotation to be counterclockwise when looking into the beam (in the ξi direction) [27]. In general, the rotation angle is the angle required to rotate the first vector v1 counterclockwise into the second vector v2 when looking into the rotation axis, which is vector v3=v1×v2. In the present case, v1=hi, v2=s, and v3=ξi is the rotation axis. In general, four cases must be considered, as shown in Fig. 5. The rotation angle depends on whether v1·v2 is positive or negative and on whether v1×v2 is parallel or antiparallel to the rotation axis. In the present case, v1×v2=hi×s is by construction parallel to ξi, in which case

αi=|cos1(v1·v2)|=|cos1(hi·s)|.
If v1×v2 is antiparallel to ξ, then
αi=2π|cos1(v1·v2)|=2π|cos1(hi·s)|.
Vector v1×v2 is parallel (antiparallel) to ξ if the z component of v1×v2 has the same (different) sign as the z component of ξ.

 figure: Fig. 5.

Fig. 5. Four cases for determining Stokes vector rotation angles. v1 is the initial direction, which is rotated into the final direction v2 by a counterclockwise rotation about direction ξ.

Download Full Size | PPT Slide | PDF

Once the rotation angle has been determined, it is applied via a rotation matrix. For the choice of a positive rotation being counterclockwise when looking into the beam, the Stokes vector rotation matrix is [27, page 25]

R̲(α)=[10000cos2αsin2α00sin2αcos2α00001].

Once a Stokes vector has been rotated into the local coordinate system of a wave facet, the 4×4 Fresnel reflection and transmission matrices are applied to obtain the reflected and transmitted (if any) Stokes vectors. There are four such Fresnel matrices, namely, for reflection from the air side of a wave facet (denoted by matrix R̲aw), transmission from air to water (T̲aw), reflection from the water side of a facet (R̲wa, which may be total internal reflection), and transmission from water to air (T̲wa). These Fresnel matrices are given in [26,28]. They depend only on the angle of incidence relative to the surface normal, θi=cos1|ξi·n|, and the index of refraction of the water.

When a ray eventually leaves the region of the sea surface, its final Stokes vector S̲f must be specified in the final meridian plane. The final horizontal and vertical directions are obtained from choosing the first direction to be the final ray direction ξf(=ξrorξt) cross z. Then, as above, vf=ξf×hf. The final rotation angle αf, which carries the Stokes vector from the most recent wave facet s and p axes to the final hf and vf axes, is obtained by Eqs. (12) or (13). It can happen, for example, that a downward reflected ray in air leaves the hexagonal grid region before intersecting the surface. Such a ray would intersect the surface at some point outside the hexagonal grid. Numerical simulation shows that, for the large spatial regions used in the present simulations, for which the initial rays are aimed toward the center of the grid, fewer than one ray per 100,000 anomalously exits the grid before intersecting the surface. Such anomalous rays are discarded for simplicity of the ray tracing algorithm, rather than being “wrapped around” the hexagon for further ray tracing.

A chain of one or more ray-surface interactions, thus, gives the final Stokes vector S̲f as

S̲f=R̲(αm,f)M̲(ψm)R̲(αm1,m)R̲(α1,2)M̲(ψ1)R̲(αi,1)S̲i,
S̲fE̲(ξiξf)S̲i,
where
  • αi,1 is the rotation angle that rotates the initial S̲i into the scattering plane of the first wave facet;
  • M̲(ψ1) is any of the four Fresnel reflectance and transmittance matrices R̲aw, R̲wa, T̲aw, or T̲wa applied to the first ray-wave interaction;
  • ψ1 represents the scattering angle for the incident and reflected or transmitted directions used in the Fresnel matrix;
  • R̲(α1,2) is the rotation matrix that takes the Stokes vector M̲(ψ1)R̲(αi)S̲i from the scattering plane of the first wave facet to the scattering plane of the second wave facet;
  • R̲(αm1,m) carries the current Stokes vector from wave facet m1 to the last facet m;
  • M̲(ψm) is the Fresnel reflection or transmission function for the last wave facet;
  • R̲(αm,f) is the rotation matrix that carries the Stokes vector from the scattering plane of the last wave facet to the final meridian plane associated with S̲f.
Note that only the incident and final Stokes vectors are referenced to a meridian plane.

Except for the case of total internal reflection of water-incident rays, each ray-surface interaction generates a reflected and a transmitted ray. The information for each (point of intersection, direction, coordinate system, and weight) is stored in a stack array. The most recently stored ray is then selected for further processing. In this way, processing continues until all rays in the stack have been traced until they have left the surface region without further interaction and their final values have been tallied.

The continuous set of all directions is discretized by partioning the directions into a set of rectangular (θ,ϕ) bins, plus two polar caps, which are collectively called quads. The quad partitioning used below has polar caps of 5 deg half angle and rectangular quads of Δθ×Δϕ=10×15deg, except for Δθ=5deg quads next to the equator of each hemisphere. The same partitioning is used in the unpolarized HydroLight radiative transfer code. Rather than perform the ray tracing de novo for particular incident Stokes vectors, quad-averaged transfer matrices E̲ are computed for all pairs of quads. There are four sets of these transfer matrices, e.g., E̲raw(QijQkl), which describe reflection for air-incident rays coming from quad Qij, corresponding to the ith θ band and jth θ band, and eventually being reflected upward after one or more ray-wave interactions into quad Qkl. For particular quads Qij and Qkl, E̲raw(QijQkl) is, thus, a 4×4 scattering (Mueller) matrix that shows how any quad-averaged incident Stokes vectors S̲i(Qij) are transformed by reflection into a final Stokes vector S̲f(Qkl). This matrix and its counterparts E̲taw, E̲rwa, and E̲twa are optical properties of the sea surface itself and do not depend on the incident Stokes vector. Thus, these matrices need be computed only once for a given wind speed (i.e., surface wave state) and water index of refraction.

The four sea-surface transfer matrices are statistically estimated as the ensemble averages of the transfer matrices computed for S independent surface realizations. The calculation for a particular ray is initialized by choosing a random direction within a quad, so that the initial ray is directed from the quad toward the middle of the realized surface. The initial E̲ is set to the 4×4 identity matrix. The random sequence of matrix multiplications seen in Eq. (14a) is then carried out. When a ray leaves the surface and crosses the hexagonal surface boundary without further interaction, the current E̲ array is tallied to the accumulating array for the pair of quads containing the incident and final ray directions. The averages of the E̲ arrays for the S surface realizations then yield the transfer matrices. Separate weight arrays are computed for singly and multiply scattered rays to enable quantification of the importance of multiple scattering.

The slope and elevation statistics of the FFT and Cox–Munk surfaces have elliptical symmetry about the along-wind direction. This means that rays need to be traced only for one quadrant of initial directions (e.g., for 0φ90deg). For the 10×15deg angular partioning (plus two polar caps) of the set of all directions used here, there are 434 quads in all. For incident rays from just one quadrant, there are 128 quads with an initial ray, for each surface realization. A run then initializes 128S rays. Multiple scattering and the proliferation of daughter rays created by reflection and transmission means that about three rays are traced to completion for each initial ray. Note that there is still Monte Carlo noise in the computed reflectances and transmittances due to tracing a finite number of rays, so these quantities are not perfectly elliptically symmetric. Numerical experimentation shows that using S=105 independent surface realizations, hence 12.8×106 initial rays, give Monte Carlo noise of only a few percent in the computed quantities for directions that are energetically important. Generating 105 surfaces and tracing roughly 36×106 rays to completion requires less than 8 h on a 2.5 GHz computer.

The E̲ matrices computed in this manner are energy-transfer matrices that give the fraction of energy in the incident quad that is reflected or transmitted into the final quad. These matrices can be converted to matrices for radiance reflectance and transmittance by incorporation of cosine and solid angle factors for the incident and final quads as shown in [25, Eq. (4.74)] and [26, Eq. (3.67)]:

R̲(QijQkl)=E̲(QijQkl)|μi|Ωij|μk|Ωkl.
Here, μ is the average of cosθ over the respective quad, and Ω is the solid angle of the quad. The notation R̲ denotes a matrix for reflection or transmission of radiance, rather than for energy. For coherent incident Stokes vectors (in the sense of [27]) with units of irradiance (Wm2), the final Stokes vector is obtained from S̲f=E̲S̲i. If the incident Stokes vector is a diffuse vector [27] with units of radiance (Wm2sr1), then the radiance version is used: S̲f=R̲S̲i. The radiance form is used if the transfer matrices are to be used as surface boundary conditions for solving the vector radiative transfer equation, e.g., as formulated in [26]. In either case, once the transfer matrices have been computed, any particular angular distribution of quad-averaged incident Stokes vectors can be transformed into reflected and transmitted Stokes vectors by simple matrix multiplications.

4. Sea Surface Reflectance and Transmittance

For an initial application of the above surface transfer functions, consider the easily understood case of an unpolarized, collimated source (e.g., the Sun in a black sky). Figures 69 show the (temporally or spatially averaged) polarized reflectance and transmittance for a fully developed sea surface and a wind speed of 10ms1. The surfaces were generated with FFTs using rescaled variance spectra to fully resolve both elevation and slope variances. The sea surface region was 200×200m with Nx=1024 and Ny=Nx/2=512, as required for mapping the rectangular FFT grid to a hexagonal grid of triangular wave facets. Ray tracing was performed for 105 surface realizations (about 7 h of computer time on a 2.4 GHz PC). The view is looking toward the point of reflection or transmission on the sea surface, i.e., looking at the glitter pattern, which is qualitatively shown by the size of the colored area. The source is, thus, located at the φ=180deg side the plotted hemisphere of directions, so the unscattered rays are traveling in the φ=0 azimuthal direction. The incident Stokes vector in each case is unpolarized light, S̲=[1,0,0,0]T. The incident direction is 50 deg from the zenith for air-incident light (Figs. 6 and 7) or 50 deg from the nadir for water incident light (Figs. 8 and 9). The upper-left panel of each figure shows the percent of the incident energy that is reflected or transmitted into each direction for quads of angular size 10×15deg in polar and azimuthal angles. Quad boundaries are shown in gray, and the specular quad is outlined. Quads receiving less than 0.001% of the total energy are not colored. The asymmetry in φ is due to Monte Carlo noise in quads receiving almost no energy. The other panels show the degree of polarization of the reflected or transmitted light.

 figure: Fig. 6.

Fig. 6. Reflected energy pattern (glitter pattern) for unpolarized air-incident light at a 50 deg incident angle from the zenith and a 10ms1 wind speed. The light source is in an otherwise black sky. The viewing direction is looking downward at the surface, facing the glitter pattern. The quad outlined in gray indicates the specular reflection quad that would receive all reflected light for a level surface. The tan color in the lower-right panel indicates values that are identically zero.

Download Full Size | PPT Slide | PDF

 figure: Fig. 7.

Fig. 7. Transmitted energy pattern for air-incident light corresponding to Fig. 6. The viewing direction is looking upward at the surface from within the water.

Download Full Size | PPT Slide | PDF

 figure: Fig. 8.

Fig. 8. Reflected energy pattern (underwater glitter pattern) for water-incident light at a 50 deg incident angle from the nadir and a 10ms1 wind speed. The light source is unpolarized in an otherwise black ocean.

Download Full Size | PPT Slide | PDF

 figure: Fig. 9.

Fig. 9. Transmitted energy pattern for water-incident light corresponding to Fig. 8. The viewing direction is looking downward at the surface from the air side.

Download Full Size | PPT Slide | PDF

For the air-incident case of Fig. 6, 3.87% of the incident energy is reflected, and the remaining 96.13% is transmitted into the water. The upper-right panel of Fig. 6 shows that the reflected light is as much as 97% horizontally polarized in the quads near θ=50deg, which is near the Brewster angle of 53 deg for a level surface. The ±45deg polarization of the lower-left panel arises from rays being reflected by wave facets tilted to the left or right relative to the incident direction. Figure 7 shows the air-incident energy that is transmitted through the surface into the water.

Figures 8 and 9 show the corresponding results for water-incident unpolarized light, i.e., light incident onto the under side of the sea surface. Figure 8 shows the “underwater glitter pattern” as seen looking upward at the bottom of the sea surface. The tracing of underwater rays, e.g., as can pass from one surface facet to another via an underwater path, assumes that the intervening water is nonabsorbing and nonscattering. For this incident angle and wind speed, 60.42% of the incident energy is reflected back into the water. For a level surface, all of the energy would be reflected because the incident angle of 50 deg is greater than the 48 deg critical angle for total internal reflection. For a windblown surface, some of the incident light intersects wave facets that are tilted, so that the incident angle with the surface is less than the angle of total internal reflection, and 39.58% of the upwelling energy leaks through the surface into the air. In the upper-right and lower-left panels of Fig. 8, the large areas of unpolarized light near the horizontal (tan color) come from total internal reflection, which cannot change unpolarized incident light into any state of polarization.

It was argued in the discussion of Fig. 2 that a value of Nx=1024 was a reasonable choice for numerical simulations. Figure 10 shows the reflected energy as a function of θ for φ=0 in Fig. 6 when values of Nx=512, 1024 and 2048 are used in the simulations and the elevation variance spectrum is rescaled as previously described. The Nx=512 curve is as much as 25% different from the Nx=1024 curve for the larger values in the central region of the glitter pattern. However, the difference in the 1024 and 2048 curves is less than 4%. Thus, Nx=1024 gives results within a few percent of the values obtained for Nx=2048 but with only one-fifth of the computer time (the ratios of (Nxlog2Nx)(Nylog2Ny) for the different choices of Nx and Ny=Nx/2). The red solid, dotted, and dashed curves in that figure are the values for Nx=1024 and three different seeds for random number generation. These curves differ by less than 4%, which shows that simulations with 105 surface realizations give statistical noise of only a few percent in the energetically important directions. Numerical experiments such as these led to the choices of Nx=1024 and 105 surface realizations for the simulations in the present work.

 figure: Fig. 10.

Fig. 10. Dependence of energy reflectance along the midline of the glitter pattern (φ=0) in Fig. 6, as a function of the number of grid points Nx used for surface generation.

Download Full Size | PPT Slide | PDF

Figure 10 also shows the reflectance for an FFT surface with Nx=1024 but without rescaling to account for unresolved slope variance. The surface is then too smooth, leading to too much reflected energy in directions near the center of the glitter pattern. Finally, the curve for a Cox–Munk surface shows that, for this wind speed, the Cox–Munk surface differs from the FFT surface by roughly ±25% in the central part of the glitter pattern. As seen in Table 1, the Cox–Munk and FFT surfaces have almost the same slope statistics for a 10ms1 wind speed. The Cox–Munk and FFT curves for the rescaled Nx=1024 or 2048 curves in Fig. 10 include multiple scattering between wave facets. The difference in these Cox–Munk and FFT curves is, thus, due to differences in resolution of the wave heights, which results in different amounts of multiple scattering and wave shadowing.

Figure 11 shows the percent of air-incident rays that undergo multiple interactions with surface wave facets. There is more multiple scattering for Cox–Munk surfaces for incident angles between 30 and 70 deg but more for the FFT surfaces at incident angles very near the horizon. The most common multiple-scatter event is one for which the incident ray is reflected into another wave facet, giving a final tally of one reflected and two refracted rays, although five rays were traced (including the intermediate one that traveled from one wave facet to another without being tallied to the final result). However, other combinations occur. For example, the incident ray may generate a reflected ray that leaves the surface without further interaction, and a refracted ray that intersects the surface again from the water side and undergoes total internal reflection. As with a single-scattering event, this multiple-scattering event gives a final tally of one reflected and one refracted ray, although four rays were traced to completion. For rays traveling under water from one facet to another, the water is assumed to be transparent. For large incident angles, rays occasionally have 10 or more interactions with the surface, but such events are rare—typically just one or two rays per 10,000 incident rays. The maximum number of multiply scattered rays is between 8% and 12% for Cox–Munk surfaces and 6% to 9% for FFT surfaces. The maximum occurs when the incident rays tend to be reflected into nearly horizontal directions, so that they are likely to intersect another wave. The drop-off for incident angles near the horizon occurs because those rays tend to intersect the sides of wave facets tilted toward the source such that the reflected ray heads upward, away from the surface.

 figure: Fig. 11.

Fig. 11. Percent of incident rays that undergo multiple interactions with the surface, as a function of wind speed and incident angle from the zenith, for FFT and Cox–Munk surfaces.

Download Full Size | PPT Slide | PDF

It is easy to show that an arbitrary sequence of Fresnel reflection and transmission matrices and rotation matrices, as seen in Eq. (14a), has the form:

E̲=[XXX0XXX0XXX0000X],
where X denotes a nonzero value. [A specific example is given in Eq. (16) below.] Thus, reflection and transmission cannot convert air-incident unpolarized or linearly polarized light into circular polarization. This is seen numerically in the lower-right panels of Figs. 6 and 7, which show that the circular polarization is identically zero. However, total internal reflection leads to matrices of the form:
E̲=[XXX0XXXXXXXX0XXX],
which still cannot convert unpolarized light into circular polarization. This is seen numerically in the lower-right panels of 8 and 9. However, if the incident upwelling underwater light is linearly polarized and undergoes total internal reflection, the reflected light can be circularly polarized. (This is, of course, the same physics used to generate circular polarization by total internal reflection in a Fresnel rhomb.) Total internal reflection is a potential source of elliptical or circular polarization in the ocean.

5. Sea Surface Reflectance and Transmittance of Polarized Sky Radiance

Simulations as just seen using a collimated unpolarized light source in an otherwise black sky or ocean are useful to illustrate the fundamental reflection and transmission features of the sea surface. However, real skies are distributed light sources and are usually partially polarized.

The sky radiance distribution is modeled as follows. Using conditions typical of a clear maritime atmosphere (80% relative humidity, 2.5 cm of precipitable water vapor, marine aerosols, ozone concentration of 300 Dobson units), the RADTRAN irradiance model [29] is run to obtain the direct (Eddir) and diffuse (Eddif) sky irradiances for the given solar zenith angle θsun and a wavelength of 550 nm. The quad containing the Sun is given a radiance of magnitude Isun=Eddir/(μsunΩsun), where μsun is the average of cosθ over the Sun’s quad, and Ωsun is the solid angle of this quad. The Sun’s direct beam is taken to be unpolarized: S̲sun=[Isun,0,0,0]T. The angular pattern of the total radiance in the remaining quads is computed using the semi-empirical, clear-sky model of Harrison and Coombes [30]. Those diffuse sky radiances are scaled so that the total irradiance from the nonsolar quad equals Eddif. This gives a calibrated sky I(θ,ϕ) that reproduces the RADTRAN direct and diffuse irradiances. The polarization pattern of the diffuse sky radiance is that of a molecular (Rayleigh) single-scattering sky with a depolarization factor of 0.0318, as given by the equations of Tilstra et al. [31]. Figure 12 shows, for a solar zenith angle of 50 deg, the resulting pattern of I, horizontal versus vertical (Q/I) and ±45deg(U/I) polarization, and the degree of total polarization DoP=(Q2+U2+V2)12/I. The circular polarization is identically zero. In the numerical calculations below, these equations are used to compute quad-averaged sky Stokes vectors.

 figure: Fig. 12.

Fig. 12. Sky radiance distribution for single scattering by atmospheric molecules according to the Rayleigh scattering equations. The upper-left panel shows the total radiance magnitude I relative to 1 in the Sun’s direction. The Sun’s location is at the left side of the plotted hemisphere of sky directions. The other panels show the horizontal versus vertical, ±45deg, and degree of total polarization in percent.

Download Full Size | PPT Slide | PDF

Multiple scattering and aerosols tend to depolarize the sky radiance. Monte Carlo simulations of atmospheric radiance patterns [32], including both molecular and aerosol components and multiple scattering, show a variety of polarization patterns, which depend on the type and concentration of the aerosols. Thus, a polarization pattern corresponding to a single-scattering Rayleigh sky serves as an upper limit to the degree of polarization that occurs for a very clear atmosphere. Use of this sky radiance distribution in simulations, thus, gives an upper bound for the effects of polarization compared to simulations using unpolarized ray tracing but only for the Rayleigh polarization pattern.

Figure 13 shows the surface-reflected sky radiance for the Sun at 50 deg and the partially polarized sky of Fig. 12. The sea surface is an FFT surface for a wind speed of 10ms1, generated with Lx=200m and Nx=1024. The increase in the radiance magnitude I going toward the horizon at ϕv=0 in the upper-left plot is due to the rapidly increasing Fresnel reflectance for the wave facet reflection angles that connect the bright region of sky near the Sun to the near-horizon directions. The level of statistical noise in this plot is generally a few percent (recall Fig. 10) but can be more for directions receiving relatively few photons. For example, the DoP values for the quad next to the horizon (θv=87.5) and ϕv=±90deg are 86.69% and 87.79%. Such differences can cause slight asymmetries in the color coding of plots for ±ϕv values. The I value for (θv,ϕv)=(87.5,0) (the orange quad in the upper-left panel) was 39.98Wm2sr1nm1 for this plot and 44.67 for an independent run starting with a different seed for random number generation. For this Sun zenith angle and atmospheric conditions, RADTRAN gives Eddir=0.6561 and Eddif=0.3509Wm2nm1. The incident sky Stokes vectors are now diffuse vectors with units of spectral radiance, so the surface reflectance function is an R̲aw, as defined in Eq. (15). All sky quads now, in principle, contribute to the reflected radiance in each quad:

S̲refl(Qkl)=ijR̲aw(QijQkl)S̲sky(Qij),
where the sum over ij indicates a sum over all sky quads.

 figure: Fig. 13.

Fig. 13. Surface-reflected radiance for a single-scattering Rayleigh sky and a wind speed of 10ms1. The Sun zenith angle is 50 deg and the viewing geometry corresponds to Fig. 6.

Download Full Size | PPT Slide | PDF

Figure 14 shows the transmitted radiance corresponding to Fig. 13. As in Fig. 7, most of the transmitted radiance is concentrated into a few quads near the refracted direction of the Sun’s direct beam, and the radiance in those quads is only weakly polarized. However, the state of linear polarization away from the direct beam direction is determined by the sky polarization in directions away from the Sun, and the degree of polarization can be large. However, those directions contain very little energy. For a level surface, all of the transmitted energy would be contained within the Snell’s cone of angle 48.2 deg from the nadir. Now, however, for the 10ms1 wind speed, the rough surface allows radiance to be transmitted into all downwelling directions, although directions near the horizon receive very little energy.

 figure: Fig. 14.

Fig. 14. Surface-transmitted radiance for a single-scattering Rayleigh sky and a wind speed of 10ms1. The Sun zenith angle is 50 deg, and the viewing geometry corresponds to Fig. 7.

Download Full Size | PPT Slide | PDF

6. Irradiance Reflectance

Surface irradiance reflectance is fundamental to energy transfer across the air-water surface. Although it is the total energy (without regard for the state of polarization) that is usually of interest, radiance reflectance by the surface depends on the state of polarization of the incident sky radiance, which, in turn, depends on the direction relative to the Sun. Likewise, the reflectance depends on the angle between the incident radiance and the normal to the sea surface. It is, therefore, worthwhile to quantify the differences that may arise from the use of Cox–Munk surfaces and/or unpolarized ray tracing, compared with the more realistic FFT surfaces and polarized ray tracing.

The sky total downwelling plane irradiance is

Ed(sky)=ijIsky(Qij)μiΩij.
For the sky used to generate Fig. 13, Ed(sky)=1.0070Wm2nm1, as expected from the input Eddir and Eddif values. The upwelling, surface-reflected irradiance computed from Irefl is Eu(refl)=0.05098Wm2nm1. Thus, for the conditions of Fig. 13, RsurfEu(refl)/Ed(sky)=0.0506. That is, about 5% of the incident energy is reflected by the sea surface, and the remaining 95% enters the water column. It should be noted that Rsurf is the irradiance reflected by the surface itself, not the albedo of the water body. The albedo is defined in the same way, but the Eu includes surface-reflected and water-leaving radiance. The albedo is, therefore, always somewhat greater than Rsurf. For the 105 surface realizations used here, the statistical noise in Rsurf is never more than 0.2% and usually much less. This variability is too small to display on the curves of Fig. 15.

 figure: Fig. 15.

Fig. 15. Comparison of surface irradiance reflectances Rsurf as functions of solar zenith angle and wind speed for a polarized sky and various combinations of sea surface model (FFT versus Cox–Munk), ray tracing (polarized versus unpolarized), single versus multiple scattering, Sun azimuthal angle, and wave age.

Download Full Size | PPT Slide | PDF

The upper-left panel of Fig. 15 shows the dependence of Rsurf on whether the ray tracing is polarized or unpolarized, for FFT surfaces and selected wind speeds. The polarized sky radiances were generated, as described above, for nominal solar zenith angles from 0 (the polar cap) to 87.5 deg (the quad next to the horizon). The corresponding first section of Table 2 shows the percent differences computed as 100[Rsurf(polarized)Rsurf(unpolarized)]/Rsurf(polarized) for Sun zenith angles of 0, 50, and 87.5 deg. For the Sun at the zenith, the difference in Rsurf for polarized versus unpolarized ray tracing is about 11%, i.e., polarized ray tracing gives a greater reflectance. For a solar zenith angle of 50 deg, the difference is only about 1%. For the near-horizon quad, the difference is 8.1% for a level surface to 11.7% for a wind speed of 15ms1, i.e., polarized ray tracing gives somewhat less reflectance. The smooth black line in this panel is the Fresnel reflectance for unpolarized light. This curve is valid only for uniform sky radiance distributions and a level sea surface. As seen here, for nonuniform skies and/or rough sea surfaces, the actual surface irradiance reflectance can be greater or less than the Fresnel reflectance. The decrease in irradiance reflectance for the Sun very near the horizon has been explained previously [33]. In brief, as the Sun goes from the zenith to the horizon, the fraction of diffuse sky irradiance increases from about 25% to over 99%. For solar zenith angles less than about 80 deg, the rapidly increasing Fresnel reflectance of the direct solar beam dominates and Rsurf increases. But for Sun zenith angles greater than 80 deg, diffuse radiance from directions nearer to the zenith dominates, and the lower Fresnel reflectance for those directions reduces Rsurf. A higher fraction of the total radiance, thus, begins to enter the water as the Sun nears the horizon and Rsurf decreases.

Tables Icon

Table 2. Differences in Surface Irradiance Reflectance Rsurf as Computed in Various Waysa

The upper-right panel of Fig. 15 compares Cox–Munk and FFT surfaces for polarized ray tracing. The corresponding second section of Table 2 gives 100[Rsurf(FFT)Rsurf(Cox-Munk)]/Rsurf(FFT). Both surface models are, of course, identical for zero wind speed, corresponding to a level sea surface. There is very little difference in Cox–Munk and FFT surfaces for solar angles less than about 50 deg. For wind speeds less than about 10ms1, the Cox–Munk surface has a greater Rsurf at large solar zenith angles, and, for wind speeds greater than 10ms1, the reverse is true. The two surfaces give almost identical Rsurf for all zenith angles when the wind speed is 10ms1.

The middle-left panel of this figure compares the commonly used unpolarized, Cox–Munk calculations with polarized, FFT computations. The corresponding section of Table 2 gives 100[Rsurf(FFT,pol)Rsurf(Cox-Munk,unpol)]/Rsurf(FFT,pol). Rsurf is 10%–12% greater for FFT surfaces than for Cox–Munk when the Sun is at the zenith. However, the Cox–Munk values are greater for large zenith angles, for which the difference is as much as 18% at low wind speeds.

The above comparisons have included all orders of multiple scattering during the ray tracing. The middle-right panel of Fig. 15 and the fourth section of Table 2 compare Rsurf computed for multiple versus single scattering, for polarized ray tracing and FFT surfaces. The tabulated differences are now 100[Rsurf(multiple)Rsurf(single)]/Rsurf(multiple). Multiple scattering always gives a greater Rsurf than single scattering, by about 2% to 8%, depending on the solar zenith angle and wind speed, with the greatest effect for nearly horizontal incident rays.

The preceding comparisons have been made with the Sun at ϕsun=0, so that the Sun’s incoming rays are parallel to the wind. The lower-left panel of Fig. 15 shows the results for FFT surfaces and polarized ray tracing when the Sun is located at ϕsun=90deg, for which the Sun’s incident rays are perpendicular to the wind speed. In this case, the surface reflectance is somewhat greater than for ϕsun=0. This is because the cross-wind wave slope variance is less than the along-wind slope, as seen in Eq. (1b) versus Eq. (1a). Thus, the incoming rays at right angles to the wind see a somewhat smoother sea surface, and the reflectance is, consequently, larger than for rays incident in the along-wind direction. The azimuthal angle is not defined for polar caps, so there is no difference for the Sun at the zenith.

The wave spectrum of [20] includes a parameter Ωc, which modifies the spectrum according to the age of the waves, that is, how long it has been since the wind began to blow. Ωc=5 gives a very young sea for which the wind has been blowing only a short time and, thus, only capillary and small gravity waves are present. Ωc decreases as the waves age and larger gravity waves develop. Ωc=1 gives a mature sea, and Ωc=0.84 gives the limiting steady-state case of a fully developed sea. The lower-right panel of the figure shows the effect of wave age for very young (Ωc=5) versus fully developed (Ωc=0.84) seas, for FFT surfaces and polarized ray tracing (and ϕsun=0). The last section of Table 2 shows the corresponding percent differences 100[Rsurf(fully devel.)Rsurf(young)]/Rsurf(fully devel.). Young seas have a somewhat smoother surface and, thus, a greater reflectance because they have not had the time to develop the slope variance that comes with a more mature sea.

In summary, the various approximations—Cox–Munk, unpolarized ray tracing, and single-scattering—and environmental conditions—Sun azimuthal angle relative to the wind direction, and wave age—all give differences of order 10% in Rsurf compared to the values computed with FFT surfaces and polarized ray tracing that includes all orders of multiple scattering. However, as is clear from Fig. 15, these differences are considerably less than the systematic decrease in Rsurf at large solar zenith angles due to wind speed. For zenith angles greater than 70 deg, Rsurf decreases by about 30% when the wind increases from W=0 to 5ms1 and by almost 45% between W=0 and W=15ms1. For solar angles less than 60 deg, Rsurf decreases by 8% between W=0 and 5ms1, and by 12% between W=0 and 15ms1.

7. Radiance Reflectance Factors

The remote-sensing reflectance Rrs is defined as the ratio of the water-leaving radiance Lw just above the sea surface to the downwelling plane irradiance Ed incident onto the sea surface: RrsLw/Ed. Rrs is the basis of most ocean-color remote sensing because this apparent optical property is strongly correlated with in-water optical properties and is minimally sensitive to environmental effects such as solar zenith angle, sky conditions, and wind speed [34]. However, any measurement of upwelling radiance Lu just above the surface includes both Lw and surface-reflected sky radiance Lsr, so that Lw=LuLsr. For a level sea surface, Lsr comes only from sky radiance in the direction that is specularly reflected into the sensor. In this case Lsr(θv,ϕv)=ρ(θv,ϕv)Lsky(θv,ϕv), where Lsky is the sky radiance and ρ is the radiance reflectance for the angle of incidence θv that connects Lsky and Lsr by specular reflection [35]. In Lsr(θv,ϕv), θv is measured from the nadir; in Lsky(θv,ϕv), θv is the same angle but measured from the zenith. For a level surface and unpolarized radiance, the radiance reflection ρ equals the unpolarized Fresnel reflectance, which is a function only of θv.

For windblown sea surfaces, Lsr arises from sky radiance that is reflected by tilted wave facets from all sky directions into the sensor. However, for clear skies, it is still reasonable [18] to define a radiance reflection factor ρ that converts a sky measurement in a particular direction into an estimate of the surface reflectance: ρ(θv,ϕv)Lsr(θv,ϕv)/Lsky(θv,ϕv). In this case, ρ no longer equals the Fresnel reflectance but rather depends on solar zenith angle, viewing direction, wave state, and sky condition. The purpose of the ρ factor is to account for all surface-reflected radiance, including direct (Sun glint) and diffuse (reflected background sky) contributions. Of course, when making measurements, a viewing direction is chosen to minimize the direct Sun glint contribution. However, some Sun glint is usually present, which is accounted for in the ρ factor. Thus, ρ is best interpreted as a scale factor that converts a sky radiance measurement in a particular direction, Lsky(θv,ϕv), into the surface reflectance in the corresponding specular direction, for given sky and sea surface conditions. The desired water-leaving radiance is then obtained from Lw=LuρLsky. Lu and Lsky are measured by the same instrument, which is assumed to be insensitive to the state of polarization. Accurate determination of the value of ρ is the key to obtaining an accurate value for Lw.

Previous studies [18,36] investigated the dependence of ρ on the Sun zenith angle, viewing direction, windspeed, and wavelength for clear sky conditions. (A table of ρ values is available [37].) These tabulated values are widely used [e.g., 3840]. However, those values were computed using Cox–Munk surfaces and unpolarized ray tracing. The next section reevaluates ρ using FFT surfaces and polarized ray tracing.

A. ρ for a Level Sea Surface

Consider first the case of a level sea surface and a partially polarized Rayleigh sky radiance distribution, as described above and illustrated in Fig. 12, for a solar zenith angle of θsun=50deg. The top panel of Fig. 16 shows the reflectance factor for a level surface and unpolarized ray tracing. The values computed by Monte Carlo ray tracing reside along the curve for unpolarized Fresnel reflectance. The computed values do not depend on the azimuthal viewing direction ϕv, which is measured relative to the Sun’s azimuthal direction. That is, ϕv=0 corresponds to looking toward the Sun, ϕv=90deg is looking at right angles to the Sun’s incident rays, and ϕv=180deg is looking away from the Sun. The bottom panel shows ρ(θv,ϕv) for azimuthal directions of ϕv=90 and 135 deg. Now, however, the ρ value depends on the azimuthal direction, even for a level sea surface. This ϕv dependence is easily understood as via a simple numerical example.

 figure: Fig. 16.

Fig. 16. Radiance reflectance factors ρ(θv,ϕv) for a level sea surface and a Sun zenith angle of 50 deg. The top panel is for unpolarized ray tracing; the bottom panel is for polarized ray tracing. In both cases, the incident radiance is that of the single-scattering Rayleigh sky shown in Fig. 12. The dashed and dotted lines show the Fresnel reflectance for incident radiance, which is linearly polarized perpendicular (R) or parallel (R) to the incident meridian plane. The solid line without symbols is the Fresnel reflectance for unpolarized incident radiance.

Download Full Size | PPT Slide | PDF

The surface radiance reflectance matrix, as defined in Eq. (15) and computed by Monte Carlo simulation for incident and reflected angles θi=θr=40deg, is

R̲(40,ϕ,40,ϕ)=[2.566×1021.963×1023.521×1070.01.963×1022.566×1027.489×1070.03.521×1077.489×1071.616×1020.00.00.00.01.616×102].
This matrix is independent of azimuthal angle ϕ for a level sea surface. The Stokes vectors of the sky radiances in directions (θv,ϕv)=(40,90) and (40, 135) are
S̲sky(40,90)=[4.931×1022.403×1021.583×1020.0],
and
S̲sky(40,135)=[3.932×1021.418×1023.262×1020.0].
These sky radiances give (DoP,Q/I,U/I)=(58.35,32.10,48.73) in percent for (40,90) and (90.45,36.06,82.95) for (40, 135), which match the values plotted at those directions in Fig. 12.

Multiplying R̲(40,ϕ,40,ϕ) times these sky radiance Stokes vectors gives the corresponding surface-reflected radiances:

S̲sr(40,90)=[9.547×1045.617×1043.883×1040.0]
and
S̲sr(40,135)=[1.287×1031.136×1035.269×1040.0].
The corresponding ρ values are
ρ(40,90)=IsrIsky=9.547×1044.931×102=0.0194ρ(40,135)=IsrIsky=1.287×1033.932×102=0.0327,
which are the values plotted for θv=40 in the bottom panel of Fig. 16.

For unpolarized ray tracing, only the (1, 1) element of R̲(40,ϕ,40,ϕ) is nonzero. In this case, for the same polarized sky radiances, the reflected radiances are

S̲sr(40,90)=[1.265×103,0,0,0]T
and
S̲sr(40,135)=[1.009×103,0,0,0]T.
These values give
ρ(40,90)=1.265×1034.931×102=0.0257
and
ρ(40,135)=1.009×1033.932×102=0.0257,
which are the values plotted in the top panel of Fig. 16.

Thus, ρ is independent of the azimuthal direction and the polarization state of the sky for unpolarized ray tracing. However, for polarized ray tracing, ρ depends on polar and azimuthal angles even for a level sea surface because the state of the sky polarization is different for different directions. The polarization state of the incident radiance determines how much total radiance is reflected by the surface, hence the value of ρ. This is an additional complication in the determination of the value of ρ to be used in processing measured radiances. Not only does ρ depend on solar zenith angle, viewing direction, and sea state, it also depends on the state of polarization of the sky radiance distribution, even though the sensor itself is not sensitive to polarization.

B. ρ for Windblown Surfaces

Figure 17 shows for reference selected ρ(θv,ϕv) values as computed by unpolarized ray tracing and Cox–Munk sea surfaces. ρ is shown a function of off-nadir viewing direction θv for Sun zenith angles of 30 and 60 deg, azimuthal viewing directions of 90 and 135 deg (relative to the Sun at ϕv=0), and wind speeds of 2, 5, 10, and 15ms1. The plot for 10ms1 reproduces the θsun=30 and 60 deg curves of [18, Fig. 6] but with less statistical noise because more rays are traced. When the Sun is high in the sky (illustrated here by the θsun=30 curves) and the sea surface is rough (wind speeds of 5ms1 or greater in the figure), the ϕv=90 curves have large ρ values at small θv values because of Sun glint seen by the sensor measuring Lu. For large θv, ρ becomes large because of the rapidly increasing Fresnel reflectance. Based on plots such as these, [18] argued that θv=30 or 40 deg and ϕv=135 is a reasonable choice for field measurements. This viewing direction is a compromise that avoids most Sun glint, keeps the instrument from looking at its own shadow, and gives minimal variability in the ρ values.

 figure: Fig. 17.

Fig. 17. Radiance reflectance factors ρ computed using Cox–Munk surfaces and unpolarized ray tracing. The curves are for Sun zenith angles θsun=30 and 60 deg, azimuthal viewing directions of ϕv=90 and 135 deg, and wind speeds of 2, 5, 10, and 15ms1; θv is the off-nadir viewing direction.

Download Full Size | PPT Slide | PDF

The solid-line curves in Fig. 18 show the corresponding values computed using FFT surfaces and polarized ray tracing. There is now a greater spread of ρ values for the given environmental conditions and viewing geometry. For θv=40deg and wind speeds of 5 and 10ms1, the unpolarized case gives tightly grouped values of ρ=0.02800.0285 and 0.0316–0.0374, respectively. However, the values determined by FFT surfaces and polarized ray tracing have a larger spread of values: 0.0179–0.0411 for 5ms1 and 0.0195–0.0458 for 10ms1. There is, thus, a factor-of-two spread in ρ values for each wind speed when computed with FFT surfaces and polarized ray tracing, versus less than 20% spread when computed with Cox–Munk surfaces and unpolarized ray tracing. This figure also shows the curves for an azimuthal viewing direction of ϕv=135, but for the Sun at a right angle to the wind speed and for a very young sea (Ωc=5 in the wave spectrum).

 figure: Fig. 18.

Fig. 18. Radiance reflectance factors ρ computed using FFT surfaces and polarized ray tracing. The solid curves correspond to the viewing and wave conditions of Fig. 17. These curves are for a fully developed sea and the Sun’s incident rays parallel to the wind speed (ϕsun=0). The dotted curves and diamond symbols show the curves for θsun=30 and 60 deg, azimuthal viewing direction of ϕv=135, but with the Sun’s azimuthal angle at ϕsun=90, so that the incoming rays are perpendicular to the wind direction. The dashed curves and box symbols show the cases of θsun=30 and 60 deg azimuthal viewing direction of ϕv=135, ϕsun=0, for a very young sea with Ωc=5.

Download Full Size | PPT Slide | PDF

Figure 19 displays ρ values in additional ways. The upper-left panel shows a plot of ρ(θv,ϕv) for a single Sun angle of θsun=40deg located at ϕv=0, as computed using FFT surfaces and polarized ray tracing. The radial distance in the plot is the off-nadir viewing angle θv; the center of the plot is looking straight down at the sea surface, and the rim of the plot is looking toward the horizon. The Sun’s specular direction for a level surface is shown by the white quad centered at (40,0). The yellow to red colors in this area of the plot show the high ρ values needed to remove Sun glint. The orange to red colors around the rim are high ρ values resulting from high Fresnel reflectances at large incident angles. The curves of Fig. 18 display the values from the center to the rim of the plot along the ϕv=90 and 135 deg directions (but for different Sun zenith angles). The black dot indicates the ρ(θv,ϕv)=(40,135) viewing direction. In the present example, ρ(40,135)=0.045. The upper-right panel of the figure shows ρ(θsun,ϕv) for an off-nadir viewing direction of θv=40deg. Now the center of the plot corresponds to the Sun at the zenith, and the rim to the Sun at the horizon. Again, the yellow-to-red colors show the high ρ resulting from Sun glint with looking in directions too near the Sun’s azimuthal direction. The blue colors along ϕv=135 indicate ρ values in the 0.02–0.05 range for all Sun zenith angles except when the Sun is at the zenith or nearly so, in which case the view includes a significant amount of direct Sun glint. The bottom two panels of the figure show the ratios of ρ computed using FFT surfaces and polarized ray tracing to the values computed using Cox–Munk surfaces and unpolarized ray tracing. Radial lines in these figures correspond to the ratios of the curves in Fig. 18 to the corresponding curves of Fig. 17 for selected Sun and viewing directions. For the Rayleigh sky and 10ms1 conditions of this figure, the lower-left panel shows that the ratio of ρ(FFT,polarized) to ρ (Cox–Munk, unpolarized) is 0.91 for θsun=40deg and ρ(θv,ϕv)=(40,90) and is 1.31 for ρ(θv,ϕv)=(40,135) (the viewing direction indicated by the black dot).

 figure: Fig. 19.

Fig. 19. ρ for a Rayleigh sky and a wind speed of 10ms1. The upper-left panel shows ρ as a function of off-nadir and azimuthal viewing directions (relative to the Sun’s azimuth at ϕv=0) for a single Sun zenith angle of θsun=40deg, computed using FFT surfaces and polarized ray tracing. The upper-right panel shows ρ as a function of the Sun’s zenith angle (radial distance in the plot) and viewing azimuth for a fixed off-nadir viewing direction of θv=40deg. The black dots indicate the ρ(θv,ϕv)=(40,135) viewing direction. The white regions indicate the specular viewing direction for θsun=40deg. The lower two panels show the ratios of ρ computed for FFT surfaces and polarized ray tracing to the values computed for Cox–Munk surfaces and unpolarized ray tracing, corresponding to the respective left and right upper panels.

Download Full Size | PPT Slide | PDF

The simulations of this section used a wavelength of 550 nm for generation of the Rayleigh single-scattering sky radiance distribution. However, ρ also depends on wavelength [36] because the degree of sky polarization and the ratio of diffuse to direct solar irradiance vary considerably from blue to red wavelengths. Similarly, in a real atmosphere, the degree and pattern of sky polarization for a given wavelength is influenced by the aerosol type and optical thickness [32]. However, Figs. 17 and 18 are sufficient to show that the previously computed ρ values based on Cox–Munk surfaces and unpolarized ray tracing [18,37] are not adequate for accurate estimation of water-leaving radiance from above-surface measurements of upwelling and sky radiances. The viewing directions of θv40deg and ϕv135deg (relative to the Sun’s azimuthal direction) remain a reasonable choice. Harmel et al. [17] studied the effects of polarization on a normalized reflectance, which is equivalent to the ρ of the present paper. Their study used an analytical representation of Cox–Munk surfaces with single scattering but included atmospheric conditions for a range aerosol optical thicknesses. They, too, reached the same conclusion as the present study: it is imperative to account for polarization effects when estimating the amount of surface-reflected radiance in a measured total above-surface radiance.

Using the correct value of ρ is critical for the estimation of Lw=LuρLsky because the measured Lsky(θv,ϕv) is typically an order of magnitude or more greater than Lw. Thus, a small error in ρ can give a large error in Lw. To address the deficiency of the ρ values previously computed using Cox–Munk surfaces and unpolarized ray tracing, a new table of ρ values has been computed using FFT surfaces and polarized ray tracing for the Rayleigh sky described above. This table is available online [37].

8. Conclusions

Wave variance spectra and fast Fourier transforms are widely used in science and the movie and video-gaming industries for the generation of random sea surfaces. The latter applications, especially, frequently ignore the scale factors needed to guarantee conservation of wave energy when going from variance spectra to sea surfaces. This paper has shown in Section 2 how the equations can be formulated to guarantee wave energy conservation. A new technique was presented in Section 2.C to account for wave elevation and wave slope variances in generated sea surfaces, without the need for large numbers of spatial grid points. Random sea surface realizations were then used in Monte Carlo simulations of polarized ray tracing of air- and water-incident light. The ray-tracing algorithm outlined in Section 3.B follows each air- or water-incident ray through any number of ray-wave interactions, during which a single incident ray can generate any number of daughter rays by reflection and transmission. Electromagnetic energy is exactly conserved photon by photon because all incident and daughter rays are traced to completion when they leave the sea surface region to be tallied.

Sea surface energy reflection and transmission were first illustrated in Section 4 for the easily understood case of an unpolarized, collimated incident light source. Section 5 then illustrated reflection and transmission of partially polarized sky radiance corresponding to the idealized case of a single-scattering, molecular (Rayleigh-scattering) atmosphere.

Section 6 studied the irradiance reflection properties of windblown sea surfaces. Irradiance reflectances Rsurf computed using Cox–Munk wind-speed wave slope statistics and/or unpolarized ray tracing were compared with the values computed using FFT surfaces and/or polarized ray tracing. It was found (Fig. 15 and Table 2) that differences in Rsurf due to ignoring polarization were typically greater than differences due to using Cox–Munk rather than FFT surfaces or differences due to using single rather than multiple scattering. The differences between Rsurf computed using Cox–Munk surfaces and unpolarized ray tracing versus FFT surfaces and polarized ray tracing range from 10% to 12% for the Sun at the zenith to as much as 18% for the Sun near the horizon. This has implications for modeling of sea surface irradiance reflectance as needed for ecosystem models [33].

The last section studied radiance reflectance factors ρ, which are used in the estimation of water-leaving radiances from above-surface measurements of sky and upwelling radiances. These factors were shown to depend on the polarization state of the sky radiance distribution even for measurements made with sensors that are not sensitive to polarization. Moreover, for a given viewing direction of the sensors, ρ as computed for FFT surfaces and polarized ray tracing, was found to vary more with the solar zenith angle and wind speed than the factors previously computed [18] using Cox–Munk surfaces and unpolarized ray tracing (Fig. 17 versus Fig. 18). The polarization dependence of ρ was studied only for the limiting cases of no sky polarization and strongly polarized sky radiance computed for a molecular, single-scattering sky, which is a model for extremely clear atmospheric conditions. The degree and pattern of atmospheric polarization depend on aerosol type and concentration. The ρ factors for atmospheres with various aerosols could be obtained by incorporating the surface modeling techniques developed here into an atmospheric radiative transfer model capable of simulating the sky polarization conditions for specific aerosol and other atmospheric conditions. The resulting ρ factors would likely differ somewhat from the values shown here. Moreover, the ρ factors depend on wavelength [36], and a detailed study of this wavelength dependence requires use of a vector atmospheric model with proper simulation of the sea surface reflectance. These observations highlight the difficulty of determining the correct value of ρ to use with particular environmental conditions of Sun zenith angle, wind speed, and atmospheric conditions.

In a similar fashion, the effects of hydrosols on the polarization state of underwater radiance distributions can be accurately simulated only via numerical models of underwater polarized radiative transfer. Detailed calculations of atmospheric and underwater polarized radiances are beyond the scope of this paper, which seeks only to show how to incorporate sea surface effects into coupled atmosphere–ocean models, so that the wave effects on transmitted and reflected polarized radiance are accurately modeled.

FUNDING INFORMATION

National Aeronautics and Space Administration (NASA) (NNH12CD06C, NNX14AQ49G).

The author thanks Yogesh Agrawal for an enlightened discussion on wave spectra, which led to the algorithm in Section 2.C. Four anonymous reviewers made a number of helpful comments. The development of surface generation and polarized ray tracing algorithms was supported by NASA contract NNH12CD06C. Publication of this paper was supported by NASA grant NNX14AQ49G.

References

1. C. Cox and W. Munk, “Measurement of the roughness of the sea surface from photographs of the Sun’s glitter,” J. Opt. Soc. Am. 44, 838–850 (1954).

2. R. W. Preisendorfer and C. D. Mobley, “Albedos and glitter patterns of a wind-roughened sea surface,” J. Phys. Oceanogr. 16, 1293–1316 (1986).

3. M. P. Levesque and D. St-Germain, “Generation of synthetic IR sea images,” Proc. SPIE 1331, 354–357 (1990).

4. M. P. Levesque, “Dynamic sea image generation,” Proc. SPIE 1486, 294–300 (1991). [CrossRef]  

5. J. W. McLean, “Modeling of ocean wave effects for LIDAR remote sensing,” Proc. SPIE 1302, 480–491 (1990). [CrossRef]  

6. J. W. McLean and J. D. Freeman, “Effects of ocean waves on airborne lidar imaging,” Appl. Opt. 35, 3261–3269 (1996). [CrossRef]  

7. S. Kay, J. Hedley, S. Lavender, and A. Nimmo-Smith, “Light transfer at the ocean surface modeled using high resolution sea surface realizations,” Opt. Express 19, 6493–6504 (2011). [CrossRef]  

8. H. M. Tulldahl and K. O. Steinvall, “Simulation of sea surface wave influence on small target detection with airborne laser depth sounding,” Appl. Opt. 43, 2462–2483 (2004). [CrossRef]  

9. G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphere-ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. 34, 1453–1472 (1989). [CrossRef]  

10. H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,” Appl. Opt. 40, 400–412 (2001). [CrossRef]  

11. M. Chami, R. Santer, and E. Dilligeard, “Radiative transfer model for the computation of radiance and polarization in an ocean-atmosphere system: polarization properties of suspended matter for remote sensing,” Appl. Opt. 40, 2398–2416 (2001). [CrossRef]  

12. P.-W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17, 2057–2079 (2009). [CrossRef]  

13. A. Tonizzo, A. Gilerson, T. Harmel, A. Ibrahim, J. Chowdhary, B. Gross, F. Moshary, and S. Ahmed, “Estimating particle composition and size distribution from polarized water-leaving radiance,” Appl. Opt. 50, 5047–5058 (2011). [CrossRef]  

14. J. Chowdhary, B. Cairns, and L. D. Travis, “Contribution of water-leaving radiances to multiangle, multispectral polarimetric observations over the open ocean: bio-optical model results for case 1 waters,” Appl. Opt. 45, 5542–5567 (2006). [CrossRef]  

15. P.-W. Zhai, G. W. Kattawar, and P. Yang, “Impulse response solution of the three-dimensional vector radiative transfer equation in atmosphere-ocean systems. I. Monte Carlo method,” Appl. Opt. 47, 1037–1047 (2008). [CrossRef]  

16. P.-W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Luker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010).

17. T. Harmel, A. Gilerson, A. Tonizzo, J. Chowdhary, A. Weidemann, R. Arnone, and S. Ahmed, “Polarization impacts on the water-leaving radiance retrieval from above-water radiometric measurements,” Appl. Opt. 51, 8324–8340 (2012). [CrossRef]  

18. C. D. Mobley, “Estimation of the remote-sensing reflectance from above-surface measurements,” Appl. Opt. 38, 7442–7455 (1999). [CrossRef]  

19. M. L. Banner, “Equlibrium spectra of wind waves,” J. Phys. Oceanogr. 20, 966–984 (1990). [CrossRef]  

20. T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, “A unified directional spectrum for long and short wind-driven waves,” J. Geophys. Res. 102, 15781–15796 (1997). [CrossRef]  

21. S. R. Massel, “On the geometry of ocean surface waves,” Oceanologia 53, 521–548 (2011).

22. C. D. Mobley, “Modeling sea surfaces: a tutorial on Fourier transform techniques,” Tech. Rep. (Sequoia Scientific, Inc., 2014).

23. J. Tessendorf, “Simulating ocean water,” Tech. Rep. (SIGGRAPH Course Notes, 2004).

24. M. S. Longuet-Higgins, D. E. Cartwright, and N. D. Smith, “Observations of the directional spectrum of sea waves using the motions of a flotation buoy,” in Ocean Wave Spectra (Prentice-Hall, 1963), pp. 111–136.

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

26. C. D. Mobley, “HydroPol mathematical documentation: invariant imbedding theory for the vector radiative transfer equation,” Tech. Rep. (Sequoia Scientific, 2014).

27. M. I. Mischenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

28. P.-W. Zhai, G. W. Kattawar, and Y. Hu, “Comment on the transmission matrix for a dielectric interface,” J. Quant. Spectrosc. Radiat. Transfer 113, 1981–1984 (2012). [CrossRef]  

29. W. W. Gregg and K. L. Carder, “A simple spectral solar irradiance model for cloudless maritime atmospheres,” Limnol. Oceanogr. 35, 1657–1675 (1990). [CrossRef]  

30. A. W. Harrison and C. A. Coombes, “Angular distribution of clear sky short wavelength radiance,” Solar Energy 40, 57–63 (1988). [CrossRef]  

31. L. G. Tilstra, N. A. J. Schutgens, and P. Stammes, “Analytical calculation of Stokes parameters Q and U of atmospheric radiation,” Tech. Rep. (Koninklijk Nederlands Meteorologisch Instituut, 2003).

32. C. Emde, R. Buras, B. Mayer, and M. Blumthaler, “The impact of aerosols on polarized sky radiance: model development, validation, and applications,” Atmos. Chem. Phys. 10, 383–396 (2010). [CrossRef]  

33. C. D. Mobley and E. S. Boss, “Improved irradiances for use in ocean heating, primary production, and photo-oxidation calculations,” Appl. Opt. 51, 6549–6560 (2012). [CrossRef]  

34. C. D. Mobley, “Reflectances,” http://www.oceanopticsbook.info/view/overview_of_optical_oceanography/reflectances (2010).

35. K. L. Carder and R. G. Steward, “A remote-sensing reflectance model of a red-tide dinoflagellate off West Florida,” Limnol. Oceanogr. 30, 286–298 (1985). [CrossRef]  

36. Z.-P. Lee, Y.-H. Ahn, C. D. Mobley, and R. Arnone, “Removal of surface-reflected light for the measurement of remote-sensing reflectance from an above-surface platform,” Opt. Express 18, 26313–26324 (2010). [CrossRef]  

37. C. D. Mobley, “Surface reflectance factors,” http://www.oceanopticsbook.info/view/remote_sensing/level_3/surface_reflectance_factors (2013).

38. S. B. Hooker, G. Zibordi, J. F. Berthon, and J. W. Brown, “Above-water radiometry in shallow coastal waters,” Appl. Opt. 43, 4254–4268 (2004). [CrossRef]  

39. G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009). [CrossRef]  

40. T. Harmel, A. Gilerson, S. Hlaing, A. Tonizzo, T. Legbandt, A. Weidemann, R. Arnone, and S. Ahmed, “Long Island Sound Coastal Observatory: assessment of above-water radiometric measurement uncertainties using collocated multi and hyperspectral systems,” Appl. Opt. 50, 5842–5860 (2011). [CrossRef]  

References

  • View by:

  1. C. Cox and W. Munk, “Measurement of the roughness of the sea surface from photographs of the Sun’s glitter,” J. Opt. Soc. Am. 44, 838–850 (1954).
  2. R. W. Preisendorfer and C. D. Mobley, “Albedos and glitter patterns of a wind-roughened sea surface,” J. Phys. Oceanogr. 16, 1293–1316 (1986).
  3. M. P. Levesque and D. St-Germain, “Generation of synthetic IR sea images,” Proc. SPIE 1331, 354–357 (1990).
  4. M. P. Levesque, “Dynamic sea image generation,” Proc. SPIE 1486, 294–300 (1991).
    [Crossref]
  5. J. W. McLean, “Modeling of ocean wave effects for LIDAR remote sensing,” Proc. SPIE 1302, 480–491 (1990).
    [Crossref]
  6. J. W. McLean and J. D. Freeman, “Effects of ocean waves on airborne lidar imaging,” Appl. Opt. 35, 3261–3269 (1996).
    [Crossref]
  7. S. Kay, J. Hedley, S. Lavender, and A. Nimmo-Smith, “Light transfer at the ocean surface modeled using high resolution sea surface realizations,” Opt. Express 19, 6493–6504 (2011).
    [Crossref]
  8. H. M. Tulldahl and K. O. Steinvall, “Simulation of sea surface wave influence on small target detection with airborne laser depth sounding,” Appl. Opt. 43, 2462–2483 (2004).
    [Crossref]
  9. G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphere-ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. 34, 1453–1472 (1989).
    [Crossref]
  10. H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,” Appl. Opt. 40, 400–412 (2001).
    [Crossref]
  11. M. Chami, R. Santer, and E. Dilligeard, “Radiative transfer model for the computation of radiance and polarization in an ocean-atmosphere system: polarization properties of suspended matter for remote sensing,” Appl. Opt. 40, 2398–2416 (2001).
    [Crossref]
  12. P.-W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17, 2057–2079 (2009).
    [Crossref]
  13. A. Tonizzo, A. Gilerson, T. Harmel, A. Ibrahim, J. Chowdhary, B. Gross, F. Moshary, and S. Ahmed, “Estimating particle composition and size distribution from polarized water-leaving radiance,” Appl. Opt. 50, 5047–5058 (2011).
    [Crossref]
  14. J. Chowdhary, B. Cairns, and L. D. Travis, “Contribution of water-leaving radiances to multiangle, multispectral polarimetric observations over the open ocean: bio-optical model results for case 1 waters,” Appl. Opt. 45, 5542–5567 (2006).
    [Crossref]
  15. P.-W. Zhai, G. W. Kattawar, and P. Yang, “Impulse response solution of the three-dimensional vector radiative transfer equation in atmosphere-ocean systems. I. Monte Carlo method,” Appl. Opt. 47, 1037–1047 (2008).
    [Crossref]
  16. P.-W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Luker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010).
  17. T. Harmel, A. Gilerson, A. Tonizzo, J. Chowdhary, A. Weidemann, R. Arnone, and S. Ahmed, “Polarization impacts on the water-leaving radiance retrieval from above-water radiometric measurements,” Appl. Opt. 51, 8324–8340 (2012).
    [Crossref]
  18. C. D. Mobley, “Estimation of the remote-sensing reflectance from above-surface measurements,” Appl. Opt. 38, 7442–7455 (1999).
    [Crossref]
  19. M. L. Banner, “Equlibrium spectra of wind waves,” J. Phys. Oceanogr. 20, 966–984 (1990).
    [Crossref]
  20. T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, “A unified directional spectrum for long and short wind-driven waves,” J. Geophys. Res. 102, 15781–15796 (1997).
    [Crossref]
  21. S. R. Massel, “On the geometry of ocean surface waves,” Oceanologia 53, 521–548 (2011).
  22. C. D. Mobley, “Modeling sea surfaces: a tutorial on Fourier transform techniques,” (Sequoia Scientific, Inc., 2014).
  23. J. Tessendorf, “Simulating ocean water,” (SIGGRAPH Course Notes, 2004).
  24. M. S. Longuet-Higgins, D. E. Cartwright, and N. D. Smith, “Observations of the directional spectrum of sea waves using the motions of a flotation buoy,” in Ocean Wave Spectra (Prentice-Hall, 1963), pp. 111–136.
  25. C. D. Mobley, Light and Water: Radiative Transfer in Natural Waters (Academic, 1994).
  26. C. D. Mobley, “HydroPol mathematical documentation: invariant imbedding theory for the vector radiative transfer equation,” (Sequoia Scientific, 2014).
  27. M. I. Mischenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).
  28. P.-W. Zhai, G. W. Kattawar, and Y. Hu, “Comment on the transmission matrix for a dielectric interface,” J. Quant. Spectrosc. Radiat. Transfer 113, 1981–1984 (2012).
    [Crossref]
  29. W. W. Gregg and K. L. Carder, “A simple spectral solar irradiance model for cloudless maritime atmospheres,” Limnol. Oceanogr. 35, 1657–1675 (1990).
    [Crossref]
  30. A. W. Harrison and C. A. Coombes, “Angular distribution of clear sky short wavelength radiance,” Solar Energy 40, 57–63 (1988).
    [Crossref]
  31. L. G. Tilstra, N. A. J. Schutgens, and P. Stammes, “Analytical calculation of Stokes parameters Q and U of atmospheric radiation,” (Koninklijk Nederlands Meteorologisch Instituut, 2003).
  32. C. Emde, R. Buras, B. Mayer, and M. Blumthaler, “The impact of aerosols on polarized sky radiance: model development, validation, and applications,” Atmos. Chem. Phys. 10, 383–396 (2010).
    [Crossref]
  33. C. D. Mobley and E. S. Boss, “Improved irradiances for use in ocean heating, primary production, and photo-oxidation calculations,” Appl. Opt. 51, 6549–6560 (2012).
    [Crossref]
  34. C. D. Mobley, “Reflectances,” http://www.oceanopticsbook.info/view/overview_of_optical_oceanography/reflectances (2010).
  35. K. L. Carder and R. G. Steward, “A remote-sensing reflectance model of a red-tide dinoflagellate off West Florida,” Limnol. Oceanogr. 30, 286–298 (1985).
    [Crossref]
  36. Z.-P. Lee, Y.-H. Ahn, C. D. Mobley, and R. Arnone, “Removal of surface-reflected light for the measurement of remote-sensing reflectance from an above-surface platform,” Opt. Express 18, 26313–26324 (2010).
    [Crossref]
  37. C. D. Mobley, “Surface reflectance factors,” http://www.oceanopticsbook.info/view/remote_sensing/level_3/surface_reflectance_factors (2013).
  38. S. B. Hooker, G. Zibordi, J. F. Berthon, and J. W. Brown, “Above-water radiometry in shallow coastal waters,” Appl. Opt. 43, 4254–4268 (2004).
    [Crossref]
  39. G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
    [Crossref]
  40. T. Harmel, A. Gilerson, S. Hlaing, A. Tonizzo, T. Legbandt, A. Weidemann, R. Arnone, and S. Ahmed, “Long Island Sound Coastal Observatory: assessment of above-water radiometric measurement uncertainties using collocated multi and hyperspectral systems,” Appl. Opt. 50, 5842–5860 (2011).
    [Crossref]

2012 (3)

2011 (4)

2010 (3)

Z.-P. Lee, Y.-H. Ahn, C. D. Mobley, and R. Arnone, “Removal of surface-reflected light for the measurement of remote-sensing reflectance from an above-surface platform,” Opt. Express 18, 26313–26324 (2010).
[Crossref]

C. Emde, R. Buras, B. Mayer, and M. Blumthaler, “The impact of aerosols on polarized sky radiance: model development, validation, and applications,” Atmos. Chem. Phys. 10, 383–396 (2010).
[Crossref]

P.-W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Luker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010).

2009 (2)

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

P.-W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17, 2057–2079 (2009).
[Crossref]

2008 (1)

2006 (1)

2004 (2)

2001 (2)

1999 (1)

1997 (1)

T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, “A unified directional spectrum for long and short wind-driven waves,” J. Geophys. Res. 102, 15781–15796 (1997).
[Crossref]

1996 (1)

1991 (1)

M. P. Levesque, “Dynamic sea image generation,” Proc. SPIE 1486, 294–300 (1991).
[Crossref]

1990 (4)

J. W. McLean, “Modeling of ocean wave effects for LIDAR remote sensing,” Proc. SPIE 1302, 480–491 (1990).
[Crossref]

M. P. Levesque and D. St-Germain, “Generation of synthetic IR sea images,” Proc. SPIE 1331, 354–357 (1990).

M. L. Banner, “Equlibrium spectra of wind waves,” J. Phys. Oceanogr. 20, 966–984 (1990).
[Crossref]

W. W. Gregg and K. L. Carder, “A simple spectral solar irradiance model for cloudless maritime atmospheres,” Limnol. Oceanogr. 35, 1657–1675 (1990).
[Crossref]

1989 (1)

G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphere-ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. 34, 1453–1472 (1989).
[Crossref]

1988 (1)

A. W. Harrison and C. A. Coombes, “Angular distribution of clear sky short wavelength radiance,” Solar Energy 40, 57–63 (1988).
[Crossref]

1986 (1)

R. W. Preisendorfer and C. D. Mobley, “Albedos and glitter patterns of a wind-roughened sea surface,” J. Phys. Oceanogr. 16, 1293–1316 (1986).

1985 (1)

K. L. Carder and R. G. Steward, “A remote-sensing reflectance model of a red-tide dinoflagellate off West Florida,” Limnol. Oceanogr. 30, 286–298 (1985).
[Crossref]

1954 (1)

Adams, C. N.

G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphere-ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. 34, 1453–1472 (1989).
[Crossref]

Ahmed, S.

Ahn, Y.-H.

Arnone, R.

Banner, M. L.

M. L. Banner, “Equlibrium spectra of wind waves,” J. Phys. Oceanogr. 20, 966–984 (1990).
[Crossref]

Berthon, J. F.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

S. B. Hooker, G. Zibordi, J. F. Berthon, and J. W. Brown, “Above-water radiometry in shallow coastal waters,” Appl. Opt. 43, 4254–4268 (2004).
[Crossref]

Blumthaler, M.

C. Emde, R. Buras, B. Mayer, and M. Blumthaler, “The impact of aerosols on polarized sky radiance: model development, validation, and applications,” Atmos. Chem. Phys. 10, 383–396 (2010).
[Crossref]

Boss, E. S.

Brown, J. W.

Buras, R.

C. Emde, R. Buras, B. Mayer, and M. Blumthaler, “The impact of aerosols on polarized sky radiance: model development, validation, and applications,” Atmos. Chem. Phys. 10, 383–396 (2010).
[Crossref]

Cairns, B.

Carder, K. L.

W. W. Gregg and K. L. Carder, “A simple spectral solar irradiance model for cloudless maritime atmospheres,” Limnol. Oceanogr. 35, 1657–1675 (1990).
[Crossref]

K. L. Carder and R. G. Steward, “A remote-sensing reflectance model of a red-tide dinoflagellate off West Florida,” Limnol. Oceanogr. 30, 286–298 (1985).
[Crossref]

Cartwright, D. E.

M. S. Longuet-Higgins, D. E. Cartwright, and N. D. Smith, “Observations of the directional spectrum of sea waves using the motions of a flotation buoy,” in Ocean Wave Spectra (Prentice-Hall, 1963), pp. 111–136.

Chaikovskaya, L. I.

Chami, M.

Chapron, B.

T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, “A unified directional spectrum for long and short wind-driven waves,” J. Geophys. Res. 102, 15781–15796 (1997).
[Crossref]

Chowdhary, J.

Coombes, C. A.

A. W. Harrison and C. A. Coombes, “Angular distribution of clear sky short wavelength radiance,” Solar Energy 40, 57–63 (1988).
[Crossref]

Cox, C.

D’Alimonte, D.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Dilligeard, E.

Elfouhaily, T.

T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, “A unified directional spectrum for long and short wind-driven waves,” J. Geophys. Res. 102, 15781–15796 (1997).
[Crossref]

Emde, C.

C. Emde, R. Buras, B. Mayer, and M. Blumthaler, “The impact of aerosols on polarized sky radiance: model development, validation, and applications,” Atmos. Chem. Phys. 10, 383–396 (2010).
[Crossref]

Fabbri, B.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Feng, H.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Freeman, J. D.

Gilerson, A.

Giles, D.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Gregg, W. W.

W. W. Gregg and K. L. Carder, “A simple spectral solar irradiance model for cloudless maritime atmospheres,” Limnol. Oceanogr. 35, 1657–1675 (1990).
[Crossref]

Gross, B.

Harmel, T.

Harrison, A. W.

A. W. Harrison and C. A. Coombes, “Angular distribution of clear sky short wavelength radiance,” Solar Energy 40, 57–63 (1988).
[Crossref]

Hedley, J.

Hlaing, S.

Holben, B. N.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Hooker, S. B.

Hu, Y.

P.-W. Zhai, G. W. Kattawar, and Y. Hu, “Comment on the transmission matrix for a dielectric interface,” J. Quant. Spectrosc. Radiat. Transfer 113, 1981–1984 (2012).
[Crossref]

P.-W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Luker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010).

P.-W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17, 2057–2079 (2009).
[Crossref]

Ibrahim, A.

Josset, D. B.

P.-W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Luker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010).

Kaitala, S.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Katsaros, K.

T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, “A unified directional spectrum for long and short wind-driven waves,” J. Geophys. Res. 102, 15781–15796 (1997).
[Crossref]

Katsev, I. L.

Kattawar, G. W.

P.-W. Zhai, G. W. Kattawar, and Y. Hu, “Comment on the transmission matrix for a dielectric interface,” J. Quant. Spectrosc. Radiat. Transfer 113, 1981–1984 (2012).
[Crossref]

P.-W. Zhai, G. W. Kattawar, and P. Yang, “Impulse response solution of the three-dimensional vector radiative transfer equation in atmosphere-ocean systems. I. Monte Carlo method,” Appl. Opt. 47, 1037–1047 (2008).
[Crossref]

H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,” Appl. Opt. 40, 400–412 (2001).
[Crossref]

G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphere-ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. 34, 1453–1472 (1989).
[Crossref]

Kay, S.

Lacis, A. A.

M. I. Mischenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

Lavender, S.

Lee, Z.-P.

Legbandt, T.

Levesque, M. P.

M. P. Levesque, “Dynamic sea image generation,” Proc. SPIE 1486, 294–300 (1991).
[Crossref]

M. P. Levesque and D. St-Germain, “Generation of synthetic IR sea images,” Proc. SPIE 1331, 354–357 (1990).

Longuet-Higgins, M. S.

M. S. Longuet-Higgins, D. E. Cartwright, and N. D. Smith, “Observations of the directional spectrum of sea waves using the motions of a flotation buoy,” in Ocean Wave Spectra (Prentice-Hall, 1963), pp. 111–136.

Lucker, P. L.

Luker, P. L.

P.-W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Luker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010).

Massel, S. R.

S. R. Massel, “On the geometry of ocean surface waves,” Oceanologia 53, 521–548 (2011).

Mayer, B.

C. Emde, R. Buras, B. Mayer, and M. Blumthaler, “The impact of aerosols on polarized sky radiance: model development, validation, and applications,” Atmos. Chem. Phys. 10, 383–396 (2010).
[Crossref]

McLean, J. W.

J. W. McLean and J. D. Freeman, “Effects of ocean waves on airborne lidar imaging,” Appl. Opt. 35, 3261–3269 (1996).
[Crossref]

J. W. McLean, “Modeling of ocean wave effects for LIDAR remote sensing,” Proc. SPIE 1302, 480–491 (1990).
[Crossref]

Mélin, F.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Mischenko, M. I.

M. I. Mischenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

Mobley, C. D.

C. D. Mobley and E. S. Boss, “Improved irradiances for use in ocean heating, primary production, and photo-oxidation calculations,” Appl. Opt. 51, 6549–6560 (2012).
[Crossref]

Z.-P. Lee, Y.-H. Ahn, C. D. Mobley, and R. Arnone, “Removal of surface-reflected light for the measurement of remote-sensing reflectance from an above-surface platform,” Opt. Express 18, 26313–26324 (2010).
[Crossref]

C. D. Mobley, “Estimation of the remote-sensing reflectance from above-surface measurements,” Appl. Opt. 38, 7442–7455 (1999).
[Crossref]

R. W. Preisendorfer and C. D. Mobley, “Albedos and glitter patterns of a wind-roughened sea surface,” J. Phys. Oceanogr. 16, 1293–1316 (1986).

C. D. Mobley, “Modeling sea surfaces: a tutorial on Fourier transform techniques,” (Sequoia Scientific, Inc., 2014).

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

C. D. Mobley, “HydroPol mathematical documentation: invariant imbedding theory for the vector radiative transfer equation,” (Sequoia Scientific, 2014).

Moshary, F.

Munk, W.

Nimmo-Smith, A.

Preisendorfer, R. W.

R. W. Preisendorfer and C. D. Mobley, “Albedos and glitter patterns of a wind-roughened sea surface,” J. Phys. Oceanogr. 16, 1293–1316 (1986).

Prikhach, A. S.

Santer, R.

Schuster, G.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Schutgens, N. A. J.

L. G. Tilstra, N. A. J. Schutgens, and P. Stammes, “Analytical calculation of Stokes parameters Q and U of atmospheric radiation,” (Koninklijk Nederlands Meteorologisch Instituut, 2003).

Seppale, J.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Slutsker, I.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

Smith, N. D.

M. S. Longuet-Higgins, D. E. Cartwright, and N. D. Smith, “Observations of the directional spectrum of sea waves using the motions of a flotation buoy,” in Ocean Wave Spectra (Prentice-Hall, 1963), pp. 111–136.

Stammes, P.

L. G. Tilstra, N. A. J. Schutgens, and P. Stammes, “Analytical calculation of Stokes parameters Q and U of atmospheric radiation,” (Koninklijk Nederlands Meteorologisch Instituut, 2003).

Steinvall, K. O.

Steward, R. G.

K. L. Carder and R. G. Steward, “A remote-sensing reflectance model of a red-tide dinoflagellate off West Florida,” Limnol. Oceanogr. 30, 286–298 (1985).
[Crossref]

St-Germain, D.

M. P. Levesque and D. St-Germain, “Generation of synthetic IR sea images,” Proc. SPIE 1331, 354–357 (1990).

Tessendorf, J.

J. Tessendorf, “Simulating ocean water,” (SIGGRAPH Course Notes, 2004).

Tilstra, L. G.

L. G. Tilstra, N. A. J. Schutgens, and P. Stammes, “Analytical calculation of Stokes parameters Q and U of atmospheric radiation,” (Koninklijk Nederlands Meteorologisch Instituut, 2003).

Tonizzo, A.

Travis, L. D.

Trepte, C. R.

P.-W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Luker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010).

P.-W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17, 2057–2079 (2009).
[Crossref]

Tulldahl, H. M.

Tynes, H. H.

Vandemark, D.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, “A unified directional spectrum for long and short wind-driven waves,” J. Geophys. Res. 102, 15781–15796 (1997).
[Crossref]

Weidemann, A.

Yang, P.

Zege, E. P.

Zhai, P.-W.

P.-W. Zhai, G. W. Kattawar, and Y. Hu, “Comment on the transmission matrix for a dielectric interface,” J. Quant. Spectrosc. Radiat. Transfer 113, 1981–1984 (2012).
[Crossref]

P.-W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Luker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010).

P.-W. Zhai, Y. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17, 2057–2079 (2009).
[Crossref]

P.-W. Zhai, G. W. Kattawar, and P. Yang, “Impulse response solution of the three-dimensional vector radiative transfer equation in atmosphere-ocean systems. I. Monte Carlo method,” Appl. Opt. 47, 1037–1047 (2008).
[Crossref]

Zibordi, G.

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

S. B. Hooker, G. Zibordi, J. F. Berthon, and J. W. Brown, “Above-water radiometry in shallow coastal waters,” Appl. Opt. 43, 4254–4268 (2004).
[Crossref]

Appl. Opt. (12)

J. W. McLean and J. D. Freeman, “Effects of ocean waves on airborne lidar imaging,” Appl. Opt. 35, 3261–3269 (1996).
[Crossref]

H. M. Tulldahl and K. O. Steinvall, “Simulation of sea surface wave influence on small target detection with airborne laser depth sounding,” Appl. Opt. 43, 2462–2483 (2004).
[Crossref]

A. Tonizzo, A. Gilerson, T. Harmel, A. Ibrahim, J. Chowdhary, B. Gross, F. Moshary, and S. Ahmed, “Estimating particle composition and size distribution from polarized water-leaving radiance,” Appl. Opt. 50, 5047–5058 (2011).
[Crossref]

J. Chowdhary, B. Cairns, and L. D. Travis, “Contribution of water-leaving radiances to multiangle, multispectral polarimetric observations over the open ocean: bio-optical model results for case 1 waters,” Appl. Opt. 45, 5542–5567 (2006).
[Crossref]

P.-W. Zhai, G. W. Kattawar, and P. Yang, “Impulse response solution of the three-dimensional vector radiative transfer equation in atmosphere-ocean systems. I. Monte Carlo method,” Appl. Opt. 47, 1037–1047 (2008).
[Crossref]

H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,” Appl. Opt. 40, 400–412 (2001).
[Crossref]

M. Chami, R. Santer, and E. Dilligeard, “Radiative transfer model for the computation of radiance and polarization in an ocean-atmosphere system: polarization properties of suspended matter for remote sensing,” Appl. Opt. 40, 2398–2416 (2001).
[Crossref]

T. Harmel, A. Gilerson, A. Tonizzo, J. Chowdhary, A. Weidemann, R. Arnone, and S. Ahmed, “Polarization impacts on the water-leaving radiance retrieval from above-water radiometric measurements,” Appl. Opt. 51, 8324–8340 (2012).
[Crossref]

C. D. Mobley, “Estimation of the remote-sensing reflectance from above-surface measurements,” Appl. Opt. 38, 7442–7455 (1999).
[Crossref]

C. D. Mobley and E. S. Boss, “Improved irradiances for use in ocean heating, primary production, and photo-oxidation calculations,” Appl. Opt. 51, 6549–6560 (2012).
[Crossref]

S. B. Hooker, G. Zibordi, J. F. Berthon, and J. W. Brown, “Above-water radiometry in shallow coastal waters,” Appl. Opt. 43, 4254–4268 (2004).
[Crossref]

T. Harmel, A. Gilerson, S. Hlaing, A. Tonizzo, T. Legbandt, A. Weidemann, R. Arnone, and S. Ahmed, “Long Island Sound Coastal Observatory: assessment of above-water radiometric measurement uncertainties using collocated multi and hyperspectral systems,” Appl. Opt. 50, 5842–5860 (2011).
[Crossref]

Atmos. Chem. Phys. (1)

C. Emde, R. Buras, B. Mayer, and M. Blumthaler, “The impact of aerosols on polarized sky radiance: model development, validation, and applications,” Atmos. Chem. Phys. 10, 383–396 (2010).
[Crossref]

J. Atmos. Ocean. Technol. (1)

G. Zibordi, B. N. Holben, I. Slutsker, D. Giles, D. D’Alimonte, F. Mélin, J. F. Berthon, D. Vandemark, H. Feng, G. Schuster, B. Fabbri, S. Kaitala, and J. Seppale, “AERONET-OC: a network for the validation of ocean color primary radiometric products,” J. Atmos. Ocean. Technol. 26, 1634–1651 (2009).
[Crossref]

J. Geophys. Res. (1)

T. Elfouhaily, B. Chapron, K. Katsaros, and D. Vandemark, “A unified directional spectrum for long and short wind-driven waves,” J. Geophys. Res. 102, 15781–15796 (1997).
[Crossref]

J. Opt. Soc. Am. (1)

J. Phys. Oceanogr. (2)

R. W. Preisendorfer and C. D. Mobley, “Albedos and glitter patterns of a wind-roughened sea surface,” J. Phys. Oceanogr. 16, 1293–1316 (1986).

M. L. Banner, “Equlibrium spectra of wind waves,” J. Phys. Oceanogr. 20, 966–984 (1990).
[Crossref]

J. Quant. Spectrosc. Radiat. Transfer (2)

P.-W. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Luker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010).

P.-W. Zhai, G. W. Kattawar, and Y. Hu, “Comment on the transmission matrix for a dielectric interface,” J. Quant. Spectrosc. Radiat. Transfer 113, 1981–1984 (2012).
[Crossref]

Limnol. Oceanogr. (3)

W. W. Gregg and K. L. Carder, “A simple spectral solar irradiance model for cloudless maritime atmospheres,” Limnol. Oceanogr. 35, 1657–1675 (1990).
[Crossref]

K. L. Carder and R. G. Steward, “A remote-sensing reflectance model of a red-tide dinoflagellate off West Florida,” Limnol. Oceanogr. 30, 286–298 (1985).
[Crossref]

G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphere-ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. 34, 1453–1472 (1989).
[Crossref]

Oceanologia (1)

S. R. Massel, “On the geometry of ocean surface waves,” Oceanologia 53, 521–548 (2011).

Opt. Express (3)

Proc. SPIE (3)

M. P. Levesque and D. St-Germain, “Generation of synthetic IR sea images,” Proc. SPIE 1331, 354–357 (1990).

M. P. Levesque, “Dynamic sea image generation,” Proc. SPIE 1486, 294–300 (1991).
[Crossref]

J. W. McLean, “Modeling of ocean wave effects for LIDAR remote sensing,” Proc. SPIE 1302, 480–491 (1990).
[Crossref]

Solar Energy (1)

A. W. Harrison and C. A. Coombes, “Angular distribution of clear sky short wavelength radiance,” Solar Energy 40, 57–63 (1988).
[Crossref]

Other (9)

L. G. Tilstra, N. A. J. Schutgens, and P. Stammes, “Analytical calculation of Stokes parameters Q and U of atmospheric radiation,” (Koninklijk Nederlands Meteorologisch Instituut, 2003).

C. D. Mobley, “Surface reflectance factors,” http://www.oceanopticsbook.info/view/remote_sensing/level_3/surface_reflectance_factors (2013).

C. D. Mobley, “Reflectances,” http://www.oceanopticsbook.info/view/overview_of_optical_oceanography/reflectances (2010).

C. D. Mobley, “Modeling sea surfaces: a tutorial on Fourier transform techniques,” (Sequoia Scientific, Inc., 2014).

J. Tessendorf, “Simulating ocean water,” (SIGGRAPH Course Notes, 2004).

M. S. Longuet-Higgins, D. E. Cartwright, and N. D. Smith, “Observations of the directional spectrum of sea waves using the motions of a flotation buoy,” in Ocean Wave Spectra (Prentice-Hall, 1963), pp. 111–136.

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

C. D. Mobley, “HydroPol mathematical documentation: invariant imbedding theory for the vector radiative transfer equation,” (Sequoia Scientific, 2014).

M. I. Mischenko, L. D. Travis, and A. A. Lacis, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (19)

Fig. 1.
Fig. 1. Example sampling of elevation and slope spectra. The red dots on the lines in the upper panels show the bounds of the elevation and slope spectra sampled using 1024 points with L = 200 m . The light red lines are the true spectra S , and the truncated heavy blue lines are the adjusted spectra S ˜ obtained from Eqs. (9 )–(11). The dashed lines at 370 rad / m indicate the conventional boundary between gravity and capillary waves. The lower panels show random realizations of short segments along the x direction of the sea surface elevation and slope for the true (light red line) and adjusted (heavy blue line) spectra.
Fig. 2.
Fig. 2. Approach to fully resolved mean square slopes (diamonds, left ordinate) and elevations (circles, right ordinate) as a function of the number of grid points N , for winds speeds of W = 5 (red) and 10 m s 1 (blue).
Fig. 3.
Fig. 3. 2D example sea surface random realization for L x = L y = 100 m and N x = N y = 512 and a wind speed of 10 m s 1 in the + x direction. Light color is high surface elevation; dark color is low elevation. Inset numbers give the significant wave height H 1 / 3 , mean square slopes (total, mss; along wind, mss x ; and cross wind, mss y ), and mean surface tilt angles from the vertical in the along-wind ( θ x ) and cross-wind ( θ y ) directions.
Fig. 4.
Fig. 4. Example mapping of a rectangular FFT grid (light blue lines) to a hexagonal grid of triangles as used in ray tracing. Red dots show the triangle vertices. FFT grid points at the right and top, shown by the dotted lines, are obtained by the inherent spatial periodicity of the surface as determined by FFT techniques. Thus, the surface elevation of point 35 is the same as that of point 27, point 57 is the same as point 1, etc. One of the triangles generated from the FFT grid is shaded in red.
Fig. 5.
Fig. 5. Four cases for determining Stokes vector rotation angles. v 1 is the initial direction, which is rotated into the final direction v 2 by a counterclockwise rotation about direction ξ .
Fig. 6.
Fig. 6. Reflected energy pattern (glitter pattern) for unpolarized air-incident light at a 50 deg incident angle from the zenith and a 10 m s 1 wind speed. The light source is in an otherwise black sky. The viewing direction is looking downward at the surface, facing the glitter pattern. The quad outlined in gray indicates the specular reflection quad that would receive all reflected light for a level surface. The tan color in the lower-right panel indicates values that are identically zero.
Fig. 7.
Fig. 7. Transmitted energy pattern for air-incident light corresponding to Fig. 6. The viewing direction is looking upward at the surface from within the water.
Fig. 8.
Fig. 8. Reflected energy pattern (underwater glitter pattern) for water-incident light at a 50 deg incident angle from the nadir and a 10 m s 1 wind speed. The light source is unpolarized in an otherwise black ocean.
Fig. 9.
Fig. 9. Transmitted energy pattern for water-incident light corresponding to Fig. 8. The viewing direction is looking downward at the surface from the air side.
Fig. 10.
Fig. 10. Dependence of energy reflectance along the midline of the glitter pattern ( φ = 0 ) in Fig. 6, as a function of the number of grid points N x used for surface generation.
Fig. 11.
Fig. 11. Percent of incident rays that undergo multiple interactions with the surface, as a function of wind speed and incident angle from the zenith, for FFT and Cox–Munk surfaces.
Fig. 12.
Fig. 12. Sky radiance distribution for single scattering by atmospheric molecules according to the Rayleigh scattering equations. The upper-left panel shows the total radiance magnitude I relative to 1 in the Sun’s direction. The Sun’s location is at the left side of the plotted hemisphere of sky directions. The other panels show the horizontal versus vertical, ± 45 deg , and degree of total polarization in percent.
Fig. 13.
Fig. 13. Surface-reflected radiance for a single-scattering Rayleigh sky and a wind speed of 10 m s 1 . The Sun zenith angle is 50 deg and the viewing geometry corresponds to Fig. 6.
Fig. 14.
Fig. 14. Surface-transmitted radiance for a single-scattering Rayleigh sky and a wind speed of 10 m s 1 . The Sun zenith angle is 50 deg, and the viewing geometry corresponds to Fig. 7.
Fig. 15.
Fig. 15. Comparison of surface irradiance reflectances R surf as functions of solar zenith angle and wind speed for a polarized sky and various combinations of sea surface model (FFT versus Cox–Munk), ray tracing (polarized versus unpolarized), single versus multiple scattering, Sun azimuthal angle, and wave age.
Fig. 16.
Fig. 16. Radiance reflectance factors ρ ( θ v , ϕ v ) for a level sea surface and a Sun zenith angle of 50 deg. The top panel is for unpolarized ray tracing; the bottom panel is for polarized ray tracing. In both cases, the incident radiance is that of the single-scattering Rayleigh sky shown in Fig. 12. The dashed and dotted lines show the Fresnel reflectance for incident radiance, which is linearly polarized perpendicular ( R ) or parallel ( R ) to the incident meridian plane. The solid line without symbols is the Fresnel reflectance for unpolarized incident radiance.
Fig. 17.
Fig. 17. Radiance reflectance factors ρ computed using Cox–Munk surfaces and unpolarized ray tracing. The curves are for Sun zenith angles θ sun = 30 and 60 deg, azimuthal viewing directions of ϕ v = 90 and 135 deg, and wind speeds of 2, 5, 10, and 15 m s 1 ; θ v is the off-nadir viewing direction.
Fig. 18.
Fig. 18. Radiance reflectance factors ρ computed using FFT surfaces and polarized ray tracing. The solid curves correspond to the viewing and wave conditions of Fig. 17. These curves are for a fully developed sea and the Sun’s incident rays parallel to the wind speed ( ϕ sun = 0 ). The dotted curves and diamond symbols show the curves for θ sun = 30 and 60 deg, azimuthal viewing direction of ϕ v = 135 , but with the Sun’s azimuthal angle at ϕ sun = 90 , so that the incoming rays are perpendicular to the wind direction. The dashed curves and box symbols show the cases of θ sun = 30 and 60 deg azimuthal viewing direction of ϕ v = 135 , ϕ sun = 0 , for a very young sea with Ω c = 5 .
Fig. 19.
Fig. 19. ρ for a Rayleigh sky and a wind speed of 10 m s 1 . The upper-left panel shows ρ as a function of off-nadir and azimuthal viewing directions (relative to the Sun’s azimuth at ϕ v = 0 ) for a single Sun zenith angle of θ sun = 40 deg , computed using FFT surfaces and polarized ray tracing. The upper-right panel shows ρ as a function of the Sun’s zenith angle (radial distance in the plot) and viewing azimuth for a fixed off-nadir viewing direction of θ v = 40 deg . The black dots indicate the ρ ( θ v , ϕ v ) = ( 40 , 135 ) viewing direction. The white regions indicate the specular viewing direction for θ sun = 40 deg . The lower two panels show the ratios of ρ computed for FFT surfaces and polarized ray tracing to the values computed for Cox–Munk surfaces and unpolarized ray tracing, corresponding to the respective left and right upper panels.

Tables (2)

Tables Icon

Table 1. Mean Square Slopes as Given by Cox–Munk Formulas in Eq. (1) and by Random Realizations of FFT Surfaces a

Tables Icon

Table 2. Differences in Surface Irradiance Reflectance R surf as Computed in Various Ways a

Equations (41)

Equations on this page are rendered with MathJax. Learn more.

σ a 2 = 0.003 + 0.00316 W 12.5 ± 0.002 ,
σ c 2 = 0.00192 W 12.5 ± 0.004 ,
σ 2 = 0.003 + 0.00512 W 12.5 ± 0.004 ,
var { z } = z 2 = Ψ ( k x , k y ) d k x d k y ,
Ψ 1 S ( k , φ ) = 1 k S ( k ) Φ ( k , φ ) .
Φ ( k , φ ) = 1 2 π [ 1 + Δ ( k ) cos ( 2 φ ) ] .
σ a 2 = k x 2 Ψ 2 S ( k x , k y ) d k x d k y = 0 π π k 2 cos 2 φ Ψ 2 S ( k , φ ) k d k d φ .
σ 2 = σ a 2 + σ c 2 = 0 k 2 S ( k ) d k .
z ^ o ( k uv ) 1 2 [ ρ ( k uv ) + i σ ( k uv ) ] × [ Ψ 1 S ( k = k uv ) 2 Δ k x Δ k y ] 1 / 2 ,
= 1 2 [ ρ ( k uv ) + i σ ( k uv ) ] Ψ 2 S ( k uv ) .
z ^ o ( u , v ) z ^ o * ( u , v ) = { 1 2 [ ρ ( u , v ) + i σ ( u , v ) ] Ψ 2 S ( u , v ) } × { 1 2 [ ρ ( u , v ) i σ ( u , v ) ] Ψ 2 S ( u , v ) } = Ψ 2 S ( u , v ) 2 [ ρ 2 ( u , v ) + σ 2 ( u , v ) ] = Ψ 2 S ( u , v ) ,
z ^ ( k uv , t ) 1 2 [ z ^ o ( + k uv ) exp ( i ω uv t ) + z ^ o * ( k uv ) exp ( + i ω uv t ) ]
z 2 = k 0 k h S ( k ) d k ,
σ 2 = k 0 k h k 2 S ( k ) d k .
f E z 2 ( N ) z 2 = 0.4219 0.4296 = 0.982 ,
S ˜ ( k ) [ 1 + δ ( k ) ] S ( k ) ,
δ ( k ) { 0 if k k p δ Ny ( k k p k Ny k p ) if k > k p ,
σ 2 k 0 k h k 2 S ( k ) d k = k 0 k Ny k 2 S ( k ) d k + k Ny k h k 2 S ( k ) d k k 0 k Ny k 2 S ˜ ( k ) d k = k 0 k Ny k 2 S ( k ) d k + δ Ny k p k Ny k 2 ( k k p k Ny k p ) S ( k ) d k .
δ Ny = k Ny k h k 2 S ( k ) d k k p k Ny k 2 ( k k p k Ny k p ) S ( k ) d k .
h i = z × ξ i | z × ξ i | , v i = ξ i × h i , ξ i = h i × v i .
s = ξ i × n | ξ i × n | , p = ξ i × s , ξ i = s × p .
α i = | cos 1 ( v 1 · v 2 ) | = | cos 1 ( h i · s ) | .
α i = 2 π | cos 1 ( v 1 · v 2 ) | = 2 π | cos 1 ( h i · s ) | .
R ̲ ( α ) = [ 1 0 0 0 0 cos 2 α sin 2 α 0 0 sin 2 α cos 2 α 0 0 0 0 1 ] .
S ̲ f = R ̲ ( α m , f ) M ̲ ( ψ m ) R ̲ ( α m 1 , m ) R ̲ ( α 1 , 2 ) M ̲ ( ψ 1 ) R ̲ ( α i , 1 ) S ̲ i ,
S ̲ f E ̲ ( ξ i ξ f ) S ̲ i ,
R ̲ ( Q i j Q k l ) = E ̲ ( Q i j Q k l ) | μ i | Ω i j | μ k | Ω k l .
E ̲ = [ X X X 0 X X X 0 X X X 0 0 0 0 X ] ,
E ̲ = [ X X X 0 X X X X X X X X 0 X X X ] ,
S ̲ refl ( Q k l ) = i j R ̲ aw ( Q i j Q k l ) S ̲ sky ( Q i j ) ,
E d ( sky ) = i j I sky ( Q i j ) μ i Ω i j .
R ̲ ( 40 , ϕ , 40 , ϕ ) = [ 2.566 × 10 2 1.963 × 10 2 3.521 × 10 7 0.0 1.963 × 10 2 2.566 × 10 2 7.489 × 10 7 0.0 3.521 × 10 7 7.489 × 10 7 1.616 × 10 2 0.0 0.0 0.0 0.0 1.616 × 10 2 ] .
S ̲ sky ( 40 , 90 ) = [ 4.931 × 10 2 2.403 × 10 2 1.583 × 10 2 0.0 ] ,
S ̲ sky ( 40 , 135 ) = [ 3.932 × 10 2 1.418 × 10 2 3.262 × 10 2 0.0 ] .
S ̲ sr ( 40 , 90 ) = [ 9.547 × 10 4 5.617 × 10 4 3.883 × 10 4 0.0 ]
S ̲ sr ( 40 , 135 ) = [ 1.287 × 10 3 1.136 × 10 3 5.269 × 10 4 0.0 ] .
ρ ( 40 , 90 ) = I sr I sky = 9.547 × 10 4 4.931 × 10 2 = 0.0194 ρ ( 40 , 135 ) = I sr I sky = 1.287 × 10 3 3.932 × 10 2 = 0.0327 ,
S ̲ sr ( 40 , 90 ) = [ 1.265 × 10 3 , 0 , 0 , 0 ] T
S ̲ sr ( 40 , 135 ) = [ 1.009 × 10 3 , 0 , 0 , 0 ] T .
ρ ( 40 , 90 ) = 1.265 × 10 3 4.931 × 10 2 = 0.0257
ρ ( 40 , 135 ) = 1.009 × 10 3 3.932 × 10 2 = 0.0257 ,

Metrics