Gerstner waves

A sine wave moves every point of the water straight up and down. A Gerstner wave moves every point in a circle. That one change sharpens the crests, flattens the troughs, and turns a rubber sheet into something that reads as swell.

František Josef Gerstner published the solution in 1802. It is an exact solution of the nonlinear equations for deep-water waves, it has a closed form, and it costs a handful of sin and cos per vertex. That combination is why it has been the workhorse of real-time water since GPU Gems Chapter 1 in 2004, and it is still the first thing to reach for before an FFT ocean.

This tutorial builds the demo above from scratch. We'll derive the displacement, get exact normals without finite differences, sum several waves safely, find foam from the Jacobian, and write the whole thing in GLSL for WebGL2 and WGSL for WebGPU. Every shader listing on this page is pulled from the source strings that are running in the header, so what you read is what you see.

Start in one dimension. The familiar travelling sine wave only moves points vertically:

Gerstner instead gives each particle a rest position x_0 and lets it orbit around that position:

At the crest (\sin\theta = 1) particles are moving forward with the wave and bunch together; in the trough they move backward and spread apart. You can see the bunching directly from how much the horizontal spacing is squeezed:

When ka reaches 1 the spacing at the crest hits zero and the crest becomes a cusp. Past 1 it goes negative and the curve loops over itself. Push the slider past 1 to see it.

Three numbers describe a wave, and choosing them well saves a lot of slider-fiddling later.

Wavelength L gives the wavenumber k. Speed should not be a free parameter at all. In deep water, gravity waves obey the dispersion relation \omega^2 = gk, so the phase speed follows from the wavelength. Long waves outrun short ones, which is a big part of why a sum of waves never looks like a looping texture.

Steepness s = ka replaces amplitude as the thing you set. It is dimensionless, it controls the crest shape directly, and the amplitude scales with wavelength for free: a 60 m swell and a 6 m chop at the same steepness look like the same wave at different scales.

Give the wave a unit direction \mathbf{D} = (D_x, D_z) on the horizontal plane. The phase now depends on how far along that direction a point sits, and the orbit's horizontal part points along \mathbf{D}:

Here \mathbf{p} = (p_x, p_z) is the vertex's rest position on a flat grid and y is up. In the shaders \mathbf{D} is a vec2, so D_z is written d.y.

Because P is a closed-form function of the grid coordinates, its tangent vectors are just its partial derivatives. No neighbour sampling, no finite-difference noise:

The order of the cross product matters: with y up and a flat sheet, T = (1,0,0) and B = (0,0,1), and B \times T = (0,1,0) points up. Notice that every a became an s: differentiating \theta brings out a factor of k, and ka = s.

One wave looks like a corrugated roof. The sea is a sum. Displacements add, and because differentiation is linear the tangent contributions add too, on top of the flat sheet's (1,0,0) and (0,0,1):

The catch: where several crests line up, their compressions add. If the total steepness stays within budget, the surface cannot fold, no matter how the crests align:

The panel above the article tracks this sum for you. Exceeding it doesn't guarantee a loop; it only means one becomes possible where enough crests meet.

Foam from the Jacobian

The horizontal part of the mapping from grid to surface has a Jacobian determinant:

It is 1 on a flat sheet, drops below 1 where water is squeezed together, and reaches 0 exactly where the surface folds. Compressed water at a sharp crest is where whitecaps form, so J is a physically motivated foam mask that comes out of numbers we've already computed. The demo whitens everything below about 0.5.

Choosing waves that look like a sea

Pick one dominant swell direction and scatter the others within roughly 30 to 60 degrees of it. Use wavelengths with awkward ratios (60, 31, 17, 9.5 in the demo) so their crests don't realign on a short period. Give longer waves more steepness than short ones.

The vertex buffer holds nothing but 2D rest positions. Height, horizontal offset, normal, and foam all come from the shader, so the same buffer serves both APIs and never has to be updated.

A uniform grid wastes vertices far from the camera and starves the near field. The demo squeezes the grid toward the centre with x = E\,u\,(0.25 + 0.75u^2) for u \in [-1, 1], giving about 0.3 m spacing under the camera and about 3 m at the 250 m edge. The rule of thumb for any wave you keep: at least four to eight vertices per wavelength wherever it's visible, or it will alias into crawling noise.

