## Abstract

A surface wave model is developed with the intention of coupling it to three-dimensional ocean circulation models. The model is based on a paper by Mellor wherein depth-dependent coupling terms were derived. To be compatible with circulation models and to be numerically economical, this model is simplified compared to popular third-generation models. However, the model does support depth and current refraction, deep and shallow water, and proper coupling with depth-variable currents.

The model is demonstrated for several simple scenarios culminating in comparisons of model calculations with buoy data during Hurricane Katrina and with calculations from the model Simulating Waves Nearshore (SWAN); for these calculations, coupling with the ocean was not activated.

## 1. Introduction

This paper follows a paper by Mellor (2003, hereafter M03), which, however, has been revised (Mellor 2008); the revisions did change Eqs. (6)–(8) below but do not affect any calculated results in this paper since coupling with an ocean has not been activated. The phase-averaged, wave–current equations of motion were extended to the third vertical dimension. In much of the literature (e.g., Phillips 1977), the wave interacting continuity and momentum equations were, a priori, vertically integrated, rendering them unsuitable for coupling with depth-dependent numerical ocean circulation models. Now, as a consequence of (the revised) M03, it is possible to couple three-dimensional circulation models with wave models; the coupling includes depth-dependent wave radiation stress terms, Stokes drift, vertical transfer of wave-generated pressure transfer to the mean momentum equation, wave dissipation as a source term in the turbulence kinetic energy equation, and mean current advection and refraction of wave energy.

There exist functional third-generation wave models such as Wave Model (WAM) (WAMDI Group 1988; Komen et al. 1994), WAVEWATCH (Tolman 1991), and Simulating Waves Nearshore (SWAN) (Booij and Holthuijsen 1999). However, for many users, these models cannot practically be coupled to circulation models. Circulation models invoke four independent variables—three horizontal and vertical coordinates and time, *x*, *y*, *z*, *t*—whereas the wave models use five independent variables—two horizontal coordinates and time, *x*, *y*, *t*, plus the wave propagation angle, *θ*, and frequency, *σ*. Consider frequency as the additional variable with, say, 30 numerically discrete, frequency bins; then add time for the computation of nonlinear wave–wave interaction and the integration of various properties including the new coupling terms (M03; Mellor 2005). Thus, one is faced with about a two orders of magnitude increase in computational effort over that required by circulation models cum sole. This is the actual experience obtained with an operational model of the New York Harbor coastal and estuarine region (www.stevens.edu/maritimeforecast/). The model builders (A. Blumberg 2006, personal communication) found that for the same horizontal grid and length of run, the SWAN model required 86 times the computer run time compared to the circulation model, in this case the Princeton Ocean Model (POM); the two models were not coupled. Circulation model applications are typically executed with marginal horizontal resolution commensurate with available computational resources so that the two orders of magnitude increased computational time will be prohibitive for many coupled circulation–wave model applications.

In the literature, one finds that Mastenbroek et al. (1993) did couple the WAM model to a surge model and was able to incorporate the vertically integrated, wave radiation stress terms into the shallow-water continuity and momentum equations; this was followed by a similar paper by Ozer et al. (2000) that included tides and surges. Zhang and Li (1997) coupled a wave equation (without frequency and directional dependence) to a three-dimensional barotropic ocean model. They did solve a separate transport equation for wavenumber and, thus, frequency. A simple eddy viscosity parameterization allowed wind stress momentum transfer into the water column. Xie et al. (2001) ran POM and WAM simultaneously and coupled the wave and circulation models only through surface and bottom stress roughness terms; they performed a series of 1-day runs to determine the effect of wave-enhanced stresses on the currents. Moon (2005) combined POM and WAVEWATCH III; and a wave stress was obtained from an integral of the wave spectrum (see relevant discussion in sections 5 and 8 below). Also, wave dissipation produced a turbulence source term for the Mellor and Yamada (1982) turbulence closure model (Mellor and Blumberg 2004). All of the above papers were for deep water, for limited regions, and for limited run times. Since vertically dependent radiation stress terms were not available, they were excluded from all wave energy and momentum equations, except for the vertically integrated surge model. Surface-weighted Doppler velocities derived from depth-dependent currents and pressure transfer of wind stress were also excluded.

Thus, there is a need for a numerically practical wave model that can be coupled to three-dimensional, time-dependent circulation models for long-term simulations and high-resolution, regional, basin, or global domains. The wave model could supply estimates of wave heights (useful, e.g., in parameterizing wind drag coefficients as well as other estimates of sea state; see Janssen et al. 2004) and depth-dependent wave radiation stress terms for the momentum equation, which would feed back depth-dependent velocity fields to the wave energy equation.

A possible candidate is the simple Great Lakes Environmental Research Laboratory (GLERL) model [Schwab et al. (1984), originally devised by Donelan (1977)], which is used by at least two regional operational models known to us. Despite its simplicity the GLERL model has been shown to be relatively consistent with observations in several studies (Schwab et al. 1984; Lin et al. 2002). However, that model precludes swell, shallow-water effects, refraction, explicit wave dissipation, and exchange of current and Stokes drift momentum.

