Numerical Resolution

The question naturally arises: “How large a region $L$ and how many points $N$ must be used in generating sea surface realizations?” Not surprisingly, the general answer is, “It depends on your application.” It is computationally possible to create 2-D surfaces with sufficiently large $N$ values in the $x$ and $y$ directions, but run times are prohibitive. Kay et al. (2011) created 2-D surfaces with $2^{16}\times 2^{16}$ = $65536\times 65536$ points in each dimension, which sampled the variance spectrum for wavelengths from 200 m gravity waves to 3 mm capillary waves. However, it took 6 hours to create just one surface on a 3 GHz computer. Note also that the storage of this surface requires $2^{32}=4.3\times 10^9$ numbers. Many ray tracing applications need to use tens of thousands of surface realizations to get good statistical averages. Thus, the specific answer for brute-force scientific calculations is usually, “A bigger $N$ than you can run on your computer.” This page shows one way to finesse certain calculations so that large-$N$ calculations can be avoided.

For creating visually appealing sea surfaces for rendering into a movie scene, $N=2^{10}=1024,2^{11}=2048,$ or perhaps at most $2^{12}=4096$ is usually adequate. This is because 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. If you are rendering an image, you don’t need to have a finer resolution for the sea surface than maps to an image pixel. For example, if you are going to simulate what a CCD with $1024\times 1024$ pixels would record, and you are looking down onto the sea surface from above at an elevation where one CCD pixel sees a $0.1m\times 0.1m$ patch of the sea surface, then you need to resolve the sea surface with a grid spacing of at most 0.1 m. Modeling a $100m\times 100m$ patch of sea surface with a spatial resolution of $1024\times 1024$ gives a grid spacing of 0.1 m, which would be adequate for the CCD image simulation.

Unfortunately, for accurate scientific calculations of sea surface reflectance, $N$ may have to be much larger. This is because the surface reflectance depends to first order on the surface wave slope, and the slope is strongly affected by the smallest waves with spatial scales down to capillary size of a few millimeters.

The $k$ value of the peak of the Pierson-Moskowitz $k$ spectrum

is easily obtained from setting $d{\mathcal{𝒮}}_{PM}\left(k\right)∕dk=0$ and solving for the value of $k$. This gives

which corresponds to a wavelength of ${\Lambda }_p=2\pi ∕k_p$. Table shows the values of $k_p$ and ${\Lambda }_p$ for a few wind speeds.

If we use $L\approx 2{\Lambda }_p$ for the size of the spatial domain, we will resolve the large gravity waves, which contain most of the wave variance. (Keep in mind that density functions do not have a unique maximum: the location of the maximum of a density function depends on which frequency variable is chosen. See the discussion of this on the page A Common Misconception. Differentiating the Pierson-Moskowitz $\omega$ spectrum

and using the dispersion relation $\omega =\sqrt{gk}$ leads to ${\Lambda }_p=2\pi {\left(1.25∕\beta \right)}^{0.25}g∕U_{19}^2$, which gives 73 m at $U_{19}=10$. However, either version of ${\Lambda }_p$ gives adequate guidance for our present purpose.)

Now consider how many points it takes to almost fully resolve, or account for, the variances of elevation and slope spectra when they are sampled for a DFT. For a specific example, take $U_{10}=10m∕s$ and $L=200m$ and use the omnidirectional, or 1-D, part of the Elfouhaily et al. spectrum presented on page Wave Variance Spectra: Examples. This $\mathcal{𝒮}\left(k\right)$ is given by Eq. (5) and the associated equations on that page. This variance spectrum and the corresponding slope spectrum $k^2\mathcal{𝒮}\left(k\right)$ are plotted as the blue curves in the next 3 figures.

As we have learned on page Wave Variance Spectra: Theory, the expected variance of the sea surface is given by

and the mean-square slope is given by

Such integrals usually must be numerically evaluated with the limits of $k=0$ and $\infty$ being replaced by some small value $k_0$ and some large but finite value $k_{\infty }$. The present example uses $k_0=0.01$ and $k_{\infty }=10^4$, with $10^6$ evaluation points in between these limits. That gives an accurate evaluation of Eqs. (3) and (4) for the $U_{10}=10m{s}^{-1}$ spectra. The results are $⟨z^2⟩=0.4296m^2$ and $mss=0.06011ra{d}^2$.

Now suppose we sample the spectrum using $N=1024$ points, which leads eventually to a sea surface realization $z\left(x_r\right)$ with 1024 points. For $L=200m$ the fundamental frequency is $k_f=2\pi ∕L=0.0314rad∕m$. This point is shown by the left-most red dot with the black center in Fig. 1. For this $L$ and $N$, the Nyquist frequency is $k_{Ny}=2\pi ∕\left(2L∕N\right)=16.085rad∕m$, which is shows by the right-most black and red dot. These two end points and the $N-2$ red vertical lines in between show the discrete points where the variance spectrum is sampled.

The elevation and slope variances that are accounted for in a sampling scheme with $N$ points are given by

and

In the present example, these integrals are $z^2\left(N=1024\right)=0.4219m^2$ and $mss\left(N=1024\right)=0.02584ra{d}^2$. Thus the finite range of the sampled frequencies includes the fraction

of the total variance of the sea surface elevation. However, the corresponding fraction of the sampled slope variance is just $f_S=0.02584∕0.06011=0.430$. Thus $N=1024$ is sampling 98% of the elevation variance but only 43% of the slope variance. This sampling is acceptable if we are interested only in creating a sea surface that looks realistic to the eye. The large gravity waves will be correctly simulated, and that is what counts for visual appearance.

