## Abstract

A new semi-Lagrangian advection scheme called multistep ray advection is proposed for solving the spectral energy balance equation of ocean surface gravity waves. Existing so-called piecewise ray methods advect wave energy over a single time step using “pieces” of ray trajectories, after which the spectrum is updated with source terms representing various physical processes. The generalized scheme presented here allows for an arbitrary number *N* of advection time steps along the same rays, thus reducing numerical diffusion, and still including source-term variations every time step. Tests are performed for alongshore uniform bottom topography, and the effects of two types of discretizations of the wave spectrum are investigated, a finite-bandwidth representation and a single frequency and direction per spectral band. In the limit of large *N*, both the accuracy and computation cost of the method increase, approaching a nondiffusive fully Lagrangian scheme. Even for *N* = 1 the semi-Lagrangian scheme test results show less numerical diffusion than predictions of the commonly used first-order upwind finite-difference scheme. Application to the refraction and shoaling of narrow swell spectra across a continental shelf illustrates the importance of controlling numerical diffusion. Numerical errors in a single-step (Δ*t* = 600 s) scheme implemented on the North Carolina continental shelf (typical swell propagation time across the shelf is about 3 h) are shown to be comparable to the angular diffusion predicted by the wave–bottom Bragg scattering theory, in particular for narrow directional spectra, suggesting that the true directional spread of swell may not always be resolved in existing wave prediction models, because of excessive numerical diffusion. This diffusion is effectively suppressed in cases presented here with a four-step semi-Lagrangian scheme, using the same value of Δ*t*.

## 1. Introduction

Phase-averaged models that aim to predict the evolution of surface gravity waves over distances much larger than the wavelength are usually based on a spectral energy balance (Gelci et al. 1957). The wave field can be represented by the spectral energy density *F* (**x**, **k**, *t*), in wavenumber vector space (**k**), as a function of geographical space (**x**) and time (*t*). Neglecting currents, the energy balance equation is given by (e.g., Whitham 1974; Willebrand 1975)

where ∇** _{x}** and ∇

**are divergence operators in geographical and wavenumber space, respectively, and**

_{k}**c**

*(the group speed) and*

_{g}**c**

*are the corresponding energy transport velocities. The source term*

_{k}*S*(

**k**,

**x**,

*t*) is the net rate of energy transfer to component

**k**resulting from interactions with the bottom and the atmosphere, and nonlinear interactions with other components of the spectrum. This transport equation, similar to Boltzmann’s equation for gas kinetics, is written in terms of local values of the wave energy, as could be measured by a fixed instrument. For each spectral component (1) can be discretized at fixed locations in physical and wavenumber vector space and integrated in time on an Eulerian grid. This type of numerical model is widely used in deep water applications where large spatial and temporal scales of wave evolution allow for relatively coarse grids (e.g., WAMDI Group 1988). Finite-difference approximations of the gradients in (1) introduce numerical diffusion in time, and in both physical and wavenumber spaces. These effects can be mitigated by using higher-order finite-difference schemes, but this enhances the “garden sprinkler effect”; that is, energy tends to flow with discrete group speeds along the discretized propagation directions, causing over large distances an artificial accumulation of energy at discrete locations rather than the smooth dispersion of a continuous spectrum. This garden sprinkler effect is usually reduced by introducing additional numerical diffusion, thereby losing part of the accuracy of the higher-order advection scheme. Tolman (2002) gave a review of this approach with new improvements.

Although source terms play a dominant role in the generation process, the local changes in wave properties are often dominated by propagation effects that are well described by linear wave theory (see, e.g., Bretherton and Garrett 1969; Mei 1989), in particular for swell-dominated sea states and in shelf regions with complex and shallow bottom topography (Lavrenov 2003). Ray theory for waves propagating through a slowly varying medium can be used to simplify (1) by writing the energy balance following a wave packet along its ray trajectory:

The conservative form of (2), that is, with *S* = 0, states that the energy density in wavenumber vector space does not change along a ray (although the wavenumber vector itself does change) and was established by Longuet-Higgins (1957). It can also be derived from more fundamental conservation laws (Lavrenov 2003, and references therein).

