Generalizing Shallow Water Simulations with Dispersive Surface Waves

This paper introduces a novel method for simulating large bodies of water as a height field. At the start of each time step, we partition the waves into a bulk flow (which approximately satisfies the assumptions of the shallow water equations) and surface waves (which approximately satisfy the assumptions of Airy wave theory). We then solve the two wave regimes separately using appropriate state-of-the-art techniques, and re-combine the resulting wave velocities at the end of each step. This strategy leads to the first heightfield wave model capable of simulating complex interactions between both deep and shallow water effects, like the waves from a boat wake sloshing up onto a beach, or a dam break producing wave interference patterns and eddies. We also analyze the numerical dispersion created by our method and derive an exact correction factor for waves at a constant water depth, giving us a numerically perfect re-creation of theoretical water wave dispersion patterns.


INTRODUCTION
The motion of water is well-described by the incompressible Euler equations where u is the water velocity, u/ is the material derivative, is water density, is pressure, and is acceleration due to gravity. Numerically approximating this equation is prohibitively expensive for the animation of large bodies of water, so researchers reduce the complexity by assuming the water surface takes the form of a height field, where the water height and velocity are both just functions of 2D spatial coordinates, instead of 3D. The most common reductions of these equations used in computer graphics are Airy wave theory and shallow water approximations. Airy wave theory [Airy 1841] assumes that the motion is a potential flow (velocity is the gradient of some potential, u = ∇ ) and that the wave amplitude is small relative to the wavelength ( ≪ , or equivalently ≪ 1 for wavenumber = 2 / ). These assumptions produce linearized water wave equations which are used extensively throughout the computer graphics literature [Canabal et al. 2016;Tessendorf 2004a]. This linear wave theory produces ACM Trans. Graph., Vol. 1, No. 1, Article 1. Publication date: January 2023. waves with an angular frequency depending on wavenumber and water depth ℎ This particular dependence of on is called the dispersion relationship, and it tells us how the frequency varies with different wavelengths. It also effectively prescribes the surface wavespeed, which is equal to / . Accurately reproducing this dispersion relation is essential for many water-related phenomena like raindrop ripples and particular interference patterns in the wake behind boats, as illustrated in Figure 1. One drawback to this approach, however, is that its ≪ 1 assumption prevents the wave heights from having a major effect on the fluid domain boundaries. For example, it cannot model water waves sloshing up a sloped beach, or spilling over a dam and filling a basin, because these scenarios require the water domain to change shape depending on the motion of the surface waves.
Another popular technique for simplifying water dynamics is to assume that the depth of the water ℎ is much smaller than the wavelength (ℎ ≪ , or equivalently ℎ ≪ 1). This assumption gives rise to the shallow water equations (SWE), which do allow waves to slosh around and spill over terrain. However, this long wavelength / shallow depth assumption is only appropriate in exceptionally limited scenarios, because SWE produces a dispersion relationship which is drastically different from Equation 2 when ℎ . Notably, SWE is unable to produce the signature water wave interference patterns described above, and instead produces ripples similar to acoustic shock waves.
Although both of these methods for animating water are theoretically sound and offer robust and efficient implementations, they both have restrictive assumptions that lead to visually obvious breakdowns in common scenarios: Airy wave simulators have to avoid flooding scenarios and gently sloped solid obstacles; shallow water solvers either add procedural textures to approximate a more appropriate dispersion relation for deeper water [Chentanez and Müller 2010], or they fail to convincingly animate short wavelengths altogether.
Interestingly, both of these simulation techniques excel where the other fails. The purpose of this paper is to introduce a principled unification of the two fluid regimes together into a single model that lacks all of the failure modes described above. We do this by decomposing the water into two regimes of water motion: a bulk flow that is best described by the shallow water assumption ℎ ≪ , and the remaining surface waves that fail the shallow water test but still obey Airy wave theory. The bulk flow naturally includes most of the fluid's mass and momentum (it includes the biggest waves with the longest wavelengths), and so we simulate it using a shallow water solver capable of simulating flooding and convective eddies. Conversely, the surface waves contain all of the surface details and high frequencies that are essential for producing convincing water wave animations, so we simulate them with a highly detailed Airy wave solver. Re-computing the decomposition at each timestep allows water motion to continuously shift between shallow and deep regimes and furthermore it permits the fluid domain to significantly change over time.
In summary this paper offers the following contributions: • The first height field method capable of simulating both Airy wave interference patterns and large bulk motions like flooding and 2D convective eddies in the same simulation; • A practical algorithm for producing a spatially-varying decomposition of shallow and surface flows; • A derivation of surface wave numerical dispersion errors, and a novel scheme for canceling them exactly in a constant depth; • Exact volume conservation and real-time performance appropriate for games and other interactive scenarios.
The remainder of this paper is organized as follows: Section 2 discusses related work, and Section 3 provides additional background from fluid dynamics and introduces the theoretical framework behind our method. Section 4 describes our algorithm in detail, including the decomposition, the shallow water solver, the modified Airy solver, and implementation details. Section 5 evaluates our method and discusses results and future work, Section 6 concludes the paper.
2 RELATED WORK 2.1 Linearized water surface waves Airy [1841] noted that the surface of an inviscid, irrotational, and incompressible fluid under the influence of gravity will exhibit waves according to the very specific dispersion relation in Equation 2. Lord Kelvin noted that the interference between these waves traveling at different speeds completely explains the visually recognizable patterns we see in wakes behind boats [Thomson 1887].
Computer graphics researchers devised numerous ways to mimic these dispersive water waves. Peachey [1986] and other researchers [Fournier and Reeves 1986;Hinsinger et al. 2002;Mastin et al. 1987] animate waves by directly computing the inverse Fourier transform of Equation 2. Tessendorf [2004a] use convolutions of speciallydesigned kernel functions on a regular grid. Subsequent research simplified the kernel computation by modifying the dispersion relation [Loviscach 2002[Loviscach , 2003, increased performance and stability with an exponential integrator [Tessendorf 2014], and increased its accuracy in the presence of boundaries and obstacles [Canabal et al. 2016].
Airy wave theory exhibits a convenient mathematical structure that benefits from numerous techniques from the analysis of partial differential equations. The irrotational and incompressible assumptions of Airy wave theory imply that the velocity is the gradient of a harmonic potential function , further explained in Section 3. Researchers exploit this harmonicity using boundary integral [Keeler and Bridson 2014] and surface-only [Da et al. 2016] models for moving vertices on the liquid surface. Other researchers use fundamental solutions of this potential [Schreck et al. 2019] to analytically compute dispersive waves in open domains. Researchers have also exploited mathematical properties of the waves themselves by computing the evolution of wavefronts [Jeschke and Wojtan 2015], moving packets of wave energy [Jeschke and Wojtan 2017], and stationary surface wavelets [Jeschke et al. 2018] to animate water waves interactively.

Non-linear water surface waves
Airy theory is sufficient for animating visually striking interference patterns, and its linearity gives rise to several computationally convenient mathematical transformations. However, this simplicity also prevents it from being useful for animating non-linear effects, like surface eddies, the fluid domain depending on wave height (waves sloshing up onto a beach), and the wave speed depending on wave height.
Direct simulation of 3D liquids is one way to achieve these nonlinear effects, but they are prohibitively complex for simple scenarios where a heightfield simulation is sufficient. Thus we consider the rich literature of 3D fluid simulation in computer graphics out of scope for this discussion, and we will discuss only theoretical simplifications of the fluid equations which give rise to 2D equations. Boussinesq [1872] developed non-linear 2D equations with a polynomial expansion of the wave height, giving rise to a non-linear dependence between wave speed and wave height. While these equations are commonly used in coastal modeling, they are only applicable for very long wavelengths and are not commonly used in computer graphics. Saint-Venant [1871] developed a 1D version of what we now consider the shallow water equations (SWE). These equations are non-dispersive, but they allow turbulence in the form of swirling eddies, and their computation easily allows for a changing domain in flooding scenarios. Presumably for these reasons and for their ease of computation, SWE are popular in computer animation. Kass & Miller [1990] introduced a linearized SWE to graphics, and Layton & Van de Panne [2002] introduced a version of SWE with non-linear convective terms. Researchers subsequently refined SWE simulation with improved numerical integration and GPU implementation [Hagen et al. 2005], simulation on surface meshes [Wang et al. 2007], increased stability and non-reflecting boundary conditions [Chentanez and Müller 2010], and enhanced convection [Pan et al. 2012]. Azencot [2018] introduced a structure-preserving integrator for the Euler-Poincaré (EPDiff) equation, a particular class of SWE that does not exhibit the linear superposition of Airy wave theory.
All of these non-linear wave equations feature more interesting wave motions than simple linear Airy theory. However, they are either too computationally expensive to solve at large scales, or they are limited to a narrow range of phenomena (like only modeling shallow water).

Extended and combined methods
We are not the first to try to extend the capabilities of the Airy and SWE models beyond their original restrictive constraints. While the original Airy theory assumes that the steepness of the wave ℎ is small, Trochoidal waves [Gerstner 1809; Rankine 1863] allow for steeper waves by allowing points on the surface to follow a circular path, and Biesel waves [Biesel 1952] approximate breaking waves by expanding these circular motions into ellipses. Both of these models have found use in computer graphics [Fournier and Reeves 1986;Tessendorf 2004b], where applications require visually interesting choppy waves. Unfortunately, while these modifications increase the visual appeal of surface waves in open water, they exacerbate the problems that Airy waves have with boundaries. Techniques like Wave Cages [Jeschke et al. 2020] can effectively prevent these waves from penetrating walls and sea floors, but until now no technique allows Airy waves to change their simulation domain (e.g., by flooding into a new region).
The work of Yu et al. [2011] is noteworthy in that it advects strips of Airy interference patterns on top of a SWE simulation. The detailed surface waves utilize closed-form Airy wave interference patterns [Fedorov and Melville 1998], but the wave motion is governed by non-dispersive SWE and is limited to the special case of standing shock waves near an obstacle.
Several researchers aimed to compensate for the restrictions of 2D wave simulations by combining them with fully 3D liquid simulations. Thuerey et al. [2006] and Chentanez et al. [2015] coupled a liquid simulation with the shallow water equations, Huang et al. [2021] coupled 3D liquid to a quasi-2D surface-only water simulation [Da et al. 2016], and Schreck et al. [2022] coupled 3D liquid to an Airy wave solver based on fundamental solutions [Schreck et al. 2019]. These methods couple the 3D and 2D solvers using domain decomposition; some regions are simulated in 3D, while different regions (presumably further away or in areas requiring less detail) are limited to 2D simulation.
More recently, researchers used a different coupling strategy to add a 2D Airy solver directly on top of another fluid simulation solved with a different technique. Kim et al. [2013] added dispersive waves to a 3D fluid solver by one-way coupling a convolution-based solver with the closest-point method, and Skrivan et al. [2020] simulated the motion of Lagrangian wave packets [Jeschke and Wojtan 2017] on the 3D surface. These 2D/3D coupling techniques work exceptionally well when one can afford the computational demand of a fully 3D solver, but they are overkill for height field-based water dynamics. In addition, their one-way coupling strategy limits surface waves essentially to high-frequency textures, instead of a means to transport water into new domains. This paper addresses precisely these limitations.
Although methods like eWave [Tessendorf 2014] make it look easy, animating wake patterns is extremely difficult in general. SWE cannot create wake patterns, because all wavelengths travel at the same speed. 3D solvers are too dissipative and low-resolution for interference patterns to appear, so specialized 2D surface wave methods are the only way to create wake patterns in practice. Thus, our work is the only practical method for animating two-way coupling effects like those in Figure 3 (wakes interacting with vortices) and Figure 1 (wakes interacting with the shore), and it does so in real time.

THEORETICAL FRAMEWORK
Here we provide additional fluid mechanics background in Sections 3.1, 3.3, and 3.2, then introduce theoretical concepts behind our method in Sections 3.4 and 3.5.

Non-linear waves in shallow water
Linear waves at any deph

Linear waves in shallow water
Shallow Water Equations Airy Wave Theory Fig. 2. SWE vs. Airy wave theory. The shallow water equations (SWE) are capable of complex non-linear wave behaviors, but are limited to constant wave speeds in shallow depths. Airy wave theory models waves at all possible depths, but they are restricted to small-amplitude linear waves. Our method combines both, significantly expanding the range of behaviors possible within a single heightfield method.

Notation
For the remainder of this paper we assume the water is a height field Wave theory often decomposes the wave height function into an average water depthh and a perturbed height functionh, such that ℎ =h +h, where most of the interesting wave dynamics happens toh, whileh is relatively static. We will use similar notation to decompose other functions, with a − over averaged quantities and ∼ over perturbed quantities. When discussing waves, k is the wavevector with wavenumber = ∥k∥, wavelength = 2 / , angular frequency , and amplitude . Airy theory assumes a velocity potential such that u = ∇ .

Linear surface waves
According to Airy wave theory [Johnson 1997], the vertical displacementh from a reference surface with water depth ℎ is with defined in Equation 2. The velocity potential is defined as where is the depth coordinate. The horizontal velocity at any water depth is defined as Integrating this velocity along the water depth gives us the horizontal volumetric flow rate in linear wave terms: Integrating along depth to 0 instead ofh simplifies the terms and makesq independent of water depth. This approximation is consistent with the ≪ 1 assumption of Airy theory.

Shallow water equations
In the shallow water regime (ℎ ≪ ), the water surface is assumed to be so close to the sea floor that the velocity is roughly constant throughout an entire column of water. The shallow water equations in conservative form are: where q = ℎu is the horizontal flow rate, and ⊗ is the dyadic product. In our implementation, q is the amount of liquid flowing from one grid cell to another, through a rectangle of height ℎ and width Δ . These equations have transport terms which are non-linear in u, q, and ℎ. If we linearize these equations, we arrive at the popular "pipe" model [Kass and Miller 1990], which is a wave equation with dispersion relation = √︁ ℎ. As an interesting aside, these linear waves are actually a subset of Airy wave theory, which prescribes = √︁ tanh( ℎ), because tanh( ℎ) ≈ ℎ for small values of ℎ. SWE and Airy theory overlap exactly in the linearized shallow water regime, but SWE is better for non-linear behaviors, and Airy theory is better for deeper water. This range of behaviors is conceptually illustrated in Figure 2. Figure 3 shows how our method can simulate both non-linear and Airy effects together.

Dynamics of decomposed waves
We generalize these wave dynamics by incorporating two different types of waves together in the same mathematical model. We first decompose the water height ℎ =h +h as mentioned earlier. As before,h is a height field with relatively high frequencies compared toh: it models small ripples on top ofh, and it obeys Airy theory.
Then, instead of assuming thath is a static average water depth, we consider it a more general low-frequency water height function modeling the "bulk" of the fluid.h is allowed to move over time according to fluid dynamics. To tractably model these dynamics, we assume that the low-frequencyh function consists of wavelengths so long that ℎ ≪ and it obeys shallow water dynamics. We then seek to know how the dynamics ofh affecth and vice-versa.
We first note that linear wave dynamics obey the superposition principle: ifh andh are both Airy waves, then so is ℎ =h +h, and the two waves have no influence on each other. In other words, it is trivial to combine these two flows together whenh is calm and obeys linearized SWE.
Although we restricth to linear wave theory, we allowh to be nonlinear. We can estimate the effect of non-linear terms by plugging ℎ =h +h and u =ū +ũ into the continuity equation (Equation 8): non-linear (10) The first term represents advection through the bulk flow, which we model directly by solving the shallow water equations forh and u. The second term is the evolution of perturbed variables, which we model with Airy wave theory. The third term is the transport of high-frequency surface displacementsh through the smooth bulk velocity fieldū, which we explicitly model with an additional advection step. Finally, the fourth term is the influence of the surface wave velocityũ across the entire water depthh; this term is zero by our assumption thatũ obeys Airy wave theory and cannot have large enough amplitudes to affect the entire depthh. We therefore do not model it in this work. In the end, we are left with 3 terms that model bulk flow, surface waves, and transport.
Although this discussion focused on the evolution of ℎ, similar reasoning applies for the transport of flow rates q (left-hand side of Equation 9), giving us three analogous terms for the update of q. Henceforth, we use − to denote "bulk" terms, ∼ to denote "surface" terms, and ≃ to denote "transport" terms.

Decomposing waves
We want to decompose ℎ and q such that their low-frequency componentsh andq satisfy the shallow water assumption ℎ ≪ . There is no exact wavelength cut-off for this property to hold, but we know that it should scale monotonically with water depth ℎ. This will allow, for example, long waves to still obey pure Airy theory in deep water, and short waves to obey SWE if the water is shallow enough. We also wish to prevent extremely steep waves from appearing inh, as they violate the small amplitude ≪ 1 assumption of Airy wave theory, and they are better treated as shock waves in the non-linear SWE viah.
Rather than filter the raw water depth function ℎ, we consider the height field = ℎ + , which combines the depth and terrain height . We do this because a smooth will give us a physically reasonable water surface, regardless of what the terrain looks like below. We use a diffusion equation as a low-pass filter [Hamming 1998]: where is a spatially-varying diffusion coefficient, and is the fictitious time over which we integrate the diffusion. Fourier analysis of this equation reveals that it filters out high frequencies proportional to the diffusion coefficient , the wavenumber squared, and the integrated time :ˆ( whereˆ( , ) is the Fourier component of corresponding to wavenumber at fictitious time . Thus diffusion acts as a filter which removes the high-frequency components˜from to get a low-frequency¯. We then recomposeh =¯− andh = ℎ −h from the filtered height field. The filter strength depends on , which should increase monotonically with water depth and decay with wave steepness. We use the same technique to compute q =q +q. We choose to decompose the ℎ and q variables (and reconstruct them again) every simulation step, because the wave frequencies may dramatically change over time. By recomputing the wave decomposition, we allow waves to adapt their behaviors as they move around, change frequencies, and change water depth.

OUR GENERALIZED WAVE MODEL
Following Chentanez and Müller and others [Chentanez and Müller 2010;Stelling and Duinmeijer 2003], we use a staggered grid discretization of the water heightfield, with ℎ stored on grid cell centers, and components of velocity and flow rate stored on cell boundaries. The water heights are time-integrated with the explicit Verlet leapfrog scheme, so ℎ is evaluated one half time step later than Algorithm 1: Overview for one time step Input q and ℎ +Δ /2 q ,q ,h +Δ /2 ,h +Δ /2 Decompose q t and ℎ +Δ /2 q +Δ Simulate bulk fluid flow q +Δ Simulate Airy waves q t+∆t Transport surface flow rateq +Δ throughū ℎ +3Δ /2 Transport surface heighth +Δ /2 throughū q +Δ Compute flow rate from q +Δ ,q t+∆t ℎ +3Δ /2 Compute wave height from q +Δ ,h +3Δ /2 q. This staggered discretization represents a finite volume scheme which enhances simulation stability and guarantees volume conservation; each substep of the algorithm accumulates the flow rates q needed to update the wave heights ℎ =h+h and the fluid momentum encoded as q =q +q. At the end of each time step, we reconstruct the final ℎ from these integrated flow rates, guaranteeing volume preservation, and resulting in complex flows like Figure 4. This section is divided into sub-sections representing major substeps of our algorithm, which is also outlined in Algorithm 1. We first decompose ℎ and q into their bulk and surface components (Section 4.1), then simulate the bulk components via the shallow water equations (Section 4.2), and simulate the surface components with an Airy wave solver (Section 4.3). We transport the surface variables to account for the third term in Equation 10 (Section 4.4), and finally integrate the overall height ℎ through time at the end of the time step (Section 4.5).

Wave decomposition
To decompose ℎ and q into bulk and surface components, we follow Section 3.5 by implementing a forward-in-time centered-in-space explicit diffusion solver on the GPU. This operation converts ℎ and q intoh andq, which we use for the bulk flow. We useh = ℎ −h andq = q −q for computing surface waves. By increasing the integration time and diffusion strength , we remove more and more high frequencies fromh andq -no diffusion produces a pure SWE simulation with ℎ =h, and diffusing all the way to a steady state produces a pure Airy wave simulation with ℎ =h.
In theory, waves should always travel slower in shallow water, but blindly applying SWE in deep water makes them travel impossibly fast. We employ a frequency filter to minimize this nonphysical behavior: for a given water depth, there exists a wavelength where the shallow water wave speed exactly equals the deep water wave speed, cutoff = 2 ℎ. Wavelengths shorter than cutoff should not be treated with SWE. We employ diffusion as a classical Gaussian filter with "cutoff wavelength" (where power is reduced by 3db) equal to cutoff , giving us a diffusion coefficient = ℎ 2 64 , as derived in Appendix C. Setting to this function of depth turns the diffusion equation into a low pass filter, preventing short wavelengths in deep water from being modeled with SWE.
To discourage very large, steep waves from being mis-classified as high-frequency surface ripples, we add another penalty term with tuning parameter = 1 100 . This heuristic discourages steep gradients in surface waves and otherwise achieves monotonic depthdependent decomposition. We find that the gradient penalty − | ∇ℎ | 2 is only needed in extreme situations in practice, like dam breaks where the height field is a step function with arbitrarily high frequencies. The effect of penalizing steep waves in the decomposition is visualized in Figure 5.
To simulate this diffusion, we sample on a staggered grid between ℎ samples. This diffusion scheme is an effective filter for our purposes, but its implementation has practical limitations. Although our system runs in real-time, this decomposition is the most expensive part of our algorithm. Explicit integration requires small time steps (128 sub-steps during each time step of wave simulation) and we explicitly clamp Δ to its maximum stable value to avoid instability in arbitrarily deep water. A more efficient decomposition (better diffusion solver or different filtering operation) should accelerate computation further.

Simulating bulk fluid flow
Given our newly decomposedq along withh −Δ /2 from the previous time step, we computeū =q /h −Δ /2 ; thish is stored in cell centers whileq is stored on cell boundaries, so we use first-order up-winding to select the most appropriateh. We then numerically integrate the shallow water equations to updateū +Δ ; for this purpose, we adopt the conservative finite-volume scheme of Stelling and Duinmeijer [2003]. This method gives explicit time stepping rules for updatingh andū, which we reproduce in Appendix A. Please see Section 5 of [Stelling and Duinmeijer 2003] for full implementation details.
At the end of this operation, we have newū +Δ values evaluated at the boundaries of grid cells. We convert back toq +Δ = u +Δh +Δ /2 , with the most recenth given by the surface decomposition in the previous section, again using first order up-winding to evaluate it on the staggered grid. Our decision to useh andq from different steps is done such that it exactly reproduces the leapfrogged SWE simulation of [Stelling and Duinmeijer 2003] in the simplified case of pure shallow water. After this operation, we have computed the flow ratesq +Δ corresponding to the "bulk" term in Equation 10.

Simulating Airy waves
We use the exponential integrator of eWave [2014] to update the flow ratesq using the derivatives of surface wave heightsh. We present pseudocode for this step in Algorithm 2, where is the fast Fourier transform computed globally over the whole domain, with as its inverse; the hatˆnotation indicates that the function exists in Fourier space.
At the start of Algorithm 2, the exponential integrator expects bothh and˜to be sampled at the same time , but the SWE solver from the previous section uses Leapfrog integration, leavingh ands taggered in time. We re-sampleh at time via linear interpolation, then compute its Fourier transform to getĥ .
In the third step of Algorithm 2, is the dispersion relation from Equation 2. We use the decomposed smooth water depthh from Section 4.1 for this purpose, so that = √︃ tanh h . However, the spatially varying water depthh is not available in Fourier space, so we approximate the behavior of waves at different water depths by calculatingq +Δ for different water depths ( = 4 in our implementation withh ∈ {1 , 4 , 16 , 64 }, refer to Figure 8 for an error analysis), evaluatingq +Δ for each value ofh , and linearly interpolating theq +Δ corresponding to the two closest depths.
The calculation ofq +Δ needs the derivativeˆh x , which we compute directly in Fourier space [Johnson 2011]. Because it is sampled on a staggered grid compared toq, we off-setˆh x by half a grid cell in Fourier space via the Fourier shift theorem (multiplyinĝ ℎ by − Δ /2 ).

Eliminating numerical dispersion.
We mentioned at the start of this section that the final wave heights ℎ will be reconstructed in a volume-preserving manner from the flow rates q at the end of the time step. This finite-volume reconstruction is based on a numerical approximation of the divergence operator in Equation 8, and it causes numerical dispersion. Consequently, our surface waves do not actually obey the desired dispersion relation if we set = √︃ tanh h . In Appendix B we derive the numerical dispersion and a correction factor , so that we achieve perfectly accurate wave speeds when we set = 1 √︃ tanh h in Algorithm 2. Figure 6 illustrates the impact of this correction where we measure wavespeed with standing waves in a rectangular pool similar to [Schreck and Wojtan 2022] at a water depth of 4 m, a gridcell size of 1 m, and a timestep Δ = 1 60 . Relative wavespeed Wavelength (m) Fig. 6. Surface wave dispersion error. We plot the measured wavespeed relative to the exact value; a value less than 1 indicates artificially slow waves, while a value greater than 1 means the waves are too fast. The uncorrected simulation (orange) exhibits numerical dispersion, while our corrected method (blue) is practically perfect even down to the Nyquist limit of = 2Δ .
At the end of this step, we have updated flow ratesq +Δ corresponding to the "surface" term in Equation 10.

Transporting surface waves
Next we transport the surface flow ratesq and wave heightsh through the bulk velocityū: We choose to apply operator splitting to these equations and numerically integrate each term separately, and we evaluate the derivatives at the midpoint of the timestep, in the style of a "leapfrog" integrator. The ∇ ·ū term amplifies or damps waves based on the convergence or divergence of the underlying bulk flow. By holdingū constant over the time step, this becomes a linear ODE that we integrate in closed form using exponential integration. In nature, this term is responsible for wave shoaling -the sharp steepening, and eventual tumbling, of surface waves when they enter shallower water with slower flows. Our model cannot handle breaking effects, so we follow previous works (which artificially damp steep waves [Jeschke and Wojtan 2015] or omit the growth term altogether [Tessendorf 2017]) by adding an optional user-tunable factor that damps this term only when it amplifies the wave in a converging background flow. We compute theū · ∇ terms using Semi-Lagrangian advection with cubic spatial interpolation [Fedkiw et al. 2001] and avoid overshooting by clamping the interpolated value to lie within the range of the four neighboring values. For convenience, we transport the most up-to-date version of( after simulating Airy waves in Algorithm 2), so the resulting term q encodes both the "surface" and "transport" terms of the flow-rate version of Equation 10. In contrast,h is not directly updated in Algorithm 2, so it is not available to be combined with the transport term. Thus,h encodes only the "transport" term in Equation 10. The detailed steps for transporting surface flow ratesq and heightsh are written in Algorithms 3 and 4.

Merging bulk and surface flow components
At this point in the time step we have computed all of the bulk, surface, and transport terms needed to update ℎ and q according to Equation 10. The integrated bulk flow rate is stored inq +Δ , and the surface wave flow rates were first updated then transported and stored inq t+∆t . The final q value is simply the sum of these terms: Next we wish to compute the final height ℎ. To guarantee volume conservation during this process, we construct ℎ from the divergence of flow rates (Equation 8) rather than updating it directly. The flow rates corresponding to the bulk and surface terms in Equation 10 have just recently been collected and stored in q +Δ , so all that remains is to compute a new flow rateq corresponding to the "transport" term for updating ℎ. We do this by multiplying the transported height by the bulk velocity:q =hū. Recall that q and u are sampled at time + Δ and are located on a staggered grid, so we need to re-sampleh +3Δ /2 in space and time to make it coincide withū +Δ . To perform this re-sampling, we re-purpose our semi-Lagrangian advection algorithm to transporth +3Δ /2 backward in time by Δ /2, to get ah +Δ which is stored on the boundaries between grid cells. We then compute˘+ Δ =h +Δū +Δ .
We plug these flow rates into Equation 8 to produce the final wave heights: The time step is now finished, and these q and ℎ variables will now feed into the next time step, where they will get decomposed into bulk and surface components and updated all over again.

Additional implementation details
To interact with dry regions, our algorithm prescribes solid boundary conditions as in [Chentanez and Müller 2010]: terrain that is higher than the water level sets the flow rate q to zero (effectively creating reflecting boundary conditions); flow rates in dry regions can otherwise be overwritten by the shallow water solver, in order to flood into new terrain. We perform this update once per time step, so the maximum flood rate is limited by the CFL condition. Using these boundary conditions allow us to compute the FFT in Algorithm 2 over the entire square domain (including both wet and dry regions), without any further special treatment for irregular water domains.
To avoid instabilities in the presence of fast flows, we use the CFL condition to clamp velocities and flow rates to their maximum stable values |u| max = Δ /(4Δ ) and |q| max = ℎΔ /(4Δ ). We apply this limiter wherever any velocity and flow rate variables are updated: after SWE integration, after transport ofh andq, and after summing flow rates in Equation 16. This limiter guarantees SWE stability at the expense of incorrect flow speeds at large time steps. This is not a problem in our experience, though it can cause two small visual artifacts: slower flooding than normal, and an artificially large wave height in places where the water would spread out quickly if it was not throttled.
For extra visual effect in Figure 1, we add a small amount of normal-mapped "texture" to the wave surfaces by computing lowamplitude deep water waves via FFT [Mastin et al. 1987], and adding the resulting surface normal to our water surface. We also use the technique of Tessendorf [2017] to add horizontal displacements to our waves in this example, in addition to the vertical wave displacement ℎ. We render the waves with a transparency that scales with depth, so the water becomes more transparent in shallow depths.
To produce boat wakes in our examples, we add point disturbances in q at every time step beneath the path of the boat. We also move the boat's vertical position to lie on the water level ℎ. We leave more thorough solid-fluid coupling to future work.

RESULTS AND DISCUSSION
We implemented our wave simulation algorithm in CUDA on an NVIDIA RTX2080 Max-Q Laptop. Our examples run in real-time (comfortably over 40fps), including the cost of rendering. Without rendering, our simulation runs around 100fps. The wave decomposition algorithm (explicit integration of a diffusion equation) takes 87% of the simulation time and is the clear bottleneck of our solver. For our results we set Δ to 1m and the grid size to 512 × 512; we set Δ to 1/60s for Figures 1 and 3 and to 1/30s for our other examples. Figure 1 shows how our algorithm computes both deep water and shallow water effects in the same scene. Two boats travel at different speeds through deeper water, producing significantly different wake patterns. An incoming wave floods onto the shore and creates a new pool of water. Figure 4 features a dam break scenario, where a wall with three open slots separates a wet domain from a dry one. The water floods through the narrow channels, creating standing wake patterns, vortices, and numerous reflected waves on the water surface. Our method uses the bulk flow to model the flooding, and the Airy waves simulate the surface details. Figure 3 models a cylindrical pillar in a steady background flow of moving at 0.45 pillar diameters per second. The pillar sheds vortices/whirlpools and forms a von Kármán vortex street, as well as a Kelvin wake pattern. The whirlpools and surface waves interact with each other. In our model, the bulk flow simulates (shallow water) vortices, and the Airy wave model creates the wake interference patterns. Finally, Figure 10 features waves spilling in from the boundaries of a square pool and interacting with complex obstacles that dynamically grow and change over time.

Stability and validation
To make our method robust enough for real-time simulation, we take care to keep each of the substeps of our algorithm as numerically stable as possible. The SWE solver enforces stability by clamping flow rates to their maximum stable value, which can be inaccurate for large time steps. The Airy solver and wave transport algorithms (Semi-Lagrangian advection and closed-form wave amplification) are each unconditionally stable when taken in isolation. The final wave height integration guarantees volume conservation exactly by construction. We merge these various methods together using operator splitting, which is generally first order accurate and not guaranteed to be stable.
We analyze the effect of varying the time step size in Figure 7, using Δ = 1/120s and 1/15s. Simulations with smaller timesteps produce smaller wavelength waves which endure longer, implying that our method exhibits some numerical damping for large Δ . Conversely, the large wavelengths can amplify in the simulations with larger time steps, because artificially clamping¯for stability can cause velocity to slow down, potentially reducing ∇ ·¯and increasingh/ . Section 4.3 uses a piecewise linear function to sample surface wave behaviors at different water depths. Figure 8 shows that the accuracy of this approach depends on how many samples we use: using more samples makes for more accurate wave speeds, at the Relative wavespeed Water depth (m) Fig. 8. Error from discretizing water depth. We plot the measured wavespeed of a wave with wavelength = 4 relative to the exact value for different water depths. The orange curve shows the error from a simulation which samples at only two depths (1m, 16m), and the blue curve at three depths, respectively (1m, 4m, 16m). expense of additional surface wave simulations. Our simulations use four depth samples, giving a maximum of 8% error (less in practice, because most surface waves do not depend much on water depth). One could further optimize these discrete values to minimize the error on all relevant water depth and wavelength combinations, but we did not deem this necessary for satisfactory visual results. Figure 9 shows the relationship between the accuracy of the surface wave speed and the wave decomposition (Section 4.1). The different curves represent different numbers of iterations in the diffusion solver each time step. Diffusing by a larger amount makes a more effective frequency filter, so we see a relationship between the number of explicit diffusion steps and the accuracy of the wave speeds. Less diffusion causes a larger proportion of waves to be treated as shallow water with a faster wave speed; a larger diffusion lowers this error, turning most of the simulation into Airy waves.
To measure the accuracy of our method's effective dispersion relation in Figures 6, 8, and 9, we measured by counting the period of standing waves at different wavelengths, just like Schreck & Wojtan [2022].
Lastly, we note that our framework generalizes both SWE and Airy wave solvers; if we force our wave decomposition to seth = ℎ andq = q, then we get a SWE simulator identical to that of [Stelling and Duinmeijer 2003]. Conversely, if we seth constant andq = 0, we get an Airy wave solver indistinguishable from [Tessendorf 2014].

Limitations and future work
To create realistic wake patterns, we take care that our method exhibits the correct wave dispersion. Our surface waves obey the exact dispersion relation of gravity waves, as discussed in Section 4.3.1 and illustrated in Figure 6. However, because our bulk flow simulator [Stelling and Duinmeijer 2003] uses a numerical discretization of spatial derivatives, it does show some unrealistic numerical dispersion in the presence of very steep waves (like in dam break scenarios).
We note that our model -and most other practical methods used in computer graphics -is based on theories that assume small wave amplitudes, so we cannot expect it to behave accurately in the presence of steep waves. Also, the choice of whether a given Relative wavespeed Wavelength (m) Fig. 9. Effect of wave decomposition on wave speed. We plot the measured wavespeed relative to the theoretical surface wave speed for simulations with varying diffusion times in our wave decomposition at a water depth of 4m. The yellow, gray, orange, and blue curves correspond to 8, 16, 32, 128 diffusion iterations, respectively.
wave is considered "bulk flow" or "surface wave" depends on our decomposition algorithm. Our model derives its efficiency by discarding the potentially complicated wave behaviors associated with the non-linear ∇ · hũ term in Equation 10. Neglecting this term implies that velocity perturbations caused by surface wavesũ do not reach all the way down to the sea floor -coupled interactions between the water surface and sea floor are associated with vortex shedding and wave breaking in nature. Linear theory suffices for interference patterns (Kelvin derived the wake pattern using linear theory), but is insufficient for these rotational effects.
Our surface wave solver operates in frequency space with an exponential integrator that assumes waves do not intersect boundaries over the course of a single time step. Consequently, depending on the time step and wave speed, waves may travel through thin boundaries. We could use the shadow kernels of Canabal et al. [2016] to prevent waves passing through boundaries in the future.
The decomposition in Section 4.1 uses an easy-to-implement explicit integration of a spatially inhomogeneous diffusion equation as a proof-of-concept. While it is already effective, this algorithm could be made more precise and efficient in the future, perhaps by using implicit integration, a spectral method, or any number of techniques from the digital signal processing literature. We look forward to refining this method in the future.
Our method assumes that the bulk and surface waves consist of low-and high-frequency height fields, respectively. This implies, because the bulk flowh interacts with the ocean floor (and because the surface wavesh do not), that the ocean floor is also a lowfrequency height field. To make our method work well with highfrequency ocean floors, we would have to take the ocean floor into account in the decomposition. Our current implementation does not do this, so extra work is needed to exhibit correct behavior with high-frequency terrain.
Equation 9 adds gravity as an external force. Although real water also experiences surface tension and viscosity forces, we did not include them in our implementation. Surface tension waves should be straightforward to add if we modify the dispersion relationship, though it might add an additional time step restriction on the bulk flow solver due to the high speed of capillary waves. Many explicit solvers add artificial damping to ensure stability, so we intentionally left it out to demonstrate robustness. In reality, viscosity damps shorter wavelengths faster than longer ones -a property that is currently missing from our simulator apart from numerical viscosity caused by Semi-Lagrangian transport of surface waves. In the future, we could add damping by approximating viscosity.
Our method also currently assumes that gravity is aligned with our height field direction. More work would be needed to simulate strong tangential forces or water droplets on vertical surfaces, perhaps by building upon the work of Wang et al. [2007].
Our shallow water solver can produce flows with significant vorticity, but our surface wave solver obeys linear wave theory, which assumes thatq is curl-free. Our decomposition in Section 4.1 does not make this distinction, and transfers high-frequency vorticity to the surface wave solver, where it is ignored. If desired, one could try to prevent the transfer of this high-frequency vorticity or add a vorticity solver to the surface wave simulation.
Finally, our algorithm for reconstructing the water surface should be improved. Our bi-cubic surface reconstruction requires a twocell-thick extrapolation into the nearby solid boundaries, and we apply a threshold on the water height to classify whether a surface is wet or dry. This combination of extrapolation and thresholding can cause the water surface to unrealistically flicker near steep walls that are close to the wet/dry threshold.

CONCLUSION
We have presented a height field water simulation algorithm that combines the benefits of shallow water simulation with those of surface wave simulation. Our results produce useful effects like flooding on variable terrain heights, wake patterns and ripples, and combinations of both surface waves and vortices in the same animation. We introduced a novel wave decomposition heuristic and a principled approach to simulating each of the terms in the decomposed fluid equations. Our method preserves surface wave speed and volume exactly, and it is well suited for parallel GPU implementation.

B CORRECTING NUMERICAL DISPERSION
Following [Marshall et al. 2004], we analyze the effect of our numerical spatial derivative operator on the shallow water wave equation: where is gravity, is the nominal depth of the water, is the flow and ℎ the free-surface displacement. and are partial derivatives with respect to time and space, respectively. The theoretical dispersion relationship of this system is = √︁ . Our implementation uses a spectral derivative (in Fourier space) to compute the first derivative, and a finite volume approximation on a staggered grid ≈ 1 Δ when reconstructing the wave heights from the divergence of flow rates in the second: The operator matrix is To discover the effective dispersion of this operation, we plug in the wave solution ( − ) and compute the derivatives: The determinant of this matrix gives us the numerical dispersion relation

C DERIVATION OF DIFFUSION COEFFICIENT
The diffusion equation (Equation 11) acts as a Gaussian filter with a cutoff wavelength cutoff = 2 ℎ. Rounding the filter's half-width at half-maximum to 1 for simplicity, this gives us a heat kernel with standard deviation equal to = cutoff 2 = 2 ℎ 2 = ℎ and thus = 2 2 = ℎ 2 2 . Fixing the diffusion timestep to the maximum stable size Δ = 0.25 and the number of diffusion iterations to a tolerable 128 iterations gives us a diffusion coefficient equal to ℎ 2 2 Δ = ℎ 2 64 .