MODELING OF WAVE-CURRENT INTERACTION USING A MULTIDIRECTIONAL WAVE- ACTION BALANCE EQUATION

This study presents an integrated numerical model to simulate wave deformation/transformation in tidal inlets or river mouths with ambient currents (e.g. tidal currents, river inflows) by carefully modeling the effect of wave-current interaction. A multidirectional wave-action balance equation is used to compute random/directional wave processes such as diffraction, refraction, shoaling, wave breaking, as well as wave-current interaction. This wave action model is coupled with a two-dimensional hydrodynamic model, the feedback effect of wave and current can be effectively simulated. This model is validated by simulating wave laboratory experiments in an inlet entrance, and waves and tidal currents in Grays Harbor, WA by using available field observation data in 1999. The capabilities of the wave model for simulating wave-current interaction and the corresponding breaking effect are confirmed in the study.


INTRODUCTION
When waves encounter currents in tidal inlets, river mouths, or nearshore zones, the properties of wave dynamics will be changed according to speed and direction of the currents (Figure 1).Strong currents, e.g., due to tide, storm surge, and river flood, interacting with waves in coastal and estuarine areas may have a drastic influence on the flooding/inundation state over the impact regions.Particularly, waves traveling in the area of strong opposite currents may become potentially hazardous for navigation as waves can be getting short and steep.Wave blocking and breaking may occur simultaneously in the case of strong opposite currents against waves.Strong interaction between currents and waves can also exacerbate morphological changes in navigation channels, and creates unusual scour holes in inlet channels and extreme shallow shoaling outside of inlets.It has been well known that the Doppler frequency shift occurs when waves take a ride on an ambient current.Early efforts on the research of wave-current interaction (e.g., the Doppler shift effect) were made in terms of laboratory experiments and analytical mathematical approaches with simplified assumptions such as steady ambient currents, monochromatic incident waves, constant water depth, small amplitude waves etc. (e.g., Uana 1942, Mei 1983).Those simplified theories and empirical relationships are unable to be applied to calculate wave fields under complex current conditions in coasts driven by complex hydrodynamic processes such as tides and river inflows.
With the advent of numerical modeling by means of digital computers, temporal/spatial variations of multiscale wave characteristics in responses to oceanographic, meteorological, and hydrological conditions can be computed by sophisticated numerical wave models.Basically, two types of wave models are commonly employed for simulations of wave fields with/without wave-current interaction: phase-resolving models and phase-averaged models.The former enables to produce variations of wave shapes when traveling in a varying bathymetry and an ambient current.An example for this type of model is the Boussinesq equations for computing temporal/spatial distributions of monochromatic/random waves.For instance, Yoon and Liu (1989) developed two-dimensional (2-D) Boussinesq equations to study weakly nonlinear water wave.Chen et al. (1998) developed a onedimensional (1-D) Boussinesq equation to investigate wave deformations, including wave blocking, in a strong steady ambient current.Nwogu and Demirbilek (2001) developed a 2-D numerical wave model, BOUSS-2D, for simulating wave dynamic processes in coastal regions and harbors.The elliptic mild-slope wave equation has also been used for studying wave-current interaction (e.g.Chen et al. 2005).Phase-resolving wave models are better suited to domains with complex bathymetric and geometric features, where the effects of wave diffraction and reflection can be important.The drawback in applications of phase-resolving models, however, is the expensive computing, due to the requirements of resolving details of wave shapes (i.e.phases).It limits the applications of those wave models within relative small scale problems.
Opposite to the phase-resolving model, a phase-averaged model is based on a wave spectral equation.Because of its efficient computing performance, phase-averaged models are suited to large scale applications of random wave growth and transformation.Examples of wave spectral models are SWAN (Booij et al. 1999), STWAVE (Resio 1993), CCHE2D-Coast (Ding and Wang 2008), and CMS-Wave (Lin et al. 2007).Zheng et al. (2008) have applied a wave-action model developed by Mase et al. (2005) to compute wave fields with the effect of wave-current interaction, and to evaluate several wave breaking models.These frequency-domain models are suited to multidirectional wave transformation over large areas of open oceans but have in recent years become increasingly popular for modeling nearshore waves.The present study focuses on this type of wave model to develop a multidirectional wave-action model coupling with a hydrodynamic model.
This study presents an integrated numerical modeling system to simulate wave deformation/transformation in tidal inlets or river mouths with ambient currents by carefully modeling the effect of wave-current interaction.A multidirectional wave-action balance equation is used to compute wave processes such as diffraction, refraction, shoaling, breaking, as well as wave-current interaction.This wave-action model is integrated with a two-dimensional hydrodynamic model, CCHE2D-Coast (Ding and Wang 2008), so that the coupling simulation of wave and current can be easily performed.The model's capabilities for simulating wave-current interaction and the corresponding breaking effect are confirmed in the model by simulating wave fields in a laboratory experimental wave flume and a real tidal inlet, Grays Harbor in the pacific coast of Washington state.

