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 cp/U10, cp/u*, or cp/Up,λ/2, where cp is the phase speed at the peak of the wave spectrum, U10 is the wind speed at 10 m above the sea surface, u* is the friction speed, and Up,λ/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 cp/U10, we will mostly follow the same practice. [However, as discussed below, cp/U10 will be converted to g/(σpU10), an equality for deep water but not for shallow water.] Alves et al. (2003) offer evidence that cp/U10 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

 
formula

and since the kinetic and potential energies are equal,

 
formula

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 37. 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,

 
formula

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, uAα, 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

 
formula
 
formula
 
formula

For the present wave model, the dispersion relation (2b) is initially solved iteratively, inverted, and cast in the form kD = f1(σ2D/g); and from (3b), n = f2(σ2D/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|/k2, where ω is obtained from (2a). Working out the vector algebra (Golding 1978) yields

 
formula

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

 
formula
 
formula
 
formula
 
formula

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

 
formula

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

 
formula

where a modified delta function, ED, is defined such that

 
formula

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

 
formula

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 ρwSθin, where ρw is the seawater density. The terms, cgα, cθ, uAα 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

 
formula
 
formula

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 α = α(Uc/cp), γ = γ(Uc/cp), σpw = σpw(Uc/cp), and Uc = U10 cos(θwθ); U10 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 Uc = U10. 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 σ4pETW/g3 as a function of U10/cp, which, for deep water, equals σpU10/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, σpU10/g ≅ constant whereas U10/cp → ∞ as kpD → 0). We define ETW 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θ, ETW = ET.

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

 
formula

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.

Fig. 1.

The inverse wave age dependence of the relation between peak frequency, total wave energy, and wind speed. The solid line is calculated from (9) whereas the dashed line is from (10) with Cσ = 0.0022.

Fig. 1.

The inverse wave age dependence of the relation between peak frequency, total wave energy, and wind speed. The solid line is calculated from (9) whereas the dashed line is from (10) with Cσ = 0.0022.

Equation (10) can also be written as

 
formula

and will be used to obtain σp. Next a small modification of (11) is possible by noting that dimensional analysis yields fcn[ETWg/U410, U10σp/g, (ETW/g)1/2/z10] = 0, where the length scale, (ETW/g)1/2 = HS/4, divided by z10 = 10 m is now added to the group of relevant nondimensional variables. Here HS = 4(ET/g)1/2 is the significant wave height (Longuet-Higgins 1952). [Note that z10 is important in the conventional velocity-dependent drag coefficient CD = CD(U10), which in dimensionless terms should be—but usually is not—written CD = CD(U10/z10g).]

The possibility that (ETW/g)1/2/z10 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(ETW/g)1/2/z10 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,

 
formula

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

 
formula

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

 
formula
 
formula
 
formula

and, similarly to (7), we have

 
formula

In Fig. 2a, cg/c and Fcθ are plotted as functions of kpD. 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 kpD is also plotted in Fig. 2a.

Fig. 2.

(a) The solid curves are spectrally weighted values of (top) cg/cp and (bottom) Fcθ as functions of kpD. The dashed curve is from (3b). The short horizontal lines are asymptotes as kpD → ∞. Wave age dependencies are small and are ignored in the model. (b) Spectrally weighted values of (left) F1 and (right) F2 for kpD = 0.2, 1.0, 2.0, and 3.0.

Fig. 2.

(a) The solid curves are spectrally weighted values of (top) cg/cp and (bottom) Fcθ as functions of kpD. The dashed curve is from (3b). The short horizontal lines are asymptotes as kpD → ∞. Wave age dependencies are small and are ignored in the model. (b) Spectrally weighted values of (left) F1 and (right) F2 for kpD = 0.2, 1.0, 2.0, and 3.0.

