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

Spectra to Surfaces: 2D

Now, again, we face the reverse process: start with a two-dimensional, one-sided variance density spectrum and generate a random 2-D realization of a sea surface. The first task is to conjure up a two-dimensional variance density spectrum, which requires some effort.


Let Ψ(k) = Ψ(kx,ky) denote a 2-D elevation variance spectrum as defined on the Wave Variance Spectra: Theory page. This spectrum has units of m2(radm)2. By definition, the integral of Ψ(kx,ky) over all frequencies gives the elevation variance:

var{z} = z2 =Ψ(k x,ky)dkxdky.

As in the 1-D case on page Spectra to Surfaces: 1D, ... indicates expectation or ensemble average over many measurements of the sea surface.

In the discrete case, the kuv values coming out of = 𝒟{z} point both “downwind” (positive kx(u) values) and “upwind” (negative kx(u) values), i.e. there are both positive and negative frequencies represented in the “two-sided” variance spectrum, denoted as Ψ2s[u,v],u = (Nx2 1),...,Ny2,v = (Ny2 1),...,Ny2 as above. In general, Ψ2s(kuv)Ψ2s(+kuv) because more energy propagates downwind than upwind at a given frequency. kx(Nx2) is the Nyquist frequency for waves propagating in the x direction; ky(Ny2) is the Nyquist frequency for waves propagating in the y direction. Ψ2s(0, 0) is the variance at zero frequencies, i.e., the variance of the mean sea surface; this term is normally set to 0 so that the mean sea surface is at height 0.

The 2-D equivalents of Eqs. (1) and (2) on page Spectra to Surfaces: 1D are

o(kuv) 1 2[ρ(kuv) + iσ(kuv)]Ψ1s (k = k uv ) 2 ΔkxΔky (1) = 1 2[ρ(kuv) + iσ(kuv)]Ψ2s (k uv ). (2)

Here, as before, ρ(kuv) and σ(kuv) are independent 𝒩(0, 1) random variables, with a different random variable drawn for each kuv value. As in the 1-D case, the expected value of o(kuv) gives back whatever variance spectrum is used for Ψ2s(kuv), but is not Hermitian. As before, the notation in these equations distinguishes between the value of the continuous spectral density function evaluated at a discrete frequency value, Ψ1s(k = kuv), and the discrete variance point function, Ψ2s(kuv).

We must define random Hermitian amplitudes for use in the inverse Fourier transform. Looking ahead to the page on Time Dependent Surfaces, which extends the results of this page to time-dependent surfaces, define the time-dependent spectral amplitude

(kuv,t) 1 2 o(kuv)eiωuvt + o(k uv)e+iωuvt . (3)

This function is clearly Hermitian, so the inverse DFT applied to (kuv,t) will give a real-valued z(xrs,t). Recall that a function of the form cos(kx ωt) gives a wave propagating in the + x direction. The corresponding o(kuv)eiωt in this equation gives a wave propagating in the + k direction, which is in the downwind half plane of all directions. In general the temporal angular frequency ωuv is a function of the spatial frequency kuv. For example, for deep-water gravity waves, ωuv2 = gk uv.

For simplicity of notation, let us momentarily drop the rs and uv subscripts on the discrete variables. The (k,t) of Eq. (3) is also consistent with the variance spectrum:

|(k,t)|2 =(k,t)(k,t) = 1 2 1 2[ρ(k) + iσ(k)]Ψ2s (k )eiωt + 1 2[ρ(k) iσ(k)]Ψ2s (k )eiωt × 1 2 1 2[ρ(k) iσ(k)]Ψ2s (k )eiωt + 1 2[ρ(k) + iσ(k)]Ψ2s (k )eiωt = 1 4 ρ2(k) iρ(k)σ(k) + iσ(k)ρ(k) + σ2(k) Ψ 2s(k)+ ρ(k)ρ(k) + iρ(k)σ(k) + iσ(k)ρ(k) σ(k)σ(k) Ψ2s (k )Ψ2s (k )ei2ωt+ ρ(k)ρ(k) iρ(k)σ(k) iσ(k)ρ(k) σ(k)σ(k) Ψ2s (k )Ψ2s (k )ei2ωt+ ρ2(k) + iρ(k)σ(k) iσ(k)ρ(k) + σ2(k) Ψ 2s(k) = 1 2[Ψ2s(k) + Ψ2s(k)].