Numerical evaluation of (2) only describes the changes in the energy of the moving wave trains; it thus also requires the computation of all ray trajectories followed by the spectral components, keeping track of changes in wavenumber vectors along the rays. All this is obtained from linear refraction theory using well-known ray-tracing algorithms (see Ardhuin et al. 2001). For spectral applications it was clearly demonstrated that backward ray tracing, from a given arrival point, is preferable to forward ray tracing, from parallel directions in deep water toward the area of interest (see, e.g., O’Reilly and Guza 1993; Bouws and Battjes 1982).

Such a “Lagrangian” method (in the sense that we follow wave groups but not water particles) is economical when rays are fixed in time, that is, in the absence of significant variable currents or temporal changes in sea levels, because rays can be precomputed once and for all. Furthermore, ray computations can be performed on a bathymetry grid different from the wave model grid, taking advantage of the full resolution of the bathymetric data available. The Lagrangian approach was initially used with *S* = 0 (e.g., Dobson 1967), providing an accurate method for the transformation of swell over complex bathymetry (e.g., O’Reilly and Guza 1991). Cavaleri and Malanotte-Rizzoli (1981) introduced a nonzero source term *S*, assuming that for a given wave component **k**, *S* (**k**) is only a function of the energy of that component, which is appropriate for commonly used parameterizations of some physical processes, such as wave generation by the wind [see also the review by Lavrenov (2003)]. Recently Ardhuin et al. (2001) presented a generalized model [Coupled Rays with Eulerian Source Terms (CREST)] that can handle arbitrary source terms, including scattering processes that couple different spectral components, such as wave scattering by the bottom topography (Ardhuin and Herbers 2002) or “quadruplet” wave–wave interactions (Ardhuin et al. 2004). The extension to arbitrary source terms is achieved by the interpolation of source terms from a fixed Eulerian grid onto the rays. The model essentially combines the well-established physical description of wave spectra evolution in the form of source terms with the accurate nondiffusive advection of finite spectral bands along rays. It also allows a stable integration of (2) with arbitrary time steps. Furthermore, the advection computations can be carried out over finite bandwidths using precomputed rays, thus eliminating the garden sprinkler effect of finite-difference schemes that use discrete spectral components (e.g., Tolman 2002).

The high computational cost of the full-Lagrangian advection scheme in CREST, associated with the variability of ray trajectories over large distances, is still prohibitive for large regional (e.g., 500 km of coastline) or global implementations on a workstation. Some alternative schemes, such as the “piecewise ray method” (Sobey 1986; Young 1988; Benoit et al. 1996) achieve numerical efficiency by propagating wave energy over only one time step along the ray trajectories, from locations (“roots”) determined by the ray calculations to the grid points (the “plant”). However, the interpolation necessary at these roots to determine the local energy from the neighboring grid points introduces numerical diffusion at each time step, similar to low-order finite-difference schemes.

Although it is well recognized that some numerical schemes commonly used in spectral wave prediction models are too diffusive to accurately propagate a localized disturbance over large distances (e.g., WAMDI Group 1988; Booij et al. 1999), the errors introduced by numerical diffusion in naturally broad wind-generated wave fields are not well known. In deep water, frequency and directional dispersion of swell causes a natural diffusion of energy in physical (**x**) space but also narrows the wave field in wavenumber (**k**) space. In shallow water, refraction can cause large spatial gradients in wave energy that enhance numerical diffusion. Wave–wave and wave–bottom scattering processes, on the other hand, provide natural diffusion mechanisms that tend to broaden the wave spectrum (Hasselmann 1966). Hindcasts of swell transformation across a wide, irregular shelf, using the wave model CREST, showed that diffusion by wave–bottom interactions (class I Bragg scattering) approximately doubled the directional spread near the shore of the North Carolina Outer Banks (Ardhuin and Herbers 2002).

The goal of the present paper is therefore to evaluate whether this type of scattering process can be resolved by numerical models that, for this purpose, may be limited by numerical diffusion. To investigate this question we propose a generalization of the Lagrangian CREST scheme, a multistep ray method, with an arbitrary number of time steps over which the wave energy is advected along the same ray, offering a clear and flexible trade-off between computational cost and numerical diffusion (section 2). The technique has similarities with INTERPOL, the semi-Lagrangian scheme described by Lavrenov and Onvlee (1995) and Lavrenov (2003, chapter 3). The present approach is illustrated in section 3 with predictions of the evolution of the directional spread of swell across the broad and shallow North Carolina–Virginia continental shelf during the Shoaling Waves Experiment (SHOWEX), comparing the effects of numerical and physical (Bragg scattering) diffusion in model results. Conclusions and perspectives for this generalized multistep propagation scheme are given in section 4.