The Fn functions in (13b) and (13c), related to the definitions (5), are explicitly defined in appendix A and are functions of ς and kpD as shown in Fig. 2b for σpU10/g = 2. The variations with respect to inverse wave age are small (in the range, 1 < σpU10/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,

 
formula

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

 
formula

where U10U(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 ρwu2*w = ρau2*a; ρw/ρa = 860 is the water to air density ratio, (15c) is from Donelan (1990), and HS = 4(ET/g)1/2 is the significant wave height.

We next obtain Sθin, which can be written

 
formula

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

 
formula

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 U10 using the law of the wall according to U(λ/2)/U10 = ln(λ/2z0)/ ln(10 m/z0), 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.

Fig. 3.

A plot of the nondimensional DHH85 spectrum, Eq. (9), and the wave growth relation, Eq. (17), as a function of frequency for U10/cp = σpU10/g = 2.0.

Fig. 3.

A plot of the nondimensional DHH85 spectrum, Eq. (9), and the wave growth relation, Eq. (17), as a function of frequency for U10/cp = σpU10/g = 2.0.

Equation (16) was integrated to obtain the so-called spreading function fspr = Sθin/STin, where STin = ∫π/2π/2 Sθin dθ; the results are plotted in Fig. 4. The dashed line is

 
formula

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.

Fig. 4.

The directional distribution of the wind energy source as functions of relative wave propagation angle, wind speed, and wave age (the closely packed curves with no labels): U10 = (top) 10, (middle) 20, and (bottom) 30 m s−1. The dependencies on wind speed and wave age are neglected in the model. The dependence on propagation angle is simply described by (18a) for β = 2.2, which is plotted as dashed curves.

Fig. 4.

The directional distribution of the wind energy source as functions of relative wave propagation angle, wind speed, and wave age (the closely packed curves with no labels): U10 = (top) 10, (middle) 20, and (bottom) 30 m s−1. The dependencies on wind speed and wave age are neglected in the model. The dependence on propagation angle is simply described by (18a) for β = 2.2, which is plotted as dashed curves.

STin can be normalized in a number of ways—by σpET or by various combinations of cp or u*w. After normalization by cpu2*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, STin/u3*w versus σpU10/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 STin/u3*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.

Fig. 5.

The normalized wind energy source term as a function of wave age and wind speed (a) normalized on u3*w; dependence on wind speed is neglected; dashed line is from (18b) integrated over all angles and is used in the model; and (b) normalized on c3p. The ascending curves are for U10 = 10, 20, and 30 m s−1, respectively. (The corresponding nondimensional parameter is U10[g(10 m)]−1/2 ≅ 1, 2 and 3)

Fig. 5.

The normalized wind energy source term as a function of wave age and wind speed (a) normalized on u3*w; dependence on wind speed is neglected; dashed line is from (18b) integrated over all angles and is used in the model; and (b) normalized on c3p. The ascending curves are for U10 = 10, 20, and 30 m s−1, respectively. (The corresponding nondimensional parameter is U10[g(10 m)]−1/2 ≅ 1, 2 and 3)

The dashed line in Fig. 5a is given by STin/u3*w = 370 exp(−0.33σpU10/g) but the calculated values abruptly decrease around σpU10/g = 0.5 and become negative for lesser values. When STin/c3p is plotted versus σpU10/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 STin/u3*w is large only because u3*w is very small. Thus, we let

 
formula

In Fig. 5, STin/u3*w is, fortuitously, a weak function of U10, whereas STin/c3p is dependent on U10. Note that a fully developed wave field is often taken to be U10/cp = σpU10/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

 
formula

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θinbEθσ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 σpU10/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

 
formula

where ub is the wave amplitude near the bottom. This formula agrees with that in Mellor (2002), where Cb is a function of nondimensional bottom roughness. Bottom roughness is generally unknown; here we suggest a typical value, Cb = 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 Hs/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

 
formula

where the fraction of breaking waves out of an ensemble, Qb, is given by the transcendental relation (Qb − 1)/ lnQb = 8(ET/g)/H2M and where HM = γ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

 
formula

where ∂σθ/∂k = cg 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 ( fspr > 0), the additional source term, ℜ, prescribed by

 
formula

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 U10 = 0. At x = 0, incoming swell is prescribed according to Eθ ∝ (β/2) sech2[β(θθ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π/10s, which according to (21) and conservation of crests, is constant; the corresponding c(x) and cg(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 ET = Σ2πθk=1Eθ(θkθ, and is normalized with its value at x = 0. The mean angle, θ = Σ2πθk=1θEθ(θkθ/E, is shown in Fig. 6e.

Fig. 6.

A simple test of refraction wherein calculations from Eq. (8) are compared with Snell’s law for a wave period of 10 s. (a) The bottom topography. (b) The phase and group speed. (c) Contours of wave energy in propagation angle and distance space calculated from Eq. (8). (d) The angle integrated (total) wave energy per unit surface area calculated (solid line) and for constant energy flux (dashed line). (e) The mean wave propagation angle as calculated (solid line) and from Snell’s law (dashed line). The grid spacing is Δx = 0.5 km and Δθ = 15°.

Fig. 6.

A simple test of refraction wherein calculations from Eq. (8) are compared with Snell’s law for a wave period of 10 s. (a) The bottom topography. (b) The phase and group speed. (c) Contours of wave energy in propagation angle and distance space calculated from Eq. (8). (d) The angle integrated (total) wave energy per unit surface area calculated (solid line) and for constant energy flux (dashed line). (e) The mean wave propagation angle as calculated (solid line) and from Snell’s law (dashed line). The grid spacing is Δx = 0.5 km and Δθ = 15°.

Calculations are obtained from Snell’s law so that

 
formula
 
formula

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

 
formula

and that is plotted as a dashed line in Fig. 7. The well-developed limit of ETg/U410 = 3.6 × 10−3 is from Pierson and Moskowitz (1964).

Fig. 7.

The calculated fetch-dependent wave energy vs nondimensional fetch distance. Calculations for wind speeds, U10 = 10 and 20 m s−1 (solid lines), are compared with Eq. (23) (dashed line) and a result from Donelan et al. (1992) (dotted–dashed line). The later two curves represent a synthesis of observational data for a range of wind speeds. The wind angle is normal to the coast located at x = 0.

Fig. 7.

The calculated fetch-dependent wave energy vs nondimensional fetch distance. Calculations for wind speeds, U10 = 10 and 20 m s−1 (solid lines), are compared with Eq. (23) (dashed line) and a result from Donelan et al. (1992) (dotted–dashed line). The later two curves represent a synthesis of observational data for a range of wind speeds. The wind angle is normal to the coast located at x = 0.

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 U10 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θ ( fspr < 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 U10, there must be some model dependence on U10 or, in dimensionless form, U10[g (10 m)]−1/2 according to dimensional analysis or from (15a)(15c), (16), and (17) together with the relation HS = 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 σpU10/g versus xg/U210 was also quite favorable.

Fig. 8.

The calculated fetch-dependent wave energy vs nondimensional fetch distance and duration for U10 = 10 m s−1. For U10 = 20 m s−1, the calculated values would be offset vertically as in Fig. 7. The “O”s are fetch-limited estimates from Hwang and Wang (2004, their Fig. 5). The “X”s are duration-limited estimates from the same source. See their paper for data scatter, which is greatest for large values of nondimensional x and t.

Fig. 8.

The calculated fetch-dependent wave energy vs nondimensional fetch distance and duration for U10 = 10 m s−1. For U10 = 20 m s−1, the calculated values would be offset vertically as in Fig. 7. The “O”s are fetch-limited estimates from Hwang and Wang (2004, their Fig. 5). The “X”s are duration-limited estimates from the same source. See their paper for data scatter, which is greatest for large values of nondimensional x and t.

Figures 9a,b demonstrate the fetch-limited behavior for U10 = 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 fspr, 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

 
formula
Fig. 9.

(a) The same as in Fig. 7, but for a wind speed of U10 = 10 m s−1 and for different mean wind angles relative to the normal to the coast at x = 0. (b) A plot of the mean wave propagation angle, defined by Eq. (24), vs fetch.

Fig. 9.

(a) The same as in Fig. 7, but for a wind speed of U10 = 10 m s−1 and for different mean wind angles relative to the normal to the coast at x = 0. (b) A plot of the mean wave propagation angle, defined by Eq. (24), vs fetch.

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

 
formula

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 U10 = 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.

Fig. 10.

The influence of a northward Gulf Stream–like jet (maximum velocity = 2 m s−1) on waves forced by 10 m s−1 winds whose angle varies from 90° (northward; wind and jet velocities in the same direction) to −90° (southward; wind and jet velocities in opposite directions). (a) The variation of wave energy. The dashed line is the background fully developed wave field in the absence of the jet. (b) The variation of the mean propagation angle for the same wind angles as in (a).

Fig. 10.

The influence of a northward Gulf Stream–like jet (maximum velocity = 2 m s−1) on waves forced by 10 m s−1 winds whose angle varies from 90° (northward; wind and jet velocities in the same direction) to −90° (southward; wind and jet velocities in opposite directions). (a) The variation of wave energy. The dashed line is the background fully developed wave field in the absence of the jet. (b) The variation of the mean propagation angle for the same wind angles as in (a).

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θuAα/∂x = sin2θ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 m3 s−2)(β/2) sech2(βθ) 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.76kh + 0.29.

Fig. 11.

Utilization of the depth-induced wave breaking formula, Eq. (20b) of Battjes and Janssen (1978), in (a) a calculation of wave breaking over (b) a shoaling bar. Data from the same authors.

Fig. 11.

Utilization of the depth-induced wave breaking formula, Eq. (20b) of Battjes and Janssen (1978), in (a) a calculation of wave breaking over (b) a shoaling bar. Data from the same authors.

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 fspr > 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 ETW = ∫Eθdθ where the integration is limited by the region where fspr > 0.

Fig. 12.

Directional pseudospectra for a sudden change of wind direction for times of 12–15 h. The contours progressively increase by the factor 2. The dashed outer circle is for a maximum frequency fp = (2π)−1σp = 0.4 Hz; the contour interval is 0.1 Hz. (a) A homogeneous steady wind of 20 m s−1 is applied for 12 h starting from rest after which the direction is changed by 90°. The significant wave heights are successively 5.16, 5.45, 5.76, and 6.02 m. (b) The wind direction is changed by 180°. The significant wave heights are successively 5.16, 5.23, 5.31, and 5.39 m.

Fig. 12.

Directional pseudospectra for a sudden change of wind direction for times of 12–15 h. The contours progressively increase by the factor 2. The dashed outer circle is for a maximum frequency fp = (2π)−1σp = 0.4 Hz; the contour interval is 0.1 Hz. (a) A homogeneous steady wind of 20 m s−1 is applied for 12 h starting from rest after which the direction is changed by 90°. The significant wave heights are successively 5.16, 5.45, 5.76, and 6.02 m. (b) The wind direction is changed by 180°. The significant wave heights are successively 5.16, 5.23, 5.31, and 5.39 m.

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 fspr 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 fspr < 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.

Fig. 13.

The Gulf of Mexico. The path of Hurricane Katrina and the location of six NDBC buoys denoted by solid squares. The arrows are the 10-m wind vectors at 12 noon on 29 Aug 2005.

Fig. 13.

The Gulf of Mexico. The path of Hurricane Katrina and the location of six NDBC buoys denoted by solid squares. The arrows are the 10-m wind vectors at 12 noon on 29 Aug 2005.

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.

Fig. 14.

The analyzed (blue lines) wind speeds and directions (measured from true north; e.g., 90° for wind toward the east) from the NOAA HRD of AOML compared to the data (red lines) obtained by the six NDBC buoys. All winds are referenced to the 10-m height.

Fig. 14.

The analyzed (blue lines) wind speeds and directions (measured from true north; e.g., 90° for wind toward the east) from the NOAA HRD of AOML compared to the data (red lines) obtained by the six NDBC buoys. All winds are referenced to the 10-m height.

To accommodate large wind velocities, drag coefficients are limited to the value CD = CDmax = .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 CD(U10) from (15b) and (15c) exceeded CDmax 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 CD(U10) for a case wherein CD 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 < σpU10/g < 1.8 so that most of the variation of CD(U10) as given by (15b) and (15c) depends on HS, which is approximately proportional to U210.

Fig. 15.

Model-calculated CD as a function of the 10-m wind speed for two of the buoys: (top) 42039 for which CD was not limited to the value 0.0025, and (bottom) 42040 for which wind speeds were greater and CD was limited. The straight line, after Garratt (1977), is CD = (0.75 + 0.067U10)10−3.

Fig. 15.

Model-calculated CD as a function of the 10-m wind speed for two of the buoys: (top) 42039 for which CD was not limited to the value 0.0025, and (bottom) 42040 for which wind speeds were greater and CD was limited. The straight line, after Garratt (1977), is CD = (0.75 + 0.067U10)10−3.

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 CD limit was not applied to the SWAN calculations and, if the limit on CD is removed, the present peak significant wave height is also 22 m.

Fig. 16.

A comparison of wave direction given by arrows and significant wave height (contours) calculated by the (a) present model and (b) SWAN on the 23rd hour of 28 Aug 2005. If the limit on CD is removed from the present model, the peak significant wave heights of both models are in agreement.

Fig. 16.

A comparison of wave direction given by arrows and significant wave height (contours) calculated by the (a) present model and (b) SWAN on the 23rd hour of 28 Aug 2005. If the limit on CD is removed from the present model, the peak significant wave heights of both models are in agreement.

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 HS, SWAN somewhat more so than the present model. The present model seems to simulate the maximum HS 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 HS for buoy 42040.

Fig. 17.

For the six buoys, significant wave heights from the buoy data (red lines) and as calculated by SWAN (blue lines) and the present model (green lines).

Fig. 17.

For the six buoys, significant wave heights from the buoy data (red lines) and as calculated by SWAN (blue lines) and the present model (green lines).

The average wave period, TaveET(∫ππ σθ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.

Fig. 18.

For the six buoys, average wave periods from the buoy data (red lines) and as calculated by SWAN (blue lines) and the present wave model (green lines).

Fig. 18.

For the six buoys, average wave periods from the buoy data (red lines) and as calculated by SWAN (blue lines) and the present wave model (green lines).

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 fspr(θ) ≥ 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.

Fig. B1. The computational stencil emphasizing the propagation angle grid and the cyclic boundary conditions. Wave energy, ei,j,m, is located at the center of each cell and the speeds, cgx,i,j,m, cgy,i,j,m, and cθ,i,j,m, are located at the edge of each cell.

Fig. B1. The computational stencil emphasizing the propagation angle grid and the cyclic boundary conditions. Wave energy, ei,j,m, is located at the center of each cell and the speeds, cgx,i,j,m, cgy,i,j,m, and cθ,i,j,m, are located at the edge of each cell.

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

REFERENCES
Alves
,
J. H. G. M.
,
M. L.
Banner
, and
I. R.
Young
,
2003
:
Revisiting the Pierson–Moskowitz asymptotic limits for fully developed wind waves.
J. Phys. Oceanogr.
,
33
,
1301
1323
.
Ardhuin
,
F.
,
T. H. C.
Herbers
,
G. P.
van Vledder
,
K. P.
Watts
,
R.
Jensen
, and
H. C.
Graber
,
2007
:
Swell and slanting fetch effects on wind growth.
J. Phys. Oceanogr.
,
37
,
908
931
.
Banner
,
M. L.
,
1990
:
Equilibrium spectra of wind waves.
J. Phys. Oceanogr.
,
20
,
966
984
.
Battjes
,
J. A.
, and
J. P. F. M.
Janssen
,
1978
:
Energy loss and set-up due to breaking of random waves. Proc. 16th Int. Conf. of Coastal Engineering, New York, NY, American Society of Civil Engineers, 569–587
.
Booij
,
N.
,
R. C.
Ris
, and
L. H.
Holthuijsen
,
1999
:
A third-generation wave model for coastal regions. 1. Model description and validation.
J. Geophys. Res.
,
104
,
7649
7666
.
Craig
,
P. D.
, and
M. L.
Banner
,
1994
:
Modeling wave-enhanced turbulence in the ocean surface layer.
J. Phys. Oceanogr.
,
24
,
2546
2559
.
Donelan
,
M.
,
M.
Skafel
,
H.
Graber
,
P.
Liu
,
D.
Schwab
, and
S.
Venkates
,
1992
:
On the growth of wind-generated waves.
Atmos.–Ocean
,
30
,
457
478
.
Donelan
,
M. A.
,
1977
:
A simple numerical model for wave and wind stress prediction. Nat. Water Res. Inst. Rep., Burlington, ON, Canada, 27 pp
.
Donelan
,
M. A.
,
1990
:
Air–sea interaction.
The Sea, B. LeMehaute and D. M. Hanes, Eds., Ocean Engineering Science, Vol. 9, John Wiley and Sons, 239–292
.
Donelan
,
M. A.
,
1999
:
Wind-induced growth and attenuation of laboratory waves. Wind-over-Wave Couplings, S. G. Sajjadi, N. H. Thomas, and J. C. R. Hunt, Eds., Clarendon Press, 183–194
.
Donelan
,
M. A.
,
J.
Hamilton
, and
W. H.
Hui
,
1985
:
Directional spectra of wind-generated waves.
Philos. Trans. Roy. Soc. London
,
A315
,
509
562
.
Donelan
,
M. A.
,
B. K.
Haus
,
N.
Reul
,
W. J.
Plant
,
M.
Stiassnie
,
H. C.
Graber
,
O. B.
Brown
, and
E. S.
Saltzman
,
2004
:
On the limiting aerodynamic roughness in very strong winds.
Geophys. Res. Lett.
,
31
.
L18206, doi:10.1029/2004GL019460
.
Emmanuel
,
K.
,
1995
:
Sensitivity of tropical cyclones to surface exchange coefficients and a revised steady state model incorporating eye dynamics.
J. Atmos. Sci.
,
52
,
3969
3976
.
Emmanuel
,
K.
,
2003
:
A similarity hypothesis for air–sea exchange of extreme wind speeds.
J. Atmos. Sci.
,
60
,
1420
1428
.
Garratt
,
J. R.
,
1977
:
Review of drag coefficients over oceans and continents.
Mon. Wea. Rev.
,
105
,
915
929
.
Golding
,
B. W.
,
1978
:
A depth dependent wave model for operational forecasting.
Turbulent Fluxes through the Sea Surface, Wave Dynamics and Prediction, A. Favre and K. Hasselmann, Eds., Plenum Press, 593–604
.
Hasselmann
,
K.
, and
Coauthors
,
1973
:
Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP). Dtsch. Hydrogr. Z., A8 (Suppl.), 1–95
.
Holthuijsen
,
L. H.
,
N.
Booij
, and
T. H. C.
Herbers
,
1989
:
A prediction model for stationary, short-crested waves in shallow water with ambient currents.
Coastal Eng.
,
13
,
23
54
.
Hwang
,
P. A.
, and
D. W.
Wang
,
2004
:
Field measurements of duration-limited growth of wind-generated ocean surface waves at young stage of development.
J. Phys. Oceanogr.
,
34
,
2316
2326
.
Janssen
,
P. A. E. M.
,
1989
:
Wave induced stress and drag of air over sea water.
J. Phys. Oceanogr.
,
19
,
1631
1642
.
Janssen
,
P. A. E. M.
,
O.
Saetra
,
C.
Wettre
,
H.
Hersbach
, and
J.
Bidlot
,
2004
:
Impact of the sea state on the atmosphere and ocean.
Ann. Hydrogr.
,
3
,
3.1
3.23
.
Kitaigorodskii
,
S. A.
,
1962
:
Application of the theory of similarity to the analysis of wind-generated water waves as a stochastic process.
Bull. Acad. Sci. USSR Ser., 1, 73 pp
.
Komen
,
G. J.
,
L.
Cavaleri
,
M.
Donelan
,
K.
Hasselmann
,
S.
Hasselmann
, and
P. A. E. M.
Janssen
,
1994
:
Dynamics and Modelling of Ocean Waves.
Cambridge University Press, 532 pp
.
Lin
,
W.
,
L. P.
Sanford
, and
S. E.
Suttles
,
2002
:
Wave measurements and modeling in Chesepeake Bay.
Cont. Shelf Res.
,
22
,
2673
2686
.
Longuet-Higgins
,
M. S.
,
1952
:
On the statistical distribution of the heights of sea waves.
J. Mar. Res.
,
11
,
245
266
.
Longuet-Higgins
,
M. S.
,
1953
:
Mass transport in water waves.
Phil. Trans. Roy. Soc.
,
A245
,
533
581
.
Mastenbroek
,
C.
,
G.
Burgers
, and
P. A. E. M.
Janssen
,
1993
:
The dynamical coupling of a wave model and a storm surge model through the atmospheric boundary layer.
J. Phys. Oceanogr.
,
23
,
1856
1866
.
Mellor
,
G. L.
,
2002
:
The oscillatory bottom boundary layer.
J. Phys. Oceanogr.
,
32
,
3075
3088
.
Mellor
,
G. L.
,
2003
:
The three-dimensional, current, and surface wave equations.
J. Phys. Oceanogr.
,
33
,
1978
1989
.
Mellor
,
G. L.
,
2005
:
Some consequences of the three-dimensional current and surface wave equations.
J. Phys. Oceanogr.
,
35
,
2291
2298
.
Mellor
,
G. L.
,
2005
:
Some consequences of the three-dimensional current and surface wave equations.
J. Phys. Oceanogr.
,
35
,
2291
2298
.
Mellor
,
G. L.
,
2008
:
The depth-dependent current and wave interaction equations: A revision.
J. Phys. Oceanogr.
,
in press
.
Mellor
,
G. L.
, and
A.
Blumberg
,
2004
:
Wave breaking and ocean surface thermal response.
J. Phys. Oceanogr.
,
34
,
693
698
.
Moon
,
I-J.
,
2005
:
Impact of a coupled ocean, wave, tide circulation system on coastal modeling.
Ocean Modell.
,
8
,
203
236
.
Ozer
,
J.
,
R.
Padilla-Hernandez
,
J.
Monbaliu
,
E. A.
Fanjul
,
J. C. C.
Albiach
,
P.
Osuna
,
J. C. S.
Yu
, and
J.
Wolf
,
2000
:
A coupling module for tides, surges and waves.
Coastal Eng.
,
41
,
95
124
.
Phillips
,
O. M.
,
1977
:
The Dynamics of the Upper Ocean.
Cambridge University Press, 336 pp
.
Pierson
,
W. J.
, and
L.
Moskowitz
,
1964
:
A proposed spectral form for fully developed wind seas based on the similarity theory of A. A. Kitaigoradskii.
J. Geophys. Res.
,
69
,
5181
5190
.
Powell
,
M. D.
,
P. J.
Vickery
, and
T. A.
Reinhold
,
2003
:
Reduced drag coefficient for high wind speeds in tropical cyclones.
Nature
,
422
,
279
283
.
Ruessink
,
B. G.
,
D. J. R.
Walstra
, and
H. N.
Southgate
,
2003
:
Calibration and verification of a parametric wave model on barred beaches.
Coastal Eng.
,
48
,
139
149
.
Schwab
,
D. J.
,
J. R.
Bennett
,
P. C.
Liu
, and
M. A.
Donelan
,
1984
:
Application of a simple numerical wave prediction model to Lake Erie.
J. Geophys. Res.
,
89
,
3586
3592
.
Smolarkiewicz
,
P. K.
,
1984
:
A fully multidimensional positive definite advection transport algorithm with small implicit diffusion.
J. Comput. Phys.
,
54
,
325
362
.
Snodgrass
,
F. E.
,
G. W.
Groves
,
K.
Hasselmann
,
G. R.
Miller
,
W. H.
Munk
, and
W. H.
Powers
,
1966
:
Propagation of swell across the Pacific.
Philos. Trans. Roy. Soc. London
,
259
,
431
497
.
Stacey
,
M. W.
,
1999
:
Simulations of the wind-forced near-surface circulation in Knight Inlet: A parameterization of the roughness length.
J. Phys. Oceanogr.
,
29
,
1363
1367
.
SWAMP Group
,
1985
:
Ocean Wave Modeling.
Plenum Press, 256 pp
.
Terray
,
E. A.
,
M. A.
Donelan
,
Y. C.
Agrawal
,
W. M.
Drennan
,
K. K.
Kahma
,
A. J. W.
Williams
III
,
P. A.
Hwang
, and
S. A.
Kitaigorodski
,
1996
:
Estimates of kinetic energy dissipation under breaking waves.
J. Phys. Oceanogr.
,
26
,
792
807
.
Tolman
,
H. L.
,
1991
:
A third-generation model for wind waves on slowly varying, unsteady, and inhomogeneous depths and currents.
J. Phys. Oceanogr.
,
21
,
782
797
.
Tolman
,
H. L.
, and
D.
Chalikov
,
1996
:
Source terms for a third-generation wind wave model.
J. Phys. Oceanogr.
,
26
,
2497
2518
.
WAMDI Group
,
1988
:
The WAM model—A third generation ocean wave prediction model.
J. Phys. Oceanogr.
,
18
,
1775
1810
.
Xie
,
L.
,
K.
Wu
,
L.
Pietrafesa
, and
C.
Zhang
,
2001
:
A numerical model of wave-current interaction through surface and bottom stresses: Wind-driven circulation in the South Atlantic Bight under uniform winds.
J. Geophys. Res.
,
106
,
C8
.
16841
16855
.
Yin
,
X-Q.
, and
L-Y.
Oey
,
2007
:
Bred-ensemble ocean forecast of Loop Current and Rings.
Ocean Modell.
,
17
,
300
326
.
Young
,
I. R.
,
2006
:
Directional spectra of hurricane wind waves.
J. Geophys. Res.
,
111
.
C08020, doi:10.1029/2006JC003540
.
Zhang
,
M. Y.
, and
Y. S.
Li
,
1997
:
The dynamic coupling of a third-generation wave model and a 3D hydrodynamic model through boundary layers.
Cont. Shelf Res.
,
17
,
1141
1170
.

APPENDIX A

Spectral Averages

The functions used in section 4 are defined by

 
formula
 
formula
 
formula
 
formula

The behavior of the various functions is useful to check the numerical calculations and to extend the numerical results for large and small kpD. First, the monochromatic functions behave as follows:

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

The constants in F1(kpD → ∞) = F2(kpD → ∞) 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(2kDς). A lookup table was prepared using numerical integrations for kpD = 0.2, 0.4,·, ·, ·, 3.0, 100.0.

For kpD ≥ 3.0, all of the functions in (5) revert to simple exponentials dependent on kpD and kpz. The result is that F1(kpD, kpz) = F1(100, kpz) and F3(kpD, kpz) = (kpD/100)F3(100, kpz). These rules are verified by the numerical integrations for kpD = 3.0.

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

 
formula

where

 
formula
 
formula

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 HS = 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

 
formula

where

 
formula

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 −π, π.

Finally the source terms, (18) and (19), are included such that

 
formula

so that the part of Sθin + Sθdis = A + Ben+1i,j,m that is dependent on ei,j,m is executed implicitly.

The discretization for (21) is

 
formula

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