Modeling Interactions and Shoaling of Solitary Waves Using a Hybrid Finite Volume and Finite Difference Solver

Modeling the propagation and shoaling transformation of solitary waves in shallow water regions is practically important to the study of impacts of nonlinear waves on coastal environments. Especially, the prediction of wave run up caused by a solitary wave like tsunami plays critical role in the planning of coastal protection. Therefore, studies of solitary waves experimentally and numerically have been performed to characterize the tsunami wave movement towards the coastal regions. In the 1870’s, Boussinesq and Rayleigh described solitary waves mathematically with the so-called Bousinessq equations (BE) [1]. Madsen et al. [2] developed a new form of BE to improve the linear dispersion characteristics of the BE in deeper water. Later, Madsen and Sorensen [3] generalized this new form of the BE for mild sloped bottom in order to model shoaling of waves from deep to shallow water. Introducing the layer-mean velocity potential, Wu [4] developed a set of generalized BE for modeling nonlinear shallow water waves. Nwogu [5] also derived a different form of BE to approximate wave propagation at much deeper environment. The finite difference (FD) schemes have been extensively used to discretize the nonlinear wave equations [2,6-8]. While modeling complex geometries, finite element (FE) method allows the computational meshes matched closely to the flow domain and obtain more accurate results. Zhong and Wang [9] stated that spurious oscillations adversely affect the numerical results when convection dominates in wave propagation. It was known that the Godunov [10] type finite volume (FV) schemes have the superiority in capturing flow discontinuities or oscillations in regular or irregular domains [11-13]. Therefore, application of FV numerical method to discretize BE has attracted researchers’ attention. The combined FVFD schemes were developed [14-17] to discretize the BE, where the FV method was applied to discretize the conservative parts and the FD scheme to discretize the dispersive terms and source terms. Stansby [14] modeled solitary wave run-up cases using FV-FD hybrid scheme. Erduran et al. [15] tested the hybrid scheme by simulating propagation of sinusoidal waves in deep and shallow water and their propagation over a submerged bar. Modeling propagation of regular waves and random waves along a sloping bottom was performed by Shiachet and Mingham [16]. Investigation of wave propagation, shoaling and runup using a finite volume solver of the Boussinesq equations was performed recently by Kazolea et al. [18]. Zhong and Wang [9] pointed out that the researchers had concentrated more on modeling monochromatic deep water waves, but nonlinear shallow-water waves such as solitary waves and their diffraction and refraction in domains with variable bottom topography had not been studied extensively. Therefore, this paper focuses on modeling cases of solitary wave transformation by discretizing the improved form of BE with a developed finite volumefinite difference solver. Nonlinear wave-wave interaction can be tested by simulating the collision of two solitary waves. During the collision of two solitary waves, waves behave like interacting elementary particles such as electrons or protons [19]. Byatt-Smith [20] derived a formulation showing that the peak of the interaction after head-on collision of two equal-amplitude waves is greater than twice of the initial wave amplitude. Yih and Wu [21] derived a corrected version of BE for long water waves and proposed a general solution for the interaction of solitary waves. Hammack et al., [22] investigated both the following and counter propagating (head-on) collisions by numerically solving Euler’s equations. As solitary waves approximately resemble steep waves on beaches, the shoaling and runup of solitary waves on a sloping bottom is a main concern in coastal regions [23]. Shoaling process can be calculated up to the breaking point based on weakly nonlinear wave theories [23]. The bed slope gradient of the cells becomes more important while discretizing flow equations on a computational slope with variable bottom topography. Centered discretization of slope term fails when the waves with sharp wave front propagate on nonflat bed. Zhou et al. [24] suggested the surface gradient method (SGM) to overcome this difficulty. In the SGM, the reconstruction of water


