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

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 =\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){\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\mathrm{m}\times 200\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\mathrm{m}{\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:

 

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

Download Full Size | PDF

Now suppose this spectrum is used to generate a surface with $L_x=200\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\mathrm{rad}{\mathrm{m}}^{-1}$ and the Nyquist frequency is $k_{\mathrm{Ny}}=16.085\mathrm{rad}{\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 $〈z^2〉(N=1024)=0.4219{\mathrm{m}}^2$ and ${\sigma }^2(N=1024)=0.02584{\mathrm{rad}}^2$. Thus, the finite range of the sampled frequencies includes the fraction,

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,

$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\mathrm{m}$ area of sea surface is simulated using an $N\times N$ FFT, for a wind speed of $W=5\mathrm{m}{\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\mathrm{m}$ and $N=2^{15}$, the surface is being resolved down to a grid spacing of $L/N=0.003\mathrm{m}$ (or a Nyquist spatial frequency of $k_{\mathrm{Ny}}=1029\mathrm{rad}{\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\mathrm{m}$ and $W=5\mathrm{m}{\mathrm{s}}^{-1}$, there is almost no difference in the elevation variances. The blue curves show the corresponding results for $L=300\mathrm{m}$ and $W=10\mathrm{m}{\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\mathrm{rad}{\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.

 

Fig. 2. Approach to fully resolved mean square slopes (diamonds, left ordinate) and elevations (circles, right ordinate) as a function of the number of grid points $N$, for winds speeds of $W=5$ (red) and $10\mathrm{m}{\mathrm{s}}^{-1}$ (blue).

Download Full Size | PDF

The surface generation techniques of this section give sea surfaces that correctly reproduce the elevation and slope statistics of a real sea surface corresponding to the chosen variance spectrum. Figure 3 shows an example of a 2D sea surface generated using these techniques. The along-wind, cross-wind, and total slope variances for this realization (computed by finite differences) are in good agreement with the respective Cox–Munk values of 0.035, 0.019, and 0.054. The significant wave height of $H_{1/3}=2.18\mathrm{m}$ for this realization with $L_x=100\mathrm{m}$ is somewhat less than the theoretical value of $H_{1/3}\equiv 4\sqrt{〈z^2〉}=2.622\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\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)$.

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

View Table | View all tables in this article

 

Fig. 3. 2D example sea surface random realization for $L_x=L_y=100\mathrm{m}$ and $N_x=N_y=512$ and a wind speed of $10\mathrm{m}{\mathrm{s}}^{-1}$ in the $+x$ direction. Light color is high surface elevation; dark color is low elevation. Inset numbers give the significant wave height $H_{1/3}$, mean square slopes (total, mss; along wind, ${\mathrm{mss}}_x$; and cross wind, ${\mathrm{mss}}_y$), and mean surface tilt angles from the vertical in the along-wind ($〈{\theta }_x〉$) and cross-wind ($〈{\theta }_y〉$) directions.

Download Full Size | PDF

3. Polarized Ray Tracing

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

A. Resampling Rectangular Grids

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

There is no unique or one best way to map a rectangular FFT grid to a grid of triangles. However, the ray tracing algorithm that determines when a ray intersects the surface does depend on the particular triangulation scheme chosen because it is necessary to determine when a ray crosses the boundaries of the facet triangles. If the FFT surface is generated using $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\mathrm{m}$). This resampling generates a hexagonal region with $3(N_y/2)(N_y/2+1)+1$ facet vertices and $3N_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].

 

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

Download Full Size | PDF

B. Ray Tracing

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

Polarized light is described by the four-component Stokes vector $\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\mathrm{deg}$ polarization, and negative $U$ is $-45\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·{\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

 

Fig. 5. Four cases for determining Stokes vector rotation angles. ${\mathbf{v}}_1$ is the initial direction, which is rotated into the final direction ${\mathbf{v}}_2$ by a counterclockwise rotation about direction $\mathit{\xi }$.

Download Full Size | PDF

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

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={\cos }^{-1}|{\mathit{\xi }}_i·\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{or}{\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

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\mathrm{deg}$, except for $\mathrm{\Delta }\theta =5\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\mathrm{deg}$). For the $10\times 15\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\mathrm{m}{\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\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\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\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.

 

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

Download Full Size | PDF

 

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

Download Full Size | PDF

 

Fig. 8. Reflected energy pattern (underwater glitter pattern) for water-incident light at a 50 deg incident angle from the nadir and a $10\mathrm{m}{\mathrm{s}}^{-1}$ wind speed. The light source is unpolarized in an otherwise black ocean.

Download Full Size | PDF

 

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

Download Full Size | PDF

For the air-incident case of Fig. 6, 3.87% of the incident energy is reflected, and the remaining 96.13% is transmitted into the water. The upper-right panel of Fig. 6 shows that the reflected light is as much as 97% horizontally polarized in the quads near $\theta =50\mathrm{deg}$, which is near the Brewster angle of 53 deg for a level surface. The $\pm 45\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{\log }_2{N}_x)(N_y{\log }_2{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.

 

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

Download Full Size | PDF

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\mathrm{m}{\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.

 

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

Download Full Size | PDF

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

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 $\cos \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\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.

 

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

Download Full Size | PDF

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

Figure 13 shows the surface-reflected sky radiance for the Sun at 50 deg and the partially polarized sky of Fig. 12. The sea surface is an FFT surface for a wind speed of $10\mathrm{m}{\mathrm{s}}^{-1}$, generated with $L_x=200\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\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\mathrm{W}{\mathrm{m}}^{-2}{\mathrm{sr}}^{-1}{\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\mathrm{W}{\mathrm{m}}^{-2}{\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:

 

Fig. 13. Surface-reflected radiance for a single-scattering Rayleigh sky and a wind speed of $10\mathrm{m}{\mathrm{s}}^{-1}$. The Sun zenith angle is 50 deg and the viewing geometry corresponds to Fig. 6.

Download Full Size | PDF

Figure 14 shows the transmitted radiance corresponding to Fig. 13. As in Fig. 7, most of the transmitted radiance is concentrated into a few quads near the refracted direction of the Sun’s direct beam, and the radiance in those quads is only weakly polarized. However, the state of linear polarization away from the direct beam direction is determined by the sky polarization in directions away from the Sun, and the degree of polarization can be large. However, those directions contain very little energy. For a level surface, all of the transmitted energy would be contained within the Snell’s cone of angle 48.2 deg from the nadir. Now, however, for the $10\mathrm{m}{\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.

 

Fig. 14. Surface-transmitted radiance for a single-scattering Rayleigh sky and a wind speed of $10\mathrm{m}{\mathrm{s}}^{-1}$. The Sun zenith angle is 50 deg, and the viewing geometry corresponds to Fig. 7.

Download Full Size | PDF

6. Irradiance Reflectance

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

The sky total downwelling plane irradiance is

 

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

Download Full Size | PDF

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\mathrm{m}{\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.

Table 2. Differences in Surface Irradiance Reflectance ${\mathsf{R}}_{\mathsf{\text{surf}}}$ as Computed in Various Ways^a

View Table | View all tables in this article

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\mathrm{m}{\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\mathrm{m}{\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\mathrm{m}{\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\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\mathrm{m}{\mathrm{s}}^{-1}$ and by almost 45% between $W=0$ and $W=15\mathrm{m}{\mathrm{s}}^{-1}$. For solar angles less than 60 deg, $R_{\text{surf}}$ decreases by 8% between $W=0$ and $5\mathrm{m}{\mathrm{s}}^{-1}$, and by 12% between $W=0$ and $15\mathrm{m}{\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$.