EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 6230 ISSN 1307-5543 – ejpam.com Published by New York Business Global Numerical Simulation of the Ripa Model Using the Upwind CE/SE Scheme Saqib Zia1,∗, Asad Rehman2, Samina Majeed1, Abuzar Ghaffari1, Shagufta Yasmeen1, Wei Sin Koh3, Ilyas Khan4,5,6 1 Department of Mathematics, COMSATS University Islamabad, Park Road Chak Shehzad, Islamabad, Pakistan 2 Department of Information Technology, Pothohar Campus Gujar Khan, University of Punjab, Pakistan 3 INTI International University, Persiaran Perdana BBN Putra, Nilai 71800, Negeri Sembilan, Malaysia 4 Department of Mathematical Sciences, Saveetha School of Engineering, SIMATS, Chennai, Tamil Nadu, India 5 Hourani Center for Applied Scientific Research, Al-Ahliyya Amman University, Amman, Jordan 6 Department of Mathematics, College of Science Al-Zulfi, Majmaah University, Al-Majmaah 11952, Saudi Arabia Abstract. This paper focuses on the numerical investigation of the shallow water equations (SWEs) with temperature gradients, commonly referred to as the Ripa system, using the finite volume Upwind Conservation Element Solution Element (CE/SE) scheme. The inclusion of source- term on the right-hand side of the model introduces additional complexity, making the equations non-conservative and challenging to solve numerically. The primary objective is to develop a numerical scheme that effectively handles non-conservative differential terms with accuracy and efficiency. To address these challenges, a finite volume upwind CE/SE technique is proposed. The method is tested on several numerical problems to evaluate its effectiveness, robustness, and sta- bility, particularly in resolving discontinuities and shocks. The results of the proposed scheme are compared with those obtained using the standard CE/SE scheme. 2020 Mathematics Subject Classifications: 65N30, 74S05, 76M10 Key Words and Phrases: Ripa model, non-conservative systems, horizontal temperature gra- dient, bottom topography, upwind CE/SE scheme, CE/SE scheme, smart grid ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.6230 Email addresses: saqibzia81@hotmail.com (S. Zia), assad013@gmail.com (A. Rehman), muqadass830@gmail.com (S. Majeed), abuzarghaffari45@gmail.com (A. Ghaffari), shagufta yasmeen@comsats.edu.pk (S. Yasmeen), weisin.koh@newinti.edu.my (W. S. Koh), i.said@mu.edu.sa (I. Khan) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 2 of 21 1. Introduction Shallow water flows (SWFs) refer to flows where the vertical scale length is significantly smaller than the horizontal scale length, a characteristic often observed in real-world sce- narios. These flows are commonly occur in nature particularly in river systems, ocean currents, and coastal waves, and play a crucial role in understanding various engineering applications and environmental processes[1, 2]. The equations governing shallow water flows are derived from the principles of conservation of momentum, mass, and energy, forming a mathematical framework to describe the behavior of inviscid and incompressible fluids in shallow regions [3] . These flows have been extensively studied across disciplines such as hydrodynamics, geophysics, and coastal engineering, reflecting their importance in research and technological development. They play a critical role in advancing our understanding of water-related phenomena and mitigating associated risks. However, the intrinsic complexity of shallow flows makes exact solution difficult, requiring precise com- putational techniques to ensure physically realistic outcomes. Saint-Venant first proposed the SWEs in 1871 to simulate flow in an open channel [4]. The model equations are a set of nonlinear hyperbolic equations derived from multi- layered models in which the velocity field, density, and horizontal pressure gradient are vertically integrated. Temperature-dependent horizontal pressure gradients cause fluid density changes within every layer [5]. These equations describe fluid flow below the sur- face of pressure (or sometimes a free surface as well). With source term and temperature gradients, SWEs are important because they allow for considering additional factors such as the effects of temperature on fluid behavior, enabling more realistic simulations of phe- nomena. In 1993, Ripa introduced SWEs with temperature gradients [6]. This model has already been approximated through several numerical techniques in the past. In their ground- breaking work [7], Savage and Hutter developed a one-dimensional shallow-water model for analyzing aerial avalanches, which was later expanded to two dimensions [8]. Numerous numerical methods have been developed over time to solve the Saint-Venant equations. These include central upwind (CUP), KFVS, CE/SE, and finite volume WENO schemes [9–13]. Based on Runge Kutta Discontinuous Galerkin (RKDG), as introduced by Cock- burn (1999), a numerical solution is presented for SWE [14]. These equations are also numerically approximated by the finite difference method [15] and space-time finite ele- ment formulation is also used for their numerical approximation [16]. A suitable numerical technique can preserve the positivity of the estimated water depth h while also capturing steady states and small perturbations. The majority of numerical approaches for solving balance laws employ upwind-type fi- nite volume techniques. These approaches need an exact or estimated Riemann-solver [17] to determine flows at cell interfaces. When accurate Riemann solutions are easily accessi- ble, such procedures are suitable, reliable, and positivity-preserving. Designing Riemann S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 3 of 21 solutions for complicated models, such as the current Ripa system with variable bottom typography and a two-phase shallow flow model, is difficult. For efficient and accurate re- sults, a central-type finite volume approach is used that does not require a Riemann solver. Chang introduced the space-time CE/SE technique [18] for the computational esti- mation of hyperbolic conservation laws. This method stands apart from conventional approaches like FDM, FEM, FVM, and spectral methods, both conceptually and method- ologically. Unlike these techniques, the CE/SE method treats space and time with equal emphasis, ensuring conservation of both local and global fluxes across space and time do- mains. It incorporates conservation and solution elements and utilizes an explicit scheme with staggered grids. A central CE/SE scheme developed for solving nonlinear hyperbolic conservation laws has the notable advantage of not requiring explicit knowledge of eigen- values. Building on this, a characteristic-based upwind CE/SE method was developed, integrating the principles of the CE/SE approach with upwind numerical fluxes derived from the Godunov-type finite volume method. This technique does not depend on the CFL number, eliminating the drawbacks associated with previous schemes. For challeng- ing computational fluid dynamics (CFD) problems like smart grid systems, multiphase flow simulations, the upwind CE/SE scheme offers improved robustness and accuracy in capturing flow field discontinuities, particularly for contact discontinuities such as at ma- terial interfaces [13]. The upwind CE/SE scheme combines the principles of upwind differencing with the CE/SE framework to deliver efficient and accurate solutions for fluid flow simulations. This variant of the CE/SE method is specifically designed to address challenges com- monly encountered in numerical fluid dynamics. By employing an upwind technique to handle advection terms in the modeling equations, it offers an effective alternative to other methods. This approach ensures that numerical fluxes are influenced by the flow direction, simplifying their computation and improving the method’s efficiency, particularly in cap- turing discontinuities and sharp gradients. By leveraging information from upwind cells, the scheme accurately resolves flow characteristics and captures dominant flow patterns. These properties make it particularly effective for modeling flows such as supersonic flows or flows around obstacles with strong directional non-homogeneity. The upwind CE/SE approach is an effective method for addressing the numerical challenges associated with Ripa systems and other complex flow behaviors. It overcomes many of the primary limitations of earlier methods, as demonstrated in [19]. This article is organized as follows: Sections 2 and 3 present the one-dimensional and two- dimensional mathematical equations of the Ripa model, while the upwind CE/SE scheme is described in Section 4. Numerical test problems are discussed in Section 5, and Section 6 concludes the paper. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 4 of 21 2. One-dimensional equations The one-dimensional equations for shallow water with a horizontal temperature vari- ation [19] i.e., the Ripa system, can be represented in the form: ∂h ∂t + ∂ ∂x (m1) = 0, (1) ∂ ∂t (m1) + ∂ ∂x ( (m1) 2 h + g 2 h(m2) ) = −gm2 ∂B ∂x , (2) ∂ ∂t (m2) + ∂ ∂x ( (m1)(m2) h ) = 0. (3) Where h stands for the height of flow, g denotes the constant of gravitational acceler- ation, m1 = hu, m2 = hθ, θ signifies the potential temperature field, u is the flow velocity in the x direction, and variable bottom topography is indicated by B = B(x), x ∈ R. The system of equations (1)-(3) allows several steady-state terms, including the following: u = 0, θ = constant, and h+B = constant, u = 0, B = constant, and p = g 2 h2θ = constant. (4) In compact form, the Eqs. (1)-(3) can be re-written as ∂qm ∂t + ∂fm ∂x = τm, m = 1, 2, 3 (5) where q1 def = h, q2 def = hu, q3 def = hθ and the fluxes are given as f1 def = m1 = q2, f2 def = (m1) 2 h + g 2 hm2 = (q2) 2 q1 + g 2 (q1)(q3), f3 def = m1m2 h = (q2)(q3) q1 . (6) Moreover, τ1 = 0, τ2 def = −gh ∂B ∂x = −gq3 ∂B ∂x , τ3 = 0. (7) From now on, conserved variables will be represented by qi, i = 1, 2, 3. The above equa- tions can be expressed in quasi-linear form as follows: ∂tq +A(q)∂x(q) = R(q) . (8) Here q = (q1, q2, q3) T , R(q) = [0,−ghθ∂xB, 0]T and A(q) indicates Jacobian matrix which has ∂fk ∂qk element at the k-th row and l-th column for k,m = 1, 2, 3 and the functions fk(q1, q2, q3) are outlined in Eq. (6) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 5 of 21 3. Two-dimensional equations The 2D Ripa system [19] is given as ∂th+ ∂x(m1) x s + ∂y(m3) y s = 0 , (9) ∂t(m1) x s + ∂x ( ((m1) x s ) 2 h + g 2 h(m2) θ s ) + ∂y ( (m1) x s (m3) y s h ) = −g(m2) θ s∂xB , (10) ∂t(m3) y s + ∂x ( (m1) x s (m3) y s h ) + ∂y ( ((m3) y s)2 h + g 2 h(m2) θ s ) = −g(m2) θ s∂yB , (11) ∂t(m2) θ s + ∂x ( (m1) x s (m2) θ s h ) + ∂y ( (m3) y s(m2) θ s h ) = 0 . (12) Here, (m1) x s = hu, (m3) y s = hv and (m2) θ s = hθ, Where h(x, y, t) represents water depth, u(x, y, t) and v(x, y, t) emphasize the fluid movement’s speed in the x and y-direction respectively. On the other hand, B(x, y) is the bottom topography, and g represents the gravitational acceleration. Moreover, θ is the field of potential temperature. The equation system above can be rewritten in a compact form as ∂qm ∂t + ∂fm ∂x + ∂gm ∂y = τm, m = 1, 2, 3, 4 . (13) The conservative variables in this case are given by q1 def = h, q2 def = (m1) x s = hu, q3 def = (m3) y s = hv, q4 def = (m2) θ s = hθ , (14) and we can express fluxes as f1 def = (m1) x s = q2, f2 def = ∂x ( ((m1) x s ) 2 h + g 2 h(m2) θ s ) = (q2) 2 q1 + g 2 q1q4, f3 def = ∂x ( (m1) x s (m3) y s h ) = q2q3 q1 , f4 def = ∂x ( (m1) x s (m2) θ s h ) = q2q4 q1 , (15) g1 def = ∂y(m3) y s = q3, g2 def = ∂y ( (m1) x s (m3) y s h ) = q2q3 q1 , g3 def = ∂y ( ((m3) y s)2 h + g 2 h(m2) θ s ) = (q3) 2 q1 + g 2 q1q4, g4 def = ∂y ( (m3) y s(m2) θ s h ) = q3q4 q1 . (16) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 6 of 21 Furthermore, τ1 = 0, τ2 def = −ghθ∂xB = −gq4∂xB, τ3 def = −ghθ∂xB = −gq4∂yB, τ4 = 0. (17) Hereafter, qm, m = 1, 2, 3 are termed as conservative variable. Now to find the nature of our system we find characteristic values for 1D and 2D model equations, where A(q) = ∂f(q) ∂q is a Jacobian matrix defined as: A(q) =  ∂f1 ∂q1 ∂f1 ∂q2 ∂f1 ∂q3 ∂f2 ∂q1 ∂f2 ∂q2 ∂f2 ∂q3 ∂f3 ∂q1 ∂f3 ∂q2 ∂f3 ∂q3  . The eigenvalues are given as, λ1 = u, λ2 = u+ √ ghθ 2 , λ3 = u− √ ghθ 2 . (18) Since eigenvalues are real, the considered 1D Ripa system is hyperbolic. The eigenvalues for the 2D Ripa system are given as, µ1 = µ2 = u, µ3 = u+ √ ghθ 2 , µ4 = u− √ ghθ 2 γ1 = γ2 = v, γ3 = v + √ ghθ 2 , γ4 = v − √ ghθ 2 As shown here, our system is not strictly hyperbolic because the two eigenvalues are the same. 4. The upwind CE/SE scheme for one-dimensional Ripa system Here, we consider the one-dimensional Ripa system’s numerical solution, ∂wi ∂t + ∂fi ∂x = τi, i = 1, 2, 3. (19) Time and space are taken as unified in the CE/SE method. With the spatial coordinate x, we first consider time t as an equal footing. Here we consider a homogeneous case for a given system to find upwind fluxes and in a divergent manner, Eq. (19) can be written as ∇ · hi = 0. (20) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 7 of 21 Figure 1: Computational mesh and the arrangement of solution points. The vector h = (fi, qi) represents the space-time flow. Differential equations can be represented as integrals using the Gauss divergence theorem.∮ S(V ) hi · dS = 0, i = 1, 2, 3 , (21) which applies to all enclosed time and space regions V . In this case, S(V ) is the boundary of V , and ds ≡ dσn, where dσ and n ≡ (nx, nt) are the corresponding boundary element’s length and unit outward normal vector on S(V ), respectively. For a simulation, every solution point has a conservation element and a solution element, which are separated in time. Fig. 2 depicts distinct meanings for the central CE/SE scheme [18] and the upwind Figure 2: Solution element associated with solution point (j,n) CE/SE scheme [7]. Although the definitions of CE are similar, they differ somewhat. CEn j has two counterparts in both definitions: rectangle AEDC represents CEn,− j and BCFB represents CEn,+ j . The scheme’s SEn j is represented as rectangle EAFB rather than solid intersecting lines at (j, n). It is assumed that SEn j , u, and fi are determined using first- order Taylor expansions (stepwise second-order linear technique). The CE/SE system uses two marching variables at each result point which are (qi) n j and (qix) n j . These variables are numerically the same as qi and ∂qi/∂x at (j, n). To construct the time marching technique for (qi) n j and (qix) n j , where i = 1, 2, 3 we apply the conservation law expressed by Eq. (19) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 8 of 21 Figure 3: Space-time flux through the boundaries of sub-CEs. on CEn,− j and CEn,+ j , in the following form,[ (qi) n j − (qix) n j ∆x 2 ] ∆x 2 = UiL ∆x 2 + (FiL − FiC) ∆t 2 , (22) [ (qi) n j + (qix) n j ∆x 2 ] ∆x 2 = UiR ∆x 2 + (FiC − FiR) ∆t 2 . (23) Where UiL, UiR, FiL, FiR, and FiC represent the average flow (u is viewed as the flow of time) via AE, EB, AD, BC, and EF. Combining Eqs. (22) and (23), the time marching approach for qi n j may be driven as, (qi) n j = 1 2 (UiL + UiR) + ∆t 2∆x (FiL − FiR). (24) And subtracting Eq. (22) from Eq. (23), one might acquire the time marching scheme for (qix) n j in the form (qix) n j ∆x 4 = 1 2 (UiR − UiL) + ∆t 2∆x (2FiC − FiL − FiR). (25) AE and AD, EB and BC are linked with SE n−1/2 j−1/2 and SE n−1/2 j+1/2 , respectively. In CE/SE and upwind CE/SE schemes, UiL and FiL and UiR and FiR are determined through the utilization of Taylor expansions in SE n−1/2 j−1/2 and SE n−1/2 j+1/2 ,i.e., UiL = (qi) n−1/2 j−1/2 + (qix) n−1/2 j−1/2 ∆x 4 , FiL = (fi) n−1/2 j−1/2 + ∆t 4 (fit) n−1/2 j−1/2 . (26) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 9 of 21 Figure 4: Conservation element and solution element for the upwind method. UiR = (qi) n−1/2 j−1/2 + (qix) n−1/2 j−1/2 ∆x 4 , (27) FiR = (fi) n−1/2 j−1/2 + ∆t 4 (fit) n−1/2 j−1/2 . (28) The distinction is only in the computation of FiC which is the flow across the internal boundary of CEn j (EF). EF is positioned between SE n−1/2 j−1/2 and SE n−1/2 j+1/2 , in the upwind conservation element and solution element techniques. The flow through EF can compute FiC and could exhibit irregular patterns. FiC = h(qi n−1/4,− j , u n−1/4,+ j ), (29) where qi n−1/4,− j = qi n−1/2 j−1/2 + ∆x 2 (qix) n−1/2 j−1/2 + ∆t 4 (qit) n−1/2 j−1/2 , qi n−1/4,+ j = qi n−1/2 j+1/2 ∆x 2 (qix) n−1/2 j+1/2 + ∆t 4 (qit) n−1/2 j+1/2 . (30) A distinct framework for producing second-order schemes is provided by equations (24) and (25). Any effective flow method in the traditional FVM can be used to calculate FiC . Since FiC does not appear directly in Eq. (25) and does not affect the automated time-stepping algorithm of qi, the only difference between the ’a’ CE/SE scheme and the upwind CE/SE scheme is that the automatic time-marching scheme of qix differs only in relation to FiC . In Eqs. (24) and (25), UiL, FiL, UiR, and FiR are determined using the Taylor expansion S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 10 of 21 method outlined earlier. However, to find FiC , we need to solve the local Riemann problem. qi(x, 0) = { WiL, if x < x0 WiR, if x > x0. (31) is calculated at x−x0 t = 0. The point x0 is situated at the inner boundary of CEn j , and the flow variable values that are known at the midpoint on both sides of this inner boundary are represented by WiL and WiR. Both states are determined using Taylor expansion in SE n−1/2 j−1/2 and SE n−1/2 j+1/2 , as shown below: WiL = (q̄i)L + (q∗ix)L + (q∗it)L, (32) WiR = (q̄i)R + (q∗ix)R + (q∗it)R. (33) These represent the average values on the left side (AC) and the right side (BC), respectively. Finally, any efficient Riemann solver is used to solve the local Riemann problem described by Eqs. (19) and (32). To prevent unwanted oscillations near dis- continuities, the spatial derivatives are reconstructed using the weighted biased averaging process (WBAP-L2) limiter as shown below.[20], (q∗ix)L = (qix) n−1/2 j−1/2W (1, θ1L, θ 2 L), (q∗ix)R = (qix) n+1/2 j−1/2W (1, θ1R, θ 2 R) (34) where θ1L = qCix (qix) n−1/2 j−1/2 , θ2L = (qix) n+1/2 j−1/2 (qix) n−1/2 j−1/2 and θ1R = qCix (qix) n+1/2 j−1/2 , θ2R = (qix) n−1/2 j−1/2 (qix) n+1/2 j−1/2 (35) with qCix = q̄iR−q̄iL ∆x 2 , and the temporal derivatives are calculated by using the chain rule as (q∗it)L = ∂fi(∂q̄i)L ∂qi (36) (q∗it)R = ∂fi(∂q̄i)R ∂qi . (37) The averaged values on (Ū)iL the left side AB and on the (Ū)iR right side BC (see Fig. (4)) are defined as (Ū)iL = (qi) n−1/2 j−1/2 + (qix) n−1/2 j−1/2 ∆x 4 , (Ū)iR = (qi) n−1/2 j+1/2 − (qix) n−1/2 j+1/2 ∆x 4 . (38) S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 11 of 21 If qix = 0, the upwind CE/SE scheme of the 2nd order decomposes into a Lax scheme. To avoid spurious oscillations in situations with strong irregularities, the derivatives in Eq. (30) need an appropriate limiter [20]. In its simplest form, it is written as WBAPL2(1, θ1, · · · , θj) =  m+ ∑j j=1 1/θj m+ ∑j j=1 1/θ 2 j , if θ > 0 0, else (39) where m ≥ 1 is the linear weight of the initial derivatives. The first derivatives have a linear weight of m ≥ 1. Since the WBAP limiter is easy to use and effective, it is also be used in the 2D case. In conclusion, we summarize the features of our proposed scheme as: (1) The majority of concepts and techniques are similar to the first ”a” scheme, includ- ing:(i) Space and time are unified and treated equally in the scheme’s formulation; (ii) The Taylor expansion coefficients qij and (qix)j are considered as independent variables, implying that (qix)j is not obtained from qij using a FDM approximation; and (iii). The technique uses the simplest stencil, updating mesh variables only based on neighboring nodes. (2) The dissipation source in the described CE/SE method comes from the upwind flux, which is essentially different from the original CE/SE approach. While the update of qi n j stays the same as in the central CE/SE method,but only (qix) n j includes the upwind flux in its time-marching scheme. Therefore, to compute (qix) n j , just a few lines of the standard CE/SE code need to be modified. (3) Although the dissipation mechanism of the present scheme is consistent with up- wind schemes, it is not consistent with traditional upwind techniques. In a typical non- staggered upwind scheme, two upwind fluxes through the cell interfaces are used to build the time-marching scheme for the cell average. The derivative (qix) n j is calculated using qi n j with a finite difference approximation. In the present scheme, however, Taylor expan- sion is used to determine the other fluxes (FiL and FiR), and only the flux through the inner boundary of CEn j (FC) is calculated using an upwind procedure. Furthermore, an independent time-marching scheme updates the derivative. (4) Even though a staggered mesh is used, the current scheme remains largely unaf- fected by the CFL number, a characteristic that holds true even when using a limiter. As a result, the upwind CE/SE schemes can be effectively applied to highly refined or adaptive meshes where the local CFL number can vary greatly (from nearly 1 to less than 10−4). 5. Numerical test problems Several numerical test problems from literature are considered here to check the per- formance of upwind CE/SE scheme. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 12 of 21 Problem 1: 1D dam break problem. The issue addressed in this problem is the dam break over the flat bottom (B ≡ 0) in the computational space ([-1,1]), as mentioned in [21]. The following initial data was used to obtain the solution: (q, u, θ) = { (5, 0, 3) if x < 0 (1, 0, 5) if x > 0 We have calculated the results q = B + h, u, and θ by the upwind CE/SE scheme at t = 0.2. Note that from Fig. 5, the results by upwind CE/SE are in close agreement with the reference solution. Clearly, the scheme offers better resolution for sharp edges. Futhermore, we have also computed the L1 error of both the schemes shown in the table 1 given below: Table 1: Comparison of L1-errors in the the schemes. Height (h) Velocity (u) N Upwind CE/SE CE/SE Upwind CE/SE CE/SE 50 0.0102 0.3020 0.0341 0.0723 100 0.0051 0.0160 0.0113 0.0375 200 0.0026 0.0081 0.0069 0.0194 400 0.0013 0.0041 0.0035 0.0102 800 6.10× 10−4 0.0016 0.0015 0.0043 1600 2.80× 10−4 5.33× 10−4 8.62× 10−4 4.82× 10−4 Problem 2: 1D symmetrical dam break. The initial data for this problem are given below [22]: (h, u, θ) = { (2, 0, 1), if |x| ≤ 0.5, (1, 0, 1.5), otherwise. The source term that shows up in the Ripa model disappears when the flat waterbed is taken into account. Fig. 6 displays the results of h, u, and θ. We can observe from the graphs that upwind CE/SE scheme resolves the sharp discontinuities effectively. Problem 3: 1D dam break over a varied bottom The problem is solved using the initial data provided in [21] for various bottom con- figurations of the Ripa model. (q, u, θ) = { (5, 0, 1), if x < 0, (1, 0, 5), if x > 0. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 13 of 21 -1 -0.5 0 0.5 1 x-axis 2.5 3 3.5 4 4.5 5 5.5 temperature CE/SE Upwind CESE reference -1 -0.5 0 0.5 1 1.5 x-axis 0 0.5 1 1.5 2 2.5 velocity CE/SE Upwind CESE reference -1 -0.5 0 0.5 1 x-axis 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 height CE/SE Upwind CESE reference Figure 5: Problem 1: Comparison of upwind CE/SE with simple CE/SE scheme at time t = 0.2. The bottom topography function, which is varying, is defined as: B(x) =  2.0(cos(10π(x+ 0.3)) + 1), if − 0.4 ≤ x ≤ −0.2, 0.5(cos(10π(x− 0.3)) + 1), if 0.2 ≤ x ≤ 0.4, 0, otherwise. This problem aims to measure upwind CE/SE performance compared to an original CE/SE scheme for non-homogeneous shallow water equations. Initially, near x = 0.3, h is dry since h = q −B(x), and h(x, 0) = 1− 0.5(cos(10π(x− 0.3)) + 1). The results are displayed in Fig.7. The solution, on the other hand, reveals that the scheme is well-balanced and keeps q, h, and θ positive. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 14 of 21 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 x-axis 1 1.2 1.4 1.6 1.8 2 h height CE/SE Upwind CESE reference -1 -0.5 0 0.5 1 x-axis 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 temperature CE/SE Upwind CESE reference -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 1.2 1.4 1.6 1.8 2 p re ss u re Pressure CE/SE central reference Figure 6: Problem 2: The breaking of dam on the level waterbed. . Problem 4: Small perturbations in a steady state condition In our next 1D experiment, we present a small disturbance of a steady-state solution [21]. Here is the initial data: (q, u, θ) = { (6, 0, 4), if x < 0, (4, 0, 9), if x > 0. Non-flat bottom topography is described as follows: B(x) =  0.85(cos(10π(x+ 0.9)) + 1), if − 1.0 ≤ x ≤ −0.8, 1.25(cos(10π(x− 0.4)) + 1), if 0.3 ≤ x ≤ 0.5, 0, otherwise. The domain [−2, 2] is divided into 200 points and the results are generated at t = 0.4. Since the problem is a smooth pulse divided into two pulses advancing in opposite directions, the subsequent families will be perturbed at various times. The first pulse exits through the secondary hump of the bottom after traveling right to the first bump of the bottom S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 15 of 21 -1 -0.5 0 0.5 1 x-axis 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 height CE/SE Upwind CESE reference -1 -0.5 0 0.5 1 x-axis 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 temperature CE/SE Upwind CESE reference -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 x-axis -0.5 0 0.5 1 1.5 2 2.5 3 3.5 velocity CE/SE Upwind CESE reference Figure 7: Problem 3: 1D Dam break over variable bottom. cover and the temperature jump. Fig. 8 shows the numerical outcomes of this problem. As shown in this picture, the solutions within the CE/SE derive, explaining that this scheme is more accurate than CE/SE scheme. Our proposed numerical scheme seems to have improved accuracy in the results. Problem 5: A one-dimensional Dam Break Scenario with a Rectangular Elevation The initial data for this problem is considered from [22]. The waterbed of the rectan- gular hump is given as: B(x) = { 8, if |x− 300| < 75, 0, otherwise, and the initial data are given by: (h, u, θ) = { (20−B(x), 0, 1), if x ≤ 300, (15−B(x), 0, 5), otherwise. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 16 of 21 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x-axis 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 height CE/SE Upwind CESE reference -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x-axis 3 4 5 6 7 8 9 10 temperature CE/SE Upwind CESE reference -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x-axis -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0.1 velocity CE/SE Upwind CESE reference Figure 8: Problem 4: Small perturbations in a steady state condition Fig.9 indicates the computational results of the numerical approximations within the com- puting region at t = 12 acquired by using the upwind CE/SE scheme. The figure displays the graphs of h, θ, u, and B. From the results, one can see that both numerical schemes can capture the sharp variations within the solution, however, the upwind CE/SE scheme resolve the sharp peaks effectively. Problem 6: Problem of Rectangular Breakage of Dam In [22], the 2-D classical rectangular dam break problem was studied. According to initial circumstances, the two reserved states are: (h, u, v, θ) = { (2, 0, 0, 1), if |x| ≤ 0.5, (1, 0, 0, 1.5), otherwise. Various case studies were performed to validate and observe our proposed scheme’s well- balancing property. There is a near match in numerical results between the two schemes. Nevertheless, we saw that, compared to the simple CE/SE scheme, the upwind C/ESE S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 17 of 21 0 100 200 300 400 500 600 x-axis 8 10 12 14 16 18 20 22 height CE/SE Upwind CESE reference 0 100 200 300 400 500 600 x-axis 4 5 6 7 8 9 10 11 temperature CE/SE Upwind CESE reference 0 100 200 300 400 500 600 x-axis -1 0 1 2 3 4 5 velocity CE/SE Upwind CESE reference Figure 9: Problem 5: One-dimensional dam break problem over the rectangular bump. Comparison of upwind CESE with simple CE/SE scheme at time t = 12. scheme captures sharp discontinuities better and has similar accuracy to other schemes from the literature. Moreover, the technique preserves the stable characteristic of the numerical problems of pressure oscillations over variable bottom topography. 6. Conclusions An upwind space-time conservation element and solution element (CESE) scheme was proposed for 1D and 2D Ripa model in rectangular coordinates, combining the advantages of the CESE method and the upwind scheme. This approach ensures strict adherence to the space-time conservation law while efficiently capturing discontinuities. Various upwind schemes can be flexibly integrated for different problems, achieving an optimal combination of the CESE method and the finite volume method (FVM). The benchmark test results from the literature demonstrate that we have successfully extended the upwind CESE scheme to the Ripa model. The outcomes of the proposed numerical scheme were validated by comparing them with those obtained from a high-resolution simple CE/SE scheme for verification.It is observed that upwind CE/SE scheme performs better than S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 18 of 21 CE/SE scheme. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 19 of 21 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x-axis 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 h Height CE/SE Upwind CESE reference 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x-axis 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5 te m p er at u re Temperature CE/SE Upwind CESE reference 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 x-axis -0.3 -0.2 -0.1 0 0.1 0.2 0.3 Velocity CE/SE Upwind CESE reference Figure 10: Problem 6: Rectangular dam break problem. 3D views of solution components (left) and their 1D comparison (right) calculated using upwind CE/SE and basic CE/SE methods at time t = 0.2 employing a uniform mesh grid ∆x =∆y= 2/100. S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 20 of 21 References [1] R. Y. Adrian, Norlizan W., Lee H. P., A. H. Z. Fariz, I. Ruqayyah, L. D. Goh, and A. Hazrina. Beam (cssb) at elevated temperature. Journal of Engineering Sciences and Technology in Physics, 19:25–39, 2024. [2] B. Al-Hadeethi, A. S. Almawla, A. H. Kamel, H. A. Afan, and A. N. Ahmed. Numeri- cal modeling of flow pattern with different spillway locations. Mathematical Modelling of Engineering Problems, 19:1219–1226, 2024. [3] A. Dehghanghadikolaei, N. Namdari, B. Mohammadian, and S. R. Ghoreishi. Deriv- ing one dimensional shallow water equations from mass and momentum balance laws. International Research Journal of Engineering and Technology, 5(6):407–419, 2018. [4] A. J. C. de Saint-Venant. Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l’introduction de marées dans leurs lits. Comptes Rendus des Séances de l’Académie des Sciences, 73(147):237–240, 1871. [5] A. Rehman, I. Ali, and S. Qamar. Exact riemann solutions of the ripa model for flat and non-flat bottom topographies. Results in Physics, 8:104–113, 2018. [6] J. Britton and Y. Xing. High-order still-water and moving-water equilibria preserve discontinuous galerkin methods for the ripa model. Journal of Scientific Computing, 82(2):30, 2020. [7] H. Shen, C. Y. Wen, and D. L. Zang. A characteristic space-time conservation ele- ment and solution element method for conservation laws. Journal of Computational Physics, 288:101–118, 2015. [8] S. B. Savage and K. Hutter. The dynamics of avalanches of granular materials from initiation to run-out. parti: Analysis. Acta Mechanica, 86:201–223, 1991. [9] S. Qamar and G. Warnecke. Application of space-time cese method to shallow-water magnetohydrodynamic equations. Journal of Computational and Applied Mathemat- ics, 196:132–149, 2005. [10] H. Z. Tang and H. M. Wu. Kinetic flux vector splitting schemes for the radiation hydrodynamical equations. Computers & Fluids, 29:917–933, 2000. [11] V. Coralic and T. Colonius. Finite-volume weno scheme for viscous compressible multi-component flows. Journal of Computational Physics, 274:95–121, 2014. [12] V. Caleffi, A. Valiani, and A. Bernini. Fourth-order balanced source-term treatment in central weno schemes for shallow water equations. Journal of Computational Physics, 218:228–245, 2006. [13] C. Y. Wen, Y. Jiang, and L. Shi. Space-time conservation element and solution element method. Advances and Applications in Engineering Sciences, 218(1):139, 2023. [14] L. Lundgren and K. Mattsson. An efficient finite difference method for shallow water equations. Journal of Computational Physics, 422:109–784, 2020. [15] F. L. B. Ribeiro, A. C. Galeao, R. G. S. Castro, and L. Landau. Finite elements for shallow water equations: stabilized formulations and computational aspects. WIT Transactions on Engineering Sciences, 29:167–179, 1970. [16] S. Chippada, C. N. Dawson, M. L. Mart́ınez, and M.F. Wheeler. A godunov-type S. Zia et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 6230 21 of 21 finite volume method for the system of shallow water equations. Computer Methods in Applied Mechanics and Engineering, 151(1–2):105–129, 1998. [17] T. Gallouët, J. M. Hérard, and N. Seguin. Some approximate godunov schemes to compute shallow-water equations with topography. Computers & Fluids, 32(4):479– 513, 2003. [18] S. C. Chang. The method of space-time conservation element and solution element—a new approach for solving the navier–stokes and euler equations. Journal of Compu- tational Physics, 119(2):295–324, 1995. [19] M. Saleem and S. Qamar. The space-time cese scheme for shallow water equations incorporates variable bottom topography and horizontal temperature gradients. Com- puters & Mathematics with Applications, 75(3):933–956, 2018. [20] W. Li, Y. X. Ren, G. Lei, and H. Luo. The multi-dimensional limiters for solving hy- perbolic conservation laws on unstructured grids. Journal of Computational Physics, 230(21):7775–7795, 2011. [21] A. Chertock, A. Kurganov, and Y. Liu. Central-upwind schemes for the system of shallow water equations with horizontal temperature gradients. Numerische Mathe- matik, 127(4):595–639, 2014. [22] R. Touma and C. Klingenberg. Well-balanced central finite volume methods for the ripa system. Applied Numerical Mathematics, 97:42–68, 2015.