## 2. Generalization of the CREST numerical scheme

CREST is a phase-averaged spectral wave prediction model based on Eq. (2) described in detail by Ardhuin et al. (2001). First, in the precomputation phase, ray trajectories for the waves are traced backward from fixed Eulerian grid points with positions **x*** _{i}* to the model boundary. Source terms are evaluated at the grid points and interpolated onto the rays. The spectrum

*F*is discretized in fixed frequency–direction bins, taking advantage of the conservation of wave frequencies along the rays. However, the spectral density

*F*is still expressed as density in wavenumber vector space because only this density is conserved along the rays in the absence of source terms. Source terms are interpolated linearly in space and direction to match the local ray position and direction.

Along each ray, arriving at **x*** _{i}* we define a Lagrangian energy density

*F*(

^{L}*t*,

*τ*) as the energy density “upstream” of

**x**

*at time*

_{i}*t*, where

*τ*is the energy advection time from the local ray position to the grid point

**x**

*. At*

_{i}**x**

*, the Lagrangian densities*

_{i}*F*(

^{L}*t*, 0) are averaged over all the rays that arrive at

**x**

*within a frequency–direction bandwidth Δ*

_{i}*f*Δ

*θ*, yielding the band-averaged discrete Eulerian spectrum

*F*(

^{E}*f*,

_{j}*θ*,

_{l}*t*) at

**x**

*. Here Δ*

_{i}*θ*is taken to be fixed, while Δ

*f*generally increases with frequency. At each grid point a source term

*S*(

*f*,

_{j}*θ*,

_{l}*t*) is determined from the full Eulerian spectrum

*F*and other local parameters such as wind stress and bottom roughness. Then

^{E}*S*is interpolated at the local ray positions to yield the source term

*S̃*(

*t*,

*τ*) along the rays, which in turn modifies

*F*(

^{L}*t*,

*τ*) (see Ardhuin et al. 2001, their Fig. 1).

Although the propagation over short distances for a finite bandwidth may be well represented by the position and direction along a single ray, the rays are scattered by complex topography. Thus several rays are needed to account for the variety of trajectories that reach a grid point with frequencies and directions in the same discrete band. As a result of this “scintillation” of rays, the cost of the CREST scheme increases rapidly with increasing model domain size. That is, the further the rays are integrated, not only must more rays be computed to describe widely scattered ray bundles, but these ray bundles also cover a larger region, thus increasing the effort to interpolate the source terms. The computational effort can be reduced by dividing the entire model domain into subdomains, at the boundaries of which the rays are terminated and the energy is interpolated from the neighboring boundary points (Ardhuin et al. 2001).

A more general and flexible approach is used here. Rays are terminated after a given number of time steps *N*. At the termination point, the wave energy density *F ^{L}* is interpolated from three neighboring grid points. This interpolation is similar to the interpolation of the source term, with the difference that for a given ray it is performed at a fixed location and time

*N*Δ

*t*, whereas the interpolated source term

*S̃*is distributed in space and time over a ray segment that corresponds to the finite advection time step Δ

*t*(Fig. 1). If no source terms are used, increasing

*N*is equivalent to increasing Δ

*t*as the energy is simply advected over larger distances along the rays. When source terms are included, the duration of the effective propagation time step

*N*Δ

*t*is essentially decoupled from the source-term integration time step Δ

*t*, and

*N*may be varied across the model domain and spectrum, for example, to better describe the propagation of swell in regions with complex bathymetry. If a ray crosses a boundary within a propagation time less than

*N*Δ

*t*, the energy is interpolated from the two adjacent boundary grid points (Ardhuin et al. 2001). Although limiting the number of advection steps

*N*is somewhat equivalent to dividing the model domain into subdomains of sizes close to

*Nc*/Δ

_{g}*t*, where

*c*is the group speed, the presented scheme is more flexible and does not introduce a bias in the propagation time due to the round-off to a lower integer number of time steps.

_{g}For small values of *N* (say, 1 or 2) the rays vary only slightly over a spectral band if the bathymetry and currents vary slowly over the distance *Nc _{g}*/Δ

*t*. Thus for low