Introduction
Modeling the propagation and shoaling transformation of solitary waves in shallow water regions is practically important to the study of impacts of nonlinear waves on coastal environments. Especially, the prediction of wave run up caused by a solitary wave like tsunami plays critical role in the planning of coastal protection. Therefore, studies of solitary waves experimentally and numerically have been performed to characterize the tsunami wave movement towards the coastal regions. In the 1870's, Boussinesq and Rayleigh described solitary waves mathematically with the so-called Bousinessq equations (BE) [1]. Madsen et al. [2] developed a new form of BE to improve the linear dispersion characteristics of the BE in deeper water. Later, Madsen and Sorensen [3] generalized this new form of the BE for mild sloped bottom in order to model shoaling of waves from deep to shallow water. Introducing the layer-mean velocity potential, Wu [4] developed a set of generalized BE for modeling nonlinear shallow water waves. Nwogu [5] also derived a different form of BE to approximate wave propagation at much deeper environment. The finite difference (FD) schemes have been extensively used to discretize the nonlinear wave equations [2,[6][7][8]. While modeling complex geometries, finite element (FE) method allows the computational meshes matched closely to the flow domain and obtain more accurate results. Zhong and Wang [9] stated that spurious oscillations adversely affect the numerical results when convection dominates in wave propagation. It was known that the Godunov [10] type finite volume (FV) schemes have the superiority in capturing flow discontinuities or oscillations in regular or irregular domains [11][12][13]. Therefore, application of FV numerical method to discretize BE has attracted researchers' attention. The combined FV-FD schemes were developed [14][15][16][17] to discretize the BE, where the FV method was applied to discretize the conservative parts and the FD scheme to discretize the dispersive terms and source terms. Stansby [14] modeled solitary wave run-up cases using FV-FD hybrid scheme. Erduran et al. [15] tested the hybrid scheme by simulating propagation of sinusoidal waves in deep and shallow water and their propagation over a submerged bar. Modeling propagation of regular waves and random waves along a sloping bottom was performed by Shiachet and Mingham [16]. Investigation of wave propagation, shoaling and runup using a finite volume solver of the Boussinesq equations was performed recently by Kazolea et al. [18]. Zhong and Wang [9] pointed out that the researchers had concentrated more on modeling monochromatic deep water waves, but nonlinear shallow-water waves such as solitary waves and their diffraction and refraction in domains with variable bottom topography had not been studied extensively. Therefore, this paper focuses on modeling cases of solitary wave transformation by discretizing the improved form of BE with a developed finite volumefinite difference solver. Nonlinear wave-wave interaction can be tested by simulating the collision of two solitary waves. During the collision of two solitary waves, waves behave like interacting elementary particles such as electrons or protons [19]. Byatt-Smith [20] derived a formulation showing that the peak of the interaction after head-on collision of two equal-amplitude waves is greater than twice of the initial wave amplitude. Yih and Wu [21] derived a corrected version of BE for long water waves and proposed a general solution for the interaction of solitary waves. Hammack et al., [22] investigated both the following and counter propagating (head-on) collisions by numerically solving Euler's equations. As solitary waves approximately resemble steep waves on beaches, the shoaling and runup of solitary waves on a sloping bottom is a main concern in coastal regions [23]. Shoaling process can be calculated up to the breaking point based on weakly nonlinear wave theories [23]. The bed slope gradient of the cells becomes more important while discretizing flow equations on a computational slope with variable bottom topography. Centered discretization of slope term fails when the waves with sharp wave front propagate on nonflat bed. Zhou et al. [24] suggested the surface gradient method (SGM) to overcome this difficulty. In the SGM, the reconstruction of water surface level rather than water depth is implemented to get rid of the spurious waves and preserve the still water condition without creating errors at the interfaces. A combined finite volume and finite difference model by solving the Madsen et al., [3] BE with a second order accurate FV scheme applied to the conservative terms and FD scheme on the dispersive source terms is developed in this study to model propagation of solitary waves, head-on collisions of solitary waves, and solitary wave shoaling over a submerged shoal. The SGM proposed by Zhou et al. [24] is included in the present study to overcome numerical errors during propagation of waves along variable topography. Zhong and Wang [9,25] developed finite element solvers to discretize fully nonlinear weakly dispersive (FNWD) and weakly nonlinear weakly dispersive (WNWD) forms of the BE. The results obtained from the present FV-FD model are compared with those from the FE solvers developed by Zhong and Wang [9,25] and other published solutions.

Governing Equations
The one-dimensional (1-D) Boussinesq equations proposed by Madsen et al., [3] for waves propagating over regions of shallow water depth ( Figure 1) are expressed as where ζ=free-surface elevation measured from the still water surface level, h=water depth measured from the bottom of the channel to the water surface level, u=velocity in x direction, S fx =friction slope in x direction, and D = dispersive term, which is given as where h=H+ ζ , q=uh, and B=dispersion coefficient. As suggested by Madsen and Sorensen [3], B=1/15. Replacing the derivatives of freesurface elevation ζ with the derivatives of water depth h, we have from Eqs. (1) and (2) the revised Boussinesq equations (BE) as