MODEL DESCRIPTION
In case of coastal flow processes in tidal inlets and river mouths, the effect of ambient currents on wave propagation cannot be neglected.Longuet-Higgins and Stewart (1961) introduced the concept of radiation stresses and showed the existence of energy transformation between waves and currents.Therefore, the coastal flow processes can be simulated considering wave and current models either in a coupled or decoupled way.In the present study two different models are used to simulate waves and currents.Both current and wave models are available in CCHE2D-Coast (Ding and Wang 2008), which is an integrated modeling system consisting of multidirectional wave-action model, hydrodynamic model, and morphodynamic model.Based on a flexible non-orthogonal mesh, this modeling system is capable of simulating combined riverine/coastal/estuarine/ocean processes driven by tides, waves, river inflows, winds, surges, and turbulence flows in a large-scale computational domain with irregular coastlines.The studies on wave-current interaction reported in this paper are done by only operating wave and current models.A detailed description of the CCHE-Coast wave model and hydrodynamic model is presented below.

Spectral Wave Model
The spectral wave model is governed by a multidirectional wave action balance equation to compute variations of wave-action density in time, space, wave directions, and frequency.The model formulation is based on the parabolic approximation equation including diffraction terms and energy dissipation due to wave breaking (Mase et al. 2005).The model can simulate unsteady/steady state spectral transformation of directional random waves, and waves can propagate from the seaward boundary toward shore and vice versa.The model takes into account the effect of an ambient horizontal current or wave-current interaction and solves the wave-action balance equation as follows: where t = time, x, y = two coordinates in two horizontal directions,   (∂/∂x, ∂/∂y) the gradient operator in the x-y plane, θ = wave angle relative to the positive x-direction, C, C g = wave celerity and group velocity, respectively; The first term on the left-hand side of (1) represents the local rate of change of action density N=N(x,y,σ,θ,t) in time, where N is the ratio E(x,y,σ,θ,t)/σ, σ is wave angular frequency (or intrinsic frequency), the spectral wave density E represents the wave energy per unit water surface area per frequency interval.The second term represents propagation of wave action density in a horizontal x-y plane (with propagation velocities c).The third term represents depthinduced and current-induced refraction (with propagation velocity c θ in θ space).The expressions for these propagation speeds are given by the linear wave theory [e.g., Holthuijsen et al. 1989], i.e., where u =depth-averaged velocity vector; k = wave number; i  = (cos, sin), a unit vector following the wave direction; h = water depth.The first term on the right-hand side of (1), introduced by Mase (2001), represents the energy dissipation due to the diffraction effect in the alongshore y-direction, which is implicitly perpendicular to wave direction;  is an empirical coefficient.Mase (2001) suggested this empirical coefficient has a possible value within a range of 2.03.0.Ding et al. (2006) found that a wider range of the  value can be used, which is dependent on problems with laboratory and field scales.To specify an incident wave spectrum in the offshore, the TMA spectrum (Bouws et al. 1985) and the Bretschneider-Mitsuyasu (B-M) spectrum (Mitsuyasu 1970) can be selected in this model.The second term in the right-hand side represents wave energy loss due to wave breaking, where ε b is a parameter for wave breaking energy dissipation.The last term in (1), Q, is source terms of wave energy due to wind forcing, bottom friction loss, nonlinear wave-wave interaction term, etc.Under the circumstance of co-existence of wave and current, wave frequency is changed with water depth and velocities.According to small-amplitude wave theory, the resulting wave number can be calculated by the following dispersion equation with the Doppler frequency shift: where g = the gravity acceleration; ω = incident (or absolute) wave angular frequency.Two mathematical roots of the equations are illustrated in Figure 2, of which the values vary with the relationship between wave and current.By excluding one unphysical root causing a negative angular frequency and a root with large wave number, it is readily found that there are three real roots for solutions of wave numbers: , wave and current traveling in the same direction;