However, if our interest is ray tracing to simulate the optical reflectance and transmission properties of the sea surface, then we must also sample almost all of the slope variance. This is because, to first order, reflectance and transmission of light are governed by the slope of the sea surface, not by its height. If we under sample the slope variance, then the generated surface will be too smooth at the smallest spatial scales, and the computed optical properties will be incorrect. The brute-force approach to sampling more of the slope spectrum is to increase $N$.

Figure 2 shows the sampling points when $N=2^{16}=65536$. Now the sampled points extend into the capillary-wave spatial scale: the smallest resolvable wave is 0.006 m. (Note that the fundamental frequency is fixed by the choice of $L$, which also fixes the spacing of the samples, $\Delta k=k_f$.) The elevation variance integral over the sampled ranges is unchanged (we are still missing a small bit of the variance from the very longest waves with frequencies below $k_f$. However, the slope integral is now $mss\left(N=65536\right)=0.05909ra{d}^2$. Thus we are now sampling a fraction $f_S=0.05909∕0.06011=0.983$ of the slope variance. The resolution of both the elevation and slope spectra are now acceptable for almost any application. Although using a very large $N$ is computationally feasible in one dimension, such a large $N$ is computationally impracticable in two dimensions, when we would need a grid of size $N\times N$. In the present example, this would would require $65536^2=4.3Gbytes$ of storage for just one real array, as well as a corresponding increase in run time for the FFT. Another approach to resolving the slope variance is therefore needed.

As we have just seen, the unsampled frequencies greater than the Nyquist frequency can account for a large part of the slope variance. These frequencies represent the smallest gravity and capillary waves. The amplitudes of these waves are small, so they contribute little to the total wave variance, but their slopes can be large. One way to account for these “missing” waves is to alias their variance into the waves of the highest frequencies. The highest 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. One way to do this is as follows.

We seek an adjusted or re-scaled elevation spectrum

such that the integral of $k^2\widetilde{\mathcal{𝒮}}\left(k\right)$ over the sampled region equals the integral of the true $k^2\mathcal{𝒮}\left(k\right)$ over the full spectral range. Then sampling the re-scaled spectrum $\widetilde{\mathcal{𝒮}}\left(k\right)$ will lead to the same slope variance as would be obtained from sampling the correct spectrum $\mathcal{𝒮}\left(k\right)$ over the entire range of frequencies.

We can choose $L$ so that the low frequencies are well sampled starting at the fundamental frequency $k_f=2\pi ∕L$. There is thus no need to modify the low-frequency part of the variance spectrum, which if done, would adversely affect the total wave variance. We want to modify only the higher frequencies of the sampled region, which contribute the most to the wave slope but little to the wave elevation. A simple approach is to take $\delta \left(k\right)$ to be a linear function of $k$ between the spectrum peak $k_p$ and the highest sampled frequency $k_{Ny}$, and zero elsewhere:

${\delta }_{Ny}$ is a parameter that depends on the spectrum (i.e., the wind speed), the size $L$ of the spatial domain, and the number of samples $N$. We must determine the value of ${\delta }_{Ny}$ so that this $\delta \left(k\right)$ “adds back in” the unresolved slope variance. The added variance will be zero at the peak frequency $k_p$ and largest at the Nyquist frequency. That is, the $\delta \left(k\right)$ function will make only a small change to the variance spectrum at the low frequencies, and the change will be largest near the highest sampled frequencies, which is consistent with the idea that the high frequency waves have the largest slopes. (There is nothing physical about the functional form of Eq. (6); it is simply a single-parameter ad hoc function that works reasonably well in practice. Nonlinear functions could be used, e.g. to concentrate more of the variance into the frequencies closest to the Nyquist frequency. However, those functions could have more than one free parameter to determine and, in any case, the end result is the same: a sea surface that reproduces the height and slope variances of the fully sampled spectrum.)

To determine ${\delta }_{Ny}$ we thus have

The right-most terms of these equations give (after recalling that $\delta \left(k\right)=0$ for $k\le k_p$)

Figure 3 shows example 1-D surface realizations created with $N=1024$ samples and both the true and scaled variance spectra. The upper two panels reproduce Fig. 1, except that the sampled points of the re-scaled spectra are shown in green. It is clear that the $\delta \left(k\right)$ function has added progressively more variance to the higher frequencies. The re-scaled variance spectrum does of course contain somewhat more elevation variance over the sampled region than does the true spectrum. As the green inset $f_E$ number shows, this re-scaling has increased the fraction of sampled/true variance from 0.982 to 1.020. However, the $f_S$ number in the upper-right panel shows that we are now sampling 99.5% of the slope variance, as opposed to 43% for the true spectrum. Having a bit too much total elevation variance is a good tradeoff for being able to model the optically important slope statistics.

The row 2 panel shows random realizations of the 1-D surfaces generated from these two spectra (with the same set of random numbers). The surface elevations $z\left(x_r\right)$ generated using the true spectrum (red line) and the re-scaled spectrum (green line) are almost indistinguishable at the scale of this figure. These significant wave heights for these two surface realizations compare well with the significant wave height of $H_{1∕3}=2.48m$ for a Pierson-Moskowitz spectrum. The row 3 panel shows the surface slopes computed from finite differences of the discrete surface heights. The Cox-Munk along-wind mean square slope is given by $ms{s}_x=0.0316U_{12.5}$, which compares well with the value of 0.032 obtained with the re-scaled spectrum. However, the value obtained from the true spectrum is only 0.022, or 70% of the Cox-Munk value. The bottom two panels replot the 0-10 m sections of the $z$ and $dz∕dx$ figures for better display of the details.

We thus see that when using the $\delta \left(k\right)$ correction to a variance spectrum, we can generate a surface that reproduces, to within a few percent of the theoretical values, both the surface elevation and slope statistics that would be obtained from the underlying true variance spectrum if it were sampled with enough points to fully resolve the elevation and slope variances.