Hybrid Finite Volume and Finite Difference Formulations
The time and spatial discretization using a predictor-corrector type scheme for the 1-D BE on a small section of the computational domain defined in Figure 2 is given as gh(s -s )-D uh u h+ gh 2 To have more accurate solutions, the Roe [26] second-order scheme is used. Different from the first-order scheme using directly the values at cell centers, the conserved variables are first reconstructed at the left and right sides of cell edges using the limiters and surface gradient method. As depicted in Figure 2, the position immediately inside of each cell edge, which is being updated, is named as the left side (L) and outside as the right side (R). The flux vectors at the cell interfaces of cell i, F i±1/ 2 , can be discretized as  Where J = Jacobian matrix and In Eq. (9), we have m eigenvector matrix, r = eigenvectors, and λ = averaged eigenvalues representing the speeds of characteristics. The average states of water depth, velocity, and eigenvectors are computed with the formulas shown below [26] =u -gh and =u + gh λ λ . It should be noted that reconstruction of the water surface level rather than the water depth at the interface provides more accurate and stable calculations when the bottom is non-flat. The water surface levels η (η =h+z) at left and right sides of the interface are calculated by using the following equations [27,28] where i ∆η , i-1 ∆η , and i+1 ∆η are respectively the limited gradients of cells i, i-1, and i+1.
are the constructed water surface elevation at the right and left sides of cell edges of i+1/2 and i-1/2. Similarly, the discharge q is reconstructed at the edges using the same formulations by replacing η by q. The Double Minmod (DB) limiter [27,29] is used and integrated in the present numerical solver to limit the oscillations. The values of i ∆η is calculated as follows: As described by van Leer [29] the Double Minmod algorithm in Eq. (13) limits the algebraic average of i-1/2 ∆η and i 1/2 ∆η + to twice of the smaller one. In comparison with the Superbee limiter [29,30], Bradford and Sanders [28] indicated that the DM limiter yields a better performance than the Superbee. Once water depths at the edges (for example, at i+1/2) are computed using The bottom slope term in Eq. (5) is computed using to preserve the balance between the flux and slope terms at stationary flow conditions. Including the finite difference schemes on the time and spatial discretization of the dispersive terms in Eq. (3) into the finite volume formulations of the conservative and source terms, the predictor formulation of the momentum equation (Eq. (6)) with unknown variable of q (q=uh) is given as In Eq. (18g), [2] 1/ 2 = ± n i F the second row of the F vector at the edge i+1/2 and i-1/2, respectively. Similarly, the numerical expression for the corrector (Eq. (7)) can also be formulated to compute The finite volume algorithm as a portion of the present hybrid FV-FD solver includes the Roe's second order scheme and an up to O(Δx 2 ) DM limiter. Here, Δx is the grid size. For the dispersive term, the central differencing scheme, which is also accurate to the order of Δx 2 , is used in the solver. Therefore, with a decrease of Δx, converged solutions are obtained. As the present solver adopts the predictorcorrector type solution procedure, Courant condition is utilized to control the stability of the numerical scheme. The condition in terms of the time step Δt is given as With the use of predictor-corrector scheme and the direct solution procedure of tri-diagonal matrix algorithm for Eq. (17), the present FV-FD solver is a cost-effective computational model for simulating propagation and transformation of nonlinear shallow-water waves.

Simulations of 1-D FV-FD Solver
The FV-FD solver is tested for cases of head-on collision of two solitary waves with equal or different amplitudes.