Here we have noted that all terms like ρ(k)ρ(k) are zero because of the independence of the random variables for different k values, as are terms like ρ(k)σ(k). The remaining term is the total variance associated with waves propagating in the downwind and upwind directions at the spatial frequency of magnitude k. It should be noted that this term is independent of time even though the waves z(x,y,t) depend on time. This is because the total variance (or energy) of the wave field is the same at all times, even though the exact shape of the sea surface varies with time.

If only a “snapshot” of the sea surface at one time is available, it is not possible to resolve how much of the total variance is associated with waves propagating in direction k compared to the opposite direction k. The forward DFT, (k,t) = 𝒟{z(x,y,t)}, and (k,t)(k,t) then gives Ψ2s(k) = Ψ2s(k), in which case the last equation reduces to

(k,t)(k,t) = Ψ 2s(k).

In any case, the amplitudes defined by Eq. (3) are Hermitian, so that the real part of the inverse Fourier transform Z(x,t) = 𝒟1{(k,t)} gives a real-valued sea surface z(x,t). That sea surface is consistent with the variance spectrum Ψ2s(k) at every time t. Wave variance (energy) is thus conserved in a round-trip calculation from variance spectrum to sea surface and back to variance spectrum. In the time-dependent case, if Ψ2s(k) > Ψ2s(k), then more variance will be contained in waves propagating in the + k direction than in the k direction. This is all that we can ask from Fourier transform techniques.

Although Eqs. (2) and (3) are compact representations of the random spectral amplitudes, the actual evaluation of these equations in a computer program warrants further examination. In particular, the Nyquist frequencies are always special cases because there is only a positive Nyquist frequency, kxNy = Nx 2 2π Lx at array element Nx 1, with a corresponding equation for the Nyquist frequency in the y direction. Writing e±iωt = cos[ω(k)t] ± i sin[ω(k)t] and expanding Eq. (3) gives

2(k,t) = ρ(k)Ψ2s (k ) + ρ(k)Ψ2s (k ) cos[ω(k)t] σ(k)Ψ2s (k ) + σ(k)Ψ2s (k ) sin[ω(k)t] + i ρ(k)Ψ2s (k ) ρ(k)Ψ2s (k ) sin[ω(k)t] + σ(k)Ψ2s (k ) σ(k)Ψ2s (k ) cos[ω(k)t]. (4)

These terms are all Nx × Ny arrays, but note that terms like ρ(k)Ψ2s (k ) represent element-by-element multiplications, not matrix multiplications.

For a particular array element (kx(u),ky(v),t) = (u,v,t), and using the indexing relation k(u) = k(N 2 u),u = 0,...,N 2 for frequencies written in math order, we can write

2 (u,v,t) = ρ(u,v)Ψ2s (u, v) + ρ(Nx 2 u,Ny 2 v)Ψ2s (Nx 2 u, Ny 2 v) cos[ω(k(u,v))t] σ(u,v)Ψ2s (u, v) + σ(Nx 2 u,Ny 2 v)Ψ2s (Nx 2 u, Ny 2 v) sin[ω(k(u,v))t] + i ρ(u,v)Ψ2s (u, v) ρ(Nx 2 u,Ny 2 v)Ψ2s (Nx 2 u, Ny 2 v) sin[ω(k(u,v))t] + σ(u,v)Ψ2s (u, v) σ(Nx 2 u,Ny 2 v)Ψ2s (Nx 2 u, Ny 2 v) cos[ω(k(u,v))t]. (5)

This equation allows for efficient computation within loops over array elements. In particular, the code can loop over the non-positive kx(u) values, u = 0,...,Nx2, and over all ky(v) values, v = 0,...,Ny 2 to evaluate the amplitudes for all non-Nyquist frequencies. The positive kx(u) values are then obtained from the negative kx(u) values by Hermitian symmetry. The Nyquist frequencies are evaluated by an equation of the same form, but with one or the other index held fixed (e.g., u = Nx 1 while v = 0,...,Ny 2). The amplitude (kx = 0,ky = 0) at array element (u,v) = (Nx2 1,Ny2 1) is usually set to zero, corresponding to the mean sea level being set to 0. For generation of time-independent surfaces, we can set t = 0 so that the cosines are one and the sines are zero, which cuts the number of terms to be evaluated in half.

If the frequencies are in the FFT order, then the last equation has the same general form, but the indexing that expresses Hermitian symmetry is different: k(u) = k(N u),u = 1,...,N 1, with k = 0 being a special case. The equation corresponding to Eq. (5) is then