Here, to include these attributes, we have developed a model that borrows a feature of the GLERL wave model and other models (SWAMP Group 1985) in that the energy distribution in frequency space is parameterized using the spectrum by Donelan et al. (1985, henceforth DHH85), and which contains elements of the Joint North Sea Wave Project (JONSWAP) spectrum (Hasselmann et al. 1973). The model avoids dealing with the wave–wave interaction process (the process is compensated by the specified spectral shape), and has the same level of complexity as circulation models whose four independent variables are *x*, *y*, *z*, *t*; the present wave model also has four independent variables, *x*, *y*, *θ*, *t*; the wave propagation angle, *θ*, will account for refraction due to bottom depth and current variations. Frequency is also dependent on *x*, *y*, *θ*, *t*; in the wind-driven portions of *θ*, the transport equation for frequency is forced by and asymptotes in time to the peak frequency of the parameterized spectrum; otherwise, the frequency is advected as swell. There is some similarity with this approach and the HISWA model (Holthuijsen et al. 1989) in that the main dependent variables are wave angle–dependent energy and frequency. The model is meant for stationary calculations for nearly fixed mean direction—the grid is oriented in the mean wave direction, say toward a coastline or a harbor—and is numerically not applicable to arbitrary basins or wind directions. Source and sink terms are derived differently than described below.

The numerical wave model code is composed of subroutines; the main subroutine can be conveniently called by a circulation model, and will approximately require the same computer resources as a circulation model.

Another motivation for the model is to directly confront unresolved research issues. For example, in M03, it was seen that some or all of the momentum transfer to the immediately underlying surface boundary layer is due to pressure transfer; surface boundary layer models currently assume that momentum transfer into the water column is entirely due to turbulence mixing.

Constraining to a specific spectral shape will, of course, not always conform to measured spectra; lost accuracy in, say, the prediction of significant wave height and other integral wave properties remains to be determined. (However, see section 9 below.)

We face a persistent and common dilemma: whether to define wave age as *c _{p}*/

*U*

_{10},

*c*/

_{p}*u*

_{*}, or

*c*/

_{p}*U*

_{p,λ/2}, where

*c*

_{p}is the phase speed at the peak of the wave spectrum,

*U*

_{10}is the wind speed at 10 m above the sea surface,

*u*

_{*}is the friction speed, and

*U*

_{p,λ/2}is the wind speed at a half-wavelength above the sea surface at the peak frequency. Since the main parameter in the spectrum according to DHH85 is

*c*/

_{p}*U*

_{10}, we will mostly follow the same practice. [However, as discussed below,

*c*/

_{p}*U*

_{10}will be converted to

*g*/(

*σ*

_{p}

*U*

_{10}), an equality for deep water but not for shallow water.] Alves et al. (2003) offer evidence that

*c*/

_{p}*U*

_{10}is the preferable wage age descriptor and from a practical point of view

*u*

_{*}is difficult to measure, but the issue does remain uncertain.

In this paper, we develop the wave model cum sole, a necessary precursor to coupling a wave model to a three-dimensional circulation model.

Let *E*_{σ,θ} = *E*_{σ,θ}(*x*, *y*, *t*, *σ*, *θ*) be the directional spectrum of kinematic energy (divided by the water density, *ρ*_{w}), a function of *x*, *y*, *t*—a point in horizontal space and time—and *σ*, *θ*, the frequency and wave direction. Henceforth, the arguments *x*, *y*, *t* will be deleted. In this paper we calculate only

and since the kinetic and potential energies are equal,

is the total kinematic wave energy per unit surface area: the product of the gravity constant and phase-averaged, squared wave elevation summed over all frequencies and directions. We restrict our attention to surface gravity waves (i.e., wavelengths in excess of about 10 cm).

In section 2, the linear wave relations are reviewed. The model is developed in sections 3–7. An array of canonical tests is executed in section 8. The model is then applied to the Gulf of Mexico during the period of Hurricane Katrina; the path of the hurricane passed near six National Data Buoy Center (NDBC) buoys to which model calculations are compared. Parallel calculations and comparisons are made using the SWAN model. Section 10 is a summary. Appendixes contain some detailed information.

## 2. Monochromatic unidirectional equations and definitions

We first define terms for a monochromatic unidirectional wave field pursuant to dealing with spectra as a function of frequency and wave propagation angle. Thus,

where *k*_{α} = *k*( cos*θ*, sin*θ*) is the wavenumber vector and *k* = |*k*_{α}|; *θ* is the wave propagation direction relative to the eastward direction; *σ* is the intrinsic frequency; the Doppler velocity, *u*_{Aα}, will be defined below; *c* is the phase speed; *g* is the gravity constant; *D* = *h* + η̂ is the water column depth where η̂ is the mean (phase averaged) surface elevation and *h* is the bottom depth. The group velocity is

For the present wave model, the dispersion relation (2b) is initially solved iteratively, inverted, and cast in the form *kD* = *f*_{1}(*σ*^{2}*D*/*g*); and from (3b), *n* = *f*_{2}(*σ*^{2}*D*/*g*). A lookup table with interpolation comprises a subroutine in the numerical code.

