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.

Theory

Let $\Psi \left(k\right)=\Psi \left(k_x,k_y\right)$ denote a 2-D elevation variance spectrum as defined on the Wave Variance Spectra: Theory page. This spectrum has units of $m^2∕{\left(rad∕m\right)}^2$. By definition, the integral of $\Psi \left(k_x,k_y\right)$ over all frequencies gives the elevation variance:

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 $k_{uv}$ values coming out of $ẑ=\mathcal{𝒟}\left\{z\right\}$ point both “downwind” (positive $k_x\left(u\right)$ values) and “upwind” (negative $k_x\left(u\right)$ values), i.e. there are both positive and negative frequencies represented in the “two-sided” variance spectrum, denoted as ${\Psi }_{2s}\left[u,v\right],u=-\left(N_x∕2-1\right),...,N_y∕2,v=-\left(N_y∕2-1\right),...,N_y∕2$ as above. In general, ${\Psi }_{2s}\left(-k_{uv}\right)\ne {\Psi }_{2s}\left(+k_{uv}\right)$ because more energy propagates downwind than upwind at a given frequency. $k_x\left(N_x∕2\right)$ is the Nyquist frequency for waves propagating in the $x$ direction; $k_y\left(N_y∕2\right)$ is the Nyquist frequency for waves propagating in the $y$ direction. ${\Psi }_{2s}\left(0,0\right)$ 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

Here, as before, $\rho \left(k_{uv}\right)$ and $\sigma \left(k_{uv}\right)$ are independent $\mathcal{𝒩}\left(0,1\right)$ random variables, with a different random variable drawn for each $k_{uv}$ value. As in the 1-D case, the expected value of ${ẑ}_o\left(k_{uv}\right)$ gives back whatever variance spectrum is used for ${\Psi }_{2s}\left(k_{uv}\right)$, 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, ${\Psi }_{1s}\left(k=k_{uv}\right)$, and the discrete variance point function, ${\Psi }_{2s}\left(k_{uv}\right)$.

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

This function is clearly Hermitian, so the inverse DFT applied to $ẑ\left(k_{uv},t\right)$ will give a real-valued $z\left(x_{rs},t\right)$. Recall that a function of the form $cos\left(kx-\omega t\right)$ gives a wave propagating in the $+x$ direction. The corresponding ${ẑ}_o\left(k_{uv}\right)e^{-i\omega 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 ${\omega }_{uv}$ is a function of the spatial frequency $k_{uv}$. For example, for deep-water gravity waves, ${\omega }_{uv}^2=g{k}_{uv}$.

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

Here we have noted that all terms like $⟨\rho \left(k\right)\rho \left(-k\right)⟩$ are zero because of the independence of the random variables for different $k$ values, as are terms like $⟨\rho \left(k\right)\sigma \left(k\right)⟩$. 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\left(x,y,t\right)$ 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, $ẑ\left(k,t\right)=\mathcal{𝒟}\left\{z\left(x,y,t\right)\right\}$, and $ẑ\left(k,t\right){ẑ}^{\ast }\left(k,t\right)$ then gives ${\Psi }_{2s}\left(-k\right)={\Psi }_{2s}\left(k\right)$, in which case the last equation reduces to

In any case, the amplitudes defined by Eq. (3) are Hermitian, so that the real part of the inverse Fourier transform $Z\left(x,t\right)={\mathcal{𝒟}}^{-1}\left\{ẑ\left(k,t\right)\right\}$ gives a real-valued sea surface $z\left(x,t\right)$. That sea surface is consistent with the variance spectrum ${\Psi }_{2s}\left(k\right)$ 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 ${\Psi }_{2s}\left(k\right)>{\Psi }_{2s}\left(-k\right)$, 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, $k_x^{Ny}=\frac{N_x}{2}\frac{2\pi }{L_x}$ at array element $N_x-1$, with a corresponding equation for the Nyquist frequency in the $y$ direction. Writing $e^{\pm i\omega t}=cos\left[\omega \left(k\right)t\right]\pm isin\left[\omega \left(k\right)t\right]$ and expanding Eq. (3) gives

These terms are all $N_x\times N_y$ arrays, but note that terms like $\rho \left(k\right)\sqrt{{\Psi }_{2s}\left(k\right)}$ represent element-by-element multiplications, not matrix multiplications.

For a particular array element $ẑ\left(k_x\left(u\right),k_y\left(v\right),t\right)=ẑ\left(u,v,t\right)$, and using the indexing relation $k\left(u\right)=-k\left(N-2-u\right),u=0,...,N-2$ for frequencies written in math order, we can write

This equation allows for efficient computation within loops over array elements. In particular, the code can loop over the non-positive $k_x\left(u\right)$ values, $u=0,...,N_x∕2$, and over all $k_y\left(v\right)$ values, $v=0,...,N_y-2$ to evaluate the amplitudes for all non-Nyquist frequencies. The positive $k_x\left(u\right)$ values are then obtained from the negative $k_x\left(u\right)$ 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=N_x-1$ while $v=0,...,N_y-2$). The amplitude $ẑ\left(k_x=0,k_y=0\right)$ at array element $\left(u,v\right)=\left(N_x∕2-1,N_y∕2-1\right)$ 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\left(u\right)=-k\left(N-u\right),u=1,...,N-1$, with $k=0$ being a special case. The equation corresponding to Eq. (5) is then