Figure 2. Solutions of nonlinear dispersion equation
If opposite current is strong enough, no solution for wave number exists for the dispersion relation by the linear wave theory.Physically, this effect is associated with wave blocking phenomenon, which prevents wave energy from traveling upstream, and trigger early wave breaking at the offshore of inlets, and causes reflection and energy dissipation.In simulations, the group velocity goes to zero in the case of wave blocked, and the wave action density is set to zero for the corresponding frequencydirection bin.In general, the solutions of wave number under the coexistence of wave and current are obtained iteratively from ( 4) and ( 5).Specifically, in the deep water, e.g., kh>3π, assuming that 1 tanh  hk , the wave number can be obtained directly from a simplified quadratic equation (Zheng et al. 2008).
According to Takayama et al. (1991), the wave energy dissipation due to wave breaking is calculated by assuming that (1) the probability distribution of breaking wave height can be represented by the Rayleigh distribution (Thornton and Guza 1983); (2) over a computational cell, the local bathymetry can be approximated as shore-parallel contours.Then, the local rate of wave breaking energy dissipation ε b in a computational cell can be calculated as the time-averaged energy loss over a time that an individual wave travels through a cell from the seaward side to the inland side, i.e. x y where Δx, Δy = the grid size in x-and y-directions, respectively.Based on the Rayleigh probability distribution of wave energy, the coefficient b   can be obtained as follows: where H bo , H bi = breaking wave heights at the offshore side and the onshore side of a computational cell, respectively; H 1/3 = the local significant wave height; α = the ratio of the significant wave height to the mean wave height, and α = 1.60 for the Rayleigh distribution of wave heights.The breaking wave height is calculated based on the breaking wave criterion proposed by Goda (1970), or the one extended by Sakai et al. (1988).
The source term Q in (1) includes the wave energy dissipation due to bottom friction and the wave energy input induced by winds.The energy loss by bottom friction is calculated by a drag law model (Collins 1972).The energy input by wind forcing can be modeled as separated sink and source terms (Lin andLin 2004a and2004b).
Boundary conditions for simulating wave deformation over a computational domain covering ocean, coast and estuary are comprised of offshore wave spectra determined by wave heights, periods, and directions.The present wave spectral model supports two kinds of offshore (deepwater) wave spectrum inputs, i.e., the TMA spectrum (Bouws et al. 1985) and the Bretschneider-Mitsuyasu (B-M) spectrum (Mitsuyasu 1970).
In CCHE2D-Coast, a general non-orthogonal quadrilateral computational grid, the wave action balance equation ( 1) is numerically discretized in a geophysical domain by means of finite difference schemes which taking into account non-uniform grid size and general quadrilateral mesh shape.In addition, the wave spectra are discretized into a number of frequency bins, based on the equal energy dividend, by which each frequency bin represents an individual wave.The bins for wave directions are also discretized to cover a half-plane wave direction (θ) domain from +π/2 to -π/2.A first-order upwinding finite difference scheme is applied to discretize the second, third, and fourth terms which represent the propagation of wave action in the horizontal plane and wave refraction due to varying water depth.The central difference scheme is applied to discretize the first term in the right-hand side for the diffraction term.A semi-implicit treatment is applied to the source terms for wave breaking and bottom friction in order to increase the stability of the numerical model.Finally, the discretized wave action balance equation is solved by means of the parabolic approximation, in which the waves are assumed to have a principal propagation direction from offshore toward onshore.An iterative solver, the Gaussian-Seidel algorithm, is used to solve the linear algebraic equations of the wave action balance equation in every y-θ plane arranged from offshore to onshore.The computed wave actions in all the frequency bins add up to the total wave energy, which can be used to calculate statistical wave parameters such as wave heights, mean wave directions, and peak/mean periods.
The radiation stresses R also can be calculated by the numerical integration of the following expressions: where R xx , R xy , R yx , and R yy are the four components of the radiation stresses R, and n =0.5+kh/sinh2kh.These radiation stresses are used to calculate wave stresses which in turn are used in the momentum equations of hydrodynamic model described in the following section.

