APPLICATION OF A TVD SOLVER IN A SUITE OF COASTAL ENGINEERING MODELS

This paper describes the recent developments in a suite of coastal engineering models using Godunov-type shockcapturing schemes. The developments include a depth-integrated, wave resolving Boussinesq model, a hydrostatic, wave-averaged circulation model, and a fully 3-D non-hydrostatic model in a surface-following σ coordinate formulation. The models implemented with the shock-capturing TVD schemes show robust performances in modeling breaking waves, nearshore circulation and coastal inundation. In this paper, we present model equations in a conservative form, MUSCLE-TVD numerical scheme and model applications. We also point out some problems caused by the TVD scheme in the recent model applications.


INTRODUCTION
Historically, most models for nearshore and coastal processes have been formulated numerically using finite difference methods.While the use of such schemes has led to both skillful models and a better understanding of nearshore processes, finite difference schemes can exhibit undesirable properties.Processes which develop discontinuities (e.g.surfzone wave breaking in wave-resolving models, or tidal bore formation in models for tidal circulation) are typically not well handled without the additional use of filtering (e.g., Kennedy et al., 2000).Additionally, demands on moving boundary algorithms have increased dramatically with the increase in focus on storm surge and tsunami inundation as well as dynamical treatment of salt marsh environments.
In view of these limitations, we have recently re-developed a suite of models using finite volume schemes of TVD (Total Variation Diminishing) form.Finite volume schemes provide a much more natural mean for treating the moving shoreline problem, and provide a robust, mass conserving treatment for shorelines translating over a wide range of grid points during tidal or individual wave excursion events.In addition, the TVD scheme allows for robust treatment of solution discontinuities through the shockcapturing mechanism.The recent publications related to model developments using the TVD-type scheme can be referred to Bradford (2005Bradford ( , 2011)), Erduran et al. (2005), Zijlema andStelling (2005, 2008), Shiach and Mingham (2009), Tonelli and Petti (2009, 2010, 2012), Roeber et al. (2010Roeber et al. ( , 2012)), Shi et al., (2012b), Ma et al. (2012), Tissier et al., (2012) and others.
Three models are described here, including a depth-integrated, wave resolving Boussinesq model, a hydrostatic, wave-averaged circulation model in 2D or quasi-3D form, and a fully 3-D non-hydrostatic model in a surface-following σ coordinate formulation.The models implemented with the shock-capturing TVD scheme show robust performances in modeling breaking waves, nearshore circulation and coastal inundation.In this paper, we present the conservative form of theoretical formulations for the three models, a multi-order MUSCLE-TVD scheme and a high-order adaptive time-stepping scheme.Model applications to large-scale wave simulation, wave-current interaction, and nearshore circulation at a complex inlet system are presented.