Head-on Collisions of Two Solitary Waves
The head-on collision of two solitary waves represents the process of two solitons colliding with each other. Similar incident wave conditions as given in Zhong and Wang [9,25] are used. In this test case, the interaction of two solitary waves with identical amplitude, i.e. a 1 =a 2 =0.1m, is investigated, where a 1 and a 2 are initial amplitudes of two propagating solitary waves. The length of the channel is 90 m (0 ≤ x ≤ 90 m) and the grid size Δx is set as 0.1 m. Zhong and Wang [9,25] developed two dimensional Boussinesq models. One solved the Fully Nonlinear Weakly Dispersive (FNWD) form of the BE and the other solved the Weakly Nonlinear Weakly Dispersive (WNWD) BE using the FE method. The present FV-FD model is simulated for head-on collision of two solitary waves and the results are compared with those computed using the models developed by Zhong and Wang [9,25]. The solitary waves are nonlinear shallow-water waves, so, nonlinear interaction results in the maximum amplitude during the collision to be greater than the sum of the two initial wave amplitudes. The peak wave amplitude during the collision process can be derived analytically and was given by Baytt-Smith [20] as where a c = peak wave amplitude during collision. The results obtained from the present solver are compared with those of WNWD and FNWD models by Zhong and Wang [9,25] in Figures 3a-3d. The plots include the time variation of the free-surface elevations before, during and after the collision. Initially, the peak of right-going solitary wave is located at x=25 m while the left-going solitary wave is situated at x=65 m ( Figure 3a). As can be seen from Figure 3b-Figure 3d, the present solver produces similar results as those from the WNWD and FNWD models. At t=4.5 sec (Figure 3b), the encountering process of the two solitary waves and partially merged waves can be noticed. The amplitudes of the two collided waves continue increasing until they merge completely into a single solitary wave with a peak amplitude at t=6.14 sec ( Figure  3c). The formulation in Eq. (19) (Baytt-Smith [20]) gives 0.205 m as the analytical peak amplitude in this test case. The maximum amplitude after collision is predicted to be at 6.14 sec by all models and the value is 0.2024 m obtained by the present solver, which is close to the analytical solution. There exists a small phase shift after the interaction of two solitary waves. As depicted in Figure 5d, the merged solitary waves are separated to recover back to their original amplitudes and shapes after colliding with each other. Overall, the results in terms of wave elevations and phases from the present solver agree reasonably well with the simulated solutions from the WNWD and FNWD models and the performance of the solver for capturing the process of head-on collision of two solitary waves is satisfactory. Head-on collision of two solitary waves with different amplitudes is also investigated to further the verification of the solver's capability of simulating interaction of  Figure  4 to provide the visual understanding of the head-on interaction of two solitary waves. According to the analytical formulation given in Eq. (19) the maximum amplitude during the collision of these two waves is calculated as 0.163 m. The peak amplitude and the occurrence time are predicted as 0.162 m and 6.2 sec, respectively by the present model. The predicted peak value is very close to the analytical solution.
From Figure 4, we notice that after the interaction, the solitary wave with small amplitude continues to propagate towards the negative x direction while the larger amplitude wave moves towards the positive x direction with both waves preserving their original shapes and amplitudes after the collision.

