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.
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
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, 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
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.
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, ED, 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 ρ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
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.
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[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,
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
and, similarly to (7), we have
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.
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,
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 U10 − 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 ρ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
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.
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.
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.
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
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
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 σ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
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).
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
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
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.
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
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.
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
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 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.
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.
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.
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.
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 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.
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.
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.
The average wave period, Tave ≡ ET(∫π−π σθ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.
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.
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.
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 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.
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
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
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 + Ben+1i,j,m that is dependent on ei,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).
Corresponding author address: Prof. George Mellor, Sayre Hall, Forrestal Campus, Princeton University, Princeton, NJ 08544-0710. Email: email@example.com