However, this equation does not explicitly show the special cases. Let zhat[u,v] be the array of $ẑ\left(k,t\right)=ẑ\left(k_x\left(u\right),k_y\left(v\right),t\right)$ values at a particular time $t$. r[u,v] and s[u,v] are the arrays of random numbers $\rho \left(u,v\right)$ and $\sigma \left(u,v\right)$, respectively. Psiroot[u,v] is $\frac{1}{2}\sqrt{{\Psi }_{2s}\left(k_x\left(u\right),k_y\left(v\right)\right)}$ (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 2 shows the pseudocode for looping over all $k_y$ values for the special case of $k_x=0$ at frequency array index $u=0$. Note that for $k_x=0$, the $\pm k_y$ values are complex conjugates.

Figure 3 shows the pseudocode for looping over all $k_x$ values for the special cases of $k_y=0$ at frequency array index $v=0$. Note that for $k_y=0$, the $\pm k_x$ values are complex conjugates.

Finally, set the $\left(k_x,k_y\right)=\left(0,0\right)$ value to 0, which sets the mean sea level to zero: zhat[0,0] = COMPLEX(0.0, 0.0).

Array zhat[u,v] = $ẑ\left(k_x\left(u\right),k_y\left(v\right)\right)$ 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 $ẑ\left(k,t\right)$ 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\left[r,s\right]=\mathcal{𝒟}\left\{zhat\left[u,v\right]\right\}$ (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 coarse grid sampling of $N_x\times N_y=64\times 64$ points, which makes it easier to see certain features in the associated plots. The physical region to be simulated is $L_x\times L_y=100\times 100m$. The two-dimensional, one-sided variance spectrum of Elfouhaily et al. (1997) is used. Thus we have

where $\mathcal{𝒮}\left(k\right)$ is the omnidirectional spectrum and $\Phi \left(k,\phi \right)$ 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 $5m{s}^{-1}$ and a fully developed sea state. The contours are of ${\Psi }_{1s}\left(k_x,k_y\right)$ evaluated at the discrete grid points; the $\frac{1}{2}\Delta k_x\Delta k_y$ factor seen in Eq. (2) has not yet been applied to create a discrete two-sided spectrum. The line along $k_y=0$ corresponds to the spectrum plotted in the third figure on page Wave Variance Spectra: Examples. We see in both plots that for a $5m{s}^{-1}$ wind, the spectrum peaks near $0.2-0.3rad∕m$.

The plots of the real and imaginary parts of $ẑ\left(k_x\left(u\right),k_y\left(v\right)\right)$ 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.