2 (u,v,t) = ρ(u,v)Ψ2s (u, v) + ρ(Nx u,Ny v)Ψ2s (Nx u, Ny v) cos[ω(k(u,v))t] σ(u,v)Ψ2s (u, v) + σ(Nx u,Ny v)Ψ2s (Nx u, Ny v) sin[ω(k(u,v))t] + i ρ(u,v)Ψ2s (u, v) ρ(Nx u,Ny v)Ψ2s (Nx u, Ny v) sin[ω(k(u,v))t] + σ(u,v)Ψ2s (u, v) σ(Nx u,Ny v)Ψ2s (Nx u, Ny v) cos[ω(k(u,v))t]. (6)

However, this equation does not explicitly show the special cases. Let zhat[u,v] be the array of (k,t) = (kx(u),ky(v),t) values at a particular time t. r[u,v] and s[u,v] are the arrays of random numbers ρ(u,v) and σ(u,v), respectively. Psiroot[u,v] is 1 2Ψ2s (kx (u), ky (v)) (incorporating the 2 seen on the left-hand side of Eq. (6)). With other obvious definitions, Fig. 1 shows the pseudo code to evaluate Eq. (6) at a particular time t.


Figure 1: Pseudocode for looping over non-zero frequencies.

Figure 2 shows the pseudocode for looping over all ky values for the special case of kx = 0 at frequency array index u = 0. Note that for kx = 0, the ± ky values are complex conjugates.


Figure 2: Pseudocode for looping over all ky values for the special case of kx = 0 at frequency array index u = 0.

Figure 3 shows the pseudocode for looping over all kx values for the special cases of ky = 0 at frequency array index v = 0. Note that for ky = 0, the ± kx values are complex conjugates.


Figure 3: Pseudocode for looping over all kx values for the special cases of ky = 0 at frequency array index v = 0.

Finally, set the (kx,ky) = (0, 0) value to 0, which sets the mean sea level to zero: zhat[0,0] = COMPLEX(0.0, 0.0).

Array zhat[u,v] = (kx(u),ky(v)) as just defined is Hermitian and has the frequencies in FFT order ready for input to an FFT routine.

Showing this level of detail may seem tedious, but it is absolutely critical that the (k,t) array be properly computed down to the last array element. Any error will show up in the generated sea surface as either a non-zero imaginary part of the complex array Z[r,s] = 𝒟{zhat[u,v]} (easy to detect) or an incorrect sea surface elevation (often much harder to detect).

Example: A Two-Dimensional Sea Surface

As a specific example of the above algorithm, consider the following. Let us use a course grid sampling of Nx × Ny = 64 × 64 points, which makes it easier to see certain features in the associated plots. The physical region to be simulated is Lx × Ly = 100 × 100m. The two-dimensional, one-sided variance spectrum of Elfouhaily et al. (1997) is used. Thus we have

Ψ(kx,ky) = 1 k𝒮(k)Φ(k,φ) Ψ(k,φ) = Ψ k = kx 2 + ky 2,φ = tan 1(k ykx) .

where 𝒮(k) is the omnidirectional spectrum and Φ(k,φ) is the nondimensional spreading function, as defined on page Wave Variance Spectra: Examples.

The upper left plot of Fig. 4 shows this spectrum for a wind speed of 5ms1 and a fully developed sea state. The contours are of Ψ1s(kx,ky) evaluated at the discrete grid points; the 1 2ΔkxΔky factor seen in Eq. (2) has not yet been applied to create a discrete two-sided spectrum. The line along ky = 0 corresponds to the spectrum plotted in the third figure on page Wave Variance Spectra: Examples. We see in both plots that for a 5ms1 wind, the spectrum peaks near 0.2 0.3radm.


Figure 4: Example two-dimensional sea surface generated from a 2-D, one-sided variance spectrum. The resolution is only 64 × 64 grid points, so as to make the features of the underlying spectrum and the Fourier amplitudes easier to see. Generated by IDL routine

The plots of the real and imaginary parts of (kx(u),ky(v)) show that most of the variance is at low frequencies and that the amplitudes have the Hermitian symmetry illustrated in the two figures of the previous page. The lower right panel of the plot shows a contour plot of the sea surface generated from the inverse FFT of the amplitudes. The significant wave height for this surface realization is 0.60 m, in good agreement with the expected value given above for the Pierson-Moskowitz spectrum, which is similar to the Elfouhaily et al. spectrum in the gravity wave region.

Comments for Spectra to Surfaces: 2D:

Loading Conversation