*N*only one or a few rays are needed to represent accurately a finite spectral band (e.g., Benoit et al. 1996). In these numerically efficient low-

*N*schemes, rays may be recomputed at regular time intervals to account for variable currents and water levels that modify the ray trajectories.

The accuracy of this semi-Lagrangian advection scheme for a small number of time steps *N*, and using many or only one ray per frequency–direction band, is investigated below through comparisons with predictions of the full-Lagrangian scheme (*N* = ∞). In that scheme, many rays are used for each spectral component, and we adjust their number so that two rays with neighboring arrival directions have close directions and positions at their upwave ends. An upper bound on this number was set to 900 rays per finite frequency–direction band, which is rarely reached, and a minimum of 30 rays. Calculations described here use in general about 50 rays per component [see Ardhuin et al. (2001) for details on the procedure for adjusting the number of rays], unless the use of a single ray is specified.

## 3. Numerical diffusion in the directional spectrum

### a. Schematic shelf

A simplified transect version of CREST for an alongshore uniform shelf was created. The interpolation grid from which rays are computed is reduced to a line of points along a cross-shelf transect. A linear cross-shore interpolation of energy spectra and source terms is performed based on the two neighboring grid points, instead of the three points of the local triangle in the full model. The time-integration scheme is described by Ardhuin et al. 2001.

The first test, presented in Fig. 2, uses a simple bottom profile consisting of a plane, gently sloping shelf (1:3000) with a relatively steep inner shelf and shelf break (Fig. 3a) that represents the general geometry of the approximately 100-km-wide U.S. east coast continental shelf off North Carolina.

Along the cross-shelf transect, grid points were regularly spaced at intervals Δ*x*, from 8 to 200 m water depth. The incident wave energy was taken to be constant in time and distributed uniformly in a single frequency band (0.07–0.074 Hz), with a directional distribution proportional to cos* ^{p}* (

*θ*−

*θ*

_{0}), including only onshore propagating components. The exponent

*p*= 12 was chosen to obtain a narrow offshore directional spread

*σ*of 15°, close to the observed values in the period 15–21 September 1999 during SHOWEX (discussed below). Here

_{θ}*σ*is based on the first-order moments and equals the standard deviation of the directional distribution in the limit of a narrow spectrum (Kuik et al. 1988). The mean incidence wave direction was set to 60° (Fig. 2b). The directional distribution was discretized with a 5° resolution in the model. The model was integrated in time, starting from energy densities set to zero across the transect and along the rays until a steady state was reached.

_{θ}Model results for different advection schemes are compared to exact analytical solutions based on the invariance of the wavenumber vector spectral densities along ray trajectories given by Snel’s law (Longuet-Higgins 1957), keeping track of wavenumber vector changes. These analytical solutions were computed with 1° resolution for a stepwise constant offshore spectrum (with a step width of 5°) that is exactly identical to the one used in the numerical model integrations. Differences between model results and Snel’s law are thus entirely the result of numerical errors in the energy and direction of waves introduced by the interpolation at the end of the rays. This error generally has a small effect on the mean wave direction (less than 0.1° and not shown) because the interpolation yields both positive and negative direction errors that tend to cancel. However, the interpolation tends to diffuse energy to neighboring rays, thus artificially broadening the wave spectra.

In the multistep ray advection scheme presented here this diffusion can be reduced by increasing the number *N* of back-refraction steps. The time step Δ*t* was kept constant at 600 s, a typical value used in realistic applications with time-varying source terms and boundary conditions. Changing the value of *N* is equivalent here to changing the value of Δ*t*, keeping the product *N*Δ*t* constant. Model runs with *N* = 20 and a finite-bandwidth representation of each spectral element (using many rays for each frequency–direction band) correspond approximately to the original CREST scheme, since the back-refracted rays generally reach the offshore boundary in less than 20 steps (with the exception of very large oblique angles). Good agreement with the analytical Snel’s law result (*σ _{θ}* errors less than 0.4°) confirms that the full-Lagrangian CREST scheme effectively eliminates numerical diffusion (Fig. 2b). As

*N*is reduced, repeated interpolations in direction and space enhance numerical diffusion. On the inner shelf the relative increase of

*σ*due to numerical diffusion is approximately 60% for

_{θ}*N*= 1, 30% for

*N*= 2, and 15% for

*N*= 4. To examine how much of this diffusion is due to spatial interpolations, results for Δ

