Page updated: March 22, 2021
Author: Curtis Mobley
View PDF

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 216 × 216 = 65536 × 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 232 = 4.3 × 109 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 = 210 = 1024, 211 = 2048, or perhaps at most 212 = 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 × 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 × 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 × 100m patch of sea surface with a spatial resolution of 1024 × 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

𝒮PM(k) = α 2k3 exp β g k2 1 U194 [m2(radm)], (1)

is easily obtained from setting d𝒮PM(k)dk = 0 and solving for the value of k. This gives

kp = 2β 3 0.5 g U192,

which corresponds to a wavelength of Λp = 2πkp. Table shows the values of kp and Λp for a few wind speeds.


5 0.276 23
10 0.069 91
15 0.031 205
20 0.017 364

Table 1: Spatial frequencies and wavelengths corresponding to peak variance for the Pierson-Moskowitz spectrum 𝒮PM(k).

If we use L 2Λ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 ω spectrum

𝒮PM(ω) = αg2 ω5 exp β g ωU19 4 [m2(rads)]. (2)

and using the dispersion relation ω = gk leads to Λp = 2π(1.25β)0.25gU 192, which gives 73 m at U19 = 10. However, either version of Λ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 U10 = 10ms 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 𝒮(k) is given by Eq. (5) and the associated equations on that page. This variance spectrum and the corresponding slope spectrum k2𝒮(k) 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

z2 =0𝒮(k)dk, (3)

and the mean-square slope is given by

σ2 = mss =0k2𝒮(k)dk. (4)

Such integrals usually must be numerically evaluated with the limits of k = 0 and being replaced by some small value k0 and some large but finite value k. The present example uses k0 = 0.01 and k = 104, with 106 evaluation points in between these limits. That gives an accurate evaluation of Eqs. (3) and (4) for the U10 = 10ms1 spectra. The results are z2 = 0.4296m2 and mss = 0.06011rad2.

Now suppose we sample the spectrum using N = 1024 points, which leads eventually to a sea surface realization z(xr) with 1024 points. For L = 200m the fundamental frequency is kf = 2πL = 0.0314radm. 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 kNy = 2π(2LN) = 16.085radm, 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.


Figure 1: Example sampling of elevation and slope spectra for N = 1024. The red dots and lines show the regions of the elevation and slope spectra sampled using 1024 points. The elevation spectrum is adequately sampled, but the slope spectrum is under sampled.

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

z2(N) =kfkNy 𝒮(k)dk,


mss(N) =kfkNy k2𝒮(k)dk.

In the present example, these integrals are z2(N = 1024) = 0.4219m2 and mss(N = 1024) = 0.02584rad2. Thus the finite range of the sampled frequencies includes the fraction

fE z2(N) z2 = 0.4219 0.4296 = 0.982

of the total variance of the sea surface elevation. However, the corresponding fraction of the sampled slope variance is just fS = 0.025840.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 = 216 = 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, Δk = kf.) 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 kf. However, the slope integral is now mss(N = 65536) = 0.05909rad2. Thus we are now sampling a fraction fS = 0.059090.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 × N. In the present example, this would would require 655362 = 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.


Figure 2: Example sampling of elevation and slope spectra for N = 65536. The red dots and lines show the sampled regions of the elevation and slope spectra. Both the elevation and the slope spectra are adequately sampled.

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

𝒮̃(k) [1 + δ(k)]𝒮(k) (5)

such that the integral of k2𝒮̃(k) over the sampled region equals the integral of the true k2𝒮(k) over the full spectral range. Then sampling the re-scaled spectrum 𝒮̃(k) will lead to the same slope variance as would be obtained from sampling the correct spectrum 𝒮(k) over the entire range of frequencies.

We can choose L so that the low frequencies are well sampled starting at the fundamental frequency kf = 2π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 δ(k) to be a linear function of k between the spectrum peak kp and the highest sampled frequency kNy, and zero elsewhere:

δ(k) 0 if k kp δNy kkp kNykpif k > kp (6)

δ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 δNy so that this δ(k) “adds back in” the unresolved slope variance. The added variance will be zero at the peak frequency kp and largest at the Nyquist frequency. That is, the δ(k) function will make only a small change to the variance spectrum at the low frequencies, and the change will be largest near the highest sampled frequencies, which is consistent with the idea that the high frequency waves have the largest slopes. (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 δNy we thus have

mss =0k2𝒮(k)dk k0k k2𝒮(k)dk =k0kNy k2𝒮(k)dk +kNyk k2𝒮(k)dk k0kNy k2𝒮̃(k)dk =k0kNy k2𝒮(k)dk + δ NykpkNy k2 k kp kNy kp 𝒮(k)dk.

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

δNy = kNykk2𝒮(k)dk kpkNyk2 kkp kNykp 𝒮(k)dk. (7)

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 δ(k) 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 fE number shows, this re-scaling has increased the fraction of sampled/true variance from 0.982 to 1.020. However, the fS 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(xr) 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 H13 = 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 mssx = 0.0316U12.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 dzdx figures for better display of the details.


Figure 3: Example 1-D surfaces and slopes created using true and adjusted variance spectra. See the preceding text for discussion.

We thus see that when using the δ(k) 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.

Comments for Numerical Resolution:

Loading Conversation