The Hydrodynamic Model
The hydrodynamic model contains the depth−averaged 2−D continuity and momentum equations to simulate flows driven by hydrological forcing such as tides, waves, river inflows, surface winds, and turbulence flows in a large-scale coastal and estuarine region, namely where  = water elevation relative to mean sea level (MSL); u =(U,V), depth-averaged velocity vector with two components in the horizontal coordinates; P s = atmospheric pressure on the free surface, which accounts for the variation of air pressure during a storm;  = water density;  t = the depthaveraged Reynolds stress;  S = wind shear stress;  b = bottom friction stress; f cor =(fV, -fU), the Coriolis force term; f = 2Ωsinφ is the Coriolis parameter, wherein Ω = angular frequency of earth rotation and φ = latitude of the study area.
For random waves, the radiation stresses ( 8), ( 9), and ( 10) can be added into the momentum equations ( 12) for simulations of wave-induced nearshore currents.However, the drawback of these expressions for calculations of radiation stresses is lack of consideration of complex nearshore flow structure such as undertow and difference in volume flux between offshore and onshore due to nonsinusoidal wave shapes.It is known that these sinusoidal formulations could not generate accurately currents inside surf zone when especially wave breaking occurs (Svendsen et al., 2003).Therefore, based on a concept for surface roller due to breaking waves, the following expression of the radiation stresses proposed by Svendsen (1984) represents the net (shortwave-averaged) force that the shortwave exert on a water column is defined as: where I = the identity matrix; Q w = the wave volume flux induced by the short wave motion; the tensor e is The scalar S m and S p , as well as the wave volume flux are calculated respectively according to the different formulations suitable for the region inside and outside the surf zone (Table 1).

Inside surf zone
Outside surf zone Note: H is wave height; A is the surface roller area, =0.9H 2 ; L is the wavelength.
Most existing hydrodynamic models only adopt the radiation stresses derived from the sinusoidal wave theory to calculate the wave-induced forcing terms over the entire coastal domain.The hydrodynamic model in CCHE2D-Coast provides several options for calculation of the radiation stresses (Ding et al. 2006).Those options are based on (1) the radiation stresses by the integration of wave energy in (8), (9), and (10), (2) the conventional radiation stresses driven by sinusoidal waves, and (3) the radiation stresses with consideration of non-sinusoidal waves and surface roller assumption.
For the third option of the radiation stresses, this hydrodynamic model needs to identify the coastal surf zone, and then uses the non-sinusoidal radiation stresses for the region inside the surf zone, and the sinusoidal wave formulation for the region outside the surf zone (deep water region), respectively.The significant wave height can be used as the representative wave height H for the case of irregular wave incidence.Table 1 summarizes the different radiation stress formulae inside and outside the surf zone.Because the radiation stresses used for the surf zone take into account the vertical variations of wave breaker structures, some 3D features of the cross-shore movement mechanisms, e.g.undertow and mass flux, are reflected accordingly in the hydrodynamic model (Ding et al. 2006).
By taking into account the influence of combined wave and current on the bottom friction stress  b , the friction law of the combined wave and current proposed by Tanaka and Thu (1994) is used to estimate the friction coefficient, namely ( ) where the overbar means the time-averaged integration over a shortwave period; C f = friction coefficient; u b =the wave orbital velocity at the bottom.The friction law of the combined wave and current is used to estimate the friction coefficient in the different flow regimes, including rough turbulent flow, smooth turbulent flow, and laminar flow.The combined friction coefficient f cw (=2C f ) is given as follows, where f c and f w are the friction coefficients due to current and wave, respectively;  is the coefficient due to nonlinear interaction of waves and currents;  is the angle between wave orthogonal and current vector.Surface wind stresses are given by the following empirical wind shear stress formulations where C d = wind drag coefficient, ρ a = density of air, W= wind speed vector measured at 10.0-m above the sea level.
As far as the depth-averaged Reynolds stress  t in Eq. ( 4) is concerned, the present hydrodynamic model provided users with two eddy viscosity turbulence models pertaining to the characteristics of waves (Ding et al. 2006), and a k-ε model (Jia et al. 2002 ).