As cited in Komen et al. (1994), the refraction speed is *c*_{θ} = −|**∇**_{h}*ω* × **k**|/*k*^{2}, where *ω* is obtained from (2a). Working out the vector algebra (Golding 1978) yields

Before defining more terms, useful combinations of hyperbolic sines and cosines are

The “sigma” variable is *ς* = (*z* − η̂)/*D* (reserving *σ* for frequency) so that *D*(1 + *ς*) = *D* + *z* − η̂.

A Doppler (or advective) velocity is required for the wave energy equation and is

where *U*_{α} = *U*_{α}(*x*, *y*, *ς*, *t*) is the ocean current plus the Stokes drift. The weight factor in the square brackets is significantly nonzero only in the wave portion of the water column. Equation (6) is from M03 (and its revision cited in section 1).

The wave radiation stress terms are

where a modified delta function, *E*_{D}, is defined such that

## 3. The wave energy equation and the specified spectrum

After integrating the full spectral equation (a function of frequency and wave angle) with respect to frequency, we arrive at

The horizontal coordinates are denoted by *x*_{α} = (*x*, *y*). The overbars represent spectral averages (and differ from the phase averaging usage in M03). The first two terms on the left of (8) determine the propagation of wave energy in time and horizontal space whereas the third term is the refraction term accounting for the change in direction of wave energy propagation. The last term on the left of (8) represents energy exchange with the mean velocity energy equation. All of the terms in (8) are functions of *θ*: *S*_{θin} is the wave energy source term dependent on wind properties; *S*_{θSdis} and *S*_{θBdis} are wave dissipation due to wave processes at the surface and bottom, respectively. All terms are kinematic (i.e., energy terms are divided by water density). Thus the atmospheric work done on the water is *ρ*_{w}*S*_{θin}, where *ρ*_{w} is the seawater density. The terms, *c*_{gα}, *c*_{θ}, *u*_{Aα} and *S*_{αβ}, differ from their counterpart terms in section 2 in that they have been spectrally averaged as described in section 4.

Following DHH85 and for the wind-driven part of the spectrum, we stipulate

Frequency is *σ* and *σ*_{p} is the peak frequency where *E*_{σ,θ} is maximum; *θ* is the wave propagation direction; *θ* is the mean wave propagatin angle [see (24)]. Deferring to DHH85 (535–538) or Komen et al. (1994, p. 187) for details, *β* = *β*(*σ*/*σ*_{p}) and the other parameters in (9) are wave age dependent such that *α* = *α*(*U*_{c}/*c*_{p}), *γ* = *γ*(*U*_{c}/*c*_{p}), *σ*_{pw} = *σ*_{pw}(*U*_{c}/*c*_{p}), and *U*_{c} = *U*_{10} cos(*θ*_{w} − *θ*); *U*_{10} is the wind speed evaluated at the 10-m height and *θ*_{w} is the wind angle. The observations in DHH85 showed that *θ* can differ from*θ*_{w}. In this model building process, we will use (9a) and (9b) as a weighting function where, however, retention or neglect of the difference between *θ*_{w} and *θ* makes little difference in the calculated results; thus we set *U*_{c} = *U*_{10}. Using (9a) and (9b), various spectrally weighted averages were obtained by numerical integration. Following Terray et al. (1996) and Banner (1990), the spectra were extended beyond 3.5*σ*_{p} by appending a *σ*^{−5} tail to 20 *σ*_{p}, whence the integration was terminated.

### a. Peak frequency

Following DHH85, to obtain the peak frequency for wind-driven waves, numerical integration of (9a) and (9b) over *σ* and *θ* yields *σ*^{4}_{p}*E*_{TW}/*g*^{3} as a function of *U*_{10}/*c*_{p}, which, for deep water, equals *σ*_{p}*U*_{10}/*g*. It is believed (Hwang and Wang 2004) that the latter term provides results that are less dependent on depth than is the former (for waves propagating onto a beach with no currents and light winds, *σ*_{p}*U*_{10}/*g* ≅ constant whereas *U*_{10}/*c*_{p} → ∞ as *k*_{p}*D* → 0). We define *E*_{TW} as the integral of the portion of *E*_{θ} that is wind driven; until we later encounter the need to separate swell from wind-driven portions of *E*_{θ}, *E*_{TW} = *E*_{T}.

The calculated integral from (9a) and (9b) is provided nondimensionally in Fig. 1 as a solid line. By analyzing data from Lake St. Clair, Donelan et al. (1992) found that

where *C*_{σ} = 0.0022. This formula also agrees quite nicely with a plot by Hwang and Wang (2004, their Fig. 5c) for the same nondimensional variables. Their study included fetch-limited and duration-limited data. Young (2006) analyzed hurricane data from a measurement site located off the northwest coast of Australia. He found good agreement with (10); citing DHH85, he noted, “the fact that both the one-dimensional and directional spectra are very similar to spectra reported under simple unidirectional winds is interpreted as being the result of the shape stabilization effects of non-linear interactions.” Refer to these papers to see data scatter around the average.

Equation (10) can also be written as

and will be used to obtain *σ*_{p}. Next a small modification of (11) is possible by noting that dimensional analysis yields fcn[*E*_{TW}*g*/*U*^{4}_{10}, *U*_{10}*σ*_{p}/*g*, (*E*_{TW}/*g*)^{1/2}/*z*_{10}] = 0, where the length scale, (*E*_{TW}/*g*)^{1/2} = *H*_{S}/4, divided by *z*_{10} = 10 m is now added to the group of relevant nondimensional variables. Here *H*_{S} = 4(*E*_{T}/*g*)^{1/2} is the significant wave height (Longuet-Higgins 1952). [Note that *z*_{10} is important in the conventional velocity-dependent drag coefficient *C*_{D} = *C*_{D}(*U*_{10}), which in dimensionless terms should be—but usually is not—written *C*_{D} = *C*_{D}(*U*_{10}/*z*_{10}*g*).]

The possibility that (*E*_{TW}/*g*)^{1/2}/*z*_{10} might be significant led to regression analyses using buoy data discussed in section 9 and model data comparisons. The results showed a weak but nevertheless discernable dependence such that replacing the constant in (11) with *C*_{σ} = 0.0015 + 0.0079(*E*_{TW}/*g*)^{1/2}/*z*_{10} resulted in some improvement in model wave periods for hurricane-scale winds in comparison with the buoy data discussed in section 9. For intermediate wind, *C*_{σ} ≅ 0.0022.

## 4. Spectrally averaged terms

Next we deal with terms that are predominantly independent of wave angle. For example,

is, strictly speaking, a function of *θ*. However, a common approximation is that *E*_{σ,θ} = fcn(*θ*)*E*_{σ}. An examination of (9) shows that the two independent variables are not exactly separable, but nevertheless trials using (12a) show that the approximation is sufficiently accurate such that

where *E*_{σ} = ∫^{π/2}_{−π/2 }*E*_{σ,θ }*d**θ*. Similarly, integrations were carried out for other terms such that, in place of (4), (6), and (7), we have

and, similarly to (7), we have

In Fig. 2a, *c*_{g}/*c* and *F*_{cθ} are plotted as functions of *k*_{p}*D*. These variables are also dependent on inverse wave age, but the dependence is weak (a few percent) and will be neglected henceforth. For comparison, Eq. (3b) where *kD* is replaced by *k*_{p}*D* is also plotted in Fig. 2a.

The *F*_{n} functions in (13b) and (13c), related to the definitions (5), are explicitly defined in appendix A and are functions of *ς* and *k*_{p}*D* as shown in Fig. 2b for *σ*_{p}*U*_{10}/*g* = 2. The variations with respect to inverse wave age are small (in the range, 1 < *σ*_{p}*U*_{10}/*g* < 5, mostly less than ±5% with a very few values near ±10%) and henceforth are neglected. Appendix A provides further information.

The wave–current interaction terms in (8) are complicated and an explicit record of the terms prior to coding is useful. Thus,

Recall that *U*_{α} = (*U*, *V*) is the ocean current plus Stokes velocities.

## 5. Wind growth source term

The friction velocity is obtained from the well-known relations

where **U**_{10} − **U**(0) is the vector difference between the 10-m wind velocity and the surface ocean current. The airside and waterside friction velocities are *u*_{*a} and *u*_{*w}, such that *ρ*_{w}*u*^{2}_{*w} = *ρ*_{a}*u*^{2}_{*a}; *ρ*_{w}/*ρ*_{a} = 860 is the water to air density ratio, (15c) is from Donelan (1990), and *H*_{S} = 4(*E*_{T}/*g*)^{1/2} is the significant wave height.

We next obtain *S*_{θin}, which can be written

*E*_{σ,θ} is given by (9a) and (9b). The wave growth according to Donelan (1999) is

where *U*(*λ*/2) is the wind speed evaluated at the half-wavelength height. Equation (17) has some similarity to an expression from Janssen (1989) although he used the airside friction velocity *u*_{*a} instead of *U*(*λ*/2); the latter can be obtained from *U*_{10} using the law of the wall according to *U*(*λ*/2)/*U*_{10} = ln(*λ*/2*z*_{0})/ ln(10 m/*z*_{0}), where *λ* = 2*π**g**σ*^{−2} tanh *kD*. Sample plots of *E*_{σ} = ∫^{π}_{−π }*E*_{σ,θ }*d**θ* using (9) and *B*_{σ} = ∫^{π}_{−π }*B*_{σ,θ }*d**θ* from (17) are shown in Fig. 3. Thus it is seen that the wind source term is shifted toward large frequencies since, for deep water, *c*^{−1} = *σ*/*g*.

Equation (16) was integrated to obtain the so-called spreading function *f*_{spr} = *S*_{θin}/*S*_{Tin}, where *S*_{Tin} = ∫^{π/2}_{−π/2 }*S*_{θin }*d**θ*; the results are plotted in Fig. 4. The dashed line is

for *β* = 2.2. The spreading function is similar to that in (9) after replacing *θ* with *θ*_{w} since the weighting function, used to find a local wave energy source term, *S*_{θin}, should depend on the local wind direction. As will be seen, the final model will produce calculations wherein *θ* differs from *θ*_{w} because of nonlocal effects. In (18a), the sech function is quite small when |*θ* − *θ*_{w}| = *π*/2; nevertheless, the cutoff for |*θ* − *θ*_{w}| > *π*/2 improved the calculations in section 7 relative to data for small fetch.

*S*_{Tin} can be normalized in a number of ways—by *σ*_{p}*E*_{T} or by various combinations of *c*_{p} or *u*_{*w}. After normalization by *c*_{p}*u*^{2}_{*w}, Terray et al. (1996) integrated *B*(*σ*, *θ*)*E*(*σ*, *θ*) over frequency and angle; they also used (17) but used observed spectra in place of (9); our results (not shown) are very similar to theirs when similarly normalized.

Normalizing by *u*_{*w}, *S*_{Tin}/*u*^{3}_{*w} versus *σ*_{p}*U*_{10}/*g* is plotted in Fig. 5a. This form is now familiar to modelers of ocean surface boundary layers. Noting that most of the wave energy source is directly dissipated into turbulence, values of *S*_{Tin}/*u*^{3}_{*w} = 100 were empirically deduced by Craig and Banner (1994) and a value of 150 by Stacey (1999) as a source term in the turbulence energy equation component of their surface boundary layer model. These values are consistent with young inverse wave ages in Fig. 5a.

The dashed line in Fig. 5a is given by *S*_{Tin}/*u*^{3}_{*w} = 370 exp(−0.33*σ*_{p}*U*_{10}/*g*) but the calculated values abruptly decrease around *σ*_{p}*U*_{10}/*g* = 0.5 and become negative for lesser values. When *S*_{Tin}/*c*^{3}_{p} is plotted versus *σ*_{p}*U*_{10}/*g* as in Fig. 5b, the transition from positive to negative values is smooth. Furthermore, the negative values, occurring when the wave speed exceeds the wind speed as manifest in Eq. (17), are small and *S*_{Tin}/*u*^{3}_{*w} is large only because *u*^{3}_{*w} is very small. Thus, we let

In Fig. 5, *S*_{Tin}/*u*^{3}_{*w} is, fortuitously, a weak function of *U*_{10}, whereas *S*_{Tin}/*c*^{3}_{p} is dependent on *U*_{10}. Note that a fully developed wave field is often taken to be *U*_{10}/*c*_{p} = *σ*_{p}*U*_{10}/*g* ≅ 0.83 and is considered to be an approximate number.

## 6. Wave dissipation

The total wave dissipation must be determined empirically. A model for white capping or wave breaking is

where *a* and *b* will be determined by reference to fetch data. The first term represents the fact that the high-frequency part of the spectrum is dissipated very nearly in situ and the second part is dissipation of the mid- (*σ* ≈ *σ*_{p}) to low-frequency part of the spectrum. This means, of course, that overall wave growth only responds to (1 − *a*)*S*_{θin} − *bE*_{θ}*σ*_{p}; nevertheless, the full dissipation will be needed as input to the turbulence kinetic energy equation when this wave model is coupled to a circulation model. [To the second term in (19), we initially added a factor involving the wave slope but finding that results were insensitive to the addition, it was withdrawn.] An estimate of *b* as a function of *a* may be obtained by equating *S*_{θSdis} = *S*_{θin} for a fully developed sea when *σ*_{p}*U*_{10}/*g* = 0.83. Thus, dominant leverage on computed results are via the parameter, *a*, as determined in section 8.

There are many prescriptions for bottom frictional dissipation, generally for turbulent flow. There are uncertainties, but following the discussion in Booij et al. (1999) we adapt to the present model such that

where *u*_{b} is the wave amplitude near the bottom. This formula agrees with that in Mellor (2002), where *C*_{b} is a function of nondimensional bottom roughness. Bottom roughness is generally unknown; here we suggest a typical value, *C*_{b} = 0.003. It should be noted that none of the suggested formulations in the wave modeling literature include the “streaming” effect described by Longuet-Higgins (1953).

The bottom can induce wave breaking if *H*_{s}/*D* is large enough. Following Battjes and Janssen 1978 [see appendix in Booij et al. (1999) for detailed discussion], we add to (20a) such that

where the fraction of breaking waves out of an ensemble, *Q*_{b}, is given by the transcendental relation (*Q*_{b} − 1)/ ln*Q*_{b} = 8(*E*_{T}/*g*)/*H*^{2}_{M} and where *H*_{M} = *γ**D*; *γ* is an empirically adjusted constant.

Finally, *S*_{θBdis} = *S*^{(1)}_{θBdis} + *S*^{(2)}_{θBdis} for use in Eq. (8).

## 7. A directionally dependent frequency

A *θ*-dependent frequency is

where ∂*σ*_{θ}/∂*k* = *c*_{g} and ∂*σ*_{θ}/∂*D* = (*σ*_{θ}/*D*)(*n* − 1/2). Equation (21a) is derived from wavenumber irrotationality, ∂*k*_{α}/∂*x*_{β} − ∂*k*_{β}/∂*x*_{α} = 0, and the conservation of crests, ∂*k*_{α}/∂*t* + ∂*ω*/∂*x*_{α} = 0. In regions of *θ* that are wind driven ( *f*_{spr} > 0), the additional source term, ℜ, prescribed by

has the effect of “nudging” (a term used, e.g., in data assimilation of various ocean properties) *σ*_{θ} toward *σ*_{p}.

## 8. Model tests

The model described above and whose independent variables are *x*, *y*, *t*, and *θ* was finite differenced and coded in FORTRAN replacing *x*, *y*, and *θ* by *i*, *j*, and *m*. Numerical details are in appendix B.

All of the tests in this section are independent of the coordinate, *y*. This was implemented by using only five rows in the *j* direction and setting cyclic boundary conditions on the bounding *j* cells. Alternately, a version of the code was created by stripping out the *j* terms. As a check on the code, the two versions gave the same results.

### a. Refraction for a monochromatic wave train

We first test the model against Snell’s Law for the case *U*_{α} = 0 and *U*_{10} = 0. At *x* = 0, incoming swell is prescribed according to *E*_{θ} ∝ (*β*/2) sech^{2}[*β*(*θ* − *θ*_{0})], where *θ*_{0} = 60° and *β* = 4.0, a rather narrow distribution with respect to *θ*. Equation (8) with a null right side is then solved to steady state. The time, spatial, and angle increments are 20 s, 500 m, and 2*π*/24 = 15°, respectively.

The bottom topography is shown in Fig. 6a. The frequency is *σ* = *σ*_{p} = 2*π*/10*s*, which according to (21) and conservation of crests, is constant; the corresponding *c*(*x*) and *c*_{g}(*x*) are plotted in Fig. 6b. The solution to Eq. (8) is given in Fig. 6c. Note that the angle domain is −*π* < *θ* < *π* but only the active portion is shown. In Fig. 6d, the total wave energy per unit surface area *E*_{T} = Σ^{2π/Δθ}_{k=1}*E*_{θ}(*θ*_{k})Δ*θ*, and is normalized with its value at *x* = 0. The mean angle, *θ* = Σ^{2π/Δθ}_{k=1}*θ**E*_{θ}(*θ*_{k})Δ*θ*/*E*, is shown in Fig. 6e.

Calculations are obtained from Snell’s law so that

and *m* is the label on each ray emanating from *x* = 0 where the initial distribution is as stated above. Averages are obtained on *m* and are plotted as the dashed lines in Figs. 6d,e. Agreement between the model and Snell’s law is improved further (the two curves are nearly indistinguishable) by decreasing the angle increment from 15° to 10°.

### b. Fetch-limited waves

The model was tested against the growth of waves for a constant offshore wind normal to a straight north–south coast located at *x* = 0. For this problem, similarity growth relations were formulated by Kitaigorodskii (1962) using dimensional analysis. As reported in Komen et al. (1994, p. 181), Kahma and Calkoen examined data from the JONSWAP experiment and data observed in the Bothnian Sea and Lake Ontario (DHH85). They first separated data into winds when the vertical density stratifications were stable or unstable. However, we deal only with the composite dataset that they represented by

and that is plotted as a dashed line in Fig. 7. The well-developed limit of *E*_{T}*g*/*U*^{4}_{10} = 3.6 × 10^{−3} is from Pierson and Moskowitz (1964).

Another fit to Lake St. Clair data by Donelan et al. (1992), using an elaborate analysis scheme to minimize inhomogeneities in the data, is shown as a dotted–dashed line in Fig. 7.

The solid lines are from the present model for two different values of *U*_{10} and for the best-fit parameters, *a* = 0.925, *b* = 0.18 × 10^{−4}; these values are hereafter held constant except for swell portions of *E*_{θ} ( *f*_{spr} < 0) as described below, where *b* is reduced by a factor of 5. The time and angle increments are 10 s and 15° respectively; at *x* = 0, the spatial step is Δ*x* = 100 m but thereafter Δ*x*(*i* + 1) = 1.10Δ*x*(*i*).

Whereas the data syntheses exclude dependence on *U*_{10}, there must be some model dependence on *U*_{10} or, in dimensionless form, *U*_{10}[*g* (10 m)]^{−1/2} according to dimensional analysis or from (15a)–(15c), (16), and (17) together with the relation *H*_{S} = 4(*E*/*g*)^{1/2}.

In Fig. 8 we show calculations of the time and distance development of wave growth together with a synthesis of data by Hwang and Wang (2004). Since duration-limited data are scarce and difficult to obtain, they deduced duration-limited growth from fetch-limited growth and compared with their own and other datasets; a synthesis of their duration-limited data is also shown in Fig. 8. A comparison (not shown) of nondimensional peak frequency *σ*_{p}*U*_{10}/*g* versus *xg*/*U*^{2}_{10} was also quite favorable.

Figures 9a,b demonstrate the fetch-limited behavior for *U*_{10} = 10 m and for different wind angles relative to the eastward direction; the coastline is north–south at *x* = 0. For wind angles not equal to zero and, because of the spread of the wind growth term in (18) denoted by *f*_{spr}, waves with mean propagation angles larger than the wind angle will propagate over a longer fetch and therefore accrue higher energies than waves with lower angles; the mean wave angle therefore differs from the wind angle until a considerable distance from the coast where the two angles coincide. The mean wave angle is given by

For small wind angles, say *θ*_{w} = *π*/6 = 30° at the coast (*x* = 0), energy input near the wind angle is small since the effective fetch is small and the flow angle is dominated by wave propagation components around *θ* = 90° corresponding to very large fetch. Consider the extreme case of wind parallel to the coast, *θ*_{w} = *π*/2 = 90°. At the coast, the wind only produces waves in the range 0 < *θ* < *π*/2 whereas, in the far field (*x* → large), the wind produced the full range 0 < *θ* < *π*; thus the mean wind angle is larger and the energy smaller at the coast than the far field. Relatively, energy propagating from the far field to the coast has dissipated before reaching the coast.

### c. Refraction due to currents

A northward Gulf Stream–like jet is prescribed according to

where *L* = 400 km is the domain zonal width. The spatial, temporal, and angle increments are 5 km, 500 s, and 15°, respectively. Since the waves are confined to the near surface, the depth dependence is not important in this application, but the code generally does account for depth dependence through the interaction terms in (8).

A background, fully developed wave field is first established in the absence of the jet with a wind speed of *U*_{10} = 10 m s^{−1} and varying directions, northward (90°) through southward (−90°); wave energy and the mean propagation angle are shown in Fig. 10 as dashed lines. Here, lateral boundary conditions are devised so that fully developed wave energy propagates through the domain boundaries. The mean propagation angles are the same as the wind direction.

With the jet in place, the solid curves in Fig. 10 show the deviations of wave energy and mean angle. Analysis of the calculations show that, in this application, the term (*k*_{α}/*k*) sin*θ* ∂*u*_{Aα}/∂*x* = sin^{2}*θ* ∂*V*/∂*x* in (13a) is mostly responsible for the current-induced deviations. However, the fact that the friction velocity in (15a) depends on the wind speed relative to the surface water speed has some effect; for example, in the −90° case (opposing wind and jet), the enhanced friction velocity—and therefore enhanced *S*_{θin}—accounts for 15% of the total increase at *x* = 200 km.

It is noted that wave models that recognize only depth-averaged currents would produce smaller deviations by the factor 4, say, for a depth of 4000 m.

### d. Depth-induced wave breaking

The model is next applied to the depth-induced wave breaking experiment of Battjes and Janssen (1978), which invokes the dissipation relation, (20b), by the same authors. Waves are produced upstream by a wavemaker and propagate into the measurement region. The frequency is 0.53 Hz and the upstream boundary condition is *E*_{θ} = (.025 m^{3} s^{−2})(*β*/2) sech^{2}(*β**θ*) for *β* = 2.2. The bottom friction term (20a) is a negligible addition to *S*_{θBdis} in this case. Figure 11 shows the significant wave height distribution for *γ* = 0.70, a best-fit value close to that chosen by Booij et al. (1999) for their third-generation model comparison with the same experiment. Their result and that in Fig. 11 shows that (20b) is a remarkable relation even though it includes an adjustable constant. By examining field data, Ruessink et al. (2003) propose a dependency on *kh* such that *γ* = 0.76*kh* + 0.29.

Since the wave breaking–induced current and setup may have some influence, this calculation should be repeated when the wave model is coupled with a circulation model.

### e. Sudden change in wind direction and swell

Figure 12a demonstrates the response of directional spectra to a sudden 90° change in a 20 m s^{−1}, homogeneous wind field after a 12-h spinup from rest. One can, if one wishes, label Fig. 12a directional pseudospectrum since the frequency distribution is parameterized using (9) consistent with the calculated values of *E*_{θ}. The portion that is wind driven is delineated by *f*_{spr} > 0.0 and the dissipation constant, as before, is *b* = 0.18 × 10^{−4}; otherwise, the dissipation constant is reduced by a factor of 5, to the value *b* = 0.036 × 10^{−4}. This result may be compared with a similar case using the third-generation model of Tolman (1991) and Tolman and Chalikov (1996). Their results are understandably more diffuse, but similarly attenuate the higher-frequency energy of the swell component compared to the present result. The swell decay rate constant is uncertain; the value, 0.036 × 10^{−4}, corresponds to decay by *e*^{−1} in about 6 days and is an estimate—with help from a spatial to temporal transformation—from Snodgrass et al. (1966) for swell dissipation not far from the source. Having introduced the possibility of swell in the model, we now note that *E*_{TW} = ∫*E*_{θ}*d**θ* where the integration is limited by the region where *f*_{spr} > 0.

Figure 12b is a case wherein the wind direction changes by 180°. Compared to the eastern portion of Fig. 12a, the swell component in the eastern portion of Fig. 12b decays faster since the winds experience a negative *S*_{θin}. This is accomplished by simple addition of a negative mirror image to *f*_{spr} in (18a) but multiplied by the factor 0.4 as suggested by the laboratory experiments of Donelan (1999), which demonstrated a rapid wave decay if opposed to the wind direction. Ardhuin et al. (2007) suggest that this factor is too high compared to field data; however, this is based on data where wind and wave directions are significantly not colinear in which case *f*_{spr} < 1.

## 9. Hurricane Katrina

During the period 25–30 August 2005, Hurricane Katrina swept through the Florida Strait and into the Gulf of Mexico with a devastating landfall in and around New Orleans during the period. The path of the hurricane is shown in Fig. 13. There were six National Data Buoy Center (http://www.ndbc.noaa.gov/) buoys in the neighborhood of Katrina’s path; they are located in Fig. 13.

Yin and Oey (2007) have applied POM to the Gulf of Mexico forced by winds analyzed and obtained from the National Oceanic and Atmospheric Administration (NOAA) Hurricane Research Division (HRD) of the Atlantic Oceanographic and Meteorological Laboratory (AOML; available online at http://www.aoml.noaa.gov/hrd/). The wind speed—referenced to 10-m height—and direction are shown in Fig. 14 and are compared with the winds obtained by the buoys. We use these winds only to force the wave model. In keeping with the purpose of this paper (slightly violated in Gulf Stream–like jet case above) to establish the wave model cum sole prior to coupling with a general circulation model, current interactions are neglected. This will also enable a direct comparison with results using the SWAN model.

To accommodate large wind velocities, drag coefficients are limited to the value *C*_{D} = *C*_{Dmax} = .0025, applicable to winds over 25 m s^{−1} (and therefore not a factor in our previous calculations); this flattening of the drag coefficient curve is a relatively recent finding by Powell et al. (2003), Donelan et al. (2004), and Emmanuel (1995, 2003). Reasons advanced by these authors for the flattening include high concentration of water droplets (spume) near the surface or high curvature of wave forms such that the airflow “skips” from crest to crest and pressure in the trough is relatively constant. At three of the buoys (42040, 42001, and 42003), the model-derived *C*_{D}(*U*_{10}) from (15b) and (15c) exceeded *C*_{Dmax} at the wind speed of 25 m s^{−1}; the other three buoys (42039, 42036, and 42038) did not exceed this limit; Fig. 15 illustrates *C*_{D}(*U*_{10}) for a case wherein *C*_{D} was not limited and a case wherein the limits were invoked. It should be noted that calculated inverse wave ages (not shown) were confined to the range 0.7 < *σ*_{p}*U*_{10}/*g* < 1.8 so that most of the variation of *C*_{D}(*U*_{10}) as given by (15b) and (15c) depends on *H*_{S}, which is approximately proportional to *U*^{2}_{10}.

Figure 16 is a plan view of the model output of significant wave height, and wave propagation direction from the present model and from SWAN on the 23rd hour of 28 August 2005. The peak value for SWAN was 22 m and for the present model 15 m. The *C _{D}* limit was not applied to the SWAN calculations and, if the limit on

*C*is removed, the present peak significant wave height is also 22 m.

_{D}Figure 17 shows comparisons of significant wave height calculated by the present model, SWAN, and the buoy data. For buoy 42038, SWAN anticipates the arrival of swell better than the present model; for buoy 42001, both have overly large maximum *H _{S}*, SWAN somewhat more so than the present model. The present model seems to simulate the maximum

*H*better for buoy 42003, but this buoy capsized at the time when both models indicate the arrival of the peak. The present model somewhat underestimates the maximum

_{S}*H*for buoy 42040.

_{S}The average wave period, *T*_{ave} ≡ *E*_{T}(∫^{π}_{−π }*σ*_{θ}*E*_{θ}*d**θ*)^{−1}, is shown in Fig. 18. For buoy 42001, the error calculated by the present model is quite small; otherwise, errors on the order of 1–2 s are incurred by both models.

## 10. Summary

A surface gravity wave model has been developed; it includes depth-dependent wave–current interaction terms in both the wave energy equation and the ocean circulation equations derived in (the revised) M03. Except for section 9, the model is presented as a stand-alone model but it is designed so that coupling with three-dimensional circulation models is possible in operational or climate applications. The stand-alone wave model comprises Eqs. (8), (11), (12b), (13), (15), (18), (19), (20), and (21) and the conditions that, if *f*_{spr}(*θ*) ≥ 0.0, *b*(*θ*) = 0.18 × 10^{−4}; otherwise, *b*(*θ*) = 0.036 × 10^{−4}.

Model calculations compare favorably with established fetch and duration “laws.” Other applications were presented to demonstrate model versatility. Hurricane Katrina model–data comparisons, which contained low to high wind forcing, demonstrated significant model generality. This simple and computationally economic model produced wave properties comparable to buoy data and comparable to those obtained with the third-generation model SWAN.

The next step will be to compare a fully coupled circulation and wave model with laboratory and field data. The wave model will then include currents in (8). The circulation model’s momentum equation will include wave parameters in the wind stress relation, incorporation of depth-dependent radiation stress terms, and incorporation of wind–wave-induced pressure momentum transfer into the water column. The model’s turbulence kinetic energy equation will include a source term equal to wave dissipation. It is to be expected that wave–ocean coupling will reduce the errors discussed above, but ultimately predictive modeling of the general circulation of either fluid will require fully coupled atmosphere–wave–ocean models.

## Acknowledgments

Some of this research was conducted while GLM was a visiting professor at the Technical University of Delft; he is indebted to the faculty and staff of the Fluid Mechanics Section for valuable consultation. GLM was supported by the NOPP surf-zone project, by ONR Grant N0014-01-1-0170, and partially by DOI/MMS Contract 1435-01-03-CT-72021. GLM and MAD are presently funded by NSF OCE0526508 and OCE0526491. LYO was supported by MMS Contract 1435-01-06-CT-39731.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**.**

**,**

### APPENDIX A

#### Spectral Averages

The functions used in section 4 are defined by

The behavior of the various functions is useful to check the numerical calculations and to extend the numerical results for large and small *k _{p}D*. First, the monochromatic functions behave as follows:

The behavior of the spectrally averaged functions obtained by numerical integration is as follows:

The constants in *F*_{1}(*k _{p}D* → ∞) =

*F*

_{2}(

*k*→ ∞) are obtained from a fit to the calculations at the two topmost points, ς = −0.0 and −0.05. For monochromatic waves, we would simply have exp(2

_{p}D*kD*ς). A lookup table was prepared using numerical integrations for

*k*= 0.2, 0.4,·, ·, ·, 3.0, 100.0.

_{p}DFor *k _{p}D* ≥ 3.0, all of the functions in (5) revert to simple exponentials dependent on

*k*and

_{p}D*k*. The result is that

_{p}z*F*

_{1}(

*k*,

_{p}D*k*) =

_{p}z*F*

_{1}(100,

*k*) and

_{p}z*F*

_{3}(

*k*,

_{p}D*k*) = (

_{p}z*k*/100)

_{p}D*F*

_{3}(100,

*k*). These rules are verified by the numerical integrations for

_{p}z*k*= 3.0.

_{p}D### APPENDIX B

#### Some Numerical Details

The following code will accommodate an orthogonal curvilinear grid in *x, y* space in the sense of a finite volume. The leapfrog tendency term is split so that

where

The orthogonal cell dimensions are Δ*x* and Δ*y* reckoned at the cell center; they may vary from cell to cell according to a particular curvilinear grid. The group velocity components and and are located at the edge of each cell. After executing (B1), *ẽ*_{i,j,m} = max(*ẽ*_{i,j,m},1 × 10^{−5}); the lower limit corresponds to *H*_{S} = 1 mm.

The terms modified by the *γ* coefficient are diffusion-like terms. For *γ* = 0, the result is pure central differencing, whereas for *γ* = 0.5, the result is pure upwind. Most calculations in the main text used *γ* = 0.2. However, at cells adjacent to boundaries we stipulate local upwinding by setting *γ* = 0.5. In general, the difference in calculated results between the two values of *γ* was quite small.

After (B1), the time step in angle space is

where

Initially *D* = 0 so that (B3) and (B4) combine for a pure upwind differencing and positive definite algorithm. However, following Smolarkiewicz (1984), the solution is iterated such that an antidiffusion *D* is calculated to reduce the diffusion incurred by upwinding; it is inserted into (B3) for each iterant (three iterations are sufficient) while maintaining positive definiteness. For the details in calculating *D*, the reader is referred to the paper by Smolarkiewicz. A grid stencil, emphasizing the angle grid, is shown in Fig. B1. The angle space is −*π* to +*π*. Cyclic boundary conditions connect the branch cut at −*π*, *π*.

so that the part of *S*_{θin} + *S*_{θdis} = *A* + *Be*^{n+1}_{i,j,m} that is dependent on *e*_{i,j,m} is executed implicitly.

The discretization for (21) is

plus similar terms for the *y* direction and terms for the right-hand side of (21a).

## Footnotes

*Corresponding author address:* Prof. George Mellor, Sayre Hall, Forrestal Campus, Princeton University, Princeton, NJ 08544-0710. Email: glmellor@princeton.edu