Ocean Modelling 28 (2009) 153–166

Contents lists available at ScienceDirect

Ocean Modelling
journal homepage: www.elsevier.com/locate/ocemod

An unstructured-grid finite-volume surface wave model (FVCOM-SWAVE):

Implementation, validations and applications
Jianhua Qi a, Changsheng Chen a,d,e,*, Robert C. Beardsley b, Will Perrie c, Geoffrey W. Cowles a, Zhigang Lai a
a
Department of Fisheries Oceanography, School for Marine Science and Technology, University of Massachusetts-Dartmouth,
706 South Rodney French Blvd, New Bedford, MA 02744, USA
b
Department of Physical Oceanography, Woods Hole Oceanographic Institution, Woods Hole, MA 02543, USA
c
Fisheries & Oceans Canada, Bedford Institute of Oceanography, Dartmouth, Nova Scotia, Canada
d
Marine Ecosystem and Environmental Laboratory, Shanghai Ocean University, 334 Jungong Road, 200090 Shanghai, PR China
e
The State Key Laboratory for Estuarine and Coastal Research, East China Normal University, Shanghai, PR China

a r t i c l e i n f o a b s t r a c t

Article history: The structured-grid surface wave model SWAN (Simulating Waves Nearshore) has been converted into an
Received 23 April 2008 unstructured-grid finite-volume version (hereafter referred to as FVCOM-SWAVE) for use in coastal ocean
Received in revised form 25 December 2008 regions with complex irregular geometry. The implementation is made using the Flux-Corrected Trans-
Accepted 13 January 2009
port (FCT) algorithm in frequency space, the implicit Crank–Nicolson method in directional space and
Available online 5 February 2009
options of explicit or implicit second-order upwind finite-volume schemes in geographic space.
FVCOM-SWAVE is validated using four idealized benchmark test problems with emphasis on numerical
Keywords:
dispersion, wave-current interactions, wave propagation over a varying-bathymetry shallow water
Unstructured grid
Finite-volume
region, and the basic wave grow curves. Results demonstrate that in the rectangular geometric domain,
Wave model the second-order finite-volume method used in FVCOM-SWAVE has the same accuracy as the third-order
FVCOM-SWAVE finite-difference method used in SWAN. FVCOM-SWAVE was then applied to simulate wind-induced sur-
face waves on the US northeast shelf with a central focus in the Gulf of Maine and New England Shelf.
Through improved geometric fitting of the complex irregular coastline, FVCOM-SWAVE was able to
robustly capture the spatial and temporal variation of surface waves in both deep and shallow regions
along the US northeast coast.
Published by Elsevier Ltd.

1. Introduction Hsu et al. (2005) recently converted SWAN into a finite-ele-

ment-based unstructured-grid version (hereafter referred to as
SWAN (Simulating WAves Nearshore) is the third-generation FE-WAVE) for use in irregular coastal settings characterized by
surface wave model developed originally by Booij et al. (1999) numerous barriers, islands, inlets and narrow navigation channels.
and improved through a team effort (SWAN Team, 2006a). This Discretizing the wave action balance equation using a non-over-
model considers the characteristics of surface waves in shallow lapping triangular mesh, the Taylor-Galerkin finite-element meth-
water by solving the wave action balance equation with inclusion od provides an accurate geometric fitting of complex coastlines,
of dissipation from bottom friction, triad and quadruplet wave– which makes FE-WAVE more suitable for nearshore applications.
wave interactions, and shallow water wave-breaking (SWAN Team, Model-model comparisons of simulated waves generated by Ty-
2006b). SWAN has become one of the most popular surface wave phoon Bilis in 2000 in the coastal region of Taiwan Strait clearly
models presently available and it is widely used for coastal ocean show that with better resolution of the complex coastal geometry,
wave simulations, engineering applications and surface wave fore- FE-WAVE provides a more realistic and accurate simulation of sig-
casts. SWAN is discretized using a curvilinear-structured grid and nificant wave heights and dominant wave periods than is possible
solved using fully implicit finite-difference algorithms. By applica- using the standard SWAN formulation (Hsu et al., 2005).
tion of a coarse-fine grid nesting approach, SWAN can be set up As an alternative unstructured-grid algorithm, the finite-vol-
with variable grids in deep and shallow ocean regions to provide ume method has recently received more attention in the coastal
high quality simulations of surface waves in the nearshore region. ocean modeling community (Casulli and Lang, 2004; Chen et al.,
2003; Fringer et al., 2006). Dividing the computational domain
* Corresponding author. Address: Department of Fisheries Oceanography, School by using a triangular mesh and solving the equations with flux-
for Marine Science and Technology, University of Massachusetts-Dartmouth, 706
based discrete algorithms, this method takes advantage of finite-
South Rodney French Blvd, New Bedford, MA 02744, USA. Tel.: +508 910 6388;
fax: +508 910 6371.
difference methods for simple code structure and computational
E-mail address: c1chen@umassd.edu (C. Chen). efficiency and finite-element methods for geometric flexibility.

1463-5003/$ - see front matter Published by Elsevier Ltd.

doi:10.1016/j.ocemod.2009.01.007
154 J. Qi et al. / Ocean Modelling 28 (2009) 153–166

