Simulation of the Shallow Water Equations with Transmission Boundary Conditions in the Bay of Bengal Region

 Corresponding author: Md. Masum Murshed Copyright © 2025 Author(s) retain the copyright of this article. This article is published under the terms of the Creative Commons Attribution License 4.0. Simulation of the Shallow Water Equations with Transmission Boundary Conditions in the Bay of Bengal Region Md. Masum Murshed * Department of Mathematics, University of Rajshahi, Rajshahi-6205, Bangladesh. World Journal of Advanced Research and Reviews, 2025, 28(01), 1802-1815 Publication history: Received on 17 September 2025; revised on 22 October 2025; accepted on 25 October 2025 Article DOI: https://doi.org/10.30574/wjarr.2025.28.1.3633 Abstract This study presents numerical simulations of the shallow water equations (SWEs) for the Bay of Bengal (BoB) using the Lagrange–Galerkin method (LGM) on a triangular mesh, with the transmission boundary conditions (TBCs). To examine the positional sensitivity of the transmission boundaries, the simulations are conducted with TBCs imposed at two distinct locations within the Bay of Bengal domain. The computed total mass and 𝐿2-norm of the surface elevation 𝜂 demonstrate that the transmission boundary condition performs efficiently and exhibits minimal dependence on its placement. These results indicate that the TBC is well-suited for modeling open-sea boundaries, ensuring smooth wave propagation without artificial reflection. The study serves as a foundational step toward developing an accurate and stable storm surge prediction framework for the Bay of Bengal using the Lagrange–Galerkin approach. Keywords: Shallow Water Equations; Transmission Boundary Conditions; Bay of Bengal; Lagrange–Galerkin Method; Numerical Simulation 1. Introduction The Shallow Water Equations (SWEs) form a coupled system consisting of a pure convection equation for the total wave height, denoted by 𝜙, and a simplified Navier–Stokes equation for the horizontal velocity 𝑢=(𝑢1,𝑢2)𝑇, obtained by depth-averaging the flow in the vertical (𝑥3) direction. These equations are widely used to model large-scale geophysical flows, such as tsunamis and storm surges, particularly in coastal and bay regions where accurate wave propagation is crucial. In numerical simulations of such phenomena, suitable boundary conditions at the open-sea boundaries are essential to prevent non-physical reflections when waves reach these boundaries (see Figure 1). To achieve this, Transmission Boundary Conditions (TBCs), as introduced in [4], are imposed on the open-sea boundaries Γ𝑇. These conditions allow outgoing waves to exit the computational domain smoothly by minimizing artificial reflections. The TBC is expressed as: 𝑢(𝑥,𝑡)=𝑐(𝑥)𝜂(𝑥,𝑡) 𝜙(𝑥,𝑡)𝑛(𝑥), …. (1) where 𝑐(𝑥) is a positive coefficient, 𝜂(𝑥,𝑡)=𝜙(𝑥,𝑡)−𝜁(𝑥) represents the free surface elevation from a reference level determined by the depth function 𝜁(𝑥), and 𝑛(𝑥) denotes the unit outward normal vector to the boundary. A number of studies [1–3, 9–15] have investigated storm surge and tidal simulations in the Bay of Bengal (BoB), which borders the coasts of Bangladesh and eastern India. Most of these models employ radiation-type boundary conditions to treat open-sea boundaries, which are mathematically similar to TBCs as described in [4]. However, as noted in [7], World Journal of Advanced Research and Reviews, 2025, 28(01), 1802-1815 1803 TBCs generally produce more stable and physically accurate results for SWEs compared to traditional radiation conditions. In [5], both theoretical and numerical stability analyses for SWEs with TBCs were performed using finite difference methods (FDM). While FDM is suitable for regular domains, such as rectangular grids, it becomes inefficient for realistic oceanic geometries that feature complex coastlines. In these cases, the finite element method (FEM), particularly with triangular meshing, offers greater flexibility and precision (see Figure 1). Figure 1 The Bay of Bengal and the coastal region of Bangladesh The Lagrange–Galerkin Method (LGM) combines the strengths of FEM with a time-discretization technique based on the material derivative, 𝜙𝑘+1(𝑥)−𝜙𝑘(𝑥−𝑢𝑘(𝑥)Δ𝑡) Δ𝑡 , Where the upwind point 𝑥−𝑢𝑘(𝑥)Δ𝑡 is used for numerical evaluation. If this point lies outside the computational domain, the nearest boundary value of 𝜙𝑘 is adopted. Unlike FDM, which can fail when 𝑢𝑘+1⋅𝑛<0 due to missing boundary data, LGM remains stable and effective under such conditions. In this study, numerical experiments are conducted to test how the performance of the transmission boundary condition depends on its spatial placement. The TBC is applied at two distinct positions within the Bay of Bengal to examine its robustness. The results confirm that the transmission boundary condition works efficiently and shows minimal sensitivity to its position. To the best of our knowledge, no previous work has implemented the Lagrange–Galerkin method for storm surge simulation in the Bay of Bengal. Therefore, this research provides a significant step toward developing a reliable and accurate storm surge prediction model for this region using LGM with transmission boundary conditions. The paper is organized as follows: Section 2 introduces the mathematical formulation of the problem. Section 3 describes the Lagrange–Galerkin numerical scheme. Section 4 presents the simulation results. Finally, Section 5 discusses the findings and concludes the study. 2. Statement of the problem Following the formulation presented in [5], the mathematical model considered in this study is based on the twodimensional Shallow Water Equations (SWEs). Let Ω⊂ℝ2 be a bounded domain, and let 𝑇>0 be a fixed time. The objective is to find (ϕ,𝑢)∶Ω×[0,𝑇]→ℝ×ℝ2 such that World Journal of Advanced Research and Reviews, 2025, 28(01), 1802-1815 1804 { 𝜕𝜙 𝜕𝑡+∇⋅(𝜙𝑢)=0 in Ω×(0,𝑇), 𝜌𝜙[𝜕𝑢 𝜕𝑡+(𝑢⋅∇)𝑢]−2𝜇∇⋅(𝜙𝐷(𝑢))+𝜌𝑔𝜙∇𝜂=0 in Ω×(0,𝑇), 𝜙=𝜂+𝜁 in Ω×(0,𝑇), (2) subject to the boundary conditions 𝑢 = 0 on Γ𝐷×(0,𝑇), (3) 𝑢 = 𝑐𝜂 𝜙𝑛 on Γ𝑇×(0,𝑇), (4) and initial conditions 𝑢=𝑢0, η=η0 in Ω, at 𝑡=0, (5) where 𝜙 denotes the total water height, and 𝑢=(𝑢1,𝑢2)𝑇 represents the horizontal velocity vector. The variable 𝜂: Ω×[0,𝑇]→ℝ corresponds to the free surface elevation relative to the reference level, while 𝜁(𝑥)>0 (𝑥∈Ω) denotes the undisturbed water depth measured from that reference level (see Figure 2). The strain-rate tensor is defined as 𝐷(𝑢)≔ ∇𝑢+(∇𝑢)𝑇 2, and 𝑛 is the outward unit normal vector on the boundary 𝜕Ω. The boundary Γ≔𝜕Ω consists of two disjoint parts: the Dirichlet boundary Γ𝐷 and the transmission boundary Γ𝑇 such that Γ=Γ𝐷∪Γ𝑇, Γ𝐷∩Γ𝑇= ∅. The Dirichlet boundary condition is imposed along the coastline, while the transmission boundary condition (TBC) is applied to the open-sea boundary to allow outgoing waves to exit the domain smoothly without artificial reflections. The physical constants are defined as follows: ρ> 0 is the density of water, μ>0 is the dynamic viscosity, and 𝑔 > 0 is the acceleration due to gravity. The coefficient 𝑐(𝑥) appearing in the transmission condition is given by 𝑐(𝑥)≔𝑐0√𝑔ζ(𝑥), where 𝑐0 is a positive constant. Throughout this study, it is assumed that ζ∈𝐶1(Ω). Figure 2 Model domain (see [5]) World Journal of Advanced Research and Reviews, 2025, 28(01), 1802-1815 1805 3. LG scheme Following the formulation presented in [5], a Lagrange–Galerkin (LG) numerical scheme is considered to approximate the shallow water equations. Let 𝒯𝒽={𝐾} denote a triangulation of the domain Ω, and 𝑀ℎ be the standard P1 (piecewise linear) finite element space. We define Ψℎ ≔𝑀ℎ or the water surface elevation η, and 𝑉ℎ(ψℎ)≔{𝑣ℎ∈𝑀ℎ2|𝑣ℎ(𝑄)=0, ∀𝑄 being a node on 𝛤𝐷 𝑣ℎ(𝑃)=𝑐(𝑃)𝜓ℎ(𝑃)−𝜁(𝑃) 𝜓ℎ(𝑃)𝑛(𝑃), ∀𝑃 being a node on Γ𝑇}, or the velocity field 𝑢. The goal of the LG scheme is to determine a sequence {(ϕℎ 𝑘,𝑢ℎ 𝑘)}𝑘=1 𝑁𝑇⊂Ψℎ×𝑉ℎ such that, for each time step 𝑘=1,…,𝑁𝑇, the following system holds: { ∫𝜙ℎ 𝑘−𝜙ℎ 𝑘−1 ∘𝑋1ℎ 𝑘−1𝛾ℎ𝑘−1 Δ𝑡 Ω𝜓ℎ𝑑𝑥=0, ∀𝜓ℎ∈Ψℎ, 𝜌∫𝜙ℎ 𝑘𝑢ℎ 𝑘−𝑢ℎ 𝑘−1 ∘𝑋1ℎ 𝑘−1 Δ𝑡 Ω⋅𝑣ℎ𝑑𝑥+2𝜇∫𝜙ℎ 𝑘𝐷(𝑢ℎ 𝑘) Ω:𝐷(𝑣ℎ)𝑑𝑥 +𝜌𝑔∫𝜙ℎ 𝑘 Ω∇𝜂ℎ 𝑘⋅𝑣ℎ𝑑𝑥=0, ∀𝑣ℎ∈𝑉ℎ, (6) 𝜙ℎ 𝑘=𝜂ℎ 𝑘+𝛱ℎFEM𝜁. Here, the mapping 𝑋1ℎ 𝑘:Ω→Ω is defined as 𝑋1ℎ 𝑘(𝑥)≔𝑥−𝑢ℎ 𝑘(𝑥)Δ𝑡, and γℎ 𝑘(𝑥)≔det (∂𝑋1ℎ 𝑘(𝑥) ∂𝑥 ) denotes the Jacobian determinant of the transformation. The composition symbol “∘” represents the functional composition, i.e. [𝑣ℎ∘𝑋1ℎ 𝑘](𝑥)≔𝑣ℎ(𝑋1ℎ 𝑘(𝑥)). The operator Πℎ FEM:𝐶(Ω)→𝑀ℎ is the standard Lagrange interpolation operator. The extension of a finite element function ψℎ to the entire plane ℝ2 is defined as ψℎ (𝑥)={ψℎ(𝑥), 𝑥∈Ω, ψℎ(𝑃𝑥), 𝑥∈ℝ2∖Ω, where 𝑃𝑥∈ Γ is the nearest nodal point to 𝑥. At each time step, the first equation in system (6) is used to determine ϕℎ 𝑘∈Ψℎ. Then, using this computed ϕℎ 𝑘, the second equation is solved to obtain the velocity 𝑢ℎ 𝑘∈𝑉ℎ. The first equation of (6) employs the concept of a mass-conservative Lagrange–Galerkin scheme as proposed in [16], ensuring numerical stability and conservation of mass throughout the simulation. World Journal of Advanced Research and Reviews, 2025, 28(01), 1802-1815 1806 Figure 3 A triangular mesh for the Bay of Bengal 4. Numerical results This section presents the numerical experiments conducted to validate the proposed Lagrange–Galerkin (LG) scheme for the shallow water equations 4.1. Simulation Results for the Bay of Bengal In this subsection, the simulation results for the Bay of Bengal region are discussed using the LG scheme introduced in Section 3. For computational simplicity, the actual geographical domain was slightly modified, and a triangular mesh was constructed as illustrated in Figure 3, following the approach in [6]. The complete setup of the numerical simulation is shown in Figure 4, and the computational domain is denoted by Ω. The horizontal extent of the domain ranges from 0 km to 1051.4 km, while the vertical extent covers 0 km to 889.59 km. Two types of boundary conditions were applied: • Zero Dirichlet boundary condition on Γ𝐷, representing the coastal boundary; and • Transmission boundary condition on Γ𝑇, representing the open-sea region. As illustrated in Figure 4, the transmission boundary Γ𝑇 is composed of three distinct segments: 𝐴𝐵=Γ𝑇1,𝐵𝐶=Γ𝑇2,𝐶𝐷=Γ𝑇3. The initial water surface elevation is prescribed as 𝜂⁰ = 0.05𝑒−0.1|𝑥 − 𝑝|², where 𝑝 denotes the nodal point nearest to the coordinates (525, 440). The simulation was performed using the following parameter values: 𝜁 = 2.0,𝜇 = 1,𝑔 = 9.8 × 10⁻³,𝜌 = 10¹², 𝑐0=0.9,ℎ =1.408044, 𝑇=5000s. The results obtained at different time levels 𝑡=0s,2800s,3120s,3240s,3740s,3940s, and 4660s are shown in the left panels of Figures 6 and 7. From these figures, it is observed that a circular wave is generated near the center of the domain and propagates outward over time. When the wavefront reaches the coastal boundary Γ𝐷, partial reflection occurs, while at the transmission boundary Γ𝑇, almost no reflection is detected. This demonstrates that the wave energy is smoothly transmitted through Γ𝑇, confirming the effectiveness of the transmission boundary condition for open-sea flow simulation. World Journal of Advanced Research and Reviews, 2025, 28(01), 1802-1815 1807 Figure 4 A figure of Bay of Bengal domain showing the setting of the transmission and the Dirichlet boundaries Figure 5 A figure of extended domain showing the setting of the transmission and the Dirichlet boundaries We also investigated the effect of slightly modifying the position of the transmission boundary. To do so, an extended domain Ω  was constructed by enlarging Ω by 100 km in the negative vertical direction, as shown in Figure 5. Consequently, Ω  extends from 0 km to 1051.4 km horizontally and from −100 km to 889.59 km vertically, with Ω⊂Ω . The Dirichlet boundary Γ𝐷 remains the same for both domains, while Γ𝑇 denotes the transmission boundary of Ω. World Journal of Advanced Research and Reviews, 2025, 28(01), 1802-1815 1808 A numerical simulation was then performed in the extended domain Ω  using the same parameter settings as before. The results at 𝑡=0s,2800s,3120s,3240s,3740s,3940s, and 4660s are presented in the right panels of Figures 6 and 7. The behavior of the wave propagation in Ω  closely resembles that observed in Ω, confirming that the transmission boundary condition remains stable and effective even when its position is slightly modified. From both sets of results (left and right panels of Figures 6 and 7), it is evident that the proposed LG scheme accurately reproduces the physical behavior of wave transmission in the open-sea region of the Bay of Bengal. 4.2. Computation of Mass and L² Norm for the Bay of Bengal This section presents the computation of the total mass and the L² norm of the surface elevation for the Bay of Bengal Figure 6 Simulation of SWEs in the Bay of Bengal at time 𝒕 = 𝟎𝒔,𝟐𝟖𝟎𝟎𝒔, and 𝟑𝟏𝟐𝟎𝒔 World Journal of Advanced Research and Reviews, 2025, 28(01), 1802-1815 1809 Figure 7 Simulation of SWEs in the Bay of Bengal at time 𝒕 =𝟑𝟐𝟒𝟎𝒔,𝟑𝟕𝟒𝟎𝒔,𝟑𝟗𝟒𝟎𝒔 and 𝟒𝟔𝟔𝟎𝒔 World Journal of Advanced Research and Reviews, 2025, 28(01), 1802-1815 1810 4.2.1. Computation of Mass The total mass M(t) at time tis defined as 𝑴(𝒕)∶=∫𝜼(𝒙,𝒕) 𝒅𝒙 ⬚ 𝜴≈∫(𝜫𝜼 𝒉 𝒌) 𝒅𝒙 ⬚ 𝜴 =: 𝑴 𝒉 𝒌, where 𝒕 = 𝒕ᵏ (𝒌 ≥ 𝟎,𝒌∈ℤ), and the integration is performed over the common part of the domains 𝛀 and 𝛀 . It should be noted that in the domain 𝛀, the boundaries AB, BC, and CD serve as transmission boundaries, whereas in the extended domain 𝛀 , the segment BC lies within the domain (see Figures 4 and 5). The computed mass variation with time is displayed in Figure 8. The time derivative of the mass is given by 𝑀′(𝑡) =∫𝜕𝜂 𝜕𝑡 ⬚ 𝛺(𝑥,𝑡) 𝑑𝑥 =∫𝜕𝜙 𝜕𝑡 ⬚ 𝛺(𝑥,𝑡) 𝑑𝑥 (from the third equation of (2) and 𝜁 = constant) =−∫𝛻·(𝜙𝑢)𝑑𝑥 ⬚ 𝛺 (from the first equation of (2)) =−∫(𝑢·𝑛)𝜙𝑑𝑠 ⬚ 𝜕𝛺 =−∫(𝑢·𝑛)𝜙𝑑𝑠 ⬚ 𝛤𝐷∪𝛤𝑇=−∫(𝑢·𝑛)𝜙𝑑𝑠 ⬚ 𝛤𝑇 =−∫𝑐𝜂𝑑𝑠 ⬚ 𝛤𝑇 (from the equation (4), here 𝑐 is a constant) Figure 8 Graphs of 𝑴 𝒉 𝒌 versus 𝒕 =𝒕𝒌 (≥ 𝟎,𝒌 ∈ ℤ) for the common part of the domains 𝛀 and 𝛀  From this relation, it follows that the variation of the total mass depends on the sign of 𝜂 along the transmission boundary 𝛤𝑇: • If 𝜂>0 on 𝛤𝑇, the total mass decreases; • If 𝜂<0 on 𝛤𝑇, the total mass increases. To better understand this behavior, the distributions of 𝜂 along the boundary segments AB, BC, and CD are plotted for both domains Ω and Ω  in Figures 9 and 10, complementing the mass evolution shown in Figure 8. The results presented in Figures 6–10 correspond to times𝑡=0s,2800s,3120s,3240s,3740s,3940s and 4660s. The following observations can be made: