## Abstract

The behavior of light at the air-sea interface has been investigated using ray tracing methods with numerically realized surfaces that incorporate features on scales from 3 millimeters to 200 meters. The directional reflection of light at the surface realizations was tested using Monte Carlo code. Estimated directionally reflected radiances were generally in good agreement with those from existing methods that model the slope statistics but not the shape of the sea surface. However, significant differences were found for some incident and exitant directions. The model was used to quantitatively estimate the pixel-to-pixel variation in ocean color images caused by spatial variation in the sea surface shape.

©2011 Optical Society of America

## 1. Introduction

Radiative transfer modeling for marine optical and infrared applications requires a detailed understanding of the behavior of light at the water surface. Applications which require this knowledge include chlorophyll estimation from ocean color measurements [1], air-water interface corrections for benthic classification in coral reef or other optically shallow water environments [2], and radiative transfer model inversion for bathymetry and other retrieval of biophysical parameters [3]. The behavior of the near-surface light field is currently of great interest, as demonstrated by the interdisciplinary observational and modeling effort of the recent RaDyO programme [4] and the development of new instrumentation targeting this component of the hydrological optics system [5–7]. The performance of the neural network-based atmospheric correction used in the CoastColour project [8] shows that improved modeling of surface reflectance is possible; a physics-based model would give more insight into the processes involved.

A commonly used state-of-the art software implementation, HydroLight [9], and many other radiative transfer models, e.g [1,10–16], use a model for the sea surface which reproduces the slope statistics of Cox and Munk [17]. However, until recently limitations on computer memory and processor speed have prevented these models from attempting to capture the elevation distribution of the water surface [18], hence the modeled surfaces are essentially flat with an elevation variance close to zero [19] (Fig. 1a ).

Slope statistics methods can be quite successful at modeling the surface albedo and spatially-averaged glint patterns [10,11], because the primary factor determining the path of light incident on the surface is the local surface slope. However, a complete model should also include elevation-dependent processes such as wave shadowing and multiple reflections between wave sides [18]. Some models include an approximation to these processes, using surfaces with small elevations. Monte Carlo models include them automatically, since each ray is traced through as many interactions as necessary for it to leave the surface region entirely. Matrix methods [12] treat the atmosphere and ocean as a series of layers, with the air-sea interface characterized by transmission and reflection functions incorporating the chosen slope statistics. These functions are commonly modified by a shadowing function, which gives the probability that a point on the surface will be illuminated [20,21]. A similar approach is taken in the discrete-ordinate method, which can also account for multiple reflection by successive application of the transmittance and reflection functions [10,13]. The successive orders of scattering method treats the sea surface as a boundary, which can include a shadowing function, and multiple scattering at the surface is built into the model [14]. However, in all of these models the surfaces used have much smaller elevations than real seas. They also cannot model spatial variation in reflectance due to large scale sea surface variation.

The work reported here compares results from radiative transfer modeling using a 3-dimensional representation of the sea surface to those obtained using a slope-statistics method. The surface realizations (Fig. 1b) incorporate both slope and elevation characteristics of observed seas and include spatial features on scales from 3 millimeters to 200 meters. They were produced by adding components of the sea wave spectrum using as wide a range of wavenumbers as possible and a fully 2-dimensional wavenumber space. The spectral synthesis technique is well established [22,23], but previously published sea surface models used for radiative transfer have either been based on a spectrum that only models the longer wavelengths [24,25] or only incorporated one spatial dimension [26]. An efficient ray-tracing algorithm was used to characterize the directional reflection and transmission across the numerically realized sea surface. The model incorporates multiple interactions at the air-water interface and allows the evaluation of bidirectional transmission and reflectance distribution functions (BRDFs) suitable for incorporation into models such as HydroLight [9].

We have used the omnidirectional spectrum of Elfouhaily *et al.* [27], which covers a wavenumber range from 10^{−2} to 10^{3} rad m^{−1} by treating the spectrum as two parts: a low wavenumber section based on the Pierson-Moskowitz and JONSWAP spectra [28,29] and a high wavenumber gravity-capillary wave section derived from theoretical and laboratory studies. The wavenumber dependence of this spectrum has been supported by the study of Hauser *et al.* [30]. However, both these authors and Ross and Dion [31] note that the spreading function used by Elfouhaily *et al.* [27] gives an upwind to crosswind slope ratio that falls with increasing wind speed, whereas the data of Cox and Munk [17] and the more recent satellite observations of Bréon and Henriot [32] show a ratio that rises with wind speed above 5 m s^{−1}. In this work the Elfouhaily *et al*. spreading function has been replaced by a slightly modified version of that proposed by Heron *et al.* [33]. This gives a better fit to the upwind to crosswind ratio observed by Cox and Munk and more recent authors (Fig. 2
). It also includes a bimodal directional distribution, as has been found in the gravity range above the spectral peak [34–36].

In Section 2 we describe how the sea surfaces are realized and validated. Section 3 outlines the ray-tracing process and presents results showing how the model output compares to that produced using a standard “slope statistics” method. Section 4 shows how the model can be used to estimate the variation in radiance in ocean color images caused by changes in the sea surface shape. Section 5 summarizes our findings and briefly discusses future work.

## 2. Creating sea surface realizations

#### 2.1. Method

Sea surface realizations were constructed from the elevation variance spectrum, based on the linear assumption that the elevation can be expressed as a harmonic series of wavenumber components:

where η(**r**,*t*) is the sea surface elevation at position **r**(*x,y)* and time *t*, **k** is the wavenumber vector (*k _{x}, k_{y}*), and

*X*(

**k**,

*t*) is the amplitude of component

**k**at time

*t*.

The component amplitudes *X*(**k**) for any point in time were calculated from the elevation variance spectrum Ψ(**k**), with a randomization element [22,23]:

where *ξ _{r}* and

*ξ*are random numbers chosen from a Gaussian distribution of mean 0 and variance 1, and Δ

_{i}*k*represent the wavenumber interval in the spectral space of Eq. (1). The spectrum value was calculated using:

_{x},Δk_{y}where *S _{Elf}*(

*k*) is the omnidirectional spectrum of Elfouhaily

*et al.*[27] and

*D*

_{Heron}(

*k, θ*) is the directional spreading function of Heron

*et al.*[33].

An inverse Fourier transform was then used to derive the surface elevation at points across a 2-d grid. Repeating these steps with different random numbers produced a set of surfaces with average slope and variance specified by the spectrum Ψ(**k**). The use of a Fourier method means that each surface has periodic boundaries; in ray-tracing the surface is treated as effectively infinite, with the pattern regularly repeating in the *x* and *y* directions.

Sea surface realizations of 65536 × 65536 vertices were created: this enabled a wavenumber range from 0.036 to 2100 rad m^{−1} to be included, allowing most features of the slope and elevation variance spectra at wind speeds up to 15 m s^{−1} to be captured (Fig. 3
). This wavenumber range produced sea surface realizations on a 200 m square grid, with a grid spacing of 3 mm (Fig. 4
). Creation of each surface took about 6 hours on a 3MHz, 8GB desktop PC.

Storing an entire surface realization in processing memory is a limiting factor on current computer hardware. Therefore, we developed a system that enables a single surface to be created and ray-traced as a set of adjacent segments, but with full cross-segment propagation of rays (see section 3). The results from this system were identical to those produced by creating and ray tracing the entire realization in one pass. The Fourier transform step was performed using FFTW3 [37], which has high numerical accuracy [38] and avoids stability problems when handling these large grids; using double precision in the FT calculation we did not experience the accuracy problems noted by Tessendorf [23] for grid sizes over 2048 points per side.

#### 2.2. Validation of surface properties

The modeled surfaces were validated by comparing their mean slopes and elevation standard deviation against empirically established values (Fig. 5 ). For mean square slope the basis of comparison was the widely used Cox and Munk model [17]:

where *U*
_{10} is the wind speed at 10 m above the sea surface. In a sample of 108 surfaces with wind speeds ranging from 1 to 15 m s^{−1}, all had mean square slopes within the range specified in Eq. (5) (Fig. 5a). The elevation standard deviation was compared to the empirical formula given by Apel [39]:

The spectrum of Elfouhaily *et al.* [27] gives a reasonable match to this approximate formula, and the realized surfaces follow the spectrum closely (Fig. 5b).

## 3. Ray tracing

#### 3.1. Method

The transfer of light across the realized surfaces was modeled using Monte Carlo ray tracing code developed specifically for this application. Rays were directed to the sea surface from randomly chosen angles from the full range of above- and below-water directions. The sphere of all directions was partitioned into 432 “quads”, each 15° horizontally by 10° vertically except at the horizon where the vertical angle was 5°, plus endcaps of 10° at the two poles (Fig. 6 ).

This is the scheme used in HydroLight [9], making it straightforward to transfer results between models. Interaction with the water surface was modeled using the Fresnel reflection formulae with wavelength-independent refractive indices of 1.00 for air and 1.34 for water [19]. The energy associated with each ray was tracked and updated at each intersection event with the surface. Upon run completion energy conservation between the incident and exitant rays was checked as an accuracy assessment step – differences were less than 0.3% for all runs. The energy loss results from rays travelling very close to the horizontal direction that are eventually discarded for time efficiency; it can be reduced at the cost of increased run time. For every wind speed ten surfaces were tested, each with five ray-tracing runs of 10^{7} rays. The average of the five runs was used to create the bidirectional quad-averaged reflectance distribution function (BRDF) for that sea surface realization. A number of sea surfaces were created and tested, enabling the variance in model output due to changing surface shape to be quantified.

The standard error for the five ray-tracing runs was below 3% for the majority of cases; calculation of the running mean for up to 25 runs indicated that five runs is enough to avoid large fluctuations in results. Twenty runs would give better convergence, but time constraints did not allow for this in the current work. Tests using ten runs gave smaller standard errors but did not change the main findings reported here.

The calculated BRDFs were used to predict the radiance reflected from the upper side of the water surface for a clear sky scenario [40], tabulated over outgoing azimuth and zenith angle. For comparison, the same ray-tracing method and analysis was applied to surfaces created according to the method described by Mobley [19], which have the mean square slopes of the Cox and Munk model [17] but with elevation variance close to zero (see Fig. 1). A shadowing function was not included in this model. All results shown are for a wavelength of 470 nm; other wavelengths were not investigated as both surface models were wavelength-independent. Note that in HydroLight and other models [1,10,14,] the provided air-water interface functions are azimuthally averaged, and therefore independent of wind direction. In our study both azimuthally-averaged and non-averaged outputs were tested.

The reflected radiance estimated by the new model and by a Cox-Munk slope statistics model (CM) was compared at all directional quads in the upper hemisphere for wind speeds 5, 10 and 15 m s^{−1}, solar zenith angles for quads centered at 0° to 80° in 10° steps and wind-sun relative azimuth angle 0° to 90° in 15° steps, giving a total of 32424 data points. Directional quads at the air-sea boundary centered at zenith 87.5° were excluded, so the largest incident and viewing zenith angle was 75°-85°. In all cases presented in this paper, results are averaged over 10 surfaces, each traced by 5x10^{7} rays, for the new model. The CM results are averaged over 10 runs, each using 2000 surfaces of 14152 facets and 4340 rays: these values were chosen to give a similar variance for repeated runs to the new model. The time for 10 surface runs was around 200 hours using the new model, about 48 hours using the CM model.

#### 3.2. Results

The two models were in generally good agreement (r^{2} = 0.9865 for the correlation between radiance estimated by the new model and the CM model), but the difference was statistically significant at 41% of the data points tested (*t*-test, unequal variances, two-tailed, p<0.05, with Simes-Hochberg correction for multiple comparisons). When azimuthally averaged models were used these figures were r^{2} = 0.9896, with significant differences at 75% of data points.

The new model predicts reduced forward scattering and increased cross- and back-scattering (Fig. 7 ), consistent with the effect of elevation-dependent processes such as wave shadowing and multiple reflections. The number of points with a significant difference was greater at lower wind speeds and higher viewing zenith angle.

Figure 8
shows sample reflected radiance distributions estimated by the two models for a wind of speed 10 m s^{−1} aligned with the sun.

The effect of wind direction and azimuthal averaging is illustrated in Fig. 9 . The non-azimuthally averaged CM model estimates the maximum radiance occurring at an angle of several degrees away from the solar plane when the wind is at 45° to the sun; the skew estimated by the new model is less strong (Fig. 9b). For the slope-statistics model, the mean square slope is only defined in the upwind and crosswind directions so a form of interpolation, fixed by the numerical construction of the surface, gives the slopes at other angles [19]. By contrast, the spreading function used in the new model [33] is defined for all angles, and has empirical support across a wide range of wavenumbers: this difference is the probable cause of the variation in asymmetry seen in Fig. 9b. Azimuthally averaging the models (Fig. 9e) removes the effect of wind direction, so the asymmetry disappears entirely. Note that the asymmetry is not seen when the wind is at 90° to the sun (Fig. 9c) because of the way surfaces are created. Waves travelling in opposite directions are given equal weighting (Eq. (2), since the realized surface is a snapshot view and has no information about wave propagation over time. This introduces an extra symmetry: the same surface would be created if the wind were blowing in the diametrically opposite direction.

## 4. Application: pixel-to-pixel variation in remote sensing images of the ocean

The model was applied to investigate the pixel-to-pixel variation caused by differing sea surface shape in images from wide swath pushbroom instruments such as the Medium Resolution Imaging Spectrometer (MERIS), on the European Space Agency’s Envisat mission, or the Ocean and Land Colour Instrument (OLCI) planned for the Sentinel 3 satellite series. The largest scale of elevation feature included in the current model (200 m) is similar to MERIS/OLCI full-resolution pixels (260 x 290 m) so the variation in results using different surfaces can give an indication of the variability in radiance due to differences in the distribution of elevation from pixel to pixel. The slope-statistics method does not model spatial distributions of elevation so cannot give this information.

Models runs were done using a clear sky scenario as in the previous section, for a range of wind speeds and sun and viewing geometries that might be encountered by ocean color instruments. The effect of the surface only was included: atmospheric correction was not considered and would need to be added to the model to obtain top-of-atmosphere values. Comparison of results obtained using 10 different sea surface realizations showed that the radiance reflected in the direction of the sensor typically varied by a few percent either side of the mean value. In some cases the variation was as high as 40-50% - the highest variations were found for higher wind speeds, larger solar zenith angles and near-nadir viewing, but there was no simple relationship between the reflected radiance and any one variable. Figure 10 gives examples for two different sun positions: the reflected radiance is plotted against viewing angle across the track of one of the above sensors, for which the field of view is typically 30°-50° from nadir. The position of the sun is fixed, as is true to within a few degrees for a single remote-sensing image.

For wind speeds around 5 m s^{−1} the upwind and crosswind slopes are nearly equal (Fig. 2) and so there is little variation with wind-sun azimuth (Fig. 10 a,d). As wind speed rises the anisotropy increases and Fig. 10 c,f suggest that at higher wind speeds the relative wind direction can have a substantial effect on the size of the estimated reflected radiance. The slope-statistics model has a weaker relation between wind speed and upwind to crosswind slope ratio, so gives a more constant azimuthal dependence. However, as noted in Section 3.2, the directional behavior of this model is less well-founded.

## 5. Conclusion and further work

This work has demonstrated that it is feasible to generate high-resolution sea surface models and use them in radiative transfer modeling of the ocean surface. Incorporating correct statistics of surface elevation as well as slope gave estimates of reflected radiance that were significantly different from a Cox-Munk slope-statistics model in 41% of cases tested, which included wind speeds up to 15 m s^{−1} and solar / instrument viewing angles in all parts of the upper hemisphere. The new model estimates reduced forward reflectance and increased sideways scattering, consistent with the effect of elevation-dependent processes such as wave shadowing and multiple reflection. The ability to implement a full directional spreading function into a 2-dimensional model of surface elevation and slope means that the new model can estimate an asymmetric radiance distribution when the wind is not aligned with the sun.

The variation in estimated reflected radiance gives an indication of the pixel-to-pixel variation in ocean color images due to differences in the shape of the sea surface: this cannot be modeled using the slope-statistics method. The variation in reflected radiance just above the sea surface was found to be typically a few per cent of the mean, but as much as 40-50% in some situations.

We believe this approach has the potential to improve modeling of radiative transfer at the air-sea interface in a number of applications where slope-statistics methods have traditionally been used. These include situations where the spatial scale is much larger than the surface features modeled: for ocean color satellites such as MERIS a surface realization is approximately the size of one pixel, but much smaller features can affect the surface reflectance. One example is the variation in reflectance due to surface shape (see Section 4); another is estimating the spatial distribution of sun glint [41]. The surface realization method could also be useful when running a model for data correction, for example to set up a training set for a neural network model [8]. The geometric optics assumption in the ray-tracing model makes it unsuitable for direct application to microwave data. However, the surface creation method could be used alongside a microwave radiative transfer model, with the required surface feature size re-evaluated.

Future work will aim to incorporate polarization. The modeled effect on total radiance reflected by the surface from an unpolarized source is likely to be small; however, if a polarized sky radiance distribution were included, or the polarization of the reflected light were of interest, then the results may be substantially different. In addition, if the model was extended to include water column scattering below the surface there would be an increased chance that rays would undergo successive polarization events, and so polarization-dependent effects may become apparent even if unpolarized illumination and measurement were assumed. Recent studies have shown that neglecting polarization can compromise the accuracy of atmospheric radiative transfer codes [42], so there may be accuracy advantages in a fully coupled vectorial model for atmosphere, air water interface and water column.

The possibility of including whitecaps and foam at higher winds speeds is also being investigated. In addition, adjustments of the wave spectrum would allow processes not produced by the local wind to be included, e.g. long-range swell or modulation of surface waves by internal waves [39,43], and the use of a non-Gaussian surface model for non-linear wave processes [44]. It would also be desirable to create surfaces that reproduce the correct spatial distribution of capillary waves, which are generated in coherent groups by gravity waves and phase-locked with their generating wave [45].

## Acknowledgments

This work was supported by ARGANS Ltd and Great Western Research (award ref. 363). J. Hedley’s contribution was partly funded by the NERC (grant number NE/E015654/1). The authors would like to thank two reviewers for their helpful comments.

## References and Links

**1. **H. R. Gordon and M. Wang, “Surface-roughness considerations for atmospheric correction of ocean color sensors. I: The Rayleigh-scattering component,” Appl. Opt. **31**(21), 4247–4260 (1992). [CrossRef] [PubMed]

**2. **E. Hochberg, S. Andrefouet, and M. Tyler, “Sea surface correction of high spatial resolution Ikonos images to improve bottom mapping in near-shore environments,” IEEE Trans. Geosci. Rem. Sens. **41**(7), 1724–1729 (2003). [CrossRef]

**3. **D. Lyzenga, N. Malinas, and F. Tanis, “Multispectral bathymetry using a simple physically based algorithm,” IEEE Trans. Geosci. Rem. Sens. **44**(8), 2251–2259 (2006). [CrossRef]

**4. **T. Dickey, “Toward the understanding and prediction of optics near the ocean surface,” presented at Ocean Optics XX, Anchorage, Alaska, 27 Sept.–1 Oct. 2010.

**5. **R. Van Dommelen, J. Wei, M. R. Lewis, and K. J. Voss, “Instrumentation, calibration and validation of a high dynamic range radiance camera,” presented at Ocean Optics XX, Anchorage, Alaska, 27 Sept.–1 Oct. 2010.

**6. **P. Bhandari, K. J. Voss, and L. Logan, “The variation of the polarized downwelling radiance distribution with depth in coastal water,” presented at Ocean Optics XX, Anchorage, Alaska, 27 Sept.–1 Oct. 2010.

**7. **M. Darecki, D. Stramski, and M. Sokolski, “An underwater porcupine radiometer system for measuring high-frequency fluctuations in light field induced by sea surface waves,” presented at Ocean Optics XIX, Barga, Italy, 6–10 Oct. 2008.

**8. **C. Brockmann, R. Doerffer, S. Sathyendranath, Z. Lee, and S. Pinnock, “Towards operational coastal ocean colour products – the CoastColour approach,” presented at Ocean Optics XX, Anchorage, Alaska, 27 Sept.–1 Oct. 2010.

**9. **C. D. Mobley and L. K. Sundman, “Hydrolight 5 Users' Guide,” (2008), http://www.sequoiasci.com/products/Hydrolight.aspx.

**10. **Z. Jin, T. P. Charlock, K. Rutledge, K. Stamnes, and Y. Wang, “Analytical solution of radiative transfer in the coupled atmosphere-ocean system with a rough surface,” Appl. Opt. **45**(28), 7443–7455 (2006). [CrossRef] [PubMed]

**11. **R. Preisendorfer and C. Mobley, “Albedos and glitter patterns of a wind-roughened sea-surface,” J. Phys. Oceanogr. **16**(7), 1293–1316 (1986). [CrossRef]

**12. **T. Nakajima and M. Tanaka, “Effect of wind-generated waves on the transfer of solar-radiation in the atmosphere ocean system,” J. Quant. Spectrosc. Radiat. Transf. **29**(6), 521–537 (1983). [CrossRef]

**13. **M. Ottaviani, R. Spurr, K. Stamnes, W. Li, W. Su, and W. Wiscombe, “Improving the description of sunglint for accurate prediction of remotely sensed radiances,” J. Quant. Spectrosc. Radiat. Transf. **109**(14), 2364–2375 (2008). [CrossRef]

**14. **P. Zhai, Y. Hu, J. Chowdhary, C. R. Trepte, P. L. Lucker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transf. **111**(7-8), 1025–1040 (2010). [CrossRef]

**15. **E. Raschke, “Multiple-scattering calculations of transfer of solar-radiation in an atmosphere-ocean system,” Bull. Am. Meteorol. Soc. **53**, 501 (1972).

**16. **G. N. Plass, G. W. Kattawar, and J. A. Guinn Jr., “Radiative transfer in the earth’s atmosphere and ocean: influence of ocean waves,” Appl. Opt. **14**(8), 1924–1936 (1975). [CrossRef] [PubMed]

**17. **C. Cox and W. Munk, “Measurement of the roughness of the sea surface from photographs of the Suns glitter,” J. Opt. Soc. Am. **44**(11), 838–850 (1954). [CrossRef]

**18. **C. Mobley, “How well does Hydrolight simulate wind-blown sea surfaces?” (Hydrolight Technical Note 2002). http://www.sequoiasci.com/products/Hydrolight.aspx.

**19. **C. Mobley, *Light and Water: Radiative Transfer in Natural Waters* (Academic Press, San Diego, Calif, 1994).

**20. **M. Sancer, “Shadow-corrected electromagnetic scattering from a randomly rough surface,” IEEE Trans. Antenn. Propag. **17**(5), 577–585 (1969). [CrossRef]

**21. **B. Smith, ““Geometrical shadowing of a random rough surface,” IEEE Trans. Antennas Propag, A **668–671**, 15 (1967).

**22. **E. Thorsos, “The validity of the Kirchhoff approximation for rough-surface scattering using a Gaussian roughness spectrum,” J. Acoust. Soc. Am. **83**(1), 78–92 (1988). [CrossRef]

**23. **J. Tessendorf, “Simulating Ocean Water,” (2004). http://tessendorf.org/reports.html*.*

**24. **Y. You, P. W. Zhai, G. W. Kattawar, and P. Yang, “Polarized radiance fields under a dynamic ocean surface: a three-dimensional radiative transfer solution,” Appl. Opt. **48**(16), 3019–3029 (2009). [CrossRef] [PubMed]

**25. **F. Schwenger and E. Repasi, “Sea surface simulation for testing of multiband imaging sensors,” Proc. SPIE **5075**, 72–84 (2003). [CrossRef]

**26. **S. Fauqueux, K. Caillault, P. Simoneau, and L. Labarre, “Multiresolution infrared optical properties for Gaussian sea surfaces: theoretical validation in the one-dimensional case,” Appl. Opt. **48**(28), 5337–5347 (2009). [CrossRef] [PubMed]

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

**28. **W. Pierson Jr and L. Moskowitz, “Proposed spectral form for fully developed wind seas based on similarity theory of S. A. Kitaigorodskii,” J. Geophys. Res. **69**(24), 5181–5190 (1964). [CrossRef]

**29. **K. Hasselmann, T. Barnett, E. Bouws, D. E. Carlson, D. E. Cartwright, K. Enke, and J. Ewing *et al.*., “Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP),” Ergnzungsheft zur Deutschen Hydrographischen Zeitschrift Reihe **8**, 95 (1973).