Numerical Approaches for Solving Hydrodynamic Equations
This integrated coastal process model (CCHE2D-Coast) has been built in a well-established numerical software package called CCHE2D (Jia et al. 2002), which is an extensively verified and validated tool to analyze 2-D shallow water flows, sediment transport, and water quality, with natural flow boundary conditions.Similar to the CCHE2D hydrodynamic model, the wave and hydrodynamic models in CCHE2D-Coast were discretized in a non-orthogonal grid system so that the models have more flexibility to simulate physical variables in complex coastal zones with irregular coastlines.A time-marching algorithm proposed by Jia et al. (2002) was used to compute tidal and wave-induced currents.A validated algorithm in CCHE2D for the treatment of wetting and drying processes was directly used for predicting the tidal flat variations and coastal inundations.This coastal process model has been extensively validated by simulating waves, wave-induced currents, and morphological changes in coastal applications in various laboratory and field scales (e.g., Ding et al. 2004, Ding and Wang 2005, Ding et al. 2006, Ding et al. 2008).

Boundary Conditions
In general, a typical computational domain in a coastal and estuarine region is surrounded by four types of boundaries: shoreline, offshore boundary, two open cross-shore boundaries, and river inflows at upstream.Inside the domain, island shorelines or offshore structure boundaries may exist.In order to compute wave fields in the computational domain, by means of the parabolic approximation for wave propagation, incident wave parameters (i.e.wave heights, directions, and periods) are needed to specify the incident wave spectra on the offshore boundary.The TMA spectrum (Bouws et al. 1985) and the Bretschneider-Mitsuyasu (B-M) spectrum (Mitsuyasu 1970) are two choices provided by the wave model.In simulation of hydrodynamic processes, the known values of velocities or discharges can be imposed on the corresponding cross-shore boundaries.The impermeable condition of currents is used on shorelines.The known tidal elevations are specified at the offshore boundary to generate tidal flows.For the initial conditions for velocities and water elevations, the cold start (starting from a static state) is generally utilized to initiate the simulations of the tidal and wave-induced currents.

MODEL VALIDATION RESULTS AND ANALYSES
Two validation cases are presented here: one is simulations of wave breaking on steady ebb currents in a laboratory idealized inlet, another wave fields in tidal currents in the Grays Harbor in the pacific coast of Washington state.