*x*= 2 km and 0.5 km are compared in Fig. 2b (squares and triangles). The nearly identical results indicate that the spatial variations of the energy densities are well resolved by the linear interpolation from the grid points to the end of the rays, even for a relatively coarse 2-km grid, and the numerical diffusion in this test case is essentially caused by the directional interpolation. For small values of

*N*, using a finite-bandwidth representation for the direction bands yields about twice as much numerical diffusion compared with results for a single ray per band (shown in Fig. 2b for

*N*= 1). Indeed, the interpolated energy at the end of a ray is an average of wave energies, previously advected by rays spanning a Δ

*θ*interval, and thus this average is clearly a source of numerical diffusion. However, if a single ray is used for each band, this interpolation at the end of the ray resolves linear variations in energy between adjacent bands and is therefore less diffusive.

The test results for full- and semi-Lagrangian schemes are compared in Fig. 2b with predictions of an explicit first-order upwind scheme that uses finite differences to represent the gradients in the Eulerian spectral energy balance [Eq. (1)]. The finite-difference scheme was implemented on a frequency ( *f* )–direction (*θ*) grid, taking advantage of the fact that *f* is conserved. On an alongshore uniform shelf, Eq. (1) reduces to an advection equation in three dimensions: *x*, *θ*, and *t*. To keep the Courant numbers, in both *x* and *θ* space, below 0.6, Δ*x* and Δ*t* were reduced to 0.25 km and 15 s, respectively. Results indicate that the semi-Lagrangian schemes, even the most diffusive *N* = 1 scheme, give less numerical diffusion than the first-order finite-difference scheme. In contrast to the semi-Lagrangian schemes, the finite-difference scheme does not account for wave refraction on subgrid scales (in space and directions) and thus introduces numerical errors in regions with significant depth changes over a grid cell. The associated strong diffusion is evident in the steeper shelfbreak and inner-shelf regions in Fig. 2b. The numerical results of the finite-difference scheme converge to the analytical result only if both Δ*x* and Δ*θ* are reduced, and the required resolution, of the order of 100 m and 0.1°, makes applications to realistic two-dimensional cases prohibitively expensive. The directional resolution Δ*θ* = 5° used here is already 3 times smaller than what is generally used in operational models, while a spatial resolution Δ*x* of the order of 1 km is a typical value for shallow-water applications.

The same test case was repeated using a more realistic cross-shelf bottom profile, taken from surveys of the actual North Carolina shelf, just south of the entrance to the Chesapeake Bay (Fig. 3a). This more irregular topography causes stronger refraction that further enhances the artificial directional smoothing of the numerical schemes. The upwind finite-difference scheme yields the largest errors, overpredicting *σ _{θ}* by a factor of 2–3 across most of the shelf (Fig. 3b). The

*N*= 1 semi-Lagrangian schemes typically overpredict

*σ*by 40%–60%, and these errors are reduced to less than 20% for the multistep

_{θ}*N*= 4 scheme.

### b. Physical and numerical diffusion in a SHOWEX hindcast

Model tests for an alongshore uniform shelf show that numerical diffusion caused by the directional discretization of the wave spectrum can broaden artificially the directional spectra predicted by wave models. For realistic bidimensional bathymetry, physical processes, in particular refraction and scattering of waves by the bottom topography, may also broaden directional wave spectra across the shelf. The importance of controlling numerical diffusion in a realistic setting with natural physical diffusion is investigated here by comparing different numerical schemes in hindcasts of observed swell evolution across the North Carolina continental shelf. During SHOWEX, a cross-shelf transect of six surface-following Datawell Directional Waverider buoys was deployed from September to December 1999. Data from three buoys, X1, X4, and X6, in 20-, 33-, and 193-m depth, respectively, are used here (Fig. 4). A description of the experiment, data processing, and hindcast procedure is given by Ardhuin (2001) and Ardhuin et al. (2003b). Wave measurements were also available from permanent instruments, including pitch-and-roll buoy 44014, maintained by the National Data Buoy Center (NDBC) in 49-m depth, and an array of pressure gauges in 8-m depth (8M in Fig. 4) at the U.S. Army Corps of Engineers Field Research Facility (FRF) in Duck, North Carolina.