Shoaling of a Solitary Wave over a Sloping Bottom
At the process of a solitary wave propagating from a region of deeper water depth to that of shallower water depth, the main wave splits into a sequence of waves with decreasing amplitude [9,25,31,32]. This is commonly called the fission phenomena of a solitary wave. In this study, two cases of solitary wave shoaling over a sloping bottom are simulated with the developed FV-FD Boussinesq model. The computational domain and changes of water depth selected as the first test case for solitary wave shoaling study are similar to the center line condition of that used by Zhong and Wang [9,25]. The dimensionless domain boundaries are from -15 to 60 in x direction and the dimensionless undisturbed water depth varies from 1 to 0.271 in the shallower water depth region. The water depth variation in dimensionless form (referenced to the deeper water depth) within the domain (also in dimensionless form) utilized is given as The dimensionless wave elevations 0 / h ζ are plotted along the x direction at different times in Figure 5 to show the shoaling process as a solitary wave propagates from a region of deeper water depth to a shallower water-depth region. Initially, a solitary wave with dimensionless amplitude of 0.08 is situated at x=0 (Figure 5a). As seen in Figure 5b when the wave enters the sloping bottom section, a wave with small amplitude is reflected back by the submerged shoal  to propagate towards the upstream boundary (Figures 5c-5d). After the solitary wave enters the submerged shoal region (e.g. at t * =20), the wave form begins to tilt forward. Here * 0 0 / t gh t h = and h 0 is the water depth in the upstream domain. The fission characteristics are well observed by the transformation of the free surface elevations shown at t * =30, 50 and 70. As soon as the main wave propagates past the sloping bottom into the shallower water depth region (t * =30), the first secondary wave starts to separate from the leading wave. With the main wave traveling further downstream, additional secondary waves with smaller amplitudes are emerged behind the first secondary wave (Figures 5d and 5e). At t * =70, the interesting wave fission phenomenon reflected by the appearance of three secondary waves with decreasing wave amplitude following the main wave can be noticed. The amplitude of primary wave increases as it propagates in the shallower water depth region towards the downstream boundary. The present model results as shown in Figure 5 match reasonably well with those obtained by Zhong and Wang [25]. Although, a small phase lag is observed at t * =50 and 70 when compared to the WNWD and the FNWD models. In general, the present model results especially the wave peaks fit better with FNWD solutions than the WNWD ones. This suggests the nonlinear effects of the present FV-FD model are similar to those in FNWD model. The shoaling process of a nonlinear solitary wave shows increased effect of nonlinearity while propagating into a region with shallower water depth. This case provides the numerical investigation in understanding the fission phenomena by presenting the deformation of the crest of the primary solitary waves and the formation of the secondary waves.
Another case studied by Madsen and Mei [33] is also selected to test the model performance in simulating the shoaling of a solitary wave over a sloping bottom. A solitary wave with amplitude a=0.914 cm is simulated within a 6-m long domain including a 1:20 sloping bottom. The schematic diagram showing the domain and sloping bottom is given in Figure 6. Four gauging stations reported in Madsen and Mei [33] for the experimental measurements are marked as 1, 2, 3, and 4 in Figure 6. They are located respectively at x=150, 276.2, 365.1 and 451.46 cm. The water depth changes from h 0 =7.62 cm at the upstream domain to h 1 =3.81cm at the downstream region. The crest of the wave is located at x=-80 cm, initially. Madsen and Mei [33] presented theoretical solutions and experimental observations for the time variation of the free-surface profiles. The results obtained from the present model are plotted in Figures 7a-7d to compare to the theoretical solutions and experimental data given by Madsen and Mei [33]. The numerical results of Yuan and Wu [32] using a finite-difference model are also included in Figure 7 for comparisons. In Figure 7a, the results at t=0 correspond to the wave condition when the wave crest reaches the    Figure 7: Comparisons of free-surface profiles among the theoretical, experimental results [33] numerical solutions by Yuan and Wu [32], and the present model solutions at gauging positions: (a) 1, (b) 2, (c) 3, (d) 4. position 1. In Figure 7, it is noticed that the numerical results including those from the present model are shown to have a better agreement with measured data than the simplified but theoretically predicted values. The present numerical model performs reasonably well when compared to the results from the experimental measurements and numerical solutions by Yuan and Wu [32] especially at the positions 1 and 2. Both numerical models predict similar wave profiles and slightly overestimate the peak of the main wave at locations 3 (Figure 7(c)). At location 4, the results are presented in Figure 7d. The present numerical model can provide slightly better predictions than those from Yuan and Wu [32] when comparing to the measured data. The formation of a secondary wave as a result of the shoaling effect is clearly shown in Figures 7c and 7d. The present solver is again demonstrated to be able to model reasonably well on the propagation and shoaling of a solitary wave over a region with a sloping bottom.

Conclusions
A one dimensional FV-FD free-surface flow solver is developed to simulate propagation of solitary waves in regions of different water depths and their head-on collisions. The extended Boussinesq equations by Madsen and Sorensen [3] are solved numerically by using the combined FV and FD methods, where FV is applied to discretize the conservative and source terms and FD is for the dispersive terms. Second order schemes and limiters are used to control unwanted small oscillations. The interaction of two solitary waves with equal-amplitude is simulated with the solver. The whole process of head-on collision including merging, interaction, and recovering back to their original shapes is presented with the plots showing the time variation of freesurface profiles. The maximum amplitude of the waves after completely merging is found to be close to the analytical results. Also, comparisons of the present model results with those of WNWD and FNWD models developed by Zhong and Wang [25] suggest the good performance received by the present model. A solitary wave propagating over a sloping bottom is also simulated. The fission process is verified with the present solver after the solitary wave interacts with the variable bottom. As a solitary wave propagates over a sloping bottom and reaches the shallower water depth region, the main wave separates into a primary wave followed with a sequence of secondary waves with decreasing wave amplitudes. The wave profiles computed with the present solver at different times agree well with those obtained by Zhong and Wang [25]. Overall, the present FV-FD solver is demonstrated to be capable of modeling interactions of solitary waves and simulating fission process over domains of variable depth.