Wave Breaking on a Current at an Idealized Inlet
To validate the wave model for simulating interaction of wave and current, a number of validations cases have been performed by simulating wave parameters in laboratory experiments conducted by Smith et al. (1998).A wave basin built in U.S. Army Corps of Engineers Waterways Experiment Station (WES) is a 150-ft (46-m) wide, 325-ft (99-m) long, and 2-ft (0.6-m) high concrete basin.The seabed parallel contours were determined by using an equilibrium profile equation from Dean (1977).
The basin floor at the offshore is 30.4 cm.The inlet throat region converges to a depth of 15.2 cm relative to a Mean Low Water (MLW) datum.The minimum width of the inlet is 244 cm across the inlet between the MLW contours.Two parallel jetties were installed in the entrance of the inlet channel which have a spacing of 366 cm and extended 550 cm offshore.A storage tank behind the wave basin supported waters into the basin to establish a current circulation.Smith et al. (1998) reported that there were 12 experiments (or runs) conducted in the basin with varying steady ebb currents and offshore incident waves.In this study, only three runs are selected to validate the present numerical model for simulating wave-current interaction.The incident wave parameters and the ebb current speed are listed in Table 2, in which H so = offshore incident significant wave height, T p = peak wave period, and U = velocity of ebb current.The offshore incident waves are generated with a TMA spectrum from using a gamma (a spectral peak enhancement factor) value of 3.3.Different from the approach by using the interpolated flow velocities in Zheng et al (2008), the flow fields over the entire computational domain were computed by the hydrodynamic model by using the steady flow velocity U as the boundary condition.The hydrodynamic model was validated by comparing computed velocities with the observed ones.Therefore, a calibrated flow field is able to provide the wave model with accurate input current data.As an example, Figure 3 shows the computed currents in Run 9 and Run 11, in which both cases share the same input currents because of the same flow boundary conditions.It indicates that the ebb currents are concentrated within the inlet channel.For validation of the flow model, in Figure 4, the profiles of crossshore velocity component along the three transects, of which the locations are shown in Figure 3, are compared with those of observations.The computed velocities are in good agreement with the measurements.For the other two runs, the computed currents have been also compared with those observed velocities.Similar agreement has been obtained.After these successful validations of the hydrodynamic model in all the three cases, the computed flow fields (i.e., velocities and water elevations/depths) are stored as input data for simulating the wave fields by considering the wavecurrent interaction.The numerically reconstructed flows have insured the best input of velocity and the good accuracy of wave field generated.
In the simulations of wave fields, the incident waves attacked the inlet propagating from the leftside boundary with the direction perpendicular to the shoreline.Through a few model tests, it was found that the extended Goda's wave breaking criterion (Sakai et al. 1988) produced better wave results than those by the original Goda's wave breaking formulation.The contour lines and the flooding color in Figure 3 represent the spatial distribution of the computed significant wave heights in Run 9 and Run 11.Due to the ebb current against the offshore wave propagation, the wave heights increased at the outside of the inlet channel.Comparing with the wave field in Run 9, the waves in Run 11 are quickly dissipated when propagating into the inlet channel, because of the smaller peak period and shorter wave length in Run 11.Due to strong opposite currents, wave blocking happened in the upstream of the inlet channel.
Detailed comparisons of wave heights in the three transects along the inlet channel for Run 5 are shown in Figure 5.As plotted in Figure 5(a), excluding the effect of wave-current interaction, the wave spectral model underestimated the wave heights in all the three sections.With the effect of wavecurrent interaction, the results of computed wave heights shown in Figure 5(b) become much closer to the observed ones.Only the wave heights along the middle line are overpredicted.For the other two cases, as shown in Figure 6 (a) and (b), the computed wave heights are in good agreement with their observations.Especially, the wave heights computed in Run 11 are reproduced very well, in which the wave blocking occurred in a wide range from the entrance of the inlet to the right boundary.This wave model has been also validated carefully for simulating waves and tidal currents at Grays Harbor which is an estuarine bay located 45 miles (72 km) north of the mouth of the Columbia River, on the southwest Pacific coast of Washington state, in the United States (Figure 7).The bay is 17 miles (27 km) long and 12 miles (19 km) wide.Two jetties are installed in the entrance of the inlet.A computational mesh with 211×271 grids covers the bay, the coast, the two jetties, and 10-km wide ocean at the offshore.The simulation period is from Sept. 19 -Oct. 15, 1999, as   First, the hydrodynamic model was validated by simulating the tidal flows during this computational period in 1999.The tidal elevations as the offshore boundary conditions are formed by typical spring-neap tides, in which the tidal range reached up to 3.0 m. Figure 8 is a snapshot of computed ebb currents.In the case, the discharge from the upstream river (i.e. the Chehalis River) was neglected.Converging from the north, south, and east bay, the ebb current at the inlet become the strongest flow in the bay. Figure 9 presents comparisons of the tidal elevations at Stations Tide #1 and Tide #3, which are located respectively at the south bank of the inlet and the North Bay (Figure 7).The computed tidal levels are in excellent agreement with the observations at the tidal gauge stations.Moreover, the tidal currents have been compared with the observed velocities, at the wave gauges stations, in which the tidal velocities have been monitored during the observation campaign in September of 1999 conducted by USACE-ERDC-CHL.The overall comparisons in all the stations showed that the computed tidal currents have been obtained with a high accuracy.As an example, Figure 10 gives the comparisons of two velocity components in the east-west direction and the northsouth direction between the computations and observations.It shows that both velocities magnitudes and phases were computed accurately.After the successful validation on the hydrodynamic model, the computations of waves and currents have been carried out by taking into account wave-current interaction.It means that the wave field is updated once over an interaction interval in the computations of the currents.Through several test runs, it was found that a one-hour interaction interval could produce enough good wave results; and more frequent wave-current interaction only made the computations expensive.
The wave parameters, including significant wave heights, mean direction, and peak periods, measured at the offshore at the CDIP03601 buoy were used as incident wave boundary conditions.TMA spectra are assumed to give the spectral shapes at the offshore, with gamma values varying from 3.3 to 4.0 during the computation period.A snapshot of wave heights and mean directions is shown in Figure 11.Comparisons of time series of wave parameters during the computational period from Sept. 19 -Oct.15, 1999 at St. 1 are plotted in Figure 12.The computed significant wave heights, mean directions, and peak periods are in good agreement with the observations at this wave gauge station.Figure 13 shows intercomparisons of wave heights at St. 4 located in the middle of the inlet close to the north bank, in which the observed wave heights are obtained from the USACE-ERDC-CHL Dataport (CHL 2009).The results clearly indicate that the wave heights computed with the inclusion of wave-current interaction follow the tidal cycles very well; without the current effect the computed wave heights could not catch this tidal-wave periodical phenomenon.So exclusion of the current effect may cause underestimation in the wave heights in ebb tides, and overprediction in flood tides.It is therefore necessary to include wave-current interaction in the simulations in order to improve prediction accuracy of waves traveling with ambient currents.