MODEL EQUATIONS FUNWAVE-TVD
FUNWAVE-TVD (Fully Nonlinear Boussinesq Model with a TVD solver) is developed based on the fully nonlinear Boussinesq equations of Chen (2006), extended to incorporate a moving reference elevation following Kennedy et al. (2001) and implemented using a TVD-type solver (Toro, 2009).A conservative form of the equations is derived in order to use a hybrid numerical method combining the finite-volume and finite-difference TVD-type schemes.The generalized conservative form of Boussinesq equations can be written as where Ψ and Θ(Ψ) are the vector of conserved variables and the flux vector function, respectively.Note that this generalized conservative form is also used in the NearCoM-TVD and NHWAVE equations described in the following sections.For the Boussinesq equations in this study, the conserved variables and flux function are given by In (2) P and Q are defined by (Pi, where u α denotes the velocity at a reference elevation z = z α , H is the total water depth, i.e., H = η + h, and ū2 is the depth averaged O(µ 2 ) contribution to the horizontal velocity field, given by in which in which In (3), where (U 1 , V 1 ), (U 1 , V 1 ), (U 2 , V 2 ), (U 3 , V 3 ) and (U 4 , V 4 ) are additional dispersive terms and can be found in Appendix A of Shi et al. (2012b).R x , R y in (3) represent the bottom stresses approximated using a quadratic friction equation.Wave breaking is modeled by switching the Boussinesq equations to NSWE as the ratio of surface elevation to water depth exceeds a certain threshold following the approach of Tonelli and Petti (2009).For modeling wave-generated nearshore currents such as rip current and alongshore current (Chen et al., 1999(Chen et al., , 2003)), a Smagorinsky (1963)-like subgrid turbulent mixing algorithm is implemented.The eddy viscosity associated with the subgrid mixing is determined by breaking-induced current field.The detailed formulations for the subgrid mixing algorithm can be found in Chen et al. (1999).

NearCoM-TVD
NearCoM-TVD is a TVD version of the Nearshore Community Model (Shi et al., 2005(Shi et al., , 2012a)).This version integrates the wave model SWAN and a modified SHORECIRC (Shi et al., 2011a) and Soulsby's (1997) sediment transport formula in the NearCoM framework.The TVD solver is implemented for solving SHORECIRC equations.In order to obtain a conservative form of governing equations required by the hybrid numerical scheme, SHORECIRC equations are re-derived using the combined contravariant and Cartesian components of a vector in the momentum equations (Shi and Sun, 1995).
The conservative form of the governing equations can be written as the same compact form as (1).The vectors of conserved variables and the flux vevtor function are written as The above equations are described in generalized curvilinear coordinates.We assume that the coordinate transformation is performed between the Cartesian coordinates x α and generalized curvilinear coordinates ξ α .The contravariant components of velocity vector can be expressed by where u β is the Cartesian component of velocity vector and In ( 11)-( 13), η is surface elevation, (u, v) represent Cartesian components of velocity.J is the Jacobian value in curvilinear coordinates, (P, Q) = (Hu 1 , Hu 2 ), the contravariant components of flow flux, (τ b x , τ b y ), (τ s x , τ s y ), ( f v, − f u), and (ROT x , ROT y ) represent the bottom friction, wind stress, Coriolis force and the rest of terms associated with the 3D dispersion effect (Putrevu and Svendsen, 1999) in x and y directions, respectively.NHWAVE NHWAVE (Non-Hydrostatic Wave Model) is a 3D non-hydrostatic wave model with a moving bottom.The governing equations are the incompressible Navier-Stokes equations in conservative form in σ coordinates.The σ coordinate is defined as where D is the total depth, i.e., D(x, y, t) = h(x, y, t) + η(x, y, t).Note that h is function of time so that the bottom movement is taken as particular application to tsunami generation by three-dimensional underwater landslides.
The compact form of the governing equations can be written in the same form as in (1).The conserved variables and numerical fluxes can be described as where (u, v) represent velocity components in (x, y) plane and ω is the vertical velocity defined in the σ coordinate image domain.The source term includes three components representing bottom slope, pressure gradient and turbulent mixing, respectively. where p here represents dynamic pressure only.Turbulent diffusion terms S τ x , S τ y , S τ z are calculated using the k − closure model.The governing equation for free surface can be written as The dynamic pressure is calculated by solving the Poisson equation in (x, y, σ) coordinate system as described in the next section.

NUMERICAL SCHEMES
A combined finite-volume and finite-difference method was applied to the three models for the spatial discretization.Following Toro (2009), two basic steps are needed to achieve the spatial scheme.The first step is to use a reconstruction technique to compute values at the cell interfaces.The second step is to use a local Riemann solver to get numerical fluxes at the cell interfaces.
The high-order MUSCL-TVD scheme can be written in a compact form including different orders of accuracy from the second to the fourth-order, according to Erduran et al. (2005) who modified Yamamoto et al.'s (1998) fourth-order approach.In x-direction, for example, the combined form of the interface construction can be written as follows: where φ L i+1/2 is the constructed value at the left-hand side of the interface i + 1 2 and φ R i−1/2 is the value at the right-hand side of the interface i − 1 2 .The values of ∆ * φ are evaluated as follows: In ( 23), minmod represents the Minmod limiter and is given by κ 1 and κ 2 in ( 21) -( 23) are control parameters for orders of the scheme in a compact form.The complete form with (κ 1 , κ 2 ) = (1/3, 1) is the fourth-order scheme given by Yamamoto et al. (1998).(κ 1 , κ 2 ) = (1/3, 0) yields a third-order scheme, while the second-order scheme can be retrieved using (κ 1 , κ 2 ) = (−1, 0).For the Boussinesq model, a higher order such as the third-order or fourth-order should be used for first-order derivatives in order to move leading order truncation errors to one order higher than dispersive terms.χ(r) in ( 21) and ( 22) is the limiter function.The original scheme introduced by Yamamoto et al. (1998) uses the Minmod limiter as used in (23).Erduran et al. (2005) found that the use of the van-Leer limiter for the third-order scheme gives more accurate results.The van-Leer limiter can be expressed as where The numerical fluxes are computed using a HLL approximate Riemann solver where The wave speeds of the Riemann solver are given by in which u s and ϕ s are estimated as and n is the normalized side vector for a cell face.
A central difference scheme is used for all terms in S.
For time stepping, a third-order Strong Stability-Preserving (SSP) Runge-Kutta scheme for nonlinear spatial discretization (Gottlieb et al., 2001) is implemented in the Boussinesq model: in which Ψ n denotes Ψ at time level n.Ψ (1) and Ψ (2) are values at intermediate stages in the Runge-Kutta integration.
A second-order SSP Runge-Kutta scheme is applied in NearCoM-TVD and NHWAVE: For NHWAVE, the two-step projection method which splits the time integration into a hydrostatic and non-hydrostatic steps is used within the SSP Runge-Kutta scheme.The Poisson equation for dynamic pressure in the σ coordinate system is solved by the preconditioned GMRES scheme.
In all three models, an adaptive time step is chosen, following the Courant-Friedrichs-Lewy (CFL) criterion.A moving shoreline is modeled using the wetting-drying scheme.The models are parallelized using the domain decomposition technique.The Message Passing Interface (MPI) with non-blocking communication is used for data communication between processors.

MODEL APPLICATIONS
To demonstrate the robustness of the TVD solver used in this study, we show examples of long-term and large-scale modeling of nearshore waves and currents.Figure 1 shows a wave simulation using FUNWAVE-TVD in a 2000× 1600 m domain at FRF, Duck, North Carolina.The computational domain includes an extension in y direction for period boundary conditions (both bathymetry and wave) and an extension in offshore direction for wave maker and sponge layer (not shown in Figure 1).The grid size is 2 m in both x and y directions, which can resolve the measurement pier (shown in white squares in Figure 1) at FRF.More applications of FUNWAVE-TVD can be found in Shi et al. (2012b), Tehranirad et al. (2011) and Kirby et al. (2012).
Figure 2 shows a result from NearCoM-TVD in a simulation of waves, tidal currents and residual flows in New River Inlet, NC, where the RIVET experiment was carried out in May, 2012.The computational domain extends approximately 60 km alongshore and 40 km offshore, and 20 km upstream from the river mouth.A curvilinear grid was adopted with the minimum grid size of about 20 m at the inlet region.In the simulation, the JONSWAP wave spectra with the significant wave height of 2.8 m and the peak period of 8 second was used.Tidal forcing is provided at the open boundaries and only M 2 was used in the simulation.Figure 2 shows the tidal averaged residual flows (Eulerian averaging).The model simulations with more realistic forcing conditions are being conducted with extensive model/data comparisons using the RIVET data.Results for the RIVET application will be reported in the near future.
NHWAVE has been validated using several laboratory measurements as presented in Ma et al. (2012).Here, we show a preliminary result from the model comparison to the RCEX field measurement by MacMahan et al. (2010) who measured surfzone waves and circulation in rip current systems.For this simulation, NHWAVE was set up with a horizontal grid resolution of 2m × 2m and the vertical resolution of 3 σ layers.Figure 3 shows nearshore currents obtained by time-and depth-averaging over modeled wave velocities.The red arrows show the wave-averaged velocities at measurement locations.An extensive model/data comparison will be carried out, emphasizing the vertical structure of rip currents using a model setup with more vertical layers.

CONCLUSIONS
A TVD-type solver was implemented in three coastal engineering models, including the fully nonlinear Boussinesq wave model FUNWAVE, the nearshore community model NearCoM and a non-hydrostatic wave model NHWAVE.The models with the TVD scheme have shown robust performance of the shockcapturing method in simulating breaking waves and moving boundaries with good numerical stability.FUNWAVE-TVD and NHWAVE have been validated and documented in several recent publications or technical reports (Shi et al., 2011b, Tehranirad et al., 2011, Shi et al., 2012b, Ma et al., 2012, Tehranirad et al., 2012).NearCoM-TVD is documented in Shi et al. (2012a) and is being validated using theoretical results and field data (Chen et al., 2012).All the three models are open source models and available at http://www.chinacat.coastal.udel.edu/programs/.
Finally, we want to point out that the TVD scheme is not problem-free scheme.Abadie et al. (2012) found that 4th-order MUSCLE-TVD scheme is unstable in the coupling of a 3D Navier-Stokes model and FUNWAVE-TVD in modeling of tsunami waves generated by the flank collapse of the Cumbre Vieja Volcano.The 3rd-order MUSCLE-TVD scheme was suggested in their case.The MUSCLE-TVD scheme with different limiters may have different degrees of numerical diffusivity as mentioned by Bradford and Sanders (2002) and Erduran et al. (2005).In FUNWAVE-TVD applications to modeling SMF tsunami waves, Grilli et al. (2012) pointed out that a simple breaking model, such as the use of Tonelli and Petti (2009) criterion, may not be sufficient for modeling the complex breaking phenomena occurring for the relatively shorter and more dispersive tsunami waves produced by a SMF.In addition, the MUSCLE-TVD scheme with the HLL approximate Riemann solver used in the present study is found to be more diffusive compared with the standard finite difference schemes.A further investigation on the numerical diffusivity caused by the TVD scheme is necessary.

Figure 1 :
Figure 1: Surface elevation at FRF, Duck NC, predicted by FUNWAVE-TVD.Grid resolution: 2m.The measurement pier is resolved by the model and shown in white rectangles.

Figure 3 :
Figure 3: Wave-averaged current field (blue arrows) predicted by NHWAVE and measured wave-averaged current velocities (red arrows).Water depth contours are black lines and locations of in situ instrumentation are marked in text (C1-C5 represent the across-shore array and L1-L5 are the alongshore array.)