The model CREST was implemented covering the entire shelf with 570 grid points extending from 34°30′N, south of Cape Hatteras, to 38°N, at the Virginia–Maryland border on Assateague Island. The model uses 29 frequency bands from 0.05 to 0.15 Hz, and 72 direction bands regularly spaced at 5° intervals. For each of these finite bands, bundles of rays were traced from all grid points. The grid was subdivided into nine subdomains (numbered 1–9 in Fig. 4) coupled at their mutual boundaries, where the ray computations were stopped. The incident wave conditions along the offshore model boundary are based on linear interpolation of spectral data from buoys X6 and 44014. The energy balance [Eq. (2)] was integrated with a 10-min time step, using source terms that represent wave–bottom Bragg scattering (Ardhuin and Herbers 2002) and bottom friction over a movable bed. The bottom friction source term is a slightly modified version of the formulation proposed by Tolman (1994) that was tuned to improve the hindcast skill of swell wave heights observed during SHOWEX and the earlier DUCK94 experiment (Ardhuin et al. 2003a). Predictions both with and without Bragg scattering were made to quantify the contribution of this process to directional diffusion. Although the model is fully spectral, comparisons presented here are restricted to the mean direction *θ _{p}* and directional spread

*σ*

_{θ,p}at the peak frequency. Analysis of spectra suggests that bottom friction does not have a strong effect on the spectral distribution of wave energy (Ardhuin et al. 2003b), while wave–bottom scattering is a linear process that does not transfer energy between wave components that have different frequencies.

The mean direction and directional spread in the previous academic examples correspond to values at the peak frequency observed on 15 September at the offshore buoy X6 when narrowband swell arrived from Hurricane Floyd at large oblique angles (150 < *θ _{p}* < 160°, where

*θ*is direction from, in nautical convention, and

_{p}*σ*

_{θ,p}= 15°; see Fig. 5a). Swell from Hurricane Gert on 17–21 September was similar, except for the mean wave direction (

*θ*= 120°) that was closer to normal incidence to the coastline.

_{p}Differences between observations at buoys X6 and 44014 (both located on the outer edge of the shelf) on 15 September suggest that errors may result from variations of the deep water wave field that are poorly represented in the model incident wave conditions, interpolated from the measured spectra at X6 and 44014, spectra from the latter being first “back refracted” to deep water. Additional uncertainty in the model predictions results from applying a uniform bottom-elevation spectrum in the Bragg scattering source term, based on limited high-resolution bathymetry data. Despite these uncertainties the predictions of *σ*_{θ,p} obtained with the full-Lagrangian scheme including Bragg scattering agree well with the observations, indicating that the dominant physical processes are well represented in the model. A more detailed analysis of over 50 days of swell-dominated conditions is given by Ardhuin (2001) and Ardhuin et al. (2003a, b). A systematic underprediction of directional spread at nearshore sites is the likely result of other scattering processes not included in the present model, such as nonlinear wave–wave interactions (see Herterich and Hasselmann 1980; Herbers and Burton 1997; Zaslavskii and Polnikov 1998).

Integrating the model with only the bottom friction for source term and the full ray advection scheme, the directional spread decreases dramatically across the shelf (Figs. 5b,c, diamonds), in particular for the large oblique incidence angles on 15 September. This directional narrowing is, however, weaker than on the alongshore uniform topography since refraction on the two-dimensional shoals tends to increase directional spread. Using the multistep ray advection scheme with *N* = 1 introduces numerical diffusion that occasionally doubles *σ*_{θ,p} at X1 (cf. squares and diamonds in Fig. 5). At X4 (Figs. 5b,c) the numerical diffusion is less severe due to the shorter propagation distance. For the most oblique incident waves (15 September) the broadening of directional spectra due to numerical diffusion for *N* = 1 is larger than that predicted by wave–bottom Bragg scattering (triangles in Fig. 5). It suggests that this source term should be used with caution when combined with numerically diffusive advection schemes, such as piecewise ray methods (*N* = 1), or first-order finite-difference schemes. However, earlier on September 14, as well as during the period 17–21 September, the physical Bragg scattering diffusion appears to be stronger than numerical diffusion, even for *N* = 1, suggesting that diffusive low-order numerical schemes may yield reasonable predictions for these more commonly observed broader wave spectra.

## 4. Summary and conclusions