**30. **D. Hauser, G. Caudal, S. Guimbard, and A. A. Mouche, “A study of the slope probability density function of the ocean waves from radar observations,” J. Geophys. Res. **113**(C2), C02006 (2008). [CrossRef]

**31. **V. Ross and D. Dion, “Sea surface slope statistics derived from Sun glint radiance measurements and their apparent dependence on sensor elevation,” J. Geophys. Res. **112**(C9), C09015 (2007). [CrossRef]

**32. **F. Bréon and N. Henriot, “Spaceborne observations of ocean glint reflectance and modeling of wave slope distributions,” J. Geophys. Res. **111**(C6), C06005 (2006). [CrossRef]

**33. **M. Heron, W. Skirving, and K. Michael, “Short-wave ocean wave slope models for use in remote sensing data analysis,” IEEE Trans. Geosci. Rem. Sens. **44**(7), 1962–1973 (2006). [CrossRef]

**34. **K. Ewans, “Observations of the directional spectrum of fetch-limited waves,” J. Phys. Oceanogr. **28**(3), 495–512 (1998). [CrossRef]

**35. **D. Wang and P. Hwang, “Evolution of the bimodal directional distribution of ocean waves,” J. Phys. Oceanogr. **31**(5), 1200–1221 (2001). [CrossRef]