FVCOM (an unstructured-grid finite-volume coastal ocean model) where n denotes the nth time step, and Dt is the time interval for
is a state-of-the-art finite-volume coastal ocean model that has the numerical integration. Eqs. (3) and (4) describe the change of
been widely used in coastal and estuarine regions (Chen et al., action density spectrum in spectral space. They are solved by the
2006a; Chen et al., 2006b; Chen et al., 2007). An integrated coastal Flux Corrected Transport method (FCT) (Boris and Book, 1973;
ocean model system has been developed around FVCOM for the Hsu et al., 2005) and the Crank–Nicolson method (Crank and Nicol-
purpose of coastal environmental prediction and management son, 1947), respectively. Eq. (5) describes the propagation of the
(described at http://fvcom.smast.umassd.edu/). Implemention of waves in geographic space. It is solved by either an explicit finite-
an unstructured-grid surface wave model within this system volume upwind advection scheme (directly adopted from FVCOM)
makes it more suitable and reliable for application to inundation or a semi-implicit finite-volume upwind advection scheme. Eq. (6)
simulations (e.g., flooding due to storm surge) and studies of fish- represents the growth, transfer and decay of the waves driven by
ery larval recruitment problems, in which the surface waves are the source terms. It is solved by a semi-implicit integration scheme
directly related to sediment resuspension at the ocean bottom. as used in the WAM model (WAMDI Group, 1988) and WAVE-
By implementing finite-volume algorithms within SWAN, we WATCH III model (Tolman, 2002). A brief description of the discrete
have converted SWAN into an unstructured-grid finite-volume algorithms used to solve Eqs. (3)–(6) is given below.
model (hereafter referred to as FVCOM-SWAVE). This model pro-
vides an alternative option for unstructured-grid wave models for 2.2.1. Action density in frequency space
the coastal ocean. FVCOM-SWAVE can also be coupled with any tri- The FCT method, proposed first by Boris and Book (1973), is a
angular mesh-based unstructured-grid ocean models for the study conservative, positive discrete algorithm suitable for steep-gradi-
and simulation of current-wave interactions. This paper describes ent problems without dispersively generated oscillations. The algo-
the finite-volume algorithms used in FVCOM-SWAVE, followed by rithm was used by Hsu et al. (2005) in FE-WAVE and also is
a series of validation experiments and an example showing the adopted here to solve Eq. (3) in FVCOM-SWAVE. The discrete ap-
application of the model to the US northeast coastal ocean. proach used in the FCT method consists of transport, anti-diffusion
and correcting stages, as given by
2. Governing equation and discrete methods
Njnþ1=4
r
¼ Njr  ðAjr þ1=2  Ajr 1=2 Þ; ð7Þ
2.1. Spectral action balance equation where jr denotes the jth frequency and the resolution is specified as
a ‘‘constant-relative-frequency” defined as D r/r. N jr represents the
The evolution of wave spectra is determined by the wave action action density at the transport stage calculated by
density spectrum balance equation expressed as
Dt 1
@N @C r N @C h N Stot Njr ¼ Nnjr  ðU  U1jr 1=2 Þ ð8Þ
þ r  ½ð~
Cg þ ~
VÞN þ þ ¼ ð1Þ Dr jr þ1=2
@t @r @h r where U is the flux defined as
where N is the wave action density spectrum; t is the time; r is the
C r;jr þ1 þ jC r;jr þ1 j C r;jr þ1  jC r;jr þ1 j
relative frequency; h is the wave direction; Crand Ch are the wave U1jr þ1=2 ¼ Nnjr þ Nnjr þ1 ð9Þ
propagation velocities in spectral space (r, h); ~ C g ¼ @ r=@~
k is the 2 2
*
~ C r;jr þ jC r;jr j n C r;jr  jC r;jr j
group velocity; k is the wave number vector; V is ambient water U1jr 1=2 n
¼ Njr 1 þ N jr ð10Þ
current vector, and r  () is the horizontal divergence operator in 2 2
geographic space. In the Cartesian coordinates, r  () = @()/@ and the superscript ‘‘1” denotes the first stage.
x + @()/@y, while in spherical coordinates, we denote k as longitude Ajr þ1=2 and Ajr 1=2 are limited anti-diffusion fluxes defined as
and u as latitude, implying r  () = @()/@k + cos1u @[cosu()]/@u. Stot n h
is the source-sink term given as Ajr þ1=2 ¼ sgnðAjr þ1=2 Þmax 0;min jAjr þ1=2 j;sgnðAjr þ1=2 ÞðN jr þ2  Njr þ1 Þ;
io
Stot ¼ Sin þ Snl3 þ Snl4 þ Sds;w þ Sds;b þ Sds;br ð2Þ sgnðAjr þ1=2 ÞðN jr  Njr 1 Þ ð11Þ
n h
where Sin is the function for the wind-induced wave growth; Snl3 is
Ajr 1=2 ¼ sgnðAjr 1=2 Þmax 0;min jAjr 1=2 j;sgnðAjr 1=2 ÞðN jr þ1  Njr Þ;
the nonlinear transfer of wave energy due to three-wave interac- io
tions; Snl4 is the nonlinear transfer of wave energy due to four-wave sgnðAjr 1=2 ÞðN jr 1  N jr 2 Þ ð12Þ
interactions; Sds,w is the wave decay due to white capping; Sds,b is
the wave decay due to bottom friction; and Sds,br is the wave decay where
due to depth-induced wave breaking. Detailed descriptions of each Dt 2 Dt 2
of these terms are given in the SWAN technical manual (SWAN Ajr þ1=2 ¼ ðU  U1jr þ1=2 Þ; Ajr 1=2 ¼ ðU  U1jr 1=2 Þ; ð13Þ
Dr jr þ1=2 Dr jr 1=2
Team, 2006b) and not included here.
C r;j þ1 þ C r;jr C r;jr 1 þ C r;jr
U2jr þ1=2 ¼ Nnjr þ1 r ; U2jr 1=2 ¼ N njr ; ð14Þ
2 2
2.2. Discretization
and superscript ‘‘2” denotes the second stage and sgnðAjr þ1=2 Þ ¼

Following the discrete approach used in FE-WAVE (Hsu et al., 1; if Ajr þ1=2 P 0
.
2005), we split Eq. (1) into four equations given as 1; if Ajr þ1=2 < 0
1
N nþ4  Nn @ðC r NÞ 2.2.2. Action density in directional space
þ ¼0 ð3Þ The action density at the (n + 2/4)th time step in wave direc-
Dt @r
2 1 tional space is calculated using a second-order accurate implicit
N nþ4  Nnþ4 @ðC h NÞ
þ ¼0 ð4Þ Crank–Nicolson scheme (Crank and Nicolson, 1947) given by
Dt @h
Dt h i
3 2
N nþ4  Nnþ4
þ r  ½ð~
Cg þ ~
VÞN ¼ 0 ð5Þ Njnþ2=4 ¼ N nþ1=4
j þ a ðC h NÞnþ2=4
j  ðC h NÞ nþ2=4
j  ð1  aÞ
Dt h h 2 Dh h 1 h þ1
3 h i
N nþ1  Nnþ4 Stot Dt nþ1=4 nþ1=4
¼ ð6Þ  ðC h NÞjh þ1  ðC h NÞjh 1 ð15Þ
Dt r 2 Dh
 J. Qi et al. / Ocean Modelling 28 (2009) 153–166 155

where jh denotes the jth direction interval and the resolution is Substituting Eq. (18) into Eq. (16) results in a 2D asymmetric and
specified by Dh. a is a weighting factor with a default value of 0.5. diagonally dominant matrix with a stencil equal to the sum of the
surrounding node points contained in a control volume. It can be
2.2.3. Action density in geographic space solved efficiently using a scalable sparse matrix solver library
Eq. (5) is solved numerically using the unstructured-grid finite- (PETSc) (Balay et al., 2007) implemented with a high performance
volume approach implemented in FVCOM (Chen et al., 2003; Chen pre-conditional HYPRE software library (HYPRE Team, 2001). This
et al., 2006b) for both Cartesian and spherical coordinates. The flux method of sparse matrix solution was implemented to solve the
form of Eq. (5) in a control volume shown in Fig. 1 can be written as Poisson equation for the non-hydrostatic version of FVCOM (Lai
et al., 2008).
Dt X
ln
Nnþ3=4 ¼ Nnþ2=4  C n;i Nli Dli : ð16Þ
X i¼1 2.2.4. Action density (Nn+1) related to source terms
Eq. (6) is solved using the same second-order, semi-implicit,
Here, X is the area of the control volume indicated by the dark
centered-difference scheme as is implemented in WAM and
shaded area in Fig. 1, li (i = 1,ln) is the perimeter of X, ln is the number
WAVEWATCH-III. This is
of edges of X and Cn,i is the component of ~ Cg þ ~
V normal to li. Eq. (16)
is solved by either the explicit upwind scheme or semi-implicit up- Dt nþ1
Nnþ1 ¼ Nnþ3=4 þ ðS þ Sn Þ ð19Þ
wind scheme. Detailed descriptions of these two solvers are given 2r
in the Appendix; the calculation of the action density at the edge of where Sn+1 and Nn+1 are nonlinearly coupled to each other. A de-
the control volume is summarized here. In the explicit approach, tailed description of this algorithm is given by the WAMDI Group
8
> lX ;n (1988) and Tolman (2002).
> Dr A P
A
>
> Nnþ2=4 þ
nþ2=4
NlX ;j DlXA ;j for C n > 0 FVCOM-SWAVE is parallelized using the MPI framework as
>
< A XA
j¼1 A
Nl ¼ ; ð17Þ implemented in FVCOM, and thus can take advantage of the sub-
>
> lX ;n stantial computing power of modern multi-processor machines
> Dr B P
B
>
> N nþ2=4
þ Nlnþ2=4 DlXA ;j for C n < 0
: B XB XA (Cowles, 2008). This version retains all options available for
j¼1
source terms in the SWAN model. The essential difference be-
where, respectively, N nA and N nB are the nth time step action densities tween the two models is in solving for wave propagation in geo-
at nodes of A and B; XA (light dashed area in Fig. 1) and XB are the graphic space by the unstructured-grid finite-volume algorithm.
total area of triangles with central nodes at A and B; lXA ;j (j=1, lXA ;n Þ However, as FE-WAVE, as developed by Hsu et al. (2005), also
and lXB ;j (j=1, lXB ;n Þ are the perimeters of XAand XB, lXA ;n , lXB ;n are the uses the same triangular mesh approach as described in
number of edges of XA and XB; and DrAand DrB are the distances FVCOM-SWAVE, these two versions can both be run using the
from node A and node B to the centroids of triangles connected to same grid, but with different spatial discretization (finite-volume
nodes A and B. Cn > 0 refers to the outward direction of the control and finite-element) methods.
volume. This is the second-order approximate advection scheme
used to solve the tracer equation in FVCOM.
In the semi-implicit approach, we have 3. Validation experiments
8 lX ;n
>
> nþ3=4 PA Four idealized benchmark tests were used to validate FVCOM-
>
> NA
> þ DXrAA Nnþ2=4
lX ;j DlXA ;j for C n > 0
< j¼1 A SWAVE with the standard SWAN model by Booij et al. (1999).
Nl ¼ lX ;n : ð18Þ These tests were designed to investigate the numerical diffusion
>
> nþ3=4 PB
nþ2=4
>
> NB þ DXrBB NlX DlXA ;j for C n < 0 of the discrete schemes, and to examine the properties of wave-
>
: j¼1 A
current interactions, and wave propagation over varying shallow
water topography. We also did standard growth curve analysis
in idealized fetch-limited cases, following the SWAMP Group
(1985), with a constant wind blowing seaward off a long straight
coastline (as shown in Booij et al. (1999)). The same procedures
are used to compare FVCOM-SWAVE and SWAN in all four ideal-
ized tests. Comparisons are made directly with the standard
SWAN model, which has been previously very carefully compared
and calibrated to growth curve data (by Booij et al. (1999)).
Inter-model or model-data comparisons in this text were made
based on both RMS and RMS (%). For the significant wave height,
Hs,the RMS (%) is defined as
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
u
u 1 X N
RMS ¼ t ½Hs;m ðiÞ  Hs;o ðiÞ2 ð20Þ
N  1 i¼1
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
u N  2
u 1 X Hs;m ðiÞ  Hs;m jmin Hs;o ðiÞ  Hs;o jmin
RMSð%Þ ¼ t  : ð21Þ
N  1 i¼1 Hs;m jmax  Hs;m jmin Hs;o jmax  Hs;o jmin

In the model-data comparison cases, Hs,m and Hs,o are the model-
predicted and observed Hs, respectively. Here, N is the sample num-
ber, and the subscripts ‘‘max” and ‘‘min” indicate maximum and
minimum. For the dominant wave period, RMS (%) is calculated
by replacing Hs by Tp (dominant peak wave periods) in expression
Fig. 1. Schematic of the unstructured triangular grid used for geographical spatial (21). In these inter-model comparison cases, the subscript ‘‘m” rep-
discretization in FVCOM-SWAVE. Definitions of variables are provided in the text. resents FVCOM-SWAVE and ‘‘o” represents SWAN.
156 J. Qi et al. / Ocean Modelling 28 (2009) 153–166

3.1. Numerical diffusion Respectively, at distances of 1 and 5 km in the x- and y- axes away
from the source, the width of the wave height field is 1.75 and
Consider a harmonic, long-crested wave propagating through a 2.5 km for the first order BSBT, 1.0 and 1.5 km for the second order
gap into a square computational domain with dimensions SORDUP, 1.5 and 1.75 km for the third-order Stelling–Lendretse
10  10 km in deep water (Fig. 2). The open gap is located in the scheme, and 0.75 and 1.75 km for FVCOM-SWAVE. At these two
lower left corner of the domain, so that the wave propagates along locations, FE-WAVE gives the same widths as FVCOM-SWAVE. If
the diagonal at an angle of 45° with respect to the positive x-axis (x only the non-stationary numerical schemes are considered, we find
is the east–west coordinate which is positive in the eastward direc- that the second-order finite-volume upwind scheme used in
tion). This harmonic wave is simulated using a Gaussian-shaped FVCOM-SWAVE can reach the same level of numerical accuracy
frequency spectrum with a peak frequency of 0.1 Hz, a standard as the third-order Stelling–Lendretse scheme, and can exceed this
deviation of 0.01 Hz and a resolution defined as 3% of the relative level of accuracy in the region close to the gap. In this special case,
frequency. The significant wave height at the gap is 1.0 m, and the effective horizontal resolution is the same for both SWAN and
the long crest of the wave is calculated using an assumed cos500(h) FVCOM-SWAVE, because the computational cells have equivalent
directional distribution. area. Because BSBT is a first-order approximation, it is not surpris-
The computational domain is tesselated with a square grid for ing that this scheme (used in SWAN) generates the largest numer-
SWAN and with right triangles for FVCOM-SWAVE. The right trian- ical diffusion. In the most current version of SWAN, this scheme is
gles are constructed by dividing each square along its diagonal line. only used in cells connected to solid boundaries for both stationary
The horizontal resolution for SWAN is 100 m, which is the same for and non-stationary waves.
FVCOM-SWAVE, where the horizontal resolution is defined using In executing this comparison experiment, we have followed the
the shortest edge of a computational cell. The resolution in direc- approaches used in Booij et al. (1999) and Hsu et al. (2005). It
tional space is 0.5°. Here, the time step in SWAN and FVCOM- should be noted that this benchmark test case only provides the
SWAVE is specified to maintain stability. relative numerical diffusion between these different finite-differ-
For this case, SWAN was run using three different discrete ence and finite-volume schemes, because the accuracy of the wave
schemes: (a) the first-order, backward space and backward time field can be affected significantly by the finite-resolution discreti-
scheme for stationary waves (BSBT), (b) the second-order upwind zation that is used in frequency and directional space. A more
iteration scheme for stationary waves (SORDUP) (Rogers et al., objective comparison should be made based on the so-called gar-
2002), and (c) the third-order Stelling–Lendrertse scheme. den sprinkler problem (SWAMP Group, 1985). This problem was
FVCOM-SWAVE was run using both explicit and implicit second- not examined in our current SWAN and FVCOM-SWAVE compari-
order finite-volume upwind schemes. son experiments.
The SWAN and FVCOM-SWAVE significant wave height distri-
butions are compared in Fig. 2. The width of the spreading of the 3.2. Wave-current interactions
significant wave height field is used as an index for numerical dif-
fusion. In general, the SORDUP result has the least numerical diffu- This idealized case considers propagation of waves in a deepwa-
sion, whereas BSBT gives the largest numerical diffusion. ter region in the presence of a uniform background current. Four

Fig. 2. Spatial distributions of the significant wave height (m) for a harmonic wave propagating along a diagonal line in a square computational domain. (a) SWAN-BSBT; (b)
SWAN-SORDUP; (c) SWAN-SL; and (d) FVCOM-SWAVE.
 J. Qi et al. / Ocean Modelling 28 (2009) 153–166 157

cases are tested: (a) wave propagation opposite to the direction of

the current; (b) wave propagation in the same direction as the cur-
rent; (c) wave propagation against the current at an angle of 30°;
and (d) wave propagation in the direction of the current at an angle
of 30°. These four cases are shown schematically in Fig. 3.
The computational domain consists of a rectangle of length
4 km and width 10 km. The assumed waves have the same distri-
bution and shape as those in the previous test case (3.1) and they
are specified along the left boundary with a direction of propaga-
tion as shown in Fig. 3. In cases (a) and (b), the speed of the current
increases in the down-wave direction, ranging from 0.0 to 2.0 m/s.
In cases (c) and (d), the current is rotated by 90°. The SWAN imple-
mentation employed a rectangular grid with a horizontal resolu-
tion of 40 m in the direction of the wave propagation (x-axis)
and 100 m in the cross-wave direction (y-axis). To maintain the
same spatial resolution as SWAN in the direction of wave propaga-
tion, a mesh of equilateral triangles with an edge length of 40 m
was used for FVCOM-SWAVE. SWAN was run using the third-order
Stelling–Leendertse scheme and FVCOM-SWAVE was run using the
second-order semi-implicit upwind scheme.
Given the same resolution in the wave propagation direction,
the distributions of significant wave height and mean wave direc-
tion in the x-direction calculated by SWAN and FVCOM-SWAVE are
almost identical (Fig. 4) and match well with the analytical solu-
tion. The root-mean-square (RMS) differences of significant wave
height (Hs) and direction (h) between SWAN and FVCOM-SWAVE
for the four cases in this experiment are calculated. The correlation
coefficients of these two models for Hs and h are essentially 1.0, and
the RMS (%) difference (see the definition given in Eq. (21)) is 0.44
Fig. 4. Cross-isobath distributions of the SWAN and FVCOM-SWAVE predicted
(or less) for Hs and 0.004 (or less) for h. This demonstrates that significant wave heights and wave propagation directions for cases (a), (b), (c), and
the second-order finite-volume advection scheme used in (d) shown in Fig. 3. Solid line: SWAN, and dashed line: FVCOM-SWAVE.
FVCOM-SWAVE has an equivalent accuracy as the third-order fi-
nite-difference scheme used in SWAN. For this idealized domain, of less than 3.0 km, but a significant bias occurs over the remainder
the fact that the wave features from structured-grid SWAN runs of the domain. In this test case, the unstructured-grid finite-vol-
can be exactly reproduced using an unstructured-grid finite-vol- ume algorithm used in FVCOM-SWAVE seems to show more prom-
ume model suggests that the finite-volume flux algorithm used ise than the unstructured-grid finite-element method used by Hsu
in FVCOM-SWAVE can ensure high numerical accuracy with little et al. (2005). However, as the FE-WAVE model is not an open
influence from the irregular grid. This suggests that the geometric source code, the comparison made here is only based on Figs. 2
grid flexibility of FVCOM-SWAVE makes it a good candidate to sim- and 3 in Hsu et al. (2005).
ulate wave-current interactions in complex coastal regions with
strong currents, headlands, islands, and irregular channels. 3.3. Wave shoaling and refraction
Hsu et al. (2005) used FE-WAVE model to conduct the same
experiment. The agreement of their solution with the original Consider surface waves propagating towards a straight beach.
SWAN varies with the x-axis. The match is good over a distance The water depth varies from the deeper outer region with a depth
of 20 m, to zero at the coast over a distance of 4.0 km. Two cases
are examined: (a) incident waves propagating in the direction per-
pendicular to the coastline and (b) incident waves propagating to-
wards to the coast with an angle of 30° (Fig. 5). These two cases are
designed to examine wave shoaling and refraction without wave-
current interaction.
The wave model boundary conditions and set-up at the deep
open boundary are the same as those given in Section 3.2. The hor-
izontal computational domain consists of a rectangular area with a
width of 20 km (y-axis) and a length of 4 km (x-axis). For SWAN,
we use a rectangular grid with a resolution of 40 m in the x-axis
(cross-isobath direction) and 200 m in the y-axis (along-isobath
direction). To have the same resolution in the cross-isobath direc-
tion, FVCOM-SWAVE uses an equilateral triangular grid with a
length of 40 m.
Running SWAN with the third-order Stelling–Leendertse scheme
produces exactly the same results as those shown in Booij et al.
(1999). The model solution matches the analytical solution from
linear wave theory, with an error of 0.1% in significant wave height
and of <0.25° in direction within the region with water depth
Fig. 3. Schematic of the surface wave propagation in a domain with an ambient >0.05 m. Given the same horizontal resolution, FVCOM-SWAVE
cg is the wave group velocity and ~
water current. ~ V is the water current vector. predicts a nearly identical solution to that of SWAN for both the
158 J. Qi et al. / Ocean Modelling 28 (2009) 153–166

Table 1
Parameters used in SWAN and FVCOM-SWAVE setups

Parameters Definition Value

cds2 Coefficient used to determine the rate of whitecapping 2.36  105
dissipation
stpm Value of the wave steepness for the Pierson–Moskowitz 3.02  105
a Proportionality coefficient in the case with inclusion of 0.0015
Cavaleri and Malanotte’swave growth term
alpha Proportionality coefficient of the rate of dissipation 1.0
gamma The ratio of maximum individual wave height over 0.73
depth
cfjon Coefficient of the JONSWAP formulation 0.067
iquad Fully explicit computation of the nonlinear transfer 2
with DIA per iteration
lambda Coefficient for quadruplet configuration in the case of 0.25
DIA
Cnl4 Proportionality coefficient for quadruplet interactions 3  107
in case of DIA
Csh1 Coefficient for shallow water scaling in case of DIA 5.5
Csh2 Coefficient for shallow water scaling in case of DIA 6/7
Csh3 Coefficient for shallow water scaling in case of DIA 1.25

Fig. 5. Schematic of the surface wave propagating towards the coast from the
consisting of a length of 400 km in the along-coast direction
deepwater region to the shallow shelf. Case (a): shoaling experiment; case (b): (y-axis) and a width of 9000 km in the cross-coast direction
refraction experiment. (x-axis). The ocean is assumed to be infinitely deep. We compared
the growth curves of total wave energy and peak frequency as esti-
mated by FVCOM-SWAVE and SWAN. Both FVCOM-SWAVE and
shoaling and refraction cases (Fig. 6). Considering the same domain SWAN are run with a cross-coastal resolution of 10 km. The param-
with water depth >0.05 m, the correlation coefficient of these two eters used in the model setups are based on the values listed in
model results is essentially 1.0 for all cases, and the RMS(%) differ-
ence is 1.84 (or less) for Hs and 0.63 (or less) for h. However, FE-
WAVE shows a relatively larger bias in comparison with SWAN
(Figs. 4 and 5 in Hsu et al. (2005)), and in fact the bias suggested
for FE-WAVE appears to be larger than that for FVCOM-SWAVE.

3.4. Growth curves for wind-generated waves

In studying growth curves for these two model systems, we

used the same source term formulations for FVCOM-SWAVE as
we used for SWAN (the so-called WAM cycle 3 physics). We as-
sumed a constant wind of 20 m/s at 10-m reference height off a
long, straight coastline in a rectangular domain with dimensions

4
SWAN
3 FVCOM-SWAVE
Hs(m)

2
(a)
(b)
1

40

30

20 (a)
θ(º)

10
(b)
0

-1 0
0 1 2 3 4
Distance (km)

Fig. 6. Cross-isobath distributions of the SWAN and FVCOM-SWAVE predicted Fig. 7. Distributions of dimensionless fetch-limited growth curves for the total
significant wave height and wave propagation directions for cases (a) and (b) shown wave energy E* (a) and peak frequency fp as a function of dimensionless fetch for
in Fig. 5. Solid line: SWAN, and dashed line: FVCOM-SWAVE. FVCOM-SWAVE, stationary and non-stationary SWAN.
 J. Qi et al. / Ocean Modelling 28 (2009) 153–166 159

Table 1, as documented by the SWAN Team (2006a,b), with the (NA) (Fig. 8). The GoM bathymetry features several deep basins,
WAM cycle 3 physics. submarine banks, and shallow shelves connected to coastal inlets,
In this case, FVCOM-SWAVE and SWAN show similar growth bays and estuaries. One of the primary objectives of developing
curves in dimensionless wave energy (E ¼ g 2 E=U 4 Þ and dimen- FVCOM-SWAVE is to include it in the FVCOM-based unstruc-
sionless peak frequency (fp ¼ U  fp =gÞ as functions of dimensionless tured-grid Northeast Coastal Ocean Forecast System (NECOFS)
fetch (X  ¼ gX=U 2 Þ (Fig. 7). Here E is the total wave energy, g is the (see http://fvcom.smast.umassd.edu/). Thus, the coupled model
gravitational acceleration, U* is the friction velocity at the sea sur- system would be able to make predictions of wind-induced surface
face, namely U 2 ¼ C D U 210 (where CD is the drag coefficient following waves and their coupling with other physical processes (e.g., wave
Wu (1982), U10 is the wind speed at the 10-m height above the sea setup, wave-current interaction, sediment transport) in this region,
surface) and X is the fetch. with its complex inner shelf bathymetry and irregular coastline. To
For large fetch, X* > 106, the dimensionless total energies of validate the accuracy and reliability of FVCOM-SWAVE for this re-
FVCOM-SWAVE and SWAN (stationary and non-stationary cases) gion, we used this model to simulate January 2007.
are within the variations shown by models considered by the In fact, the GoM FVCOM-SWAVE grid covers the entire US east
SWAMP Group (1985). For smaller fetch, X* < 106, the growth curve coast region with an open boundary running from a land point at
for FVCOM-SWAVE is very close to that shown for the stationary about 65°W, 10°N towards the northeast to about 20°W, 50°N
SWAN run; however these curves are significantly lower than that and then turning towards the east coast of Greenland (Fig. 9).
resulting from the non-stationary SWAN run. Similarly, both The computational domain consists of the unstructured triangular
FVCOM-SWAVE and non-stationary SWAN-predicted peak frequen- grid, with a horizontal resolution varying from 25 to 50 km in the
cies are within the range of variation examined by the SWAMP interior of the NA to 0.5–1.0 km along the coast of the GoM, New
Group (1985), except for very small fetches of X*  105. In this case, England Shelf, Long Island Sound and over Georges Bank
the SWAN results are identical to those shown by Booij et al. (Fig. 10). The total number of nodes is 30760 and total number
(1999). The results of this comparison are consistent with the find- of cells is 57013. Although our main interest is in the GoM and
ings of Hsu et al. (2005). The latter pointed out that in its non-station- adjacent coastal regions, running FVCOM-SWAVE in this large do-
ary mode, SWAN significantly overestimates the total wave energy main can minimize open boundary issues. The flexibility in the
growth curve as a function of fetch where X* < 106. Here, we show spatial discretization enabled through the use of unstructured
that this bias is resolved by the non-stationary mode of FVCOM- meshes allows coarse resolution in the open ocean and fine resolu-
SWAVE. tion in the coastal region, with limited reduction in computational
efficiency. In January 2007, westerly wind prevailed over the GoM
region. The dominant wind patterns blew from the land to the
4. Application to the Gulf of Maine ocean, so that essentially no wave forcing was required to specify
the remote open boundary. In this application experiment, no
4.1. Design of numerical experiments ambient currents were included.
The GoM FVCOM-SWAVE was driven by the wind fields derived
The Gulf of Maine (GoM) is located along the US northeast coast. from the GoM-WRF (Weather Research and Forecast) and the NA-
It is a semi-enclosed basin opening to the North Atlantic Ocean WRF models. The GoM-WRF model was established as a part of NE-

Fig. 8. Locations of NDBC environmental buoys in the US northeast coastal region used for the model-data comparison. GoM, Gulf of Maine; BF, Bay of Fundy; GB, Georges
Bank; LIS, Long Island Sound; CB, Chesapeake Bay. Filled triangles: buoys located in zone 1 (water depths <50 m); filled diamonds: buoys located in zone 2 (water depths in
50–200 m range); and filled circles: buoys located in zone 3 (water depths >200 m). The number listed on each buoy is its name.
160 J. Qi et al. / Ocean Modelling 28 (2009) 153–166

Table 2
Correlations and RMS (%) errors of assimilated wind speed via observations at buoy
stations.

Buoy station Sample Correlation RMS error

identification number coefficient (%)
44004 734 0.84 1.25
44005 742 0.95 1.11
44007 743 0.87 1.48
44009 744 0.87 1.29
44013 739 0.92 0.99
44014 741 0.84 1.46
44017 743 0.87 1.29
44018 698 0.94 0.83
44024 692 0.92 0.98
44025 741 0.87 1.32
44027 743 0.94 0.85
44029 740 0.88 1.38
44030 740 0.90 1.17
44031 698 0.87 1.45
44032 728 0.90 1.14
44033 724 0.85 1.71
44034 704 0.92 1.07
44035 600 0.74 2.73
44037 668 0.94 0.82
44038 728 0.93 0.89
Mean 0.89 1.26

Fig. 9. Unstructured grid of the FVCOM-SWAVE configured for the US northeast in the NA-WRF region, it causes little influence on the wave simu-
coastal region. lation in the GoM and adjacent coastal regions where only the
hourly GoM-WRF winds were used.
The source terms used in GoM FVCOM-SWAVE were the same
COFS and uses a 9-km horizontal resolution covering the entire as those used in the SWAN simulations. The wave growth due to
GoM and adjacent New England Shelf and Scotian Shelf. Hourly wind input (Sin) is determined by the sum of linear and exponen-
wind velocity data from all meteorological buoys in these regions tial growths. The linear term is specified by an empirical formula
were assimilated into this weather model to produce a reliable of Cavaleri and Malanotte-Rizzoli (1981) with a filter for the trun-
and accurate wind field. Since most buoys were located near the cation at the Pierson–Moskowitz frequency (Pierson and Mosko-
western NA coast, the assimilation helped reduce the RMS errors witz, 1964). The exponential wave growth term is specified
of the WRF winds in this region (Table 2). In the region outside using Komen et al. (1984)’s empirical expression that is a function
of the GoM-WRF domain, we used the reanalysis 3-hourly 32-km of inverse wave age U*/cph (where cph is the phase speed). The dis-
resolution NA-WRF-predicted wind fields (available at http://no- sipation includes white-capping, Sds,w, bottom friction, Sds,b, and
mads.ncdc.noaa.gov/). Linear interpolation was used to merge depth-induced breaking, Sds,br. Sds,w is derived from Hasselmann
NA-WRF and GoM-WRF outputs to form an hourly wind-forcing (1974)’s pulse-based model, Sds,b follows from the ‘‘JONSWAP”
field covering the entire computational domain. Although this ap- empirical function with Collins (1972)’s drag law expression
proach tends to ‘‘filter” out wind variations with time scales < 3 h and Madsen et al. (1988)’s eddy-viscosity formulation. Sds,br is

46

45

44
Latitude ( N)

43

42

41

40

39
75 74 73 72 71 70 69 68 67 66 65 64 63 62 61 60
Longitude ( W)

Fig. 10. An enlarged view of the unstructured triangular grid in the Gulf of Maine and New England Shelf region. The horizontal resolution of the grid varies from 25 km off
Georges Bank to 0.5–1.0 km near the coast.
 J. Qi et al. / Ocean Modelling 28 (2009) 153–166 161

12

Wind velocity (m/s)

Buoy 44004
8
4
0
-4
-8
-12

12
Wind velocity (m/s)

Buoy 44032
8
4
0
-4
-8
-12

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30
Days (January 2007)

Fig. 11. Time series of wind velocity at the 10-m height above the sea surface at buoy 44004 and buoy 44032 for January 1–31, 2007.

estimated by the bore-based function derived by Battjes and Jans- of water depth, we divided the computational domain into three
sen (1978). The nonlinear transfer of wave energy due to wave– zones. Letting h be the mean water depth, zone 1 refers to the re-
wave interactions is calculated by a sum of Snl3 (triad wave inter- gion in which h < 50 m; zone 2 is the region with
actions) and Snl4 (quadruplet wave interactions). Snl3 is estimated 50 m 6 h 6 200 m, and zone 3 is the region with h > 200 m. There
from the Lumped Triad Approximation (LTA) derived by Elde- are 8 buoys in zone 1, 11 buoys in zone 2 and 2 buoys in zone 3.
berky (1996) and Snl4 is specified by the Discrete Interaction
Approximation (DIA) by Hasselmann et al. (1985). The parameters
used in this experiment are shown in Table 1; detailed definitions 4.2. Model-data comparisons
are given by the SWAN Team (2006b).
The GoM FVCOM-SWAVE results are compared with significant During January 2007, the GoM surface weather was character-
wave height and period data recorded by 21 NOAA buoys. The buoy ized by cold-air frontal passages with a typical time scale of 3–5
locations are shown in Fig. 8 and the data are available at the Na- days before January 20; thereafter, the area was dominated by
tional Data Buoy Center (http://www.ndbc.noaa.gov/hmd.shtml). southeastward to southward winds that were comparatively more
To evaluate the performance of GoM FVCOM-SWAVE as a function steady (Fig. 11). The maximum 10-m wind speed was 10–13 m/s.

Table 3 Table 4
Mean and standard deviations (SD) of observed and calculated significant wave Mean and standard deviations (SD) of observed and calculated peak wave periods,
heights, correlation coefficients, and RMS errors at buoy stations. correlation coefficients, and RMS errors at buoy stations.

Station Mean ± SD Mean ± SD Correlation RMS RMS Station Mean ± SD Mean ± SD Correlation RMS RMS
(sample (observed) (m) (calculated) coefficient (m) (%) (Sample (observed) (calculated) coefficient (sec) (%)
number) (m) number) (sec) (sec)
Zone 1 (h < 50 m) Zone 1 (h < 50 m)
44033 (701) 0.73 ± 0.36 1.15 ± 0.34 0.60 0.53 24.9 44033 (659) 5.84 ± 2.71 6.27 ± 2.28 0.22 3.17 33.0
44031 (693) 1.06 ± 0.52 0.93 ± 0.37 0.75 0.37 21.2 44031 (679) 6.30 ± 2.65 6.82 ± 2.38 0.27 3.09 32.0
44007 (733) 0.95 ± 0.50 0.87 ± 0.35 0.73 0.35 22.4 44007 (733) 7.36 ± 3.12 7.10 ± 2.23 0.41 3.00 26.7
44017 (733) 1.51 ± 0.64 1.41 ± 0.50 0.84 0.36 17.6 44017 (733) 6.64 ± 1.93 6.87 ± 1.64 0.28 2.17 27.8
44025 (726) 1.42 ± 0.62 1.39 ± 0.49 0.83 0.35 16.4 44025 (726) 6.36 ± 1.96 6.50 ± 1.76 0.38 2.08 23.4
44009 (739) 1.30 ± 0.52 1.30 ± 0.42 0.75 0.35 17.0 44009 (739) 6.01 ± 2.13 6.51 ± 2.26 0.34 2.58 23.4
44014 (741) 1.62 ± 0.74 1.53 ± 0.47 0.86 0.41 11.9 44014 (741) 7.27 ± 1.87 7.05 ± 2.03 0.17 2.52 23.5
Mean 1.23 ± 0.56 1.23 ± 0.42 0.77 0.39 18.8 Mean 6.54 ± 2.34 6.73 ± 2.08 0.30 2.66 27.1
Zone 2 (50 m 6 h 6 200 m) Zone 2 (50 m 6 h 6 200 m)
44027 (742) 1.77 ± 0.92 1.62 ± 0.69 0.83 0.54 13.0 44027 (742) 6.85 ± 2.15 6.78 ± 2.02 0.49 2.10 23.2
44034 (677) 1.62 ± 0.85 1.46 ± 0.62 0.83 0.51 20.6 44034 (677) 6.92 ± 2.62 7.09 ± 2.02 0.34 2.71 24.9
44032 (726) 1.40 ± 0.68 1.27 ± 0.49 0.71 0.50 24.3 44032 (723) 6.59 ± 2.60 6.55 ± 2.30 0.39 2.72 30.4
44030 (736) 0.99 ± 0.44 0.98 ± 0.36 0.67 0.33 20.3 44030 (735) 6.13 ± 2.79 6.90 ± 2.78 0.25 3.50 33.6
44029 (734) 0.99 ± 0.44 0.95 ± 0.36 0.64 0.35 36.4 44029 (734) 5.54 ± 2.42 5.97 ± 2.79 0.41 2.89 28.6
44013 (737) 0.93 ± 0.46 0.88 ± 0.40 0.72 0.33 20.5 44013 (737) 5.31 ± 2.88 6.10 ± 3.04 0.46 3.18 32.5
44018 (558) 1.78 ± 0.70 1.85 ± 0.67 0.86 0.36 15.3 44018 (558) 7.32 ± 1.69 7.19 ± 1.56 0.38 1.83 21.1
44008 (729) 2.33 ± 0.90 1.99 ± 0.74 0.87 0.57 13.2 44008 (729) 7.54 ± 1.64 7.21 ± 1.48 0.40 1.75 20.5
44038 (719) 2.18 ± 1.06 1.99 ± 0.88 0.86 0.57 11.2 44038 (719) 7.24 ± 2.00 7.09 ± 1.55 0.44 1.92 24.9
44005 (479) 1.81 ± 0.80 1.52 ± 0.60 0.83 0.50 12.8 44005 (479) 6.51 ± 1.98 6.89 ± 2.19 0.46 2.42 25.6
44024 (700) 2.71 ± 1.30 2.21 ± 0.87 0.86 0.88 13.0 44024 (700) 8.09 ± 2.06 7.75 ± 1.78 0.44 2.09 29.6
Mean 1.68 ± 0.78 1.16 ± 0.61 0.79 0.49 18.2 Mean 6.73 ± 2.26 6.86 ± 2.14 0.41 2.46 26.8
Zone 3 (h> 200 m) Zone 3 (h > 200 m)
44037 (690) 2.03 ± 0.97 1.88 ± 0.80 0.83 0.56 12.7 44037 (690) 6.80 ± 1.91 6.73 ± 1.57 0.43 1.88 23.7
44004 (734) 2.70 ± 1.26 2.23 ± 0.89 0.91 0.76 17.4 44004 (734) 7.81 ± 1.65 7.35 ± 1.61 0.58 1.57 15.3
Mean 2.37 ± 1.11 2.06 ± 0.85 0.87 066 15.1 Mean 7.30 ± 1.78 7.04 ± 1.59 0.51 1.73 19.5
162 J. Qi et al. / Ocean Modelling 28 (2009) 153–166

As shown in Fig. 11, the wind speed and direction varied signifi- evident at buoy 44032 (50-m isobath) and buoy 44004 (3182-
cantly from the near-coastal regions to the open ocean, which is m isobath).

Observed FVCOM-SWAVE
20
6 Buoy 44007 (23.7m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44009 (28m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44014 (47.5m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44017 (44.8m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30

20
6 Buoy 44025 (36.3m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44031 (20.1m)
15
TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44033 (<50m)
15
TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44035 (10.1m)
15
TD (sec)
Hs (m)

4
10
2 5

0 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Days (January 2007) Days (January 2007)

Fig. 12. Comparisons of the GoM FVCOM-SWAVE predicted (solid line) and observed (dashed line) significant wave heights (left) and dominant wave periods (right) at buoys
in zone 1 (water depth <50 m).
 J. Qi et al. / Ocean Modelling 28 (2009) 153–166 163

The GoM FVCOM-SWAVE simulation performed reasonably SWAVE reasonably captured the temporal variation in wave
well in capturing the spatial and temporal variability of wind- heights that were observed (Fig. 13). Cr was 0.64–0.87, with
induced surface waves along the coast. Tables 3 and 4 show RMS (%) error varying from 11.2 to 36.4 with a mean RMS (%)
the means, standard deviations (SD), correlation coefficients, error of 18.2. In zone 3, the peaks of the observed Hs values
RMS and RMS (%) errors of Hs and Tp at buoy locations. The are similar to those found in zone 2. The GoM FVCOM-SWAVE
mean and SD values are calculated based on the measured data reliably simulated the Hs peaks (Fig. 14), with Cr in the range of
used in this comparison. In zone 1 (in water depth <50 m), the 0.83–0.91, RMS(%) errors of 12.7 and 17.4 and a mean RMS(%)
observations show that the peak Hs values are 1–3 m in the error of 15.1.
near-coastal regions with mean values of 0.7–1.6 m and stan- The model-predicted Tp values have a low correlation with
dard deviations of 0.4–0.7 m. We neglect buoy 44035 in this observations in all three zones, although the model clearly cap-
comparison because it was located in a sheltered area (in Pas- tures the wind-induced temporal variation pattern for Tp. The
samaquoddy Bay) not open to wind-generated waves of the model-predicted means and standard deviations of Tp are in
GoM and had high errors in the wind fields and recorded very agreement with the range of the observations. The mean RMS
low values for Hs 0.1–0.3 m. These spatial and temporal dis- (%) error for zone 1 is 27.1, which is comparable with our value
tributions were well captured by the GoM FVCOM-SWAVE of 26.8 for zone 2. A close examination of the differences in the
(Fig. 12). The model-predicted Hs estimates were in good agree- model and observed Tp time series in Figs. 12–14 shows that the
ment with 7 near-shore coastal buoys ranging from the Bay of observations exhibit more notable variation with time than the
Fundy to the Chesapeake Bay. The model predicted means and model. In particular, the GoM FVCOM-SWAVE estimates show
standard deviations of Hs show the same range as observations: a tendency to under-predict observed Tp peaks. Although this
the RMS (%) error varying from 11.9 to 24.9 with a mean RMS can be attributed to biases in the wind fields and model inade-
(%) error of 18.8. The largest observed Hs values occurred in quacies, it results from the common tendency of spectra to have
zone 2, in the range of 1–7 m with maximum mean and stan- multiple peaks and thus Tp estimates vary among the competing
dard deviation values of 2.7 and 1.3 m. The GoM FVCOM- spectral peaks.

Observed FVCOM-SWAVE
20
6 Buoy 44029 (65.4m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44030 (>50m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44032 (>50m)
15
TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44034 (>50m)
15
TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44038 (>50m)
15
TD (sec)
Hs (m)

4
10
2 5

0 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Days (January 2007) Days (January 2007)

Fig. 13. Comparisons of the GoM FVCOM-SWAVE predicted (solid line) and observed (dashed line) significant wave heights (left) and dominant wave periods (right) at buoys
in zone 2 (water depth in the range 50–200 m).
164 J. Qi et al. / Ocean Modelling 28 (2009) 153–166

Observed FVCOM-SWAVE
20
6 Buoy 44005 (195.7m)
15

TD (sec)
Hs (m)
4
10
2 5

0 0
20
6 Buoy 44008 (62.5m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44013 (55m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44018 (74.4m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44024 (180.1m)
15

TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44027 (182m)
15 TD (sec)
Hs (m)

4
10
2 5

0 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Days (January 2007) Days (January 2007)

Fig. 13 (continued)

Observed FVCOM-SWAVE
20
6 Buoy 44004 (3182.1m)
15
TD (sec)
Hs (m)

4
10
2 5

0 0
20
6 Buoy 44037 (251.2m)
15
TD (sec)
Hs (m)

4
10
2 5

0 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30
Days (January 2007) Days (January 2007)

Fig. 14. Comparisons of the GoM FVCOM-SWAVE predicted (solid line) and observed (dashed line) significant wave heights (left) and dominant wave periods (right) at buoys
in zone 3 (water depth >200 m).
 J. Qi et al. / Ocean Modelling 28 (2009) 153–166 165

5. Conclusions OOS Program, NSF projects (OCE0234545, OCE0606928,

OCE0712903, OCE0732084 and OCE0726851) and MIT Sea Grant
Compared to the conventional SWAN formulation, FVCOM- project (NA06OAR41700019). Support also comes from the Cana-
SWAVE provides an alternative version of the wave model based dian Panel on Energy Research and Development (PERD) and Go-
on an unstructured-grid finite-volume approach. In particular the MOOS – the Gulf of Maine Ocean Observing System.
latter is more suitable for application in coastal and estuarine re-
gions characterized with irregular coastal geometry. Four idealized
benchmark test problems in rectangular domains are used to dem- Appendix A. Discretization of the action density equation in
onstrate that the second-order unstructured-grid finite-volume geographic space
method used in FVCOM-SWAVE has the same level of accuracy
as the third-order finite-difference method used in SWAN. At a The time integration equation of action density in geographic
fourth test, fetch-limited growth curve analysis following SWAMP space (Eq. (5) in the text) is rewritten here as
Group (1985) shows that FVCOM-SWAVE behaves acceptably in 3
Nnþ4  Nnþ4
2
*
comparison to SWAN growth curves. Finally, an application of þ r  ðC NÞ ¼ 0 ðA1Þ
Dt
the models to the Northwest Atlantic, particularly the Gulf of * * *
Maine and related US northeast coastal waters, suggests that where C ¼ C þ V . Integrating (A1) in the control area X (see Fig. 1)
g
FVCOM-SWAVE is robust and can capture the temporal and spatial yields
variation of waves generated by Nor’easters and other high-wind
I
3 2 Dt * 2 Dt
events over both continental shelf and near-shore regions. Nnþ4 ¼ Nnþ4  t r  ðC NÞdX ¼ Nnþ4  C n Nl dl: ðA2Þ
X X X l
Our experiments indicate that the default parameterizations of
SWAN (which is adopted by FVCOM-SWAVE) tend to underestimate In the second-order approximated upwind scheme,
8
the significant wave heights during strong wind events. Similar
< NlA
> for C n > 0
problems were also reported in SWAN’s application by Rogers
Nl ¼ N lB for C n < 0 ðA3Þ
et al. (2002), who found that SWAN tended to underpredict low- >
:
and medium-frequency energy in a wind-sea portion of the spec-
trum. The accuracy of the SWAN simulation is directly affected by where
parameterization of whitecapping dissipation, as well as nonlinear
wave–wave interactions, which are central to the growth and devel- NlA ¼ NA þ r  NA Dr A ðA4Þ
opment of wind-generated waves (Zhang et al., 2006; Zhang and Per- NlB ¼ NB þ r  N B Dr B : ðA5Þ
rie, 2008). The drag coefficient parameterization for surface wind
Integrating (A4) in area XA (a sum of the areas of the light gray tri-
stress is somewhat old, dating from more than 25 years ago. In the
angles), we have
current version of SWAN, the surface wind stress is calculated by I
an empirical function. This function does not be capture the wave- Dr A Dr A
NlA ¼ NA þ t r  NA dXA ¼ NA þ NlX dlXA : ðA6Þ
age dependence advocated in modern parameterizations for wave XA XA XA lX
A
A
drag (Zhang and Perrie, 2008). It is well know that physical processes
Similarly, integrating (A5) in area XB(a sum of the areas of triangles
in wave models have definite limitations, particularly the DIA for-
with edges shown in dashed line, we have,
mulation for nonlinear wave–wave interactions. Resio and Perrie
I
(2008) pointed out very definite problems in simulating shallow Dr B Dr B
NlB ¼ NB þ t r  NB dXB ¼ NB þ NlXB dlXB ðA7Þ
water waves, which are the focus for our present paper. Wind input XB XB XB lX
B
(Sin) and dissipation (Sds) are usually tuned to compensate whatever
DIA tries to give. This approach, however, cannot compensate for The finite-volume discrete expressions of (A2), (A6) and (A7) can be
deficiencies in shallow water simulations. Wave-current interaction written in the flux form as
also may be a factor that can influence wave evolution. These issues
3 2 Dt X
ln
are being addressed in our ongoing validation of FVCOM-SWAVE. Nnþ4 ¼ Nnþ4  C n;i Nli Dli ðA8Þ
The development of advanced inundation and wave forecast sys- X i¼1
tem has received increased attention in recent years due to high- lX ;n
Dr A X A

wave events combined with massive flooding (e.g., New Orleans NlA ¼ NA þ N DlX ;j ðA9Þ
XA j¼1 lXA ;j A
during Hurricane Katrina in 2005) and other recent major coastal
flooding events where waves were high. Accurate simulation of lX ;n
Dr B X B

the wind-driven surface waves is required for the reliable prediction NlB ¼ NB þ N Dl : ðA10Þ
XB j¼1 lXB ;j XB ;j
of storm surge processes. In order to set up an accurate reliable fore-
cast system for surface waves using FVCOM-SWAVE, new formula-
In the explicit scheme, Nl is given by its value at the (n + 2/4)th time
tions for the source terms must be developed and implemented in
step, while in the implicit scheme, it is an unknown variable at the
the model. Candidates for these new parameterizations are emerg-
(n + 3/4)th time step.
ing in recent years (Resio and Perrie, 2008; Perrie and Resio, 2008)
and we hope to make appropriate tests with a large number of severe
storms. These validation experiments should include a large selec- References
tion of different wind regimes to validate the accuracy, reliability
Battjes, J.A. Janssen, J.P.F.M. 1978. Energy loss and set-up due to breaking of random
and limitation of the model for coastal applications.
waves, Proc. 16th Int. Conf. Coastal Engineering, ASCE, pp. 569–587.
Balay, S., Buschelman, K., Eijkhout, V., Gropp, W., Kaushik, D., Knepley, M., McInnes,
Acknowledgement L.C., Smith, B., Zhang, H. 2007. PETSc Users Manual. Technical Report ANL-95/
11-Revision 2.3.3, Argonne National Laboratory.
Booij, N., Ris, R.C., Holthuijsen, L.H., 1999. A third-generation wave model for coastal
The development of FVCOM-SWAVE was supported by the Mas- regions, Part I, Model description and validation. J. Geophys. Res. 104 (C4),
sachusetts Marine Fisheries Institute (MFI) (DOC/NOAA/NA04NMF 7649–7666.
4720332 and DOC/NOAA/NA05NMF4721131), the NOAA NERAC- Boris, J.P., Book, D.L., 1973. Flux corrected transport I, SHASTA, a fluid transport
algorithm that works. J. Comp. Phys. 11, 38–69.
166 J. Qi et al. / Ocean Modelling 28 (2009) 153–166

Casulli, V., Lang, G. 2004. Mathematical model UnTRIM, validation document, Komen, G.J., Hasselmann, S., Hasselmann, K., 1984. On the existence of a fully
version June 2004 (1.0). Technical report, Federal Waterways Engineering and developed wind-sea spectrum. J. Phys. Oceanogr. 14, 1271–1285.
Research Institute, Wedeler Landstr. 157, D-22559 Hamburg, pp. 78. Lai, Z., Chen, C., Cowles, G., Beardsley, R.C. 2008. A non-hydrostatic version of
Cavaleri, L., Malanotte-Rizzoli, P., 1981. Wind wave prediction in shallow water: FVCOM-validation experiment I: surface standing and solitary waves. J.
Theory and application. J. Geophys. Res. 86 (C11), 10961–10973. Geophys. Res., submitted.
Chen, C., Liu, H., Beardsley, R.C., 2003. An unstructured, finite volume, three- Madsen, O.S., Poon, Y.-K., Graber, H.C. 1988. Spectral wave attenuation by bottom
dimensional, primitive equation ocean model: Application to coastal ocean and friction: Theory. Proc. 21th Int. Conf. Coastal Engineering, ASCE, 492–504.
estuaries. J. Atmos. Oceanic Technol. 20, 159–186. Perrie, W., Resio, D. 2008. A two-scale approximation for efficient representation of
Chen, C., Beardsley, R.C., Cowles, G. 2006a. An unstructured grid, finite-volume nonlinear energy transfers in a wind wave spectrum. Part 2 Application to
coastal ocean model (FVCOM) system. Special Issue entitled ‘‘Advances in observed wave spectra. Submitted to J. Phys. Oceanogr.
Computational Oceanography”. Oceanography, 19(1), 78–89. Pierson, W.J., Moskowitz, L., 1964. A proposed spectral form for fully developed
Chen, C., Beardsley, R.C., Cowles, G. 2006b. An unstructured grid, finite-volume wind seas based on the similarity theory of S.A. Kitaigorodskii. J. Geophys. Res.
coastal ocean model FVCOM user manual. SMAST/UMASSD, 2006. 69 (24), 5181–5190.
Chen, C., Huang, H., Beardsley, R.C., Liu, H., Xu, Q., Cowles, G., 2007. A finite-volume Resio, D., Perrie, W., 2008. A two-scale approximation for efficient representation of
numerical approach for coastal ocean circulation studies: comparisons with finite nonlinear energy transfers in a wind wave spectrum. Part 1: Theoretical
difference models. J. Geophys. Res. 112, C03018. doi:10.1029/2006JC003485. Development. J. Phys. Oceanogr. 38 (12), 2801–2816.
Collins, J.I., 1972. Prediction of shallow water spectra. J. Geophys. Res. 77 (15), Rogers, W.E., Kaihatu, J.M., Petit, H.A.H., Booij, N., Holthuijsen, L.H., 2002. Diffusion
2693–2707. reduction in an arbitrary scale third generation wind wave model. Ocean Engng.
Cowles, G., 2008. Parallelization of the FVCOM Coastal Ocean Model. International 29, 1357–1390.
Journal of High Performance Computing Applications 22 (2), 177–193. SWAMP Group. 1985. Ocean Wave Modelling. Plenum, New York, 256 pp.
Crank, J., Nicolson, P., 1947. A practical method for numerical evaluation of partial SWAN Team. 2006a. SWAN Cycle III version 40.51 Technical documentation. Delft
differential equations of the heat conduction type. Proceedings of the University of Technology, Faculty of Civil Engineering and Geosciences,
Cambridge Philosophical Society 43, 50–64. Environmental Fluid Mechanics Section, P.O.Box 5048, 2600 GA Delft, The
Eldeberky, Y. 1996. Nonlinear transformation of wave spectra in the nearshore Netherlands.
zone. Ph.D. thesis, Delft University of Technology, Department of Civil SWAN Team. 2006b. SWAN Cycle III version 40.51 user manual. Delft University of
Engineering, The Netherlands. Technology, Faculty of Civil Engineering and Geosciences, Environmental Fluid
Fringer, O.B., Gerritsen, M., Street, R.L., 2006. An unstructured-grid, finite-volume, Mechanics Section, P.O. Box 5048, 2600 GA Delft, The Netherlands.
nonhydrostatic, parallel coastal ocean simulator. Ocean Modelling 14, 139–173. Tolman, H.L. 2002. User manual and system documentation of WAVEWATCH-III
Hasselmann, K., 1974. On the spectral dissipation of ocean waves due to version 2.22. Environmental Modeling Center, Marine Modeling and Analysis
whitecapping. Bound.-layer Meteor. 6 (1–2), 107–127. Branch, NOAA.
Hasselmann, S., Hasselmann, K., Allender, J.H., Barnett, T.P., 1985. Computations and WAMDI Group. 1988. The WAM model – A third generation ocean wave prediction
parameterizations of the nonlinear energy transfer in a gravity wave spectrum. model, J. Phys. Oceanogr. 18, 1775–1810.
Part II: Parameterizations of the nonlinear transfer for application in wave Wu, J., 1982. Wind-stress coefficients over sea surface from breeze to hurricane. J.
models. J. Phys. Oceanogr. 15, 1378–1391. Geophys. Res. C87, 9704–9706. doi:10.1029/JC087iC12p09704.
Hsu, T., Ou, S., Liau, J., 2005. Hindcasting nearshore wind waves using a FEM code for Zhang, W., Perrie, W., 2008. The influence of air-sea roughness, sea spray and storm
SWAN. Coastal Engineering 52, 177–195. translation speed on waves. J. Physical Oceanogr. 38 (4), 817–839.
HYPRE Team (2001), HYPRE-high performance preconditions, user’s manual. Center Zhang, W., Perrie, W., Li, W., 2006. Impacts of waves and sea spray on mid-latitude
for Applied Scientific Computing, Lawrence Livermore National Laboratory, 48pp. storm structure and intensity. Monthly Weather Review 134 (9), 2418–2442.