1. Introduction
Internal solitary waves (ISW) are water waves which are generated and propagated within the interior of density stratified water on the Earth. The solitary wave phenomenon was first observed by John Scott Russell (1808-1882), who observed a solitary wave (known as a soliton) in the Union Canal in Scotland [1] . At first, Russel’s observations could not be understood by the existing water wave theories of that time until Lord Rayleigh submitted an article [2] to support Russel’s experimental observation based on Boussinesq’s research [3] . Korteweg and de Vries published an equation to deal with weakly nonlinear waves, such as those observed by Russell [4] . ISW (either elevation or depression) are nonlinear and in theory keep their amplitude and speed while propagating in the fluid due to a balance in nonlinearity and dispersion. They propagate along a pycnocline (a sharp transition between the upper mixed layer and the deep layers, where the density gradient is strongest) in the water body. ISWs account for up to one-third of the energy dissipation required to explain the slowing of the rotation of the Earth [5] . These waves can vary in amplitude from tens to hundreds of meters [6] [7] [8] and may travel thousands of kilometers from their sources [9] [10] before breaking on the continental slope or shore.
ISW interaction phenomena are able to carry particles from one level to another level easily. From a biological/environmental view point, ISWs have a significant effect on pollutant and nutrient transport in the seas and oceans. Shoaling can also result in energy dissipation. Furthermore, shoaling and run-up of ISWs cause external forces on deep water construction and hydraulic engineering issues. In addition, ISWs play a role in erosion of coastal land. For all of these reasons, it is important to study ISWs, particularly their effect on coastal regions. Helfrich and Melville [6] worked on solitary waves propagating in coastal areas and noticed that this type of wave is a ubiquitous feature around the Earth’s coastal regions. The shoaling of internal waves and their interactions in the Earth’s oceans have been also studied by Wunsch and Ferrari [11] and Ferrari and Wunsch [12] . The latter concluded that internal waves appear to be a major source of wave breaking and mixing within the ocean’s interior.
As a solitary wave propagates, it encounters the slope of the coast and breaks. While several researchers investigated the propagation of the solitary waves [13] [14] , others investigated the breaking of a solitary wave on an inclined surface. Kao, Pan and Renouard [15] considered shoaling and breaking of an internal solitary wave of elevation for a sloped shelf and concluded that shear instability can cause wave breaking. Saffarinia and Kao [16] studied the breaking of a soliton on the stratified pycnocline in a two-layer system using the Navier-Stokes and diffusion equations. They defined a criterion for the breaking of a solitary wave and concluded that breaking occurs if the particle velocity exceeds the wave celerity in the fluid domain.
Pede [17] conducted research about an internal wave of elevation encountering a shoal and a shelf. Li et al. [18] studied a numerical tank using OpenFOAM and captured the interface of fluids using the volume of fluid (VOF) method. They provided comparisons between several theoretical solutions and the numerical results. Li et al. [19] simulated three-dimensional internal solitary waves using Internal Solitary Wave Open Foam (ISWFoam). They compared the numerical results of ISWFoam with laboratory results such as the ISW profile and the wave breaking location and found good similarity concluding that ISWFoam could effectively simulate the ISW breaking phenomenon and the interaction between ISWs and a complex topography.
Vieira et al. [20] used a hyperbolic tangent profile to model the density stratification to better mimic the continuous form of the actual pycnoclines in ocean. The bolus (a sub-surface dense fluid vortex formed as ISW shoals on slopes as stated in [21] ) size and displacement upslope were examined as a function of the pycnocline thickness, incoming wave energy, density change across the pycnocline, and topographic slope. They concluded that boluses in continuous stratifications would tend to be bigger and transport material farther than in corresponding two-layers of uniform density. Arthur and Fringer [22] investigated transport phenomena due to breaking internal waves on slopes using particle-tracking in three-dimensional simulations. They demonstrated onshore and offshore transport within the bolus and lateral particle transport away from the slope due to turbulence developed by the breaking. Masunaga et al. [23] modeled sediment resuspension due to boluses with a transport equation for suspended sediment, obtaining resuspension processes that are in good agreement with observational data and investigating the effect of varying the topographic slope.
This paper presents an investigation into the interaction between internal solitary waves and a shoal. The novelty of this study lies in estimating the critical angle of the shoal, above which bolus formation is not expected to occur. As previously mentioned, when the wave encounters the shoal, a bolus, characterized by a region of high vorticity of dense fluid, forms. Boluses have the potential to transport substances from deep to shallow levels, which could lead to corrosion and damage to man-made equipment submerged in oceans. The initial section of this paper outlines simulations focusing on the interaction between waves and slopes, presenting diverse findings including the maximum altitude reached by the bolus. Additionally, measurements of the vorticity of the pycnocline post-interaction with the slope were conducted, revealing that vorticity escalates during the bolus formation process as it encounters the shoal.
The subsequent segment of this paper investigates the existence of a critical angle above which bolus formation does not occur, while also exploring how parameters such as the Reynolds number influence this process. Various scenarios of shoals are analyzed to showcase instances where bolus formation is absent. Subsequently, a highly shallow shoal is scrutinized, affirming the hypothesis that shallower shoals lead to larger boluses. This data holds potential for predicting the impact of internal solitary waves on transport in coastal areas.
As far as the authors are aware, no prior study has attempted to predict bolus formation based on shoal slope, as demonstrated in the present research.
2. Mathematical Calculation
2.1. Numerical Domain
In this study, the numerical domain consists of two fluids with different densities (
${\rho}_{1}$ and
${\rho}_{2}$ ), initially positioned with the denser fluid at the bottom of the domain. A density step is introduced to initiate the formation of a solitary wave. Figure 1 illustrates the schematic of the numerical domain, where the tank’s
depth is denoted as d, and its length before the slope is defined as 5d.
As detailed in [24] , the length of 5d is chosen to allow the leading solitary wave to travel and separate from the initial disturbance before encountering the shoal. The slope is adjusted to determine a critical angle for bolus formation by altering the value of L, where the initial value for L is set to 5d. The angle of the tank slope, denoted as
$\phi $ , is calculated as
$\phi =\text{atan}\left(\frac{d}{L}\right)$ .
${h}_{1}$ and
${h}_{2}$ represent the initial heights of the upper and lower fluids, respectively. The vertical coordinate (Y) is referenced from the upper surface of the domain, while the horizontal coordinate (X) originates at the beginning of the slope.
2.2. Governing Equations and Flow Configuration
A two-dimensional, incompressible, Newtonian, viscous, and unsteady flow is under consideration. The governing equations encompass the continuity equation, the Navier-Stokes equation incorporating the Boussinesq approximation, and the convection-diffusion equation. According to the Boussinesq approximation, density differences can be disregarded except where they are multiplied by the acceleration due to gravity. Consequently, the governing equations are modified as follows:
$\rho \nabla \cdot V=0$ (1)
$\rho \frac{\partial V}{\partial t}+\rho \left(V\cdot \nabla \right)V=\nabla \cdot \left[-p+\mu \left(\nabla V+\nabla {V}^{\text{T}}\right)\right]+\Delta \rho g\theta $ (2)
$\frac{\partial \theta}{\partial t}+\nabla \cdot \left(-D\nabla \theta \right)+V\cdot \nabla \theta =0$ (3)
In the provided equations,
$V$ represents velocity,
$\rho $ denotes the density of fluid 1,
$\Delta \rho ={\rho}_{2}-{\rho}_{1}$ , and
$\mu $ stands for the viscosity of fluid 1. The density difference ratio, denoted as
$\theta $ , is defined as
$\theta =\frac{\rho -{\rho}_{1}}{{\rho}_{2}-{\rho}_{1}}$ , ranging from 1 at the bottom of the tank to 0 at the upper level.
To non-dimensionalize Equations (1)-(3), all lengths are scaled by the total depth of the tank, d. Additionally, velocities are scaled using the wave celerity, c_{0}.
${c}_{0}=\sqrt{\frac{g\left|\frac{{\rho}_{1}-{\rho}_{2}}{{\rho}_{2}}\right|{h}_{1}{h}_{2}}{d}}$ (4)
where
${\rho}_{2}$ is used for the scaling of density.
With this nondimensionalization, the nondimensional numbers governing the problem are
$Sc=\frac{\mu}{\rho D}\mathrm{,}\text{\hspace{1em}}Re=\frac{\rho {c}_{0}d}{\mu}\mathrm{,}\text{\hspace{1em}}F{r}_{d}^{2}=\frac{{h}_{1}{h}_{2}}{{d}^{2}}\frac{\left({\rho}_{1}-{\rho}_{2}\right)}{{\rho}_{2}}=\frac{{c}_{0}^{2}}{gd}$ (5)
The Schmidt number, Sc, represents the ratio of viscous diffusion to molecular diffusion. The Reynolds number, symbolized as Re, signifies the ratio of inertial force to viscous force. Lastly, the Froude number, designated as Fr, denotes the ratio of inertia to gravitational force.
2.3. Boundary and Initial Conditions
In Figure 1, the left, bottom, and shoal boundaries are modeled as no-slip and non-diffusive walls. The upper boundary of the domain is represented as a pure-slip free surface (rigid lid). According to [15] , this configuration is appropriate for the upper level of a tank because the Froude number (Fr) is small, reflecting the negligible density difference relative to the bulk density.
Initially, the tank is stratified, containing two layers of fluid with distinct densities. At time = 0, the tank remains at rest. The disturbance is initiated through the initial density distribution which is a smoothed step as illustrated in Figure 2. This distribution is generated by:
$\theta =0.5\left(1-\text{erf}\left(y-\stackrel{\xaf}{y}\left(x\right)\right)/\eta \right)$ (6)
where
$\stackrel{\xaf}{y}\left(x\right)=-{h}_{2}+\frac{a}{2}\left(1-\text{erf}\left(\frac{x-\stackrel{\xaf}{x}}{\eta}\right)\right)$ and
$\stackrel{\xaf}{x}=-5d+w$
In the provided formula,
$\eta $ specifies the thickness of the density transition. Table 1 displays the parameters employed for all subsequent simulations, unless explicitly stated otherwise in the text.
As the bump is released, the difference in densities and the tendency of the system to restore itself, cause a flow in the tank. The initial energy is converted to kinetic energy, and an internal solitary wave is formed.
2.4. Numerical Method
The simulations were conducted using COMSOL Multiphysics. A mixed mesh comprising quadrilateral and triangular elements was employed. An illustration of the unstructured mesh is depicted in Figure 3. Figures 3(b)-(c) provide close-up views of the selected region (highlighted by the black oval) in Figure 3(a), offering a clearer insight into the mesh distribution within the domain. Details regarding the grid resolution are provided in Table 2. Element counts for different configurations are specified in the results.
COMSOL employs an adaptive time-stepping scheme regulated by a user-specified tolerance. Throughout all simulations presented here, this tolerance was set to 0.001, and the free adaptation scheme was selected utilizing the backwards differentiation formula, as described in [25] .
3. Results and Discussion
To begin, simulations of an internal solitary wave on three distinct slopes are
Figure 3. Unstructured mesh distribution.
outlined, with subsequent discussions on the results of each case. Subsequently, the impact of varying Reynolds numbers (Re) on a particular case is investigated.
3.1. Parameter Values
Table 3 provides a summary of the length-to-depth ratio (L/D), the corresponding slope angle, and the number of elements utilized for the three cases. It is noteworthy that different mesh sizes are employed for the scenario featuring an 11˚ slope, as the domain is wider compared to the other cases.
3.2. ISW Hitting the Shoal with Slope of 11˚
Figure 4 shows density contours at times t = 0, 4, 6.75, 7.74, and 8.89 s for the case of θ = 11˚. Time = 4 corresponds to the onset of the wave onto the shoal and t = 6.75 s corresponds to the initial formation of the bolus. Figure 4(d) shows the run-up of the bolus and Figure 4(e) shows a point just before the bolus reverses direction. Figure 5 shows a close view of the last three times in Figure 4. These figures show that a bolus can transport the high density fluid
Figure 4. Density-difference contour of the wave, (a) time = 0; (b) time = 4; (c) time = 6.75; (d) time = 7.74; (e) time = 8.89, Slope = 11˚.
significant distances up the shoal. The maximum height the bolus reaches is approximately x = 2.55, y = −0.47 which occurs at time = 9.5 s. This corresponds to a vertical transport of 0.5 d.
Figure 5. Closer view at density-difference contour, (a) time = 6.75; (b) time = 7.74; (c) time = 8.89, Slope = 11˚.
Figure 6 displays velocity vector plots overlaid on density contours. Inside the bolus as it ascends the shoal, the fluid exhibits a counter-clockwise rotation. At time t = 7.74, the velocity vector indicates fluid drainage down the shoal from the rear of the bolus (notice the reversed direction of the velocity vector arrows); nevertheless, the fluid rotation within the bolus remains counter-clockwise. By time t = 8.89, more fluid drains down onto the shoal, causing the bolus to spread out.
Vorticity contours of the internal solitary wave at the same time instances are depicted in Figure 7. In Figure 7(a), the bolus is identified as a region characterized by high positive vorticity. As time progresses, illustrated in Figure 7(b), the vorticity gradually diminishes. Employing the method of image vortices, it can be deduced that the combination of positive vorticity interacting with the solid wall leads to an induced velocity of the vorticity in the direction of the bolus’s motion.
3.3. Bolus Properties
To ascertain the characteristics of the bolus, both its location and vorticity values are recorded and plotted as it ascends the shoal at three distinct angles (11˚, 21˚, 27˚). The location refers to where the maximum positive vorticity of the bolus is observed.
In Figure 8(a), it is evident that as the wave interacts with the shoal at a 11˚ angle, vorticity undergoes a rapid escalation, reaching its peak before gradually diminishing as the bolus progresses up the slope over time. This accumulation of
Figure 6. Velocity vector with corresponding density contour, (a) time = 6.75; (b) time = 7.74; (c) time = 8.89, Slope = 11˚.
Figure 7. Close view at vorticity contour, (a) time = 7.74; (b) time = 8.89, Slope = 11˚.
Figure 8. Comparison of three cases, (a) Vorticity of the bolus; (b) Speed of the bolus.
vorticity is closely associated with the bolus breaking away from the leading edge of the wave. Due to the higher density of fluid within the bolus compared to its surroundings, fluid drains down the slope, leading to a reduction in vorticity to zero. A similar trend is observed in the vorticity patterns for angles of 21˚ and 27˚. Vorticity exhibits an increase as the wave encounters the shoal, followed by a decline as the bolus advances up the slope. Notably, the peak vorticity value at a 27˚ angle is lower compared to the other two cases (Figure 8(a)).
Figure 8(b) shows the speed of the bolus as time passes. This figure shows that speed of the movement decreased until the bolus spread out. The speed of the onset wave is about 1. Hence, the bolus moves up the shoal at nearly half of the onset wave speed for angle of 11˚. The maximum speed for case with angle of 21˚ case is almost one third of the onset wave speed and less than of the previous case. In addition, the bolus spreads out in a shorter time in comparison with the previous case (11˚). It can be concluded that the case with higher value of vorticity, has higher value of speed. The maximum speed of for the last case (angle of 27˚) is less than the two previous cases and the bolus stops in a shorter time range compared to the previous cases. The peak value of vorticity is the lowest one in these three cases (Figure 8(b)).
Table 4 presents the maximum values of the x location, y location, vorticity, speed, and total length of the bolus during the run-up process for the three cases.
On the 21˚ shoal, the bolus ascends near the coordinates (x = 1.30 and y = −0.48) before spreading out. Consequently, the bolus traverses a total length of 1.4 units along the shoal, which is shorter than the distance covered in the case with a 11˚ angle (x = 2.55, y = −0.47). In the third case, the bolus advances up the shoal close to the point (x = 0.99 and y = −0.49) before dispersing. It travels a total length of 1.11 units along the shoal.
The length of the run-up in the 27˚ angle case is shorter compared to the two previous cases, but the maximum run-up in the y direction remains nearly the same. The magnitudes of the variables mentioned above for the 27˚ angle case are the lowest among all cases. This suggests that as the angle decreases, the bolus travels further along the slope. However, it reaches almost the same depth in all three cases. Additionally, both the maximum vorticity and the maximum speed of movement on the shoal decrease as the angle increases.
3.4. Determination of the Critical Angle
To ascertain the critical angle (or a range of angles) for bolus formation, multiple cases need to be analyzed. It has been observed that when the angle is 90˚, bolus formation cannot be observed [25] . This suggests that as the slope increases from 11˚ to 90˚, bolus formation does not occur. In this study, a single-density contour is utilized to identify the critical angle. For angles exceeding the critical angle, this density contour remains predominantly connected, whereas for angles below it, the contour becomes disconnected. The slope of the shoal is gradually adjusted to determine the angle at which a separation in the density contour occurs. The chosen density levels for examination are 0.65 and 0.85.
Figure 9 and Figure 10 depict an angle slightly below the critical angle for a
Table 4. Maximum values of variables.
Figure 9. Density contour level of 0.85 - Angle 21˚, (a) time = 6.67; (b) time = 6.70.
Figure 10. Density contour level of 0.85 - Angle of 22˚, (a) time = 6.5; (b) time = 6.6; (c) time = 6.7.
density contour of 0.85. At this angle, the density contour separates before ascending the shoal. In Figure 10, it can be observed that the density contour remains continuous as it advances and retreats on a steeper shoal. Consequently, the critical slope for bolus formation at a density contour level of 0.85 is determined to be 21˚.
The choice of density contour has a degree of sensitivity in determining the critical angle. For a density contour of 0.65 (results not presented here to avoid redundancy), the critical angle is established to be 27˚, which is 6˚ higher than that determined for the contour level of 0.85.
3.5. Effect of Reynolds Number for a Slope of 27˚
In this section, the results concerning the calculation of four different Reynolds numbers—1250, 2500, 5000, and 10,000—on a slope of 27˚ are presented. All other variables are kept constant, including the mesh resolution (446,596 elements) with Sc = 10. As illustrated in Figure 11(a), an increase in the Reynolds number results in the bolus moving further along the slope in the y-direction. Moreover, as the Reynolds number increases, the time at which maximum vorticity occurs also increases (Figure 11(b)).
Figure 11(c) illustrates the variations in bolus speed values as it ascends the shoal. A higher Reynolds number corresponds to a slight increase in maximum
Figure 11. Constant Sc = 10 and different Reynolds numbers; 1250, 2500, 5000, 10,000 (a) y location of the bolus in the run-up; (b) Vorticity of the bolus in run-up; (c) Speed of the bolus in run-up.
speed. It can be inferred that although the maximum peak vorticity values differ significantly, they are distributed over a smaller area. Consequently, the circulation value is not proportionally high, resulting in a smaller difference between speed values. It is noteworthy that for Re = 10,000, the bolus was not detected after time = 6.9 s, which renders the data invalid at time = 7.1 s, unlike the other Reynolds numbers.
3.6. Mesh Resolution Study
To ensure the robustness of the aforementioned solutions against mesh resolution, the bolus location is tracked over time using five different mesh resolutions (27,747, 55,743, 111,882, 222,014, 446,596 elements). This analysis is conducted for the case with a slope of 27˚ and Re = 2500, while keeping other parameters consistent with previous cases.
The results depicted in Figures 12(a)-(c) showcase the evolution of the x-location, y-location, and vorticity of the bolus over time. It is evident from the data that the position and vorticity of the bolus remain nearly identical for the third, fourth, and fifth mesh resolutions.
Table 5 provides a summary of the maximum difference observed between
Figure 12. (a) x-location of the bolus in run-up; (b) y-location of the bolus in run-up; (c) Vorticity of the bolusup.
Table 5. Angle 27˚-Mesh resolution study.
the last two mesh resolutions. The greatest disparity is found in the vorticity, which amounts to only 3%.
Based on the observations from Figure 12 and Table 5, it becomes apparent that as the number of elements increases to 222,014, the discrepancy between the case with 446,596 elements and the case with 222,014 elements exhibits no substantial variance. Nonetheless, it is important to note that in this paper, the mesh grid size corresponding to 446,596 elements was utilized for the preceding cases.
4. Conclusion
The research delved into the propagation of an internal solitary wave of elevation on the continuously stratified pycnocline within a two-layer system, and its interaction with various slopes. It was determined that the critical slope angle for bolus formation falls within the range of 21˚ to 27˚ and that the vertical height to which the bolus rises is nearly independent of the slope. As the Reynolds number increases, the maximum vorticity, the bolus speed, and the height obtained by the bolus all increase. This paper provided a good sense of the 2D behavior of bolus propagation, but further studies are needed at higher Reynolds numbers to understand how turbulent and 3D effects change the behavior.