CONCLUSIONS
This study presents an integrated numerical model to simulate wave deformation/transformation in tidal inlets or river mouths with ambient currents (e.g.tidal currents, river inflows) by carefully modeling the effect of wave-current interaction.A multidirectional wave-action balance equation is used to simulate temporal and spatial variations of wave spectra driven by wave processes such as refraction, diffraction, shoaling, wave breaking, bottom friction, winds, as well as the Doppler frequency shift effect of waves generated by ambient currents.This wave action model is integrated with a two-dimensional hydrodynamic model, CCHE2D, which is capable of simulating multiple-scale flows existing in coasts and estuaries, such as tidal currents, nearshore currents, storm surge flooding flows, river flows, turbulence flows, driven by various hydrological forcing such as tides, waves, atmospheric pressures, river inflows, winds, turbulence, and the Colioris force.Therefore, this integrated coastal model has capabilities in simulating wave-current interaction by different physical processes with multiple-scale flow motions.
This integrated model is carefully validated by simulating the wave-current interaction at an inlet in a laboratory wave basin driven by steady ebb current which can be viewed as steady river inflows.Before the computation of wave field, the flow field over the whole computational domain was reconstructed by the hydrodynamic model, and the computed velocities were compared with the measurements.Using the calibrated and reconstructed flow field as input currents, the simulations of three experimental cases reproduced the wave fields in which the effect of the waves interacting with the steady currents is well captured.A detailed comparison of wave heights along three transects between simulations and measurements shows that the computed results are in good agreement with the laboratory observations.
To confirm the model's capabilities in simulating wave fields in a field scale, this wave model was further validated by reproducing the waves and tidal currents in Grays Harbor, WA over a monthlong period starting on September 19, 1999.Using available field observation data for tidal elevations, currents, wave parameters in the period provided in the USACE-ERDC-CHL Dataport (CHL Data 2009), the wave and hydrodynamic modules in the integrated coastal model were validated throughout.First, the validation of the hydrodynamic model was achieved by reproducing the tidal currents over the monthlong period.And then the wave model was run by interacting with the hydrodynamic model in which the wave forcing was included.In comparison with the observed wave data, it is indicated that the wave-action model has well captured the periodical variations of wave heights at the stations in the inlet channel, and this fluctuation effect in waves clearly follows tidal cycles in the bay.In contrast, without consideration of wave-current interaction, the integrated could not catch this long periodical phenomenon in waves in the inlet; but it caused underestimation in the wave heights in ebb tides, and overprediction in flood tides.
This study concludes that it is important to include wave-current interaction in the simulations in order to improve the accuracy of predicted waves in the case of ambient currents (e.g.tidal flows and river inflows).It also emphasizes that a careful model validation in various spatial/temporal scales is the only way to ensure the high quality of numerical models in predicting waves and currents in coasts and estuaries.In the future, this integrated coastal model called CCHE2D-Coast will be applied to investigate morphodynamic response to waves and currents in inlets and river mouths.

Figure 1 .
Figure 1.Wave-Current Interaction in a River Mouth: Sketch of Coastal Processes

Figure 3 .
Figure 3. Computed currents and wave heights in Run 11

Figure 4 .
Figure 4. Comparisons of cross-shore velocities in three transects

Figure 5 .Figure 6 .
Figure 5. Comparisons of wave heights (Run 5) the observations over this period have been conducted by U.S. Army Corps of Engineers (USACE), Engineering Research and Development Center (ERDC)'s Coastal & Hydraulics Laboratory (CHL), and the data can be assessed in their website at http://sandbar.wes.army.mil/.

Figure 10 .
Figure 8. Computed ebb current at Grays Harbor in September 1999

Figure 11 .-
Figure 11.Computed significant wave heights and mean directions in September, 1999

Figure 12 .
Figure 12.Computed wave parameters (significant wave heights, mean direction, and peak period) at St. 1

Figure 13 .
Figure 13.Comparison of significant wave heights at St. 4