The vertex shader is the heart of the tutorial. gerstner() returns one wave's displacement and accumulates its tangent contributions through inout parameters; main() sums four of them, packed as vec4(direction.x, direction.y, steepness, wavelength).

The fragment shader shades with the normal, the view direction, and the foam mask. Section 9 walks through it.

On the JavaScript side there is very little to do: upload the grid once, then each frame set five uniforms and issue one drawElements. The grid has over 148,000 vertices, so the index buffer needs Uint32Array and gl.UNSIGNED_INT, which WebGL2 supports without an extension.

The WGSL is a line-by-line port. The one structural difference is that WGSL has no inout, so the tangent accumulators are passed as ptr<function, vec3f> and updated with *tangent += …. Uniforms live in a single struct bound at group 0.

WGSL's uniform layout rules are strict, so it's worth writing the byte offsets down once. Packing time into camPos.w keeps everything in 16-byte rows with no padding surprises:

The other trap is depth range. WebGL clip space runs z from −1 to 1; WebGPU runs it from 0 to 1. Reusing a GL projection matrix in WebGPU clips away the near half of your scene. The demo builds the matrix with a flag:

The renderer creates a pipeline with a single float32x2 vertex attribute, a depth buffer, and 4× MSAA, then records one render pass per frame.

The shading model is deliberately small, and every term earns its place:

• Fresnel. Water reflects about 2% of light head-on and nearly all of it at grazing angles. Schlick's approximation with F_0 = 0.02 blends between the water's body colour and the reflected sky. This one term does most of the work of making it read as water.
• Sky reflection. A two-colour gradient indexed by the reflected ray's height. The ray's y is mirrored with abs() so steep back-facing slopes don't reflect a nonexistent sky below the horizon.
• Body colour. Deep blue-green in troughs, brighter teal toward crests, plus a scatter term when the sun is behind a crest. It's a cheap stand-in for light passing through thin water.
• Sun glint. A very tight Blinn-Phong highlight. The exact normals from section 4 are what make it sparkle along crest edges instead of smearing.
• Foam and fog. The Jacobian mask from section 5, and squared-exponential fog toward the horizon colour. The page behind the canvas is a CSS gradient that ends in that same colour, which hides the edge of the finite grid.

Sooner or later you want a boat to float. On the CPU this is less obvious than it looks: because points move sideways, the vertex that ends up above world position \mathbf{q} did not start at \mathbf{q}. You need the rest position \mathbf{p} whose displaced position lands on \mathbf{q}.

Fixed-point iteration solves it: guess \mathbf{p} = \mathbf{q}, subtract the horizontal displacement at the guess, repeat. When the steepness budget holds, each step shrinks the error by a factor of at least \sum s_i, so a few iterations are plenty.

Because the inverse is iterative, check it the way you'd check any inverse: run the result forward again and measure how far it lands from where you asked. The readout below does that at a point 12 m in front of the origin, live, against whatever waves you've dialled in.

Surface height there:

…

. Round-trip miss after 4 iterations:

…

.

Push the total steepness above 1 and watch the miss grow: the iteration stops being a contraction exactly where the surface can fold.

Keep the CPU wave list identical to the GPU one, including any distance-based fading you add later, or your boat will float on a different ocean from the one on screen.

• Time precision. After a few hours, k * c * t in 32-bit floats loses enough precision that waves visibly stutter. Compute each wave's phase k c t \bmod 2\pi on the CPU in double precision and pass that instead of raw time.
• Zero-length directions. normalize(vec2(0)) is NaN, and one NaN wave blanks the whole mesh. Validate directions before upload.
• Distant aliasing. Short waves become shimmering noise far from the camera where the grid is coarse. Fade their steepness to zero with distance (in the vertex shader and in your CPU query).
• More waves. Four is a teaching number. Move them into a uniform block or a WebGPU storage buffer and sum 16 to 64 with directions and wavelengths drawn from a spectrum.
• Normal maps. Geometry carries the swell; small ripples belong in a scrolling normal map blended with the analytic normal.
• FFT oceans. Tessendorf's method uses this same horizontal displacement idea, but drives thousands of waves from a statistical spectrum with an inverse FFT. Everything in this tutorial (steepness, the Jacobian foam, the height query) carries over.

For further reading, Mark Finch's "Effective Water Simulation from Physical Models" (GPU Gems, 2004) is the classic real-time treatment, and Jerry Tessendorf's "Simulating Ocean Water" course notes are the next step up.