## 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 ${10}^{5}$.

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]:

relate the along wind (${\sigma}_{a}^{2}$), cross wind (${\sigma}_{c}^{2}$), and total (${\sigma}^{2}={\sigma}_{a}^{2}+{\sigma}_{c}^{2}$) variances of the nondimensional sea surface slopes to the wind speed ${W}_{12.5}$ in meters per second at 12.5 m above sea level (the coefficients of ${W}_{12.5}$ have units of $\mathrm{s}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}$). 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–$14\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$. 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 [3
–7] 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\times 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 [9 –13] or Cox–Munk surfaces [14 –17]. 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(\mathbf{x})=z(x,y)$ be the sea surface elevation in meters at point $\mathbf{x}=(x,y)$ at a particular time. The spatial extent of the sea surface is $0\le x<{L}_{x}$ and $0\le y<{L}_{y}$, with ${L}_{x}$ and ${L}_{y}$ in meters. A wind-centered coordinate system is used, with the $+x$ direction chosen to be downwind; $-x$ is then upwind, and $\pm y$ are the cross-wind directions. This surface is sampled on a rectangular grid of ${N}_{x}$ by ${N}_{y}$ points, where both ${N}_{x}$ and ${N}_{y}$ are powers of 2 to enable use of FFTs. The spatial sampling points are then $x(r)=r\mathrm{\Delta}x,r=0,\dots ,{N}_{x}-1$ with $\mathrm{\Delta}x={L}_{x}/{N}_{x}$. A similar equation holds for the $y(s)=s\mathrm{\Delta}y$ points. A discrete sample of $z(x,y)$ is denoted $z({\mathbf{x}}_{\mathrm{rs}})$ 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 $\widehat{z}(u,v)={\mathrm{FFT}}_{2D}\{z(r,s)\}$. These amplitudes are defined at positive and negative discrete angular spatial frequencies ${k}_{x}(u)=u\mathrm{\Delta}{k}_{x},u=-({N}_{x}/2+1),\dots ,{N}_{x}/2$, with $\mathrm{\Delta}{k}_{x}=2\pi /{L}_{x}$. The fundamental frequency of magnitude $\mathrm{\Delta}{k}_{x}$ corresponds to the longest resolvable wave. The positive frequency ${k}_{x}({N}_{x}/2)=2\pi /(2\mathrm{\Delta}x)$ is the Nyquist frequency, which corresponds to the shortest resolvable wavelength of $2\mathrm{\Delta}x$. Corresponding results hold for the ${k}_{y}(v)$ spatial frequencies. These spatial frequencies have units of radians per meter. A discrete frequency pair ${\mathbf{k}}_{\mathrm{uv}}=[{k}_{x}(u),{k}_{y}(v)]$ is labeled by $(u,v)$. The sea surface $z({\mathbf{x}}_{\mathrm{rs}})$ is real valued, so the amplitudes are Hermitian: ${\widehat{z}}^{*}(-{\mathbf{k}}_{\mathrm{uv}})=\widehat{z}(+{\mathbf{k}}_{\mathrm{uv}})$, where ${\widehat{z}}^{*}$ denotes the complex conjugate of $\widehat{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 $\mathrm{\Psi}(u,v)={|\widehat{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: $\mathrm{\Psi}({k}_{x},{k}_{y})=\mathrm{\Psi}(u,v)/(\mathrm{\Delta}{k}_{x}\mathrm{\Delta}{k}_{y})$. 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 ${z}^{2}(r,s)$, ${|\widehat{z}(u,v)|}^{2}$, and $\mathrm{\Psi}(u,v)$ are all point functions with units of ${\mathrm{m}}^{2}$, whereas the continuous function $\mathrm{\Psi}({k}_{x},{k}_{y})$ is a spectral density with units of ${\mathrm{m}}^{2}/(\mathrm{rad}/\mathrm{m}{)}^{2}$. The integral of $\mathrm{\Psi}({k}_{x},{k}_{y})$ over all spatial frequencies gives the variance of the zero-mean sea surface elevation:

The FFT of $z(r,s)$ leads to a discrete two-sided (2S) spectrum, denoted here by ${\mathrm{\Psi}}_{2S}({\mathbf{k}}_{\mathrm{uv}})$. “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 ${k}_{x}$ value correspond to waves propagating generally downwind ($|\phi |<\pi /2$, where $\phi =\mathrm{tan}({k}_{y}/{k}_{x})$); negative ${k}_{x}$ corresponds to waves propagating upwind ($|\phi |>\pi /2$).

The resulting ${\mathrm{\Psi}}_{2S}(\mathbf{k})={\mathrm{\Psi}}_{2S}({k}_{x},{k}_{y})$ is equal for directions $+\mathbf{k}$ and $-\mathbf{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=|\mathbf{k}|={({k}_{x}^{2}+{k}_{y}^{2})}^{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: ${\mathrm{\Psi}}_{1S}(k)=2{\mathrm{\Psi}}_{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, ${\mathrm{\Psi}}_{2S}(\mathbf{k})\gg {\mathrm{\Psi}}_{2S}(-\mathbf{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,\phi )$ as

The $\mathrm{\Psi}$ spectra give the variance of the sea surface elevations. The corresponding variance spectrum of the sea surface slopes is ${k}^{2}\mathrm{\Psi}$. Thus, the variance of the wave slopes (the mean square slope) in the along-wind ($\pm x$) direction is given by either of

#### 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

The amplitudes ${\widehat{z}}_{o}({\mathbf{k}}_{\mathrm{uv}})={\widehat{z}}_{o}(u,v)$ are random variables. The statistical expectation of ${|{\widehat{z}}_{o}({\mathbf{k}}_{\mathrm{uv}})|}^{2}$ is

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 ${\omega}_{\mathrm{uv}}$ from the spatial angular frequency ${k}_{\mathrm{uv}}$. For deep-water gravity waves, the dispersion relation is ${\omega}_{\mathrm{uv}}^{2}=g{k}_{\mathrm{uv}}$, 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 $\mathrm{\Phi}(k,\phi )=C(s){\mathrm{cos}}^{2s}(\phi /2)$, where $s$ is a spreading parameter that depends on $k$, and $C$ is a normalization coefficient. These functions are asymmetric about $\phi =\pi /2$ and give much stronger propagation downwind than upwind. In this case, essentially all of the variance is contained in the downwind half of ${\mathrm{\Psi}}_{2S}(\mathbf{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)={\mathrm{FFT}}_{2D}^{-1}\{\widehat{z}({\mathbf{k}}_{\mathrm{uv}},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/\sqrt{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 ${\mathrm{\Psi}}_{1S}$ is not divided by 2 in Eq. (4) and the $1/\sqrt{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 ${L}_{x}$, ${L}_{y}$, ${N}_{x}$, and ${N}_{y}$ 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, ${L}_{x}$ and ${L}_{y}$ can be chosen large enough to resolve the gravity waves with the greatest elevation variance (estimated from the peak of the $\mathcal{S}(k)$ spectrum). The number of grid points can be small; ${2}^{9}=512$ or ${2}^{10}=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 $200\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\times 200\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ surfaces with ${2}^{16}=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 $\mathcal{S}(k)$ of Elfouhaily *et al.* [20] for a wind speed of $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$. The upper-right panel shows the corresponding slope spectrum ${k}^{2}\mathcal{S}(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 ${k}_{0}$ to a very high frequency ${k}_{h}$, which cover the frequency range where the spectra are non-negligible in magnitude:

Now suppose this spectrum is used to generate a surface with ${L}_{x}=200\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ using ${N}_{x}=1024$ sampling points. The fundamental spatial frequency corresponding to the longest resolvable wavelength is then ${k}_{f}=2\pi /{L}_{x}=0.0314\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{rad}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}$ and the Nyquist frequency is ${k}_{\mathrm{Ny}}=16.085\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{rad}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}$; these frequencies are shown by the two red dots in the upper panels of Fig. 1. These two frequencies and the ${N}_{x}-2$ evenly spaced frequencies in between are the frequencies at which the elevation variance spectrum is sampled. Using ${k}_{f}$ and ${k}_{\mathrm{Ny}}$ as the lower and upper limits in Eqs. (7) and (8) gives $\u3008{z}^{2}\u3009(N=1024)=0.4219\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{m}}^{2}$ and ${\sigma}^{2}(N=1024)=0.02584\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\mathrm{rad}}^{2}$. Thus, the finite range of the sampled frequencies includes the fraction,

of the total variance of the sea surface elevation. However, the corresponding fraction of the sampled slope variance is just ${f}_{S}=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,

such that the integral of ${k}^{2}\tilde{\mathcal{S}}(k)$ over the sampled region ${k}_{f}$ to ${k}_{\mathrm{Ny}}$ equals the integral of the true ${k}^{2}\mathcal{S}(k)$ over the entire spectral range. Sampling the rescaled spectrum ${k}^{2}\tilde{\mathcal{S}}(k)$ over the ${k}_{f}$ to ${k}_{\mathrm{Ny}}$ frequency range then gives the same slope variance, as would be obtained by sampling the true spectrum ${k}^{2}\mathcal{S}(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 ${k}_{f}$. 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 $\delta (k)$ to be a linear function of $k$ between the spectrum peak ${k}_{p}$ and the highest sampled frequency ${k}_{\mathrm{Ny}}$, and zero elsewhere:

The parameter ${\delta}_{\mathrm{Ny}}$ is determined by noting that

The right-most terms of these equations give (after recalling that $\delta (k)=0$ for $k\le {k}_{p}$)

The heavy blue lines in the upper panels of Fig. 1 show the $\tilde{\mathcal{S}}(k)$ and ${k}^{2}\tilde{\mathcal{S}}(k)$ spectra for the present example. It is clear that the $\delta (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 ${f}_{E}$ value shows, this rescaling has increased the fraction of sampled/true variance from 0.982 to 1.020. However, the upper ${f}_{S}$ 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({x}_{r})$ 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 $\delta (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\times L=100\times 100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ area of sea surface is simulated using an $N\times N$ FFT, for a wind speed of $W=5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$. The open red diamonds show the corresponding mean square slopes for the rescaled spectra. For $N={2}^{7}=128$ points, the mss is only 0.0119, versus a value of $\text{mss}=0.0328$ for the rescaled spectrum. Not until $N={2}^{15}=\mathrm{32,768}$ does the mss for an unadjusted spectrum equal the value for the adjusted spectrum. For $L=100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ and $N={2}^{15}$, the surface is being resolved down to a grid spacing of $L/N=0.003\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ (or a Nyquist spatial frequency of ${k}_{\mathrm{Ny}}=1029\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{rad}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}$), which fully resolves the capillary waves. At this value, the rescaled and unscaled spectra are the same, i.e., $\delta =0$ in Eq. (9). Larger values of $N$ then make no difference in the mss values. The red dots show the elevation variances ${z}^{2}$ for unadjusted (solid dots) and adjusted (open dots) spectra. For $L=100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ and $W=5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$, there is almost no difference in the elevation variances. The blue curves show the corresponding results for $L=300\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ and $W=10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$. In that case, the unadjusted spectrum does not fully resolve the mss until $N={2}^{17}=\mathrm{131,072}$ points are used, which corresponds to a grid point spacing of 0.0023 m or a Nyquist frequency of $1373\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{rad}\text{\hspace{0.17em}}{\mathrm{m}}^{-1}$. Now, however, the small $N$ values give a substantial increase to the elevation variance, as seen by the open blue circles. For $N=128$, ${z}^{2}$ 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.

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 ${H}_{1/3}=2.18\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ for this realization with ${L}_{x}=100\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ is somewhat less than the theoretical value of ${H}_{1/3}\equiv 4\sqrt{\u3008{z}^{2}\u3009}=2.622\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ obtained by numerical integration of the elevation variance spectrum, as in Eq. (2). However, increasing ${L}_{x}$ to 400 m accounts for the very longest waves and gives ${H}_{1/3}=2.613\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$. 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 ${\sigma}_{a}^{2}/{\sigma}_{c}^{2}$) depend on choice of spreading function $\mathrm{\Phi}(k,\phi )$, but the total variance depends only on the omnidirectional spectrum $\mathcal{S}(k)$.

## 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 ${N}_{y}={N}_{x}/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 ${N}_{x}=16$ and ${N}_{y}=8$ (and ${L}_{x}={L}_{y}=8\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$). This resampling generates a hexagonal region with $3({N}_{y}/2)({N}_{y}/2+1)+1$ facet vertices and $3{N}_{y}^{2}/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].

#### 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 $\underline{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 $+45\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ polarization, and negative $U$ is $-45\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$. 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 $\mathbf{z}$ axis (outward normal to the mean sea surface) and the direction ${\mathit{\xi}}_{i}$ of the incident ray. The associated Stokes vector ${\underline{S}}_{i}$ is specified relative to horizontal and vertical directions, ${\mathbf{h}}_{i}$ and ${\mathbf{v}}_{i}$, which are, respectively, perpendicular and parallel to the incident meridian plane. Unit vectors ${\mathbf{h}}_{i}$ and ${\mathbf{v}}_{i}$ are determined so that “perpendicular” cross “parallel” is in the direction of light propagation. For an incident ray, this rule gives

An incident ray always intersects a wave facet. The incident direction ${\mathit{\xi}}_{i}$ and the wave facet normal $\mathbf{n}$ define the scattering plane containing the incident, reflected (${\mathit{\xi}}_{r}$), and transmitted (${\mathit{\xi}}_{t}$) rays. The incident Stokes vector must be “rotated” (transformed) from the incident meridian plane into a coordinate system whose axes are perpendicular ($\mathbf{s}$) and parallel ($\mathbf{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

The rotation angle $0\le {\alpha}_{i}<2\pi $ that rotates the incident Stokes vector into the scattering plane is the angle that rotates the initial perpendicular direction ${\mathbf{h}}_{i}$ into the direction perpendicular to the scattering plane. (This same angle rotates ${\mathbf{v}}_{i}$ into $\mathbf{p}$.) The rotation angle can, therefore, be obtained from the dot product of ${\mathbf{h}}_{i}$ and $\mathbf{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 $-{\mathit{\xi}}_{i}$ direction) [27]. In general, the rotation angle is the angle required to rotate the first vector ${\mathbf{v}}_{1}$ counterclockwise into the second vector ${\mathbf{v}}_{2}$ when looking into the rotation axis, which is vector ${\mathbf{v}}_{3}={\mathbf{v}}_{1}\times {\mathbf{v}}_{2}$. In the present case, ${\mathbf{v}}_{1}={\mathbf{h}}_{i}$, ${\mathbf{v}}_{2}=\mathbf{s}$, and ${\mathbf{v}}_{3}={\mathit{\xi}}_{i}$ is the rotation axis. In general, four cases must be considered, as shown in Fig. 5. The rotation angle depends on whether ${\mathbf{v}}_{1}\xb7{\mathbf{v}}_{2}$ is positive or negative and on whether ${\mathbf{v}}_{1}\times {\mathbf{v}}_{2}$ is parallel or antiparallel to the rotation axis. In the present case, ${\mathbf{v}}_{1}\times {\mathbf{v}}_{2}={\mathbf{h}}_{i}\times \mathbf{s}$ is by construction parallel to ${\mathit{\xi}}_{i}$, in which case

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]

Once a Stokes vector has been rotated into the local coordinate system of a wave facet, the $4\times 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 ${\underline{R}}_{\mathrm{aw}}$), transmission from air to water (${\underline{T}}_{\mathrm{aw}}$), reflection from the water side of a facet (${\underline{R}}_{\mathrm{wa}}$, which may be total internal reflection), and transmission from water to air (${\underline{T}}_{\mathrm{wa}}$). These Fresnel matrices are given in [26,28]. They depend only on the angle of incidence relative to the surface normal, ${\theta}_{i}={\mathrm{cos}}^{-1}|{\mathit{\xi}}_{i}\xb7\mathbf{n}|$, and the index of refraction of the water.

When a ray eventually leaves the region of the sea surface, its final Stokes vector ${\underline{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 ${\mathit{\xi}}_{f}(={\mathit{\xi}}_{r}\text{\hspace{0.17em}}\text{or}\text{\hspace{0.17em}}{\mathit{\xi}}_{t})$ cross $\mathbf{z}$. Then, as above, ${\mathbf{v}}_{f}={\mathit{\xi}}_{f}\times {\mathbf{h}}_{f}$. The final rotation angle ${\alpha}_{f}$, which carries the Stokes vector from the most recent wave facet $\mathbf{s}$ and $\mathbf{p}$ axes to the final ${\mathbf{h}}_{f}$ and ${\mathbf{v}}_{f}$ 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 ${\underline{S}}_{f}$ as

- • ${\alpha}_{i,1}$ is the rotation angle that rotates the initial ${\underline{S}}_{i}$ into the scattering plane of the first wave facet;
- • $\underline{M}({\psi}_{1})$ is any of the four Fresnel reflectance and transmittance matrices ${\underline{R}}_{\mathrm{aw}}$, ${\underline{R}}_{\mathrm{wa}}$, ${\underline{T}}_{\mathrm{aw}}$, or ${\underline{T}}_{\mathrm{wa}}$ applied to the first ray-wave interaction;
- • ${\psi}_{1}$ represents the scattering angle for the incident and reflected or transmitted directions used in the Fresnel matrix;
- • $\underline{R}({\alpha}_{1,2})$ is the rotation matrix that takes the Stokes vector $\underline{M}({\psi}_{1})\underline{R}({\alpha}_{i}){\underline{S}}_{i}$ from the scattering plane of the first wave facet to the scattering plane of the second wave facet;
- • $\underline{R}({\alpha}_{m-1,m})$ carries the current Stokes vector from wave facet $m-1$ to the last facet $m$;
- • $\underline{M}({\psi}_{m})$ is the Fresnel reflection or transmission function for the last wave facet;
- • $\underline{R}({\alpha}_{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 ${\underline{S}}_{f}$.

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 $(\theta ,\varphi )$ 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 $\mathrm{\Delta}\theta \times \mathrm{\Delta}\varphi =10\times 15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$, except for $\mathrm{\Delta}\theta =5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ 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 $\underline{\mathcal{E}}$ are computed for all pairs of quads. There are four sets of these transfer matrices, e.g., ${\underline{\mathcal{E}}}_{\mathrm{raw}}({Q}_{ij}\to {Q}_{kl})$, which describe reflection for air-incident rays coming from quad ${Q}_{ij}$, corresponding to the $i$th $\theta $ band and $j$th $\theta $ band, and eventually being reflected upward after one or more ray-wave interactions into quad ${Q}_{kl}$. For particular quads ${Q}_{ij}$ and ${Q}_{kl}$, ${\underline{\mathcal{E}}}_{\mathrm{raw}}({Q}_{ij}\to {Q}_{kl})$ is, thus, a $4\times 4$ scattering (Mueller) matrix that shows how any quad-averaged incident Stokes vectors ${\underline{S}}_{i}({Q}_{ij})$ are transformed by reflection into a final Stokes vector ${\underline{S}}_{f}({Q}_{kl})$. This matrix and its counterparts ${\underline{\mathcal{E}}}_{\mathrm{taw}}$, ${\underline{\mathcal{E}}}_{\mathrm{rwa}}$, and ${\underline{\mathcal{E}}}_{\mathrm{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 $\mathfrak{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 $\underline{\mathcal{E}}$ is set to the $4\times 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 $\underline{\mathcal{E}}$ array is tallied to the accumulating array for the pair of quads containing the incident and final ray directions. The averages of the $\underline{\mathcal{E}}$ arrays for the $\mathfrak{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\le \phi \le 90\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$). For the $10\times 15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ 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 $128\mathfrak{S}$ 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 $\mathfrak{S}={10}^{5}$ independent surface realizations, hence $12.8\times {10}^{6}$ initial rays, give Monte Carlo noise of only a few percent in the computed quantities for directions that are energetically important. Generating ${10}^{5}$ surfaces and tracing roughly $36\times {10}^{6}$ rays to completion requires less than 8 h on a 2.5 GHz computer.

The $\underline{\mathcal{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)]:

## 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 6 –9 show the (temporally or spatially averaged) polarized reflectance and transmittance for a fully developed sea surface and a wind speed of $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$. The surfaces were generated with FFTs using rescaled variance spectra to fully resolve both elevation and slope variances. The sea surface region was $200\times 200\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ with ${N}_{x}=1024$ and ${N}_{y}={N}_{x}/2=512$, as required for mapping the rectangular FFT grid to a hexagonal grid of triangular wave facets. Ray tracing was performed for ${10}^{5}$ 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 $\phi =180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ side the plotted hemisphere of directions, so the unscattered rays are traveling in the $\phi =0$ azimuthal direction. The incident Stokes vector in each case is unpolarized light, $\underline{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\times 15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ 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 $\phi $ 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.

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 $\theta =50\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$, which is near the Brewster angle of 53 deg for a level surface. The $\pm 45\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ 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 ${N}_{x}=1024$ was a reasonable choice for numerical simulations. Figure 10 shows the reflected energy as a function of $\theta $ for $\phi =0$ in Fig. 6 when values of ${N}_{x}=512$, 1024 and 2048 are used in the simulations and the elevation variance spectrum is rescaled as previously described. The ${N}_{x}=512$ curve is as much as 25% different from the ${N}_{x}=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, ${N}_{x}=1024$ gives results within a few percent of the values obtained for ${N}_{x}=2048$ but with only one-fifth of the computer time (the ratios of $({N}_{x}{\text{\hspace{0.17em}}\mathrm{log}}_{2}\text{\hspace{0.17em}}{N}_{x})({N}_{y}\text{\hspace{0.17em}}{\mathrm{log}}_{2}\text{\hspace{0.17em}}{N}_{y})$ for the different choices of ${N}_{x}$ and ${N}_{y}={N}_{x}/2$). The red solid, dotted, and dashed curves in that figure are the values for ${N}_{x}=1024$ and three different seeds for random number generation. These curves differ by less than 4%, which shows that simulations with ${10}^{5}$ 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 ${N}_{x}=1024$ and ${10}^{5}$ surface realizations for the simulations in the present work.

Figure 10 also shows the reflectance for an FFT surface with ${N}_{x}=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 $\pm 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 $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$ wind speed. The Cox–Munk and FFT curves for the rescaled ${N}_{x}=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.

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:

## 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 (${E}_{d}^{\text{dir}}$) and diffuse (${E}_{d}^{\text{dif}}$) sky irradiances for the given solar zenith angle ${\theta}_{\text{sun}}$ and a wavelength of 550 nm. The quad containing the Sun is given a radiance of magnitude ${I}_{\text{sun}}={E}_{d}^{\text{dir}}/({\mu}_{\text{sun}}{\mathrm{\Omega}}_{\text{sun}})$, where ${\mu}_{\text{sun}}$ is the average of $\mathrm{cos}\text{\hspace{0.17em}}\theta $ over the Sun’s quad, and ${\mathrm{\Omega}}_{\text{sun}}$ is the solid angle of this quad. The Sun’s direct beam is taken to be unpolarized: ${\underline{S}}_{\text{sun}}={[{I}_{\text{sun}},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 ${E}_{d}^{\text{dif}}$. This gives a calibrated sky $I(\theta ,\varphi )$ 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 $\pm 45\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}(U/I)$ polarization, and the degree of total polarization $\mathrm{DoP}={({Q}^{2}+{U}^{2}+{V}^{2})}^{\frac{1}{2}}/I$. The circular polarization is identically zero. In the numerical calculations below, these equations are used to compute quad-averaged sky Stokes vectors.

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 $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$, generated with ${L}_{x}=200\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}$ and ${N}_{x}=1024$. The increase in the radiance magnitude $I$ going toward the horizon at ${\varphi}_{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 (${\theta}_{v}=87.5$) and ${\varphi}_{v}=\pm 90\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ are 86.69% and 87.79%. Such differences can cause slight asymmetries in the color coding of plots for $\pm {\varphi}_{v}$ values. The $I$ value for $({\theta}_{v},{\varphi}_{v})=(87.5,0)$ (the orange quad in the upper-left panel) was $39.98\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}\text{\hspace{0.17em}}{\mathrm{m}}^{-2}\text{\hspace{0.17em}}{\mathrm{sr}}^{-1}\text{\hspace{0.17em}}{\mathrm{nm}}^{-1}$ 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 ${E}_{d}^{\text{dir}}=0.6561$ and ${E}_{d}^{\text{dif}}=0.3509\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}\text{\hspace{0.17em}}{\mathrm{m}}^{-2}\text{\hspace{0.17em}}{\mathrm{nm}}^{-1}$. The incident sky Stokes vectors are now diffuse vectors with units of spectral radiance, so the surface reflectance function is an ${\underline{\mathcal{R}}}_{\mathrm{aw}}$, as defined in Eq. (15). All sky quads now, in principle, contribute to the reflected radiance in each quad:

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 $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$ wind speed, the rough surface allows radiance to be transmitted into all downwelling directions, although directions near the horizon receive very little energy.

## 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

For the sky used to generate Fig. 13, ${E}_{d}(\mathrm{sky})=1.0070\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}\text{\hspace{0.17em}}{\mathrm{m}}^{-2}\text{\hspace{0.17em}}{\mathrm{nm}}^{-1}$, as expected from the input ${E}_{d}^{\text{dir}}$ and ${E}_{d}^{\text{dif}}$ values. The upwelling, surface-reflected irradiance computed from ${I}_{\text{refl}}$ is ${E}_{u}(\text{refl})=0.05098\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{W}\text{\hspace{0.17em}}{\mathrm{m}}^{-2}\text{\hspace{0.17em}}{\mathrm{nm}}^{-1}$. Thus, for the conditions of Fig. 13, ${R}_{\text{surf}}\equiv {E}_{u}(\text{refl})/{E}_{d}(\mathrm{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 ${R}_{\text{surf}}$ 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 ${E}_{u}$ includes surface-reflected and water-leaving radiance. The albedo is, therefore, always somewhat greater than ${R}_{\text{surf}}$. For the ${10}^{5}$ surface realizations used here, the statistical noise in ${R}_{\text{surf}}$ is never more than 0.2% and usually much less. This variability is too small to display on the curves of Fig. 15.The upper-left panel of Fig. 15 shows the dependence of ${R}_{\text{surf}}$ 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[{R}_{\text{surf}}(\text{polarized})-{R}_{\text{surf}}(\text{unpolarized})]/{R}_{\text{surf}}(\text{polarized})$ for Sun zenith angles of 0, 50, and 87.5 deg. For the Sun at the zenith, the difference in ${R}_{\text{surf}}$ 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 $15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$, 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 ${R}_{\text{surf}}$ 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 ${R}_{\text{surf}}$. A higher fraction of the total radiance, thus, begins to enter the water as the Sun nears the horizon and ${R}_{\text{surf}}$ decreases.

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[{R}_{\text{surf}}(\mathrm{FFT})-{R}_{\text{surf}}(\text{Cox-Munk})]/{R}_{\text{surf}}(\mathrm{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 $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$, the Cox–Munk surface has a greater ${R}_{\text{surf}}$ at large solar zenith angles, and, for wind speeds greater than $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$, the reverse is true. The two surfaces give almost identical ${R}_{\text{surf}}$ for all zenith angles when the wind speed is $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$.

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[{R}_{\text{surf}}(\mathrm{FFT},\text{pol})-{R}_{\text{surf}}(\text{Cox-Munk},\text{unpol})]/{R}_{\text{surf}}(\mathrm{FFT},\text{pol})$. ${R}_{\text{surf}}$ 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 ${R}_{\text{surf}}$ computed for multiple versus single scattering, for polarized ray tracing and FFT surfaces. The tabulated differences are now $100[{R}_{\text{surf}}(\text{multiple})-{R}_{\text{surf}}(\text{single})]/{R}_{\text{surf}}(\text{multiple})$. Multiple scattering always gives a greater ${R}_{\text{surf}}$ 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 ${\varphi}_{\text{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 ${\varphi}_{\text{sun}}=90\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$, for which the Sun’s incident rays are perpendicular to the wind speed. In this case, the surface reflectance is somewhat greater than for ${\varphi}_{\text{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 ${\mathrm{\Omega}}_{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. ${\mathrm{\Omega}}_{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. ${\mathrm{\Omega}}_{c}$ decreases as the waves age and larger gravity waves develop. ${\mathrm{\Omega}}_{c}=1$ gives a mature sea, and ${\mathrm{\Omega}}_{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 (${\mathrm{\Omega}}_{c}=5$) versus fully developed (${\mathrm{\Omega}}_{c}=0.84$) seas, for FFT surfaces and polarized ray tracing (and ${\varphi}_{\text{sun}}=0$). The last section of Table 2 shows the corresponding percent differences $100[{R}_{\text{surf}}(\text{fully devel.})-{R}_{\text{surf}}(\text{young})]/{R}_{\text{surf}}(\text{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 ${R}_{\text{surf}}$ 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 ${R}_{\text{surf}}$ at large solar zenith angles due to wind speed. For zenith angles greater than 70 deg, ${R}_{\text{surf}}$ decreases by about 30% when the wind increases from $W=0$ to $5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$ and by almost 45% between $W=0$ and $W=15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$. For solar angles less than 60 deg, ${R}_{\text{surf}}$ decreases by 8% between $W=0$ and $5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$, and by 12% between $W=0$ and $15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$.

## 7. Radiance Reflectance Factors

The remote-sensing reflectance ${R}_{\mathrm{rs}}$ is defined as the ratio of the water-leaving radiance ${L}_{w}$ just above the sea surface to the downwelling plane irradiance ${E}_{d}$ incident onto the sea surface: ${R}_{\mathrm{rs}}\equiv {L}_{w}/{E}_{d}$. ${R}_{\mathrm{rs}}$ 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 ${L}_{u}$ just above the surface includes both ${L}_{w}$ and surface-reflected sky radiance ${L}_{\mathrm{sr}}$, so that ${L}_{w}={L}_{u}-{L}_{\mathrm{sr}}$. For a level sea surface, ${L}_{\mathrm{sr}}$ comes only from sky radiance in the direction that is specularly reflected into the sensor. In this case ${L}_{\mathrm{sr}}({\theta}_{v},{\varphi}_{v})=\rho ({\theta}_{v},{\varphi}_{v}){L}_{\mathrm{sky}}({\theta}_{v},{\varphi}_{v})$, where ${L}_{\mathrm{sky}}$ is the sky radiance and $\rho $ is the radiance reflectance for the angle of incidence ${\theta}_{v}$ that connects ${L}_{\mathrm{sky}}$ and ${L}_{\mathrm{sr}}$ by specular reflection [35]. In ${L}_{\mathrm{sr}}({\theta}_{v},{\varphi}_{v})$, ${\theta}_{v}$ is measured from the nadir; in ${L}_{\mathrm{sky}}({\theta}_{v},{\varphi}_{v})$, ${\theta}_{v}$ is the same angle but measured from the zenith. For a level surface and unpolarized radiance, the radiance reflection $\rho $ equals the unpolarized Fresnel reflectance, which is a function only of ${\theta}_{v}$.

For windblown sea surfaces, ${L}_{\mathrm{sr}}$ 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 $\rho $ that converts a sky measurement in a particular direction into an estimate of the surface reflectance: $\rho ({\theta}_{v},{\varphi}_{v})\equiv {L}_{\mathrm{sr}}({\theta}_{v},{\varphi}_{v})/{L}_{\mathrm{sky}}({\theta}_{v},{\varphi}_{v})$. In this case, $\rho $ no longer equals the Fresnel reflectance but rather depends on solar zenith angle, viewing direction, wave state, and sky condition. The purpose of the $\rho $ 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 $\rho $ factor. Thus, $\rho $ is best interpreted as a scale factor that converts a sky radiance measurement in a particular direction, ${L}_{\mathrm{sky}}({\theta}_{v},{\varphi}_{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 ${L}_{w}={L}_{u}-\rho {L}_{\mathrm{sky}}$. ${L}_{u}$ and ${L}_{\mathrm{sky}}$ are measured by the same instrument, which is assumed to be insensitive to the state of polarization. Accurate determination of the value of $\rho $ is the key to obtaining an accurate value for ${L}_{w}$.

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

#### A. $\rho $ 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 ${\theta}_{\text{sun}}=50\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$. 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 ${\varphi}_{v}$, which is measured relative to the Sun’s azimuthal direction. That is, ${\varphi}_{v}=0$ corresponds to looking toward the Sun, ${\varphi}_{v}=90\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ is looking at right angles to the Sun’s incident rays, and ${\varphi}_{v}=180\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ is looking away from the Sun. The bottom panel shows $\rho ({\theta}_{v},{\varphi}_{v})$ for azimuthal directions of ${\varphi}_{v}=90$ and 135 deg. Now, however, the $\rho $ value depends on the azimuthal direction, even for a level sea surface. This ${\varphi}_{v}$ dependence is easily understood as via a simple numerical example.

The surface radiance reflectance matrix, as defined in Eq. (15) and computed by Monte Carlo simulation for incident and reflected angles ${\theta}_{i}={\theta}_{r}=40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$, is

Multiplying $\underline{\mathcal{R}}(40,\varphi ,40,\varphi )$ times these sky radiance Stokes vectors gives the corresponding surface-reflected radiances:

For unpolarized ray tracing, only the (1, 1) element of $\underline{\mathcal{R}}(40,\varphi ,40,\varphi )$ is nonzero. In this case, for the same polarized sky radiances, the reflected radiances are

and These values give and which are the values plotted in the top panel of Fig. 16.Thus, $\rho $ is independent of the azimuthal direction and the polarization state of the sky for unpolarized ray tracing. However, for polarized ray tracing, $\rho $ 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 $\rho $. This is an additional complication in the determination of the value of $\rho $ to be used in processing measured radiances. Not only does $\rho $ 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. $\mathsf{\rho}$ for Windblown Surfaces

Figure 17 shows for reference selected $\rho ({\theta}_{v},{\varphi}_{v})$ values as computed by unpolarized ray tracing and Cox–Munk sea surfaces. $\rho $ is shown a function of off-nadir viewing direction ${\theta}_{v}$ for Sun zenith angles of 30 and 60 deg, azimuthal viewing directions of 90 and 135 deg (relative to the Sun at ${\varphi}_{v}=0$), and wind speeds of 2, 5, 10, and $15\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$. The plot for $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$ reproduces the ${\theta}_{\text{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 ${\theta}_{\text{sun}}=30$ curves) and the sea surface is rough (wind speeds of $5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$ or greater in the figure), the ${\varphi}_{v}=90$ curves have large $\rho $ values at small ${\theta}_{v}$ values because of Sun glint seen by the sensor measuring ${L}_{u}$. For large ${\theta}_{v}$, $\rho $ becomes large because of the rapidly increasing Fresnel reflectance. Based on plots such as these, [18] argued that ${\theta}_{v}=30$ or 40 deg and ${\varphi}_{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 $\rho $ values.

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 $\rho $ values for the given environmental conditions and viewing geometry. For ${\theta}_{v}=40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ and wind speeds of 5 and $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$, the unpolarized case gives tightly grouped values of $\rho =0.0280-0.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 $5\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$ and 0.0195–0.0458 for $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$. There is, thus, a factor-of-two spread in $\rho $ 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 ${\varphi}_{v}=135$, but for the Sun at a right angle to the wind speed and for a very young sea (${\mathrm{\Omega}}_{c}=5$ in the wave spectrum).

Figure 19 displays $\rho $ values in additional ways. The upper-left panel shows a plot of $\rho ({\theta}_{v},{\varphi}_{v})$ for a single Sun angle of ${\theta}_{\text{sun}}=40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ located at ${\varphi}_{v}=0$, as computed using FFT surfaces and polarized ray tracing. The radial distance in the plot is the off-nadir viewing angle ${\theta}_{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 $\rho $ values needed to remove Sun glint. The orange to red colors around the rim are high $\rho $ 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 ${\varphi}_{v}=90$ and 135 deg directions (but for different Sun zenith angles). The black dot indicates the $\rho ({\theta}_{v},{\varphi}_{v})=(40,135)$ viewing direction. In the present example, $\rho (40,135)=0.045$. The upper-right panel of the figure shows $\rho ({\theta}_{\text{sun}},{\varphi}_{v})$ for an off-nadir viewing direction of ${\theta}_{v}=40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$. 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 $\rho $ resulting from Sun glint with looking in directions too near the Sun’s azimuthal direction. The blue colors along ${\varphi}_{v}=135$ indicate $\rho $ 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 $\rho $ 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 $10\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{m}\text{\hspace{0.17em}}{\mathrm{s}}^{-1}$ conditions of this figure, the lower-left panel shows that the ratio of $\rho (\mathrm{FFT},\text{polarized})$ to $\rho $ (Cox–Munk, unpolarized) is 0.91 for ${\theta}_{\text{sun}}=40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ and $\rho ({\theta}_{v},{\varphi}_{v})=(40,90)$ and is 1.31 for $\rho ({\theta}_{v},{\varphi}_{v})=(40,135)$ (the viewing direction indicated by the black dot).

The simulations of this section used a wavelength of 550 nm for generation of the Rayleigh single-scattering sky radiance distribution. However, $\rho $ 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 $\rho $ 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 ${\theta}_{v}\approx 40\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ and ${\varphi}_{v}\approx 135\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{deg}$ (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 $\rho $ 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 $\rho $ is critical for the estimation of ${L}_{w}={L}_{u}-\rho {L}_{\mathrm{sky}}$ because the measured ${L}_{\mathrm{sky}}({\theta}_{v},{\varphi}_{v})$ is typically an order of magnitude or more greater than ${L}_{w}$. Thus, a small error in $\rho $ can give a large error in ${L}_{w}$. To address the deficiency of the $\rho $ values previously computed using Cox–Munk surfaces and unpolarized ray tracing, a new table of $\rho $ 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 ${R}_{\text{surf}}$ 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 ${R}_{\text{surf}}$ 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 ${R}_{\text{surf}}$ 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 $\rho $, 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, $\rho $ 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 $\rho $ 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 $\rho $ 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 $\rho $ factors would likely differ somewhat from the values shown here. Moreover, the $\rho $ 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 $\rho $ 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]