**36. **C. E. Long and D. T. Resio, “Wind wave spectral observations in Currituck Sound, North Carolina,” J. Geophys. Res. **112**(C5), C05001 (2007). [CrossRef]

**37. **M. Frigo and S. G. Johnson, “The Design and Implementation of FFTW3,” Proc. IEEE **93**(2), 216–231 (2005). [CrossRef]

**38. **M. Frigo and S. G. Johnson, “FFT Accuracy Benchmark Results,” http://www.fftw.org/accuracy/.

**39. **J. Apel, “An Improved Model of the Ocean Surface-Wave Vector Spectrum and Its Effects on Radar Backscatter,” J. Geophys. Res. **99**(C8), 16269–16291 (1994). [CrossRef]

**40. **R. H. Grant, G. M. Heisler, and W. Gao, “Photosynthetically-active radiation: sky radiance distributions under clear and overcast conditions,” Agric. For. Meteorol. **82**(1-4), 267–292 (1996). [CrossRef]

**41. **S. Kay, J. D. Hedley, and S. Lavender, “Sun glint correction of high and low spatial resolution images of aquatic scenes: a review of methods for visible and near-infrared wavelengths,” Remote Sens. **1**(4), 697–730 (2009). [CrossRef]

**42. **S. Y. Kotchenova, E. F. Vermote, R. Levy, and A. Lyapustin, “Radiative transfer codes for atmospheric correction and aerosol retrieval: intercomparison study,” Appl. Opt. **47**(13), 2215–2226 (2008). [CrossRef] [PubMed]

**43. **W. J. Plant, W. C. Keller, K. Hayes, G. Chatham, and N. Lederer, “Normalized radar cross section of the sea for backscatter: 2. Modulation by internal waves,” J. Geophys. Res. **115**(C9), C09033 (2010). [CrossRef]

**44. **F. Nouguier, C. Guérin, and B. Chapron, “““Choppy wave” model for nonlinear gravity waves,” J. Geophys. Res. **114**(C9), C09012 (2009). [CrossRef]

**45. **J. R. Zaneveld, E. Boss, and P. Hwang, “The influence of coherent waves on the remotely sensed reflectance,” Opt. Express **9**(6), 260–266 (2001). [CrossRef] [PubMed]