**Page updated:**
April 16, 2020 **Author:** Curtis Mobley

View PDF

# Time-Dependent Surfaces

One ﬁnal step remains: the addition of time dependence to generate a sequence of the sea surface realizations. Many scientiﬁc applications do not need time dependence. This is the case when many independent random realizations of sea surfaces are used for Monte Carlo ray tracing to compute the average optical reﬂectance and transmittance properties of wind-blown sea surfaces. However, time dependence is crucial for applications such as computer-generated imagery for movies.

We already have the needed theory in hand from the page on Spectra to Surfaces: 2D. The fundamental are equations from that page are

$$\begin{array}{llll}\hfill {\u1e91}_{o}\left({k}_{uv}\right)\equiv & \phantom{\rule{1em}{0ex}}\frac{1}{\sqrt{2}}\left[\rho \left({k}_{uv}\right)+i\sigma \left({k}_{uv}\right)\right]\sqrt{\frac{{\Psi}_{1s}\left(k={k}_{uv}\right)}{2}\Delta {k}_{x}\Delta {k}_{y}}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill =& \phantom{\rule{1em}{0ex}}\frac{1}{\sqrt{2}}\left[\rho \left({k}_{uv}\right)+i\sigma \left({k}_{uv}\right)\right]\sqrt{{\Psi}_{2s}\left({k}_{uv}\right)}\phantom{\rule{0.3em}{0ex}},)\phantom{\rule{2em}{0ex}}& \hfill \text{(1)}\end{array}$$ and

These quantities are evaluated as described on that page, but with the evaluations done at particular times $t$. Most importantly, these equations make no simplifying assumptions about the $\pm k$ symmetry of the two-sided variance spectrum ${\Psi}_{2s}\left(k\right)$.

As we have learned, the Fourier transform of a snapshot of the sea surface gives a two-sided variance spectrum with identical values for $-k$ and $+k$. This represents equal amounts of energy propagating in the $-k$ and $+k$ directions, i.e., equal amounts of energy in waves propagating upwind and downwind. Such a situation in nature gives standing waves. Here also, if ${\Psi}_{2s}\left(-k\right)={\Psi}_{2s}\left(k\right)$, the time-dependent surface will be standing waves composed of waves of equal energy propagating the $+k$ and $-k$ directions. To obtain waves propagating downwind, as is the case in a real ocean, we must use an asymmetric two-sided spectrum with ${\Psi}_{2s}\left(-k\right)<<{\Psi}_{2s}\left(+k\right)$, so that almost all of the energy is propagating downwind. Note, however, that we cannot simply set ${\Psi}_{2s}\left(-k\right)=0$, which represents no energy at all propagating upwind. This is because ${\Psi}_{2s}\left(-k\right)=0$ destroys the Hermitian property of Eq. (2). Thus we must conjure up an asymmetric variance spectrum that allows a nonzero (although perhaps negligible) amount of energy to propagate upwind.

An asymmetric two-sided variance spectrum can be constructed using an asymmetric spreading function $\Phi \left(k,\phi \right)$ as described on the previous page, Spreading Function Eﬀects . Spreading functions of the form

$$\Phi \left(k,\phi \right)={C}_{s}{cos}^{2S\left(k\right)}\left(\phi \u22152\right)$$ | (3) |

described there allow some energy to propagate in $-k$ directions, i.e. when $\left|\phi \right|>\pi \u22152$. With this choice of a spreading function, we can let the magnitude of ${\Psi}_{2s}\left(+k\right)$ equal the magnitude of the one-sided variance function ${\Psi}_{1s}\left(+k\right)$, which gives the total variance, because we assume that a negligible amount of the total energy propagates in $-k$ directions. With this assumption, there is no division of ${\Psi}_{1s}\left(+k\right)$ by 2 in the ﬁrst line of Eq. (1). That is, we are starting with a two-sided, asymmetric spectrum, not with a one-sided spectrum based on the assumption of upwind-downwind symmetry.

Once an asymmetric ${\Psi}_{2s}\left(\pm k\right)$ has been deﬁned, we can evaluate Eq. (1) for a set of random numbers $\rho \left({k}_{uv}\right)$ and $\sigma \left({k}_{uv}\right)$. This is done only once. Then to generate a sequence of sea surface realizations for times $t=0,\Delta t,2\Delta t,...$, we multiply the time-independent amplitudes by the time-dependent cosines and sines as was shown previously in