A multistep ray advection scheme is proposed here for phase-averaged wave models. It is based on the Lagrangian energy balance [Eq. (2)] and generalizes the fully Lagrangian advection scheme of the CREST wave model (Ardhuin et al. 2001) by reducing the time interval over which energy is advected along a given ray trajectory. This time interval, a multiple *N* of the predefined time step Δ*t* used for the source-term integration, can be adjusted between the original fully Lagrangian scheme (*N* = ∞, in practice limited by the model domain size) and the efficient but diffusive piecewise ray methods (*N* = 1; Sobey 1986). This scheme belongs to the class of semi-Lagrangian schemes with the earlier INTERPOL scheme of Lavrenov and Onvlee (1995), and it is unconditionally stable, can be implemented on an arbitrary irregular spatial grid, and effectively eliminates the garden sprinkler effects of finite-difference schemes by using finite bandwidths.

The finite-bandwidth representation used in CREST, computing as many rays as necessary to resolve their variability, is probably more accurate (and slower) than the angular smoothing technique used in INTERPOL for the shallow-water cases described here with strong refraction effects. However, in deeper water the latter technique is probably more efficient. For large values of *N*, a large number of rays is needed to describe a finite bandwidth, resulting in an expensive precomputation effort. The effort of both computing rays and integrating the energy balance equation is reduced drastically in a low-*N* scheme requiring only a single or a few rays to describe a band. This efficiency makes possible a frequent update of the rays to represent the variable refraction due to unsteady currents and water levels, already performed in the TELEMAC-based Operational Model Addressing Wave Action Computation (TOMAWAC) piecewise ray model (M. Benoit 2001, personal communication). In the academic test, presented in Fig. 2, the *N* = 1 single-ray advection scheme was 5 and 2.5 times faster than *N* = 20 and 4 multistep schemes with finite-bandwidth representations, yielding a maximum error of *σ _{θ}* of 35%, compared with 4% and 15% for

*N*= 20 and 4, respectively (Fig. 6).

The *N* = 1 finite-bandwidth scheme is both less accurate (55% error) and slower (by a factor of 1.5) than the *N* = 1 single-ray scheme but has the benefit of eliminating the garden sprinkler effect that may appear in large domains with regular spatial grids. On the real shelf bathymetry the *N* = 1 finite-bandwidth ray advection scheme was already quite efficient, making the entire computation (including source-term computation and interpolation) 4 times faster than the original CREST scheme.

The effectiveness of the present semi-Lagrangian scheme in reducing numerical diffusion was demonstrated with academic tests and observations of the evolution of the directional spread *σ _{θ}* of swell across the shelf during the Shoaling Waves Experiment in 1999. Although this parameter is rarely given much attention, it has important implications for nearshore morphodynamics (Reniers et al. 2004) and forces on structures. A single-step ray advection scheme may artificially broaden the directional spectrum, sometimes more so than the predicted natural diffusion resulting from wave scattering by the bottom topography. These errors are effectively suppressed by increasing

*N*(e.g., a four-step scheme with a 600-s time step was sufficient for the North Carolina shelf). These results demonstrate the usefulness of semi-Lagrangian schemes for studying natural scattering processes that may be buried in the numerical noise when other schemes are used.

Further work is needed, including model intercomparisons for realistic wind–sea cases and verification of predicted full spectra and other wave parameters (e.g., Lavrenov and Onvlee 1995), to evaluate the potential benefits of the scheme presented here, in terms of both accuracy and cost, for routine wave forecasting or hindcasting.

## Acknowledgments

This research is supported by the Coastal Dynamics Program of the U.S. Office of Naval Research. Bill O’Reilly, Paul Jessen, and Mark Orzech contributed significantly to the field experiment planning and deployment, and the wave and bathymetry data analysis, and made fruitful comments on the present work. Buoys X1–X6 were deployed by staff from the Naval Postgraduate School and the Scripps Institution of Oceanography Center for Coastal Studies. Wave data from the 8-m-depth array and buoy 44014 were provided by the Field Research Facility of the U.S. Army Corps of Engineers Waterways Experiment Station’s Coastal Engineering Research Center, and the National Data Buoy Center, respectively. Permission to use these data is appreciated. Anonymous reviewers made helpful suggestions.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Dr. Fabrice Ardhuin, SHOM/CMO/RED, Vagues et domaine littoral, 13 Rue du Chatellier, Brest 29609, France. Email: ardhuin@shom.fr