$$\begin{array}{llll}\hfill 2& \phantom{\rule{0.3em}{0ex}}\u1e91\left(u,v,t\right)\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill =& \phantom{\rule{1em}{0ex}}\left[\rho \left(u,v\right)\sqrt{{\Psi}_{2s}\left(u,v\right)}+\rho \left({N}_{x}\phantom{\rule{0.3em}{0ex}}-u,{N}_{y}\phantom{\rule{0.3em}{0ex}}-v\right)\sqrt{{\Psi}_{2s}\left({N}_{x}\phantom{\rule{0.3em}{0ex}}-u,{N}_{y}\phantom{\rule{0.3em}{0ex}}-v\right)}\right]cos\left[\omega \left(k\left(u,v\right)\right)t\right]\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill -& \phantom{\rule{1em}{0ex}}\left[\sigma \left(u,v\right)\sqrt{{\Psi}_{2s}\left(u,v\right)}+\sigma \left({N}_{x}\phantom{\rule{0.3em}{0ex}}-u,{N}_{y}\phantom{\rule{0.3em}{0ex}}-v\right)\sqrt{{\Psi}_{2s}\left({N}_{x}\phantom{\rule{0.3em}{0ex}}-u,{N}_{y}\phantom{\rule{0.3em}{0ex}}-v\right)}\right]sin\left[\omega \left(k\left(u,v\right)\right)t\right]\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill +& \phantom{\rule{1em}{0ex}}i\left\{\right.\left[\rho \left(u,v\right)\sqrt{{\Psi}_{2s}\left(u,v\right)}-\rho \left({N}_{x}\phantom{\rule{0.3em}{0ex}}-u,{N}_{y}\phantom{\rule{0.3em}{0ex}}-v\right)\sqrt{{\Psi}_{2s}\left({N}_{x}\phantom{\rule{0.3em}{0ex}}-u,{N}_{y}\phantom{\rule{0.3em}{0ex}}-v\right)}\right]sin\left[\omega \left(k\left(u,v\right)\right)t\right]\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill +& \phantom{\rule{1em}{0ex}}\left[\sigma \left(u,v\right)\sqrt{{\Psi}_{2s}\left(u,v\right)}-\sigma \left({N}_{x}\phantom{\rule{0.3em}{0ex}}-u,{N}_{y}\phantom{\rule{0.3em}{0ex}}-v\right)\sqrt{{\Psi}_{2s}\left({N}_{x}\phantom{\rule{0.3em}{0ex}}-u,{N}_{y}\phantom{\rule{0.3em}{0ex}}-v\right)}\right]cos\left[\omega \left(k\left(u,v\right)\right)t\right]\left\}\right.\phantom{\rule{0.3em}{0ex}}.\phantom{\rule{2em}{0ex}}& \hfill \text{(4)}\end{array}$$

We thus obtain the amplitudes $\u1e91\left({k}_{uv},t\right)$ at the current time $t$. The inverse Fourier transform of this $\u1e91\left({k}_{uv},t\right)$ gives the sea surface $z\left({x}_{rs}\right)$ at time $t$.

The physics (or lack thereof) of this process is simple. We start with a realization of the sea surface at time zero. This surface contains waves of many discrete frequencies ${k}_{uv}$ traveling in all directions ${\phi}_{uv}$. Then to get the surface at any later time $t$, we simply propagate the sinusoidal waves of each frequency ${k}_{uv}$ in their original direction of travel through a phase change determined by the time step and the dispersion relation $\omega \left({k}_{uv}\right)$. For deep-water gravity waves, the dispersion relation is $\omega \left({k}_{uv}\right)=\sqrt{g{k}_{uv}}$. It thus visually appears that the waves are propagating and the sea surface proﬁle is changing with time. However, this Fourier transform technique is really just moving a collection of independent, non-interacting sinusoids through frequency-dependent phase angles to create a new surface realization from the sum of the phase-shifted sinusoids. In a real ocean, waves of diﬀerence frequencies can interact with each other (redistribute energy among frequencies) in highly complex and nonlinear ways, so that the sea surface statistics may not be time-independent. This is, in particular, how little waves grow to big waves when the wind begins to blow over a level surface. Modeling the nonlinear evolution of a sea surface is beyond the abilities of Fourier transform techniques which are, at heart, just a mathematical artiﬁce based on a time-independent directional variance density spectrum.

Figure 1 shows a sequence of surface realizations for a $10\phantom{\rule{1em}{0ex}}m\phantom{\rule{1em}{0ex}}{s}^{-1}$ wind speed and a spatial grid of size ${L}_{x}\times {L}_{y}=100\times 100\phantom{\rule{1em}{0ex}}m$, and a grid resolution of ${N}_{x}\times {N}_{y}=256\times 256$. A time series of images made with such a course grid could be useful for some non-scientiﬁc purposes.

We have seen that the spatial pattern of a sea surface generated by a Fourier transform is periodic. This is convenient for tiling a small computed region into a visually acceptable larger region. When time dependence is included and a ﬁnite-length time series of images is created, the sequence of images is not periodic in time because the various sinusoids comprising the surface do not have a common period. As pointed out by Tessendorf (2004), this can be remedied by ”quantizing” the temporal frequency $\omega \left({k}_{uv}\right)$ as follows.

Let ${T}_{r}$ be the length of time over which the time sequence of surface realizations should exactly repeat. The number of frames ${N}_{t}$ in the video loop determines the time step between frames, $\Delta t={T}_{r}\u2215{N}_{t}$. For a smoothly moving sea surface, ${N}_{t}$ must be large enough that $\Delta t$ is less than about 0.1 s. Deﬁne ${\omega}_{o}\equiv 2\pi \u2215{T}_{r}$. For deep-water gravity waves, the true temporal frequency $\omega \left({k}_{uv}\right)=\sqrt{g{k}_{uv}}$ is replaced by

$$\stackrel{\u0303}{\omega}\left({k}_{uv}\right)=\u230a\frac{\sqrt{g{k}_{uv}}}{{\omega}_{o}}\u230b{\omega}_{o}\phantom{\rule{0.3em}{0ex}},$$ | (5) |

where $\lfloor x\rfloor $ converts a real number $x$ into its integer part (e.g., 15.23 is converted to 15). This operation slightly alters the temporal frequency of each wave component so that each component returns to exactly its initial shape after time ${T}_{r}$. A video loop can then be created from the sequence of surfaces, and the loop will match perfectly when the frame for time ${T}_{r}-\Delta t$ goes to the frame for time ${T}_{r}$, which is the same surface as time 0, after which the surfaces repeat. This adjustment to $\omega $ is greatest for the lowest frequencies, but the adjustment becomes smaller and smaller as the repeat time becomes larger and larger. It is thus easy to create a time-dependent small area of sea surface that can be tiled in both space and time to create an image of a larger sea surface over a longer time. This is good enough to fool movie-goers.

In order to employ a re-scaled variance spectrum as described on page Numerical Resolution, determine the value of the ${\delta}_{Ny}$ parameter using the omnidirectional variance spectrum, as seen in Eq. (7) of that page. Then apply the $\delta \left(k\right)$ correction to the directional spectrum $\Psi \left({k}_{x},{k}_{y}\right)$ with $k={k}_{uv}$ for each $\left({k}_{x}\left(u\right),{k}_{y}\left(v\right)\right)$ point of the rectangular grid.

An example of a 20-second (simulated time) video loop created using all of these tricks can be seen at video loop. This shows a time series of 2D sea surface realizations generated using the variance spectrum of Elfouhaily et al. (1997) with a frequency-dependent cosine-2S spreading function. This omnidirectional elevation variance spectrum was adjusted as described on page Numerical Resolution at the higher spatial frequencies so that unresolved slope variance (from frequencies higher than the ${k}_{x}$ and ${k}_{y}$ Nyquist frequencies resolvable by a 512 by 512 DFT grid) is fully resolved. The true dispersion relation $\omega =\sqrt{gk}$ was quantized for each discrete spatial frequency so that the surface is exactly periodic after 20 seconds. Note in particular that the signiﬁcant wave heights ${H}_{1\u22153}$ are slightly greater on average than the value of 2.25 m predicted by a Pierson-Moskowitz spectrum, which has somewhat less energy than the Elfouhaily et al. spectrum used here. Note also that the along-wind ($ms{s}_{x}$) and cross-wind ($ms{s}_{y}$) mean square slopes are very close to what is expected by Cox-Munk statistics: 0.031 and 0.019 respectively. The sequence of plots used to create the video loop was created by IDL routine cgAnimate˙2D˙SeaSurface.pro. The plots created by the IDL routine were then combined into the video using VideoMach.