Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2167 
https://internationalpubls.com 

Artificial Viscosity Stabilization of Meshless Methods for Shallow 

Water Flows in Channels 

Mariam Ijerch1, Mohamed Sadik2 

1Engineering Sciences Laboratory, Faculty of Sciences, Agadir, Ibn Zohr University, Morocco 

(Email: mariam.ijerch@edu.uiz.ac.ma) 

2Engineering Sciences Laboratory, Faculty of Sciences, Agadir, Ibn Zohr University, Morocco 

(Email: m.sadik@uiz.ac.ma) 

 

Article History: 

Received: 12-01-2025 

Revised: 15-02-2025 

Accepted: 01-03-2025 

Abstract:  

The purpose of this paper is to present an efficient localized meshless methods based on radial 

basis functions to accurately analyze the shallow water equations (SWEs) in open channels. 

It is a hyperbolic system of first-order nonlinear partial differential equations. Due to their 

classification as advection equations, they allow for discontinuous solutions, often 

characterized by shock behaviour in simulations. Therefore, the development of an efficient 

and accurate numerical model for analyzing the SWEs holds critical importance in scientific 

research. To address non-physical oscillations near discontinuities, we developed a technique 

involving artificial viscosity, which is integrated with the local meshless methods for spatial 

discretization of the SWEs, while temporal discretization is achieved using the fourth-order 

Runge-Kutta method. A series of experiments has been conducted to evaluate the 

effectiveness and accuracy of the suggested methods. These experiments included examining 

flow at a steady state with and without friction, as well as studying the dam-break scenario 

on a wet bed. The results compared with analytical solutions, demonstrate that the artificial 

viscosity combined with the proposed methods proficiently capture shocks and effectively 

handles discontinuous flow by introducing an appropriate viscosity coefficient into the 

equations. Overall, the results are satisfactory and show good agreement. 

Keywords: Shallow water equations, Open channel, Radial basis function partition of unity 

method (RBF-PUM), Radial basis function finite difference (RBF-FD) method, QR 

factorization. 

 

1. Introduction

The shallow water equations (SWEs) were first proposed in 1871 by Saint-Venant [1] to describe 

various applications in ocean and hydraulic engineering, such as coastal flow, open-channel flow in 

natural rivers and prismatic channels, flood waves, failure of dams, etc. Under some essential 

assumptions [2], the SWEs are derived by three-dimensional Navier-Stokes equations where the 

horizontal wave length is much larger than the depth of the fluid. SWEs form a system of non-linear 

hyperbolic partial differential equations (PDEs) governed by conservation of mass and balance of 

momentum. The critical problem in solving this system of equations is the spurious oscillation happens 

in the discontinuous areas such as shock waves, dam break, or hydraulic jump problems.   The 

analytical solution for these problems is only available for a limited number of special cases [3]. 

Therefore it is significant to develop a reliable and accurate numerical solver.  

mailto:mariam.ijerch@edu.uiz.ac.ma
mailto:m.sadik@uiz.ac.ma


Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2168 
https://internationalpubls.com 

In the past decades, plenty of mesh-based numerical models have been developed for solving the 

SWEs, such as the Finite Difference Method (FDM) [4], Finite Element Method (FEM) [5], and Finite 

Volume Method (FVM) [6, 7]. However, for multidimensional problems defined on a complex 

geometries, mesh generation have a common difficulty that it is a time-consuming and their 

performance depends on the quality of the mesh. To overcome these difficulties, many researchers 

have paid attention to use the newly developed meshless method as the promising alternative to 

traditional methods which have been employed to analyse various problems. The meshless methods 

were also adopted to solve SWEs by many scholars. Among them, the development of meshless 

methods based on Radial Basis Functions (RBFs), such as multiquadric method [8], Compactly 

Supported RBF (CSRBF) [9], Local Radial Basis Function based Differential Quadrature Method 

(LRBFDQM) [10], Extrapolated Local Radial Basis Functions Collocation Method (ELRBFCM) [11], 

and others. The main advantages of RBF-based methods for solving partial differential equations lie 

in their simplicity, flexibility, and applicability to various PDEs. They are particularly effective in 

handling high-dimensional problems with complex geometries. 

Since the shallow water equations are a hyperbolic non-linear partial differential equations, the 

information of characteristic and the correct directions of wave transmission are quite important to 

numerical simulation. The first author who introduced the radial basis functions method in the field of 

PDEs was Kansa in 1990 [12, 13]. He developed the global RBF method which is a numerical method 

that gives an approximate solution for PDEs. It can be considered as one of the collocation methods 

which are efficient and expected to produce spectral accuracy for the numerical solution. However, the 

infinitely smooth RBFs (Table 1) produce dense matrices thus ill-conditioned matrices, especially for 

the parameter values that lead to the highest accuracy. To overcome this issue, it is essential to 

introduce certain locality in the system to be solved. Not only the conditioning will be improved, but 

also the accuracy of the numerical solution will be guaranteed. In this context, we conducted our 

studies using local and stable meshless methods in addition to Kansa's method to solve SWEs. One 

such method is the RBF Finite Difference method (RBF-FD),  first described in [14] in 2000.  The idea 

of this method is to form computational stencils, similar to those in finite difference method, by 

constructing derivative discretizations based on RBF interpolants on localized node sets [15].  Another 

localized method is the RBF Partition of Unity Method (RBF-PUM) which is based on approximating 

the solution of PDEs in a regular covering of overlapping subdomains and then patching up the local 

solutions using the partition of unity principle [16, 17]. The RBF-PUM, increasing the sparsity of linear 

systems, enables us to offer significant savings in the computational time, but in order to guarantee the 

stability of the solution, as 𝜀 ⟶ 0  (called flat limit regime), a stable evaluation algorithm is employed. 

The stable RBF-QR algorithm bypasses some limitations, including complications related to the 

selection of 𝜀, ill-conditioning of RBF-PUM matrices. The main logic behind the RBF-QR method is 

to replace the bad basis (ill-conditioned) with a good one that can span the same space [18, 19]. 

The RBF methods that we investigate are: the RBF generated finite difference(RBF-FD) method, the 

RBF partition of unity method (RBF-PUM) and Kansa's method. The first two methods are localized 

such that the final system of equations is sparse, and the third method is global and gives a dense 

system of equations. It is known, based on numerical experiments, that all three methods in general 

are unstable when computing collocated solutions to time-dependent conservation laws without 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2169 
https://internationalpubls.com 

additional stabilization [20]. In [21], the author suggests activating artificial viscosity in the time 

domain rather than in localized regions of space and using hyperviscosity provides an elegant approach 

to spatially selecting the diffusion effect based on supplementing the numerical scheme by a 

hyperviscosity term. Hyperviscosity was first introduced to spectral methods in 1998 by Ma under the 

name (super)spectral viscosity [22]. Despite the promising analytical results demonstrated by Ma, the 

method was never widely adopted and has only been mentioned in a few research papers where it was 

applied to the Navier-Stokes equations. The meshless community first adopted spectral viscosity, now 

called hyperviscosity, in 2011 when Fornberg [20] demonstrated that hyperviscosity could stabilize 

time-stepping schemes by shifting spurious eigenvalues into the stable region. Subsequently, 

hyperviscosity was applied to a variety of fluid dynamics cases [23, 24]. Initially, the parameters of 

hyperviscosity were only briefly mentioned and included a large number of empirical approximations. 

 Recently, new artificial hyperviscosity formulations have been proposed in the literature [25], which 

are very effective for RBF-FD methods for solving surface convection-diffusion equations and for 

nonlinear conservation laws [26]. additionally, meshless schemes combined with hyperviscosity 

technique have been proposed to the authors' knowledge such in [27] used the RBF-based method for 

solving the shallow water equations to capture shocks. Moreover, hyperviscosity and its connection to 

the stability of the RBF Partition of Unity (RBF-PU) method was only well established by the work of 

Liu et al. [28], who showed that the hyperviscosity improve the stability of the RBF-PU method for 

solving convection-diffusion equations on surfaces. Therefore, we suggested to combine the RBF-FD 

(using Inverse Multiquadratic function) and RBF-PUM-QR with hyperviscosty technique in order to 

establish a localized numerical meshless model which can accurately and effectively solve the SWEs. 

This paper is organized as follows. In Section 2, the governing equations and some properties are 

explained. The Kansa's method, RBF-FD method and RBF-PUM with QR factorization with 

Hyperviscosity technique is presented in Section 3. Then, two challenging tests are adopted to evaluate 

the accuracy and stability of the numerical methods in Section 4. Finally, some conclusions are given 

in Section 5. 

2. Governing equations and properties 

Saint-Venant [1] introduced the mathematical model for shallow water, known as Saint-Venant's 

Equations or Shallow Water Equations (SWEs). These equations are based on fundamental physical 

principles, specifically the laws of mass conservation and momentum conservation. The model is 

applicable to free surface flows when the horizontal scales of water mass significantly exceed the 

vertical scale and the flow in the vertical direction is negligible. The water density in our model is 

considered to be constant and Coriolis force is neglected. As the source terms the bed slope term and 

bed friction term are considered. The equations of interest consist of the following system:  

                                         {
𝜕𝑡ℎ + 𝜕𝑥𝑞 = 0,                                        

𝜕𝑡𝑞 + 𝜕𝑥 (
𝑞2

ℎ
+

1

2
𝑔ℎ2) = 𝑔ℎ𝑆𝑏 − 𝑆𝑓 ,

                                               (1)                                                                                                                                                           

where ℎ(𝑥, 𝑡) is the height of the water,  𝑢(𝑥, 𝑡) is the flow velocity, 𝑞(𝑥, 𝑡) = ℎ(𝑥, 𝑡)𝑢(𝑥, 𝑡) is the 

flow discharge and 𝑔  is the acceleration due to gravity. Concerning the source terms, 𝑆𝑏 = −𝜕𝑥𝑏 is 

the bed slope where 𝑏: ℝ → ℝ+ is a given smooth function which describes the topography and 𝑆𝑓 is 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2170 
https://internationalpubls.com 

the friction term, it is a given function of ℎ  and 𝑢 , two examples widely used  in hydrology are the 

Manning and the Darcy-Weisbach friction laws, which are given respectively by 

                                                        𝑆𝑓 = 𝐶𝑓𝑔
𝑢|𝑢|

ℎ
1
3

= 𝐶𝑓𝑔
𝑞|𝑞|

ℎ
7
3

,                                                                (2) 

where 𝐶𝑓  =  𝑛2,  𝑛 is Manning's coefficient. 

                                                        𝑆𝑓 = 𝐶𝑓𝑔𝑢|𝑢| = 𝐶𝑓𝑔
𝑞|𝑞|

ℎ2 ,                                                                      (3) 

where 𝐶𝑓 =
𝑓

8𝑔
, 𝑓 is dimensionless coefficient. In order to emphasize the properties of (1) without 

friction, we rewrite the one-dimensional equations in quasilinear form 

                                                        ∂t𝐐 + A(𝐐) ∂xQ = S(𝐐),                                                                       (4) 

with 

             𝑸 = (ℎ
𝑞

),            A(𝐐) = (
0 1

−𝑢2 + 𝑔ℎ 2𝑢
),            𝑆(𝑸) = ( 0

−gh ∂xb
).                       (5) 

When ℎ > 0, the system is strict hyperbolic with the eigenvalues of 𝐴 are given by 

                                       𝜆1 =  𝑢 − √𝑔ℎ        𝑎𝑛𝑑        𝜆2 =  𝑢 + √𝑔ℎ.                                            (6) 

These eigenvalues represent the wave speeds  which are basic characteristics of the flow. When ℎ = 0 

(dry zone) the eigenvalues coincide and the system looses hyperbolicity. For considering the well-

posed problem in 1D described by the SWEs, the number of boundary conditions should be confirmed 

by the Froude number 

                                                                    Fr =  
|𝑢|

√𝑔ℎ
.                                                                          (7) 

If Fr < 1 the flow is considered subcritical, indicating the presence of one positive and one negative 

eigenvalue. In such cases, it is necessary to define one boundary value on the left and another on the 

right. On the other hand, when Fr > 1 then the flow is supercritical which means that two values have 

to be specified at one of the boundaries depending on the sign of the velocity. A transcritical regime 

can exist in the solution if parts of the flow are subcritical and other parts are supercritical. 

2.1 Hyperviscosity technique 

Although the global stability theory of RBF-FD and RBF-PUM methods is largely undeveloped, 

Fornberg and Lehto [20] proposed a practical method to correct the spectra of RBF-FD differentiation 

matrices for hyperbolic operators on the sphere by introducing artificial hyperviscosity. Their approach 

effectively shifted rogue eigenvalues to the left half of the complex plane. This issue is also caused by 

using large RBF-FD stencils, where increasing the stencil size 𝑛𝑙𝑜𝑐  enhances the accuracy of the 

approximation. Based on the above analysis, we will use the Kansa, RBF-FD, and RBF-PUM methods 

combined with an artificial hyperviscosity formulation to solve the shallow water equations, where the 

convection term dominates and leads to non-physical oscillations that would otherwise amplify with 

time. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2171 
https://internationalpubls.com 

Stabilization of the methods is achieved by applying a hyperviscosity as a filter to the right hand side 

of the model, the equation solved takes the form 

                                                        

                                                   𝜕𝑡𝑸 =  −𝐹(𝑸) + 𝑆(𝑸) + 𝐻(𝑸),                                                        (8) 

 

where 

 

              𝑸 = (ℎ
𝑞

),               F(𝐐) = (
𝜕𝑥𝑞

𝜕𝑥(
𝑞2

ℎ
+

1

2
𝑔ℎ2)

),              𝑆(𝑸) = ( 0
𝑔ℎ𝑆𝑏−𝑆𝑓

),                            (9) 

 

and 𝐻 is the hyperviscosity filter defined by 

 

                                                             𝐻 = 𝜇
𝜕2𝑘

𝜕𝑥2𝑘 ,                                                                        (10) 

 

where 𝜇 ∈ 𝑅 and 𝑘 ∈ 𝑁∗ small numbers that must be selected. The constant μ has to be selected 

carefully to provide enough stabilisation without ruining the solution [20, 23]. There is no broad 

consensus in the literature with an abundance of proposed estimates and scalings for the constant μ. 

However, there are many studies suggesting that the hyperviscosity parameter should be chosen 

through trial and error. Fornberg and Lehto [20] provide suggestions for selecting μ in the RBF-FD 

method, they found scaling 𝜇 ~ 𝑁−2𝑘 worked well, where N is the total number of nodes in the domain. 

Then, in [23], Flyer et al. empirically computed μ and k for the shallow water equations on the sphere. 

They found experimentally that for given values of N, 𝑛𝑙𝑜𝑐 and 𝜀 , choosing 𝜇 = 𝜇𝑐𝑁−𝑘  where 

𝜇𝑐 ranging from O(1) to O(10−2) provides stability with good accuracy for PDEs. 

 

In this paper, we experimentally determine the hyperviscosity parameter through trial and error, the 

order k set to one, drawing on the work of Dehghan and Abbaszadeh [29], who used a value of 10−3 

to stabilize their method for the dam break problem. We also reference [27], where values ranging 

from 0.01 to 0.5 were employed with an order k set to one, and [30], where a small value around 10−5 

was deliberately chosen with the same order k. 

 

3. Numerical methods 

 

3.1 Kansa method 

 

First, we introduce the well-known Kansa method [12, 13]. It is also referred to as Radial Basis 

Function Collocation Method (RBFCM). The main idea of the method is approximating the solution 

in terms of linear combination of infinitely differentiable Radial Basis Functions (RBFs) 𝜑 which are 

listed in Table 1, this functions depends only on the distance to a center point and shape parameter 𝜀. 

This method is particularly useful for solving PDEs in complex geometries with irregularly distributed 

points, offering an efficient and flexible numerical technique. However, The main disadvantage of the 

method is that it represents the discretization of the PDE as a dense matrix. These matrices are 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2172 
https://internationalpubls.com 

extremely sensitive to the choice of the shape parameter for the RBFs, which makes it difficult to solve 

problems with a large number of unknowns. This occurs because in the use of RBF interpolation and 

add more nodes, the matrices involved become less stable. This is especially true when a poor choice 

has been made when designating the RBF centers, and also when infinitely smooth basis functions are 

used with associated shape parameters set to extreme values. To have more details on the theoretical 

background, we refer the reader to [35, 36]. 

  

Table 1.  Infinitely smooth RBFs 

 

Multiquadric               Inverse multiquadratic         Inverse quadratic         Gaussian 

(𝟏 + 𝝐𝟐𝒓𝟐)
𝟏
𝟐

                       (𝟏 + 𝝐𝟐𝒓𝟐)
−𝟏
𝟐

                       (𝟏 + 𝝐𝟐𝒓𝟐)−𝟏                    𝒆−𝝐𝟐𝒓𝟐
                              

 

3.1.1  Method’s principle 

 

In this section, we introduce the general notation and quantities we need for RBF approximation of 

SWEs. Assume that X =  {𝑥1, 𝑥2, . . . , 𝑥𝑁} is a set of distinct points in Ω ⊂  ℝ𝑑. We start from given an 

interpolant 𝐼𝑢 that approximates a function 𝑢: 

 

                                         𝐼𝑢(𝑥) =  ∑ 𝜆𝑖
𝑛
𝑖=1 𝜑(‖𝑥 − 𝑥𝑖‖),               𝑥 ∈  𝛺,                            (11) 

 

where 𝜆𝑖  ∈  ℝ are the unknown coefficient to be determined and 𝜑 can be any radial basis function 

such in Table 1, the shape parameter 𝜀 appearing in the RBFs dictates the flatness of the radial basis 

function and plays a key role in the convergence rate of the approximations and the condition number 

of the coefficient matrices. 

By applying the interpolation criteria 𝐼𝑢(𝑥𝑗)  =  𝑢(𝑥𝑗), the coefficient 𝜆𝑗 , 𝑗 =  1, . . . , 𝑁 is discovered. 

And we obtain a linear system as follow: 

 

A�̅� = 𝑈,                                                          (12) 

where 

        𝐴 =  (
𝜑(‖𝑥1 − 𝑥1‖) ⋯ 𝜑(‖𝑥1 − 𝑥𝑁‖)

⋮ ⋱ ⋮
𝜑(‖𝑥𝑁 − 𝑥1‖) ⋯ 𝜑(‖𝑥𝑁 − 𝑥𝑁‖)

),                       (13) 

 

                        �̅� =  [𝜆1, 𝜆2, … , 𝜆𝑁]𝑇,              U = [𝑢(𝑥1), 𝑢(𝑥2), … , 𝑢(𝑥𝑁)]𝑇.           (14) 

 

When solving PDE, we prefer to work with the discrete approximate instead of the coefficients. From 

(11) and (12) we can write  

 

                                    𝐼𝑢(𝑥) = �̅�(𝑥)𝐴−1𝑈 = 𝐷(𝑥)𝑈,                                       (15) 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2173 
https://internationalpubls.com 

where �̅�(𝑥)  =  [𝜑 (∥ 𝑥 −  𝑥1 ∥) , 𝜑 (∥ 𝑥 −  𝑥2 ∥) , … , 𝜑 (∥ 𝑥 −  𝑥𝑁 ∥)]  and 𝐷 = �̅�𝐴−1.  Micchlli 

[37] proved that the strict positive definiteness of an RBF listed in Table 1 guarantees that the 

interpolation matrix A in (13) is also positive definite, and thus invertible. Because our final target is 

to solve PDEs, we apply a linear differential operator 𝐿 to the RBF approximation, we get 

 

     (𝐼𝑢)𝐿(𝑥) = �̅�𝐿(𝑥)𝐴−1𝑈 = 𝐷𝐿(𝑥)𝑈,                                 (16) 

 

where 𝐷𝐿 =  �̅�𝐿𝐴−1 is known as a differentiation matrix. 

 

3.1.2  Spatial discretization of SWEs using Kansa method 

In this section, we will describe the implementation of Kansa method for shallow water equations. Let 

{𝑥1, 𝑥2, . . . , 𝑥𝑁} be a set of scattered nodes considered in the study area. Using (16), the first and second 

spatial derivatives of a function 𝑸(𝑥, 𝑡) can be expressed as follows: 

 

𝜕𝑸(𝑥,𝑡)

𝜕𝑥
= 𝐷𝑥𝑸(𝑥, 𝑡),          

𝜕2𝑸(𝑥,𝑡)

𝜕𝑥2 = 𝐷𝑥𝑥𝑸(𝑥, 𝑡),                        (17) 

 

where 𝐷𝑥  and 𝐷𝑥𝑥  are first and second derivative matrices with respect to 𝑥 . The system (8) is 

approximated by the following equivalent system: 

 

𝜕𝑡𝑸 =  −�̅�(𝑸) + 𝑆̅(𝑸) + �̅�(𝑸),                                        (18) 

 

where �̅� is the convective flux, �̅� contains the hyperviscosity and 𝑆̅ is the source term due to bottom 

elevator. Using the Kansa method, we obtain 

 

          𝐹 ̅(𝑸) = (
𝐷𝑥𝑞

𝐷𝑥(
𝑞2

ℎ
+

1

2
𝑔ℎ2)

),        𝑆 ̅(𝑸) = ( 0
−𝑔ℎ𝐷𝑥𝑏

),         𝐻 ̅̅ ̅(𝑸) = (𝜇𝐷𝑥𝑥ℎ
𝜇𝐷𝑥𝑥𝑞

),             (19) 

 

Equivalent to 

 

         𝜕𝑡𝑸 = 𝑅𝐻𝑆𝑘𝑎𝑛𝑠𝑎(𝑸)      where     𝑅𝐻𝑆𝑘𝑎𝑛𝑠𝑎(𝑸) = −�̅�(𝑸) + 𝑆̅(𝑸) + �̅�(𝑸).     (20) 

 

3.2  Radial Basis Function Finite Difference Method (RBF-FD) 

 

First introduced by Tolstykh [14], RBF-FD formulas are derived through RBF interpolation over local 

sets of nodes on the surface. This type of method is conceptually similar to the standard Finite 

Difference (FD). The fundamental principle underlying RBF-FD revolves around the approximation 

of the differential operator for solutions at interior nodes. The approximation is achieved by means of 

a linear combination of function values at the neighbouring node locations and then determining the 

RBF-FD weights by assuming that the approximations become exact for all RBFs that are centred at 

the neighbouring nodes, as opposed to the standard FD, which focuses on approximating polynomials 

for the same nodes sets. After performing calculations at all interior nodes, the approximate solution 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2174 
https://internationalpubls.com 

can be computed from the associated linear system of equations, for which the RBF-FD system matrix 

is sparse, thus can be effectively inverted. 

 

3.2.1 Method’s principle 

 

We will next present an outline of the RBF along with finite-difference (RBF-FD) formulation. Let 

X =  {𝑥1, 𝑥2, . . . , 𝑥𝑁} is a set of distinct points in Ω ⊂  ℝ𝑑 and let 𝐼(𝑥𝑗) = {𝑥1
𝑗 , 𝑥2

𝑗 , … , 𝑥𝑛𝑙𝑜𝑐

𝑗 } be a set of 

points to form a stencil weight at 𝑥𝑗  . Note that the number of points 𝑛𝑙𝑜𝑐  in each stencil is either 

constant or vary with 𝑗, without loss of generality we suppose that 𝑛𝑙𝑜𝑐 is constant. In RBF-FD, any 

linear differential operator 𝐿  that acts on 𝑢(𝑥)  evaluated at 𝑥𝑗  can be approximated by a linear 

weighted combination of the function values at the points of 𝐼(𝑥𝑗), i.e. 

 

              𝐿𝑢(𝑥𝑗) ≈  ∑ 𝜔𝑘
(𝑗)𝑛𝑙𝑜𝑐

𝑘=1 𝑢(𝑥𝑘
(𝑗)

).                                      (21) 

 

The RBF-FD weights  𝜔𝑘
(𝑗)

, 𝑘 =  1, … , 𝑛𝑙𝑜𝑐  are computed by assuming the approximations (21) 

become exact for all RBFs 𝜑 (Table 1) that are centred at the nodes 𝑥𝑖
(𝑗)

, 𝑖 =  1, … , 𝑛𝑙𝑜𝑐 , i.e., we 

assume that (21) becomes exact for function 𝑢(𝑥𝑗)  =  𝜑(∥ 𝑥 − 𝑥𝑖
(𝑗)

∥)|
𝑥=𝑥𝑗

 , 𝑖 =  1, … , 𝑛𝑙𝑜𝑐 , this 

assumption leads to 𝑛𝑙𝑜𝑐 × 𝑛𝑙𝑜𝑐 linear system: 

 

(

𝜑 (‖𝑥1
(𝑗)

− 𝑥1
(𝑗)‖) ⋯ 𝜑 (‖𝑥1

(𝑗)
− 𝑥𝑛𝑙𝑜𝑐

(𝑗) ‖)

⋮ ⋱ ⋮

𝜑 (‖𝑥𝑛𝑙𝑜𝑐

(𝑗)
− 𝑥1

(𝑗)‖) ⋯ 𝜑 (‖𝑥𝑛𝑙𝑜𝑐

(𝑗)
− 𝑥𝑛𝑙𝑜𝑐

(𝑗) ‖)

) (
𝜔1

(𝑗)

⋮

𝜔𝑛𝑙𝑜𝑐

(𝑗)
) = (

 𝐿𝜑(∥𝑥 − 𝑥1
(𝑗)

∥)|
𝑥=𝑥𝑗

 

⋮

 𝐿𝜑(∥𝑥 − 𝑥𝑛𝑙𝑜𝑐

(𝑗)
∥)|

𝑥=𝑥𝑗

)(22) 

 

Note that the system matrix is invertible due to the properties of the RBFs mentioned before. Therefore, 

the RBF-FD weights can always be calculated. To obtain the RBF-FD weights, we have to solve the 

system (22) for each stencil center 𝑥𝑗  , 𝑗 =  1, … , 𝑁  to form 𝑁  rows of the differentiation matrix 

denoted by 𝑊𝐿 where 𝐿 refer to a linear differential operator that acts on solutions, 𝑊𝐿 is a sparse 

matrix and contains 𝑛𝑙𝑜𝑐 non-zero elements per row. When solving PDE, we apply a similar procedure 

of what has done for Kansa method. After obtaining the weights, we use the approximation equation 

(21) to conclude an approximate solution for the PDE. 

 

3.2.2  Spatial discretization of SWEs using RBF-FD method 

 

In this section, we will present implementation details on the computational of RBF-FD formulation 

for solving equations (1). Let {𝑥1, 𝑥2, . . . , 𝑥𝑁} be a set of scattered nodes. Using (21), the first and 

second spatial derivatives of a function 𝑸(𝑥, 𝑡)  can be expressed as follows: 

 

𝜕𝑸(𝑥𝑗,𝑡)

𝜕𝑥
= ∑

𝜕 𝜔𝑘
(𝑗)

𝜕𝑥

𝑛𝑙𝑜𝑐
𝑘=1 𝑸(𝑥𝑘

(𝑗)
, 𝑡),   

𝜕2𝑸(𝑥𝑗,𝑡)

𝜕𝑥2 = ∑
 𝜕2𝜔𝑘

(𝑗)

𝜕𝑥2

𝑛𝑙𝑜𝑐
𝑘=1 𝑸(𝑥𝑘

(𝑗)
, 𝑡), 𝑗 = 1, … , 𝑁  (23) 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2175 
https://internationalpubls.com 

 

After collocating the points, we obtain 

 

                                   
𝜕𝑸(𝑥,𝑡)

𝜕𝑥
= 𝑊𝑥𝑸(𝑥, 𝑡),          

𝜕2𝑸(𝑥,𝑡)

𝜕𝑥2 = 𝑊𝑥𝑥𝑸(𝑥, 𝑡),                           (24) 

 

where 𝑊𝑥  and 𝑊𝑥𝑥  are first and second derivative matrices with respect to 𝑥 . The system (8) is 

approximated by the following equivalent system: 

 

                                     𝜕𝑡𝑸 =  −�̅�(𝑸) + 𝑆̅(𝑸) + �̅�(𝑸),                                        (25) 

 

Using the RBF-FD method to perform numerical evaluation, we obtain 

 

          𝐹 ̅(𝑸) = (
𝑊𝑥𝑞

𝑊𝑥(
𝑞2

ℎ
+

1

2
𝑔ℎ2)

),        𝑆 ̅(𝑸) = ( 0
−𝑔ℎ𝑊𝑥𝑏

),         𝐻 ̅̅ ̅(𝑸) = (𝜇𝑊𝑥𝑥ℎ
𝜇𝑊𝑥𝑥𝑞

),             (26) 

 

Equivalent to 

 

       𝜕𝑡𝑸 = 𝑅𝐻𝑆𝑅𝐵𝐹−𝐹𝐷(𝑸),  where  𝑅𝐻𝑆𝑅𝐵𝐹−𝐹𝐷(𝑸) = −�̅�(𝑸) + 𝑆̅(𝑸) + �̅�(𝑸).     (27) 

3.3 Radial Basis Function Partition of Unity Method (RBF-PUM) 

 

The Partition of Unity Method (PUM) was first proposed by Babuška and Melenk [31] and applied the 

method to solve PDEs. The main idea of RBF-PUM is to subdivide the domain into overlapping 

subdomains or patches in which we construct the local RBF approximation on each patches, then the 

global solution is given by the sum of these local approximations multiplied by partition of unity 

weight functions. 

 

3.3.1 Method’s principle 

 

Let {Ω𝑖}
𝑖=1

𝑛𝑝
 be an open covering of an open set 𝛺 i.e., 𝛺 ⊆ ⋃ 𝛺𝑖

𝑛𝑝

𝑖=1
. The RBF-PUM is a localized 

method based on subdividing the domain 𝛺 on np patches 𝛺1, … , 𝛺𝑛𝑝
. First, we define a partition of 

unity functions  {𝜔𝑖}𝑖=1

𝑛𝑝
 subordinated to the covering  {Ω𝑖}

𝑖=1

𝑛𝑝  such that 

 

                                                ∑ 𝜔𝑖(𝑥) = 1,               𝑥 ∈ Ω,
𝑛𝑝

𝑖=1
                                         (28) 

 

where the weight function 𝜔𝑖 ∶  𝛺𝑖  →  ℝ is compactly supported, non-negative and continuous. For 

each patch we may thus construct a local RBF interpolant 𝐼𝑢
𝑖 ∶  𝛺𝑖  →  ℝ of the form: 

 

                                                 𝐼𝑢
𝑖 (𝑥) =  ∑ 𝜆𝑗

𝑖𝑛𝑖
𝑗=1 𝜑(‖𝑥 − 𝑥𝑗

𝑖‖),                                               (29)  

 

where 𝑛𝑖 is the number of nodes in Ωi. Therefore, the global RBF-PUM interpolant is defined as 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2176 
https://internationalpubls.com 

 

              𝐼𝑢(𝑥) =  ∑ 𝜔𝑖(𝑥)𝐼𝑢
𝑖 (𝑥)

𝑛𝑝

𝑖=1
= ∑ 𝜔𝑖(𝑥)

𝑛𝑝

𝑖=1
∑ 𝜆𝑗

𝑖𝑛𝑖
𝑗=1 𝜑(‖𝑥 − 𝑥𝑗

𝑖‖),     𝑥 ∈ Ω,     (30) 

 

By enforcing the interpolation condition, we have the following form of the global interpolant [32] 

 

                                               𝑈 = ∑ 𝑅𝑖𝑊𝑖
𝑛𝑝

𝑖=1
𝐴𝑖𝜆𝑖,̅̅ ̅                                                         (31) 

 

where 𝑈 is a vector containing the global solution at𝑥𝑘, 𝑘 =  1, … , 𝑁, 𝑅𝑖 is a permutation projection 

operator which maps the local index set into the global one, 𝑊𝑖 is a diagonal matrix with element 

𝜔𝑖(𝑥) on it and 𝐴𝑖 is the local RBF matrix. The weight function 𝜔𝑖 are constructed using Shepard’s 

method [33] 

 

                                       𝜔𝑖(𝑥) =
𝜓𝑖(𝑥)

∑ 𝜓𝑙(𝑥)
𝑛𝑝
𝑙=1

 ,              𝑖 = 1, … , 𝑛𝑝                                               (32) 

 

where 𝜓(𝑥)  is compactly supported function with support on 𝛺𝑖 . We have used the following 

compactly supported Wendland function [34] 

 

                                      𝜓(𝑟) = {
(1 − 𝑟)2(32𝑟3 + 25𝑟2 + 8𝑟 + 1)           0 ≤ 𝑟 ≤ 1,
0                                                                       𝑟 > 1       

                     

 

In this paper, we will use circular patches due to their flexibility in practice. Thus, the Wendland 

functions will be scaled to get 

 

𝜓𝑖(𝑥) = 𝜓 (
‖𝑥−𝜉𝑖‖

𝜌𝑖
) ,              𝑖 = 1, … , 𝑛𝑝                            (33) 

 

where 𝜉𝑖 and 𝜌𝑖 are, respectively, the centres and the radii of patches 𝛺𝑖, 𝑖 = 1, … , 𝑛𝑝. Let 𝑢𝑘
𝑖  be the 

value of local solution at the node 𝑥𝑘 located in 𝛺𝑖. Because our final target is to solve PDEs, the 

problem here is that there would be more unknown than equations. This can be fixed by requiring the 

local solutions 𝑢𝑘
𝑖  to coincide with the global one. The interpolation property implies that 

 

𝑢𝑖 = 𝐴𝑖𝜆�̅�      ⇒       𝜆�̅� = 𝐴𝑖
−1𝑢𝑖                                     (34) 

 

Therefore, the approximation of any linear differential operator 𝐿 can be derived 

 

𝐿𝑈 = ∑ 𝑅𝑖𝐿(𝑊𝑖
𝑛𝑝

𝑖=1
𝐴𝑖)𝜆𝑖̅̅ ̅̅                                            (35) 

  

              =  ∑ 𝑅𝑖𝐿(𝑊𝑖
𝑛𝑝

𝑖=1
𝐴𝑖)𝐴𝑖

−1𝑢𝑖                                   (36) 

 

                                                        =  ∑ 𝑅𝑖𝐷𝐿
𝑖𝑛𝑝

𝑖=1
𝑢𝑖                                                    (37) 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2177 
https://internationalpubls.com 

 

where 𝐷𝐿
𝑖  is a local differentiation matrix which defined as 𝐷𝐿

𝑖 = 𝐿(𝑊𝑖𝐴𝑖)𝐴𝑖
−1. 

 

3.2.2  Spatial discretization of SWEs using RBF-PUM  

 

We intend to illustrate the discretization of RBF–PU method for the shallow water equations. Let 

{𝑥1, 𝑥2, . . . , 𝑥𝑁} be a set of 𝑁 distinct points. Using (35), the first and second spatial derivatives of a 

function 𝑸(𝑥, 𝑡)   can be expressed as follows: 

 

                                
𝜕𝑸(𝑥,𝑡)

𝜕𝑥
= ∑ 𝑅𝑖𝐷𝑥

𝑖𝑛𝑝

𝑖=1
𝑸𝑖 ,   

𝜕2𝑸(𝑥,𝑡)

𝜕𝑥2 = ∑ 𝑅𝑖𝐷𝑥𝑥
𝑖𝑛𝑝

𝑖=1
𝑸𝑖 ,                        (38) 

 

Hence, the global differentiation matrices 𝐷𝑥 
̅̅ ̅̅  and 𝐷𝑥𝑥 

̅̅ ̅̅ ̅ is computed by assembling the above local 

matrices into the global one, (38) becomes 

 

                                   
𝜕𝑸(𝑥,𝑡)

𝜕𝑥
= 𝐷𝑥 

̅̅ ̅̅ 𝑸(𝑥, 𝑡),          
𝜕2𝑸(𝑥,𝑡)

𝜕𝑥2 = 𝐷𝑥𝑥 
̅̅ ̅̅ ̅𝑸(𝑥, 𝑡),                           (39) 

 

The system (8) is approximated by the following equivalent system 

 

                                     𝜕𝑡𝑸 =  −�̅�(𝑸) + 𝑆̅(𝑸) + �̅�(𝑸).                                        (40) 

 

Using the RBF-PUM to perform numerical evaluation, we obtain 

 

          𝐹 ̅(𝑸) = (
𝐷𝑥 ̅̅ ̅̅ 𝑞

𝐷𝑥 ̅̅ ̅̅ (
𝑞2

ℎ
+

1

2
𝑔ℎ2)

),        𝑆 ̅(𝑸) = ( 0
−𝑔ℎ𝐷𝑥 ̅̅ ̅̅ 𝑏

),         𝐻 ̅̅ ̅(𝑸) = (𝜇𝐷𝑥𝑥 ̅̅ ̅̅ ̅̅ ℎ

𝜇𝐷𝑥𝑥 ̅̅ ̅̅ ̅̅ 𝑞
),             (41) 

 

Equivalent to 

 

  𝜕𝑡𝑸 = 𝑅𝐻𝑆𝑅𝐵𝐹−𝑃𝑈𝑀(𝑸),  where  𝑅𝐻𝑆𝑅𝐵𝐹−𝑃𝑈𝑀(𝑸) = −�̅�(𝑸) + 𝑆̅(𝑸) + �̅�(𝑸).     (42) 

 

3.4  RBF-PUM with QR factorization (RBF-PUM-QR) 

 

It is well known that the use of infinitely smooth RBFs (Table 1) can provide spectral accuracy for 

solving PDEs. As mentioned previously, this type of RBFs is usually formulated by including a free 

shape parameter 𝜀, which is generally applied to control the fatness of the functions. Decreasing 𝜀 

often improves the accuracy of approximation. However, mathematical investigations show that when 

𝜀 →  0 (flat limit region), the RBF methods suffers from severe ill-conditioning. Therefore, RBF-

PUM requires a stable evaluation method to converge as 𝜀 →  0 [38, 39]. The obvious question 

becomes how to devise algorithms that balance numerical stability and computational efficiency when 

dealing with small values of ε. This was the subject of many papers published during the last two 

decades (e.g. [40, 41, 42]). 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2178 
https://internationalpubls.com 

In order to use small values of the shape parameter in RBF-PUM, we consider the RBF-QR algorithm 

[40]. The essential concept behind the RBF-QR algorithm is finding a way to form a well conditioned 

basis in the same function space spanned by a finite set of nearly flat RBFs [41]. The change of basis 

will allow to remove the ill-conditioning issue while maintaining the same accuracy or even better one. 

The construction of the new basis was done by expanding the Gaussian RBF (Table (1)) on Chebyshev 

polynomials in one dimension, on Chebyshev polynomials and trigonometric functions in two 

dimensions and on a combination between Chebyshev polynomials and spherical harmonics in three 

dimensions. The new basis derived from RBF-QR algorithm spans the same approximation space as 

the Gaussian basis which produce an accurate numerical solutions and bypasses some limitations, 

including complications related to the selection of ε, ill-conditioning of RBF-PUM matrices. To have 

more details on the construction of the new basis, we refer the reader to [40]. 

 

3.5 Temporal discretization of the numerical methods  

 

To achieve a higher order of accuracy, it is sensible to use a high order time discretization as well. The 

ODE systems (20), (21), (42) are integrated in time using the fourth order Runge-Kutta scheme. The 

procedure to advance the solution from the time 𝑡𝑛 to the next time  𝑡𝑛+1 is carried out as 

 

𝑘1 =  ∆𝑡 𝑅𝐻𝑆(𝑸𝑛), 

 

              𝑘2 =  ∆𝑡 𝑅𝐻𝑆 (𝑸𝑛 +
1

2
𝑘1), 

 

              𝑘3 =  ∆𝑡 𝑅𝐻𝑆 (𝑸𝑛 +
1

2
𝑘2), 

 

         𝑘4 =  ∆𝑡 𝑅𝐻𝑆(𝑸𝑛 + 𝑘3), 

 

      𝑸𝑛+1 = 𝑸𝑛 +
1

6
(𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4),                       (43) 

 

where 𝑅𝐻𝑆 can be 𝑅𝐻𝑆𝑘𝑎𝑛𝑠𝑎 or 𝑅𝐻𝑆𝑅𝐵𝐹−𝐹𝐷 or𝑅𝐻𝑆𝑅𝐵𝐹−𝑃𝑈𝑀, 𝑛 represents the time level and 𝛥𝑡 is the 

time step, this needs to be chosen carefully to guarantee the stability of the scheme. In all our 

simulations, we chose 𝛥𝑡 using the following formula: 

 

∆𝑡 = 𝐶𝐹𝐿
𝑑𝑚𝑖𝑛

max (|𝑢|+√𝑔ℎ)
 ,                                       (44) 

 

where 𝑑𝑚𝑖𝑛 denotes the smallest nodal distance between collocation points and 𝐶𝐹𝐿 is the Courant 

number such that 0 <  𝐶𝐹𝐿 <  1. In the following, Algorithm 1 and Algorithm 2 illustrates the full 

discretization of the SWEs using RBF-FD method and RBF-PUM-QR method, respectively. 

 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2179 
https://internationalpubls.com 

Algorithm 1 Full discretization of 1D SWEs using RBF-FD method 

1: Enter the required simulation parameters, i.e 𝑵, 𝝁, 𝒏𝒍𝒐𝒄, 𝑻, 𝒏𝒕 . .. 

2: Construct the RBF-FD differentiation matrices : 

a: Specify the computational domain 𝜴. 

b: Generate a set of collocation points (nodes) {𝒙𝒊}𝒊=𝟏
𝑵  in the domain 𝜴. 

c: Choose an RBF 𝝋 and determine the shape parameter 𝜺. 

d: For each node 𝒙𝒊, find a stencil of neighboring nodes {𝒙𝟏
𝒊 , 𝒙𝟐

𝒊 , … , 𝒙𝒏𝒍𝒐𝒄
𝒊 }. 

e: Form a local RBF interpolation matrix using the stencil nodes. 

f: Construct the differentiation matrices 𝑾𝒙 and 𝑾𝒙𝒙 by differentiating the RBFs and solving for          

the weights the system (22). 

3: Make the right hand side of ODE obtained : 

 

𝑹𝑯𝑺𝒉 = @(𝒉, 𝒒) − 𝑾𝒙𝒒 + 𝝁𝑾𝒙𝒙𝒉                                                 (45.a) 

𝑹𝑯𝑺𝒒 = @(𝒉, 𝒒) − 𝑾𝒙 (
𝒒𝟐

𝒉
+

𝟏

𝟐
𝒈𝒉𝟐) −𝒈(𝒉. 𝑾𝒙𝒃) + 𝝁𝑾𝒙𝒙𝒒          (45.b) 

 

4: Enter initial condition. 

5: for 𝒋 = 𝟏: 𝒏𝒕 do 

6:       𝒕 = 𝒋. 𝒅𝒕  

7:       𝒌𝟏 = 𝒅𝒕. 𝑹𝑯𝑺𝒉(𝒉, 𝒒)                                 𝒌𝟏 = 𝒅𝒕. 𝑹𝑯𝑺𝒒(𝒉, 𝒒) 

8:       𝒌𝟐 = 𝒅𝒕. 𝑹𝑯𝑺𝒉(𝒉 +
𝒌𝟏

𝟐
, 𝒒)                         𝒌𝟐 = 𝒅𝒕. 𝑹𝑯𝑺𝒒(𝒉, 𝒒 +

𝒌𝟏

𝟐
) 

9:       𝒌𝟑 = 𝒅𝒕. 𝑹𝑯𝑺𝒉(𝒉 +
𝒌𝟐

𝟐
, 𝒒)                         𝒌𝟑 = 𝒅𝒕. 𝑹𝑯𝑺𝒒(𝒉, 𝒒 +

𝒌𝟐

𝟐
) 

10:     𝒌𝟒 = 𝒅𝒕. 𝑹𝑯𝑺𝒉(𝒉 + 𝒌𝟑, 𝒒)                        𝒌𝟒 = 𝒅𝒕. 𝑹𝑯𝑺𝒒(𝒉, 𝒒 + 𝒌𝟑) 

11:     𝒉 = 𝒉 +
𝟏

𝟔
(𝒌𝟏 + 𝟐𝒌𝟐 + 𝟐𝒌𝟑 + 𝒌𝟒)           𝒒 = 𝒒 +

𝟏

𝟔
(𝒌𝟏 + 𝟐𝒌𝟐 + 𝟐𝒌𝟑 + 𝒌𝟒) 

12:      Apply boundary condition. 

13: end 

14: Compute error. 

 

 

Algorithm 1 Full discretization of 1D SWEs using RBF-PUM-QR method 

1: Enter the required simulation parameters, i.e 𝑵, 𝝁, 𝒏𝒑, 𝑻, 𝒏𝒕 . .. 

2: Construct the RBF-PUM-QR differentiation matrices : 

a: Specify the computational domain 𝜴. 

b: Generate a set of collocation points (nodes) {𝒙𝒊}𝒊=𝟏
𝑵  in the domain 𝜴. 

    c: Divide the domain into overlapping subdomains {Ω𝒊}𝒊=𝟏

𝒏𝒑
 such that 𝜴 ⊆ ⋃ 𝜴𝒊

𝒏𝒑

𝒊=𝟏
. 

    d: For each subdomain Ω𝒊, select a set of nodes {𝒙𝒋
𝒊}𝒋=𝟏

𝒏𝒊  within  Ω𝒊 . 

    e: For each point {𝒙𝒋
𝒊}𝒋=𝟏

𝒏𝒊 , define the partition of unity using Shepard’s method (32). 

    f: Construct the differentiation matrices 𝑫𝒙 
̅̅ ̅̅  and 𝑫𝒙𝒙 

̅̅ ̅̅ ̅̅  using the differentiation weights obtained 

from the QR factorization of the local RBF-PUM systems, specifically using Gaussian RBFs [40]. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2180 
https://internationalpubls.com 

3: Make the right hand side of ODE obtained : 

 

𝑹𝑯𝑺𝒉 = @(𝒉, 𝒒) − 𝑫𝒙 
̅̅ ̅̅ 𝒒 + 𝝁𝑫𝒙𝒙 

̅̅ ̅̅ ̅̅ 𝒉                                                 (46.a) 

𝑹𝑯𝑺𝒒 = @(𝒉, 𝒒) − 𝑫𝒙 
̅̅ ̅̅ (

𝒒𝟐

𝒉
+

𝟏

𝟐
𝒈𝒉𝟐) − 𝒈(𝒉. 𝑫𝒙 

̅̅ ̅̅ 𝒃) + 𝝁𝑫𝒙𝒙 
̅̅ ̅̅ ̅̅ 𝒒        (46.b) 

 

4: Enter initial condition. 

5: for 𝒋 = 𝟏: 𝒏𝒕 do 

6:       𝒕 = 𝒋. 𝒅𝒕  

7:       𝒌𝟏 = 𝒅𝒕. 𝑹𝑯𝑺𝒉(𝒉, 𝒒)                                 𝒌𝟏 = 𝒅𝒕. 𝑹𝑯𝑺𝒒(𝒉, 𝒒) 

8:       𝒌𝟐 = 𝒅𝒕. 𝑹𝑯𝑺𝒉(𝒉 +
𝒌𝟏

𝟐
, 𝒒)                         𝒌𝟐 = 𝒅𝒕. 𝑹𝑯𝑺𝒒(𝒉, 𝒒 +

𝒌𝟏

𝟐
) 

9:       𝒌𝟑 = 𝒅𝒕. 𝑹𝑯𝑺𝒉(𝒉 +
𝒌𝟐

𝟐
, 𝒒)                         𝒌𝟑 = 𝒅𝒕. 𝑹𝑯𝑺𝒒(𝒉, 𝒒 +

𝒌𝟐

𝟐
) 

10:     𝒌𝟒 = 𝒅𝒕. 𝑹𝑯𝑺𝒉(𝒉 + 𝒌𝟑, 𝒒)                        𝒌𝟒 = 𝒅𝒕. 𝑹𝑯𝑺𝒒(𝒉, 𝒒 + 𝒌𝟑) 

11:     𝒉 = 𝒉 +
𝟏

𝟔
(𝒌𝟏 + 𝟐𝒌𝟐 + 𝟐𝒌𝟑 + 𝒌𝟒)           𝒒 = 𝒒 +

𝟏

𝟔
(𝒌𝟏 + 𝟐𝒌𝟐 + 𝟐𝒌𝟑 + 𝒌𝟒) 

12:      Apply boundary condition. 

13: end 

14: Compute error. 

 

 

3.6 Friction source term discretization 

 

The non-linear nature of the friction terms (2), (3) and their interactions with other source terms 

remains to be a challenge for developing numerically accurate and stable schemes to solve the SWEs 

for simulating very shallow flows as found in the applications involving overland flows and wet/dry 

fronts. Indeed, the friction term actually dominates the stability of a numerical scheme. In particular, 

when the water depth becomes very small, the friction formulations may lead to an exaggerated force 

that can even reverse the flow, which is obviously physically incorrect [43]. To deal with this problem, 

at the beginning of each step of the Runge-Kutta scheme, the friction effect is evaluated and used 

implicitly by the splitting method described by [27, 43], and it is equivalent to solve the following 

ordinary differential equations 

 
𝜕𝑞

𝜕𝑡
=  𝑆𝑓 .                                                        (47) 

 

Since the friction term is only involved in the momentum equation, only the flow rates need to be 

evaluated. The above equation is then discretized by a full implicit method as: 

 

𝑞𝑛+1−𝑞𝑛

∆𝑡
= 𝑆𝑓

𝑛+1,                                              (48) 

 

 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2181 
https://internationalpubls.com 

where the friction term 𝑆𝑓
𝑛+1 may be expressed using a Taylor series as 

 

𝑆𝑓
𝑛+1 = 𝑆𝑓

𝑛 + (
𝜕𝑆𝑓

𝜕𝑞
)

𝑛

∆𝑞 + 𝑜(∆𝑞2),                            (49) 

 

Where ∆𝑞 = 𝑞𝑛+1 − 𝑞𝑛. Ignoring the higher-order terms and substituting the above equation into (48) 

lead to the following formula for updating water discharge 𝑞 at the new time step: 

 

𝑞𝑛+1 = 𝑞𝑛 + ∆𝑡
𝑆𝑓

𝑛

1−∆𝑡(
𝜕𝑆𝑓

𝜕𝑞
)

𝑛                                         (50) 

= 𝑞𝑛 + ∆𝑡 (
𝑆𝑓

𝐷
)

𝑛

= 𝑞𝑛 + ∆𝑡𝐹,                              (51) 

where 

𝐷 = {
1 + 2∆𝑡 (

𝑔𝑛2|𝑢|

ℎ
4
3

)                           

1 + 2∆𝑡 (
𝑓|𝑢|

8ℎ
)                                

                         (52) 

 

𝐷 is the coefficient derived for a full implicit scheme where the first line is for manning law and the 

second line is for Darcy-Weisbach law and 𝐹  is the friction source term including the implicit 

coefficient. This updated water discharge is used as an initial condition for the operators in equation 

(43). 

 

4. Numerical results 

In order to verify the feasibility and the ability of the proposed methods combined with artificial 

viscosity to solve the problem with strong discontinuity appears in SWEs, we present sets of numerical 

tests, including various frictionless steady-state solutions, transient solutions and steady-state solutions 

with friction to address different numerical challenges. These problems are widely used to test 

numerical algorithms for the SWEs (e.g. [44, 45, 51]). The goal of the first numerical test is to assess 

the ability of the suggested schemes to capture steady state solutions, in addition to preserving them. 

The second numerical example is performed to test if the numerical methods catches the transitory 

behaviour and shock behaviour of the solution properly. The third numerical test verifies the accuracy 

of the proposed methods in preserving moving steady states that involve both topography and friction. 

The following parameters are defined: 𝑁 is the total number of nodes in the whole computational 

domain, 𝑛𝑙𝑜𝑐  is the stencil number of the RBF-FD method based on the inverse multiquadric radial 

basis function, the optimal shape parameter 𝜀 is selected using the computed error of the approximate 

solutions, 𝑛𝑝 is the number of patch of the RBF-PUM-QR method. For Kansa method based on the 

multiquadric radial basis function, the shape parameter 𝜀 =
𝜀0√𝑁

𝑑𝑚𝑖𝑛
  in which 𝑑𝑚𝑖𝑛 is the minimum 

distance between two nodes and 10−4 ≤ 𝜀0 ≤ 10−1 [46] and 𝑇 is the terminal time of the simulation. 

To show the advantages of the methods in terms of accuracy, we compute the relative error as follows: 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2182 
https://internationalpubls.com 

𝐸𝑟𝑟 =  √1

𝑁
∑ (

𝑸𝑛𝑢𝑚(𝑥𝑖,𝑡)−𝑸𝑎𝑛𝑎(𝑥𝑖,𝑡)

𝑸𝑎𝑛𝑎(𝑥𝑖,𝑡)
)

2
𝑁
𝑖=1 ,                                    (53) 

 

where 𝑸𝑛𝑢𝑚 and 𝑸𝑎𝑛𝑎 are respectively, the numerical and analytical solutions. In all our computations 

a fixed courant number 𝐶𝐹𝐿 =  0.5. The numerical tests are performed on a core CPU i5 2.11GH 

computer in a MATLAB 2018a tools. 

 

4.1 Test 1: Steady flow over a bump 

 

According to [47], the source term appears in (8) is a crucial point in preserving steady states. In order 

to prove the ability of the proposed methods to catch these states, we apply the schemes to a series of 

benchmark cases of different boundary and initial conditions [48]. These benchmarks have been 

derived by considering the Bernoulli equation that governs the steady state solutions of the shallow-

water equations with non-flat topography and a vanishing friction contribution, The experiments from 

[48] are called the subcritical flow, the transcritical flow without shock and the transcritical flow with 

shock. The computational domain is a rectangular channel with 𝐿 =  25 𝑚 in length. The inflow and 

outflow boundary conditions are imposed along the left and the right boundary segments, respectively. 

The bed elevation of a bump is defined as follows: 

 

𝑏(𝑥) = {0.2 − 0.05(𝑥 − 10)2                     8𝑚 < 𝑥 < 12𝑚
0                                                                𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒.

                          (54) 

 

The boundary conditions are defined in terms of two quantities, 𝑞0  and ℎ𝐿 , whose values vary 

depending on the specific experiment being analysed: 

 

• on the left boundary, the water height satisfies a homogeneous Neumann condition and the 

discharge is set to a specific 𝑞0. 

• on the right boundary, the water height is set to ℎ𝐿  when the flow is subcritical (and a 

homogeneous Neumann boundary condition is prescribed otherwise), and the discharge follows a 

homogeneous Neumann boundary condition. 

 

In addition, the initial conditions set to ℎ(𝑥) + 𝑏(𝑥) = ℎ𝐿𝑚  and 𝑞(𝑥) = 0𝑚2/𝑠  throughout the 

domain.  

All these results are frictionless and displayed at 𝑇 = 200𝑠 using𝑁 = 1001, 𝑛𝑙𝑜𝑐 = 100 and 𝑛𝑝 =

300, the analytical solutions are also plotted within the obtained numerical results. 

4.1.1 Subcritical flow 

For this test, a discharge of 𝑞0 = 4.42𝑚2/𝑠 and the constant water surface level of the water level 

ℎ𝐿 = 2𝑚  are set as upstream and downstream boundary conditions, respectively. The steady-state 

numerical solutions of water height ℎ and discharge 𝑞 compared with the analytical solution [48] are 

plotted in Figure 1. The numerical examples are performed for Kansa method using 𝜇 = 8 × 10−3, for 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2183 
https://internationalpubls.com 

RBF-PUM-QR method using 𝜇 = 8 × 10−3, and for RBF-FD method using 𝜇 = 2 × 10−3 with the 

optimal shape parameter 𝜀 = 5 . The selection of an optimal shape parameter 𝜀  and the viscosity 

coefficient 𝜇 in all tests in this paper was made using trial and error strategy. The strategy is to perform 

a series of experiments by varying the values of the shape parameter or viscosity coefficient, and then 

pick the best one corresponding to the smallest error. In Figure 1, it is evident that the impact of a 

bump in subcritical flow is clearly observable and the model can well capture the drop of the water 

surface elevation at the bump. A more detailed analysis of the numerical error is listed in Table 2. 

Figure 2  shows the computed 𝐿2 error for different values of the viscosity coefficient 𝜇. From Figure 

2, we remark that for Kansa method the minimum error is for obtained around 𝜇 = 10−3, for RBF-FD 

method 𝜇 = 10−2  and for RBF-PUM-QR method 𝜇 = 10−3 . We make a error analysis for the 

subcritical flow problem using RBF-FD method (see Figure 3 (a)). Four different stencil sizes are used, 

𝑛𝑙𝑜𝑐 = 20, 𝑛𝑙𝑜𝑐 = 100 and 𝑛𝑙𝑜𝑐 = 150. The result seen in Figure 3 (a) plotted in logarithmic scales , 

where the 𝐿2-norm of the numerical error is plotted against the grid size, we display in Figure 3 (b) the 

error of the numerical solution at 𝑇 =  200 seconds in terms of the shape parameter in order to show 

that stable computations were achieved for 𝜀 = 5. The computed solutions are oscillation-free which 

demonstrates the ability of the proposed methods to accurately capture steady state solutions. 

 

Table 2. Error comparison of water level and discharge for different flow conditions (subcritical, 

transcritical without shock and transcritical with shock) using Kansa, RBF-FD and RBF-PUM-QR 

methods.  

 

                                                                                    Kansa                    RBF-FD                 RBF-PUM-QR 

                                              Water height ℎ[𝑚]     6.73 × 10−6            3.94 × 10−6              1.38 × 10−6 

         Subcritical 

                                              Discharge 𝑞[𝑚2/𝑠]    7. 02 × 10−6            4.21 × 10−6              1.05 × 10−6 

                                              Water height ℎ[𝑚]     3.99 × 10−5            1.95 × 10−5              1.26 × 10−5 

Transcritical without shock 

                                              Discharge 𝑞[𝑚2/𝑠]    1. 88 × 10−5            1.06 × 10−5              6.90 × 10−6 

                                              Water height ℎ[𝑚]     4.70 × 10−4            4.51 × 10−4              4.40 × 10−4 

  Transcritical with shock 

                                              Discharge 𝑞[𝑚2/𝑠]    3.50× 10−4            2.50 × 10−4              3.10 × 10−4 

 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2184 
https://internationalpubls.com 

                             (a) Kansa method                                                      (b) RBF-FD method 

 

 

 

 

 

 

 

 

 

 

(c) RBF-PUM-QR method 

Figure 1. Subcritical flow over a bump using (a) Kansa, (b) RBF-FD, and (c) RBF-PUM-QR at T = 

200 seconds 

4.1.2 Transcritical flow without shock 

In this case, the upstream inflow of the discharge 𝑞0 = 1.53𝑚2/𝑠  is imposed and a downstream 

condition for the height ℎ𝐿 = 0.40𝑚 only if the flow is subcritical. If the flow is supercritical, no 

condition is imposed. The numerical examples are performed for Kansa method using 𝜇 = 7 × 10−3, 

for RBF-FD method using 𝜇 = 4 × 10−3, the optimal shape parameter 𝜀 = 7 and for RBF-PUM-QR 

method using 𝜇 = 10−3. As shown in Figure 4, the water height and discharge are compared with the 

analytical solution [48], which can be observed that the water surface drops significantly as the flow 

passes through the bump. The transcritical flow happened nearby the bump, since the subcritical flow 

is at the upstream and the supercritical flow is at the downstream. The errors of water height ℎ and 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2185 
https://internationalpubls.com 

discharge 𝑞 are listed in Table 2. It can be clearly seen that the numerical and analytical solutions are 

in good agreement. 

 

 

 

 

 

 

 

 

 

 

 

 

                          (a) Kansa method                                                            (b) RBF-FD method 

 

 

 

 

 

 

 

 

 

 

 

 

                                                               (c) RBF-PUM-QR method 

Figure 2. 𝐿2 error as a function of the viscosity coefficient for the ”Subcritical flow” problem for (a) 

Kansa, (b) RBF-FD, and (c) RBF-PUM-QR at T = 200 seconds 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2186 
https://internationalpubls.com 

 

 

 

 

 

 

 

 

 

 

(a) (b)  

Figure 3. 

(a) 𝐿2 Errors of the RBF-FD method for solving ”Subcritical flow” problem using 𝜇 = 8 ×

10−3. (b) 𝐿2 error in terms of the shape parameter 𝜀. 

 

4.1.3 Transcritical flow with shock 

This case is similar to the previous one but with different boundary conditions. Here, a discharge of 

𝑞0 = 0.18𝑚2/𝑠 is imposed as the upstream boundary condition and a water level of ℎ𝐿 = 0.33𝑚 as 

the downstream boundary condition; The flow regime changed from subcritical flow to supercritical 

flow and back to subcritical flow with a hydraulic jump over the bump. For this numerical test, we 

used the following parameters 𝜇 = 8 × 10−3,  𝜇 = 6 × 10−3 and 𝜇 = 9 × 10−3, for Kansa, RBF-FD, 

RBF-PUM-QR methods respectively. In our simulation, we set the shape parameter 𝜀 = 5 . Figure 5, 

shows the numerical solutions concerning the water height and discharge, compared with the analytical 

solution. However, near the jump discontinuity we can observe differences that are even more clear 

when we compare analytical and numerical discharge. A more quantitative analysis of the accuracy is 

listed in Table 1, we observe that capturing the discharge 𝑞 correctly in the Kansa method is more 

difficult than in the other methods. From the results and comparisons in this example, the merits of the 

proposed meshless methods are verified. 

 

4.2 Test 2: Dam break problem 

The dam-break problem is the most commonly investigated problem and have also become the 

standard test scenario for numerical methods for shallow water equations. Toro [49] has the dam-break  

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2187 
https://internationalpubls.com 

 

                             (a) Kansa method                                                      (b) RBF-FD method 

 

 

 

 

 

 

 

 

 

 

(c) RBF-PUM-QR method 

Figure 4. Transcritical flow without shock over a bump using (a) Kansa, (b) RBF-FD, and (c) RBF-

PUM-QR at T = 200 seconds 

in high regard, because its solution includes both continuous and discontinuous solutions at the same 

time, and several researches have shown that the shallow water equations are suitable for the 

representation of dam-break flows e.g. [3, 50, 51] and thus, in order to examine the capability of the 

proposed methods considering artificial viscosity term to solve the problem of discontinuous initial 

condition, we consider dam-break test is adopted as the second test. The one-dimensional dam break 

is an initial value problem consisting of two still bodies of water of different heights. The upstream 

body (the reservoir) is separated from the downstream body of water by a partition that is 

instantaneously removed at time 𝑡 =  0 seconds. The two bodies of water are allowed to interact under 

the force of gravity. We consider the dam-break problem in a rectangular channel with flat bottom, 

𝑏(𝑥)  =  0 and frictionless. The channel is of length 10𝑚 and the initial conditions are given by 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2188 
https://internationalpubls.com 

 

                             (a) Kansa method                                                      (b) RBF-FD method 

 

 

 

 

 

 

 

 

 

 

(c) RBF-PUM-QR method 

Figure 5. Transcritical flow with shock over a bump using (a) Kansa, (b) RBF-FD, and (c) RBF-

PUM-QR at T = 200 seconds 

 

              𝑞(𝑥, 0) = 0         and          ℎ(𝑥, 0) = {
0.005                 𝑥 < 5𝑚                   
0.001                𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒              

 (55) 

 

As boundary conditions, a zero discharge and a free boundary are considered at the left and right ends 

of the channel. The analytical solution for this simple dam break test consists of a backward 

propagating rarefaction and a forward-moving shock wave [49]. The comparison of the numerical 

results at 𝑇 =  2, 4, 6 𝑠 are displayed in Figure 6 for RBF-FD method using 𝜇 = 2 × 10−4 and the 

optimal shape parameter 𝜀 =  12. Table 3 depicts the error comparison of water height ℎ and discharg   



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2189 
https://internationalpubls.com 

Table 2. Error comparison of water level and discharge for Dam break problem over a wet bed. 

 

                                                               Kansa                    RBF-FD                 RBF-PUM-QR 

    Water height ℎ[𝑚]                        3.07 × 10−4            1.71 × 10−4              2.10 × 10−4 

        

    Discharge 𝑞[𝑚2/𝑠]                       1.48× 10−3            9.76 × 10−4              8.85 × 10−4 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6: Dam-break on wet bed using RBF-FD method at 𝑇 =  2, 4, 6 𝑠, top: water height, bottom: 

discharge. 

for all RBF-based methods. The results show that the proposed meshless numerical methods with 

hyperviscosity has ability of shock capturing and provides stability in dealing with problems with 

discontinuous flow. Additionally, Figure 7 shows the computed 𝐿2 error for different values of the 

viscosity coefficient 𝜇, this coefficient controls the accuracy of the method where small values are 

expected to provide accurate solutions, and should be a limit value of around 10−4 to avoid instability 

issues. Another consequence is that the accuracy of RBF-FD method does not indefinitely increase 

with stencil size 𝑛𝑙𝑜𝑐. As shown in Figure 8 (left) the stencil size has no bearing on accuracy, hence 

𝑛𝑙𝑜𝑐  =  100 is the stencil size chosen for this test case. With 𝜇 =  10−4, we display in Figure 8 (right) 

the curve of the 𝐿2 in terms of the shape parameter ε, we remark that the minimum error is obtained 

for 𝜀 =  12. From this study, we can note that increasing the number of points per stencil leads to 

enhanced accuracy. 

 

 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2190 
https://internationalpubls.com 

4.3 Test 3: Steady flow with friction 

 

 

 

 

 

 

 

 

 

 

 

 

                          (a) Kansa method                                                            (b) RBF-FD method 

 

 

 

 

 

 

 

 

 

 

 

                                                               (c) RBF-PUM-QR method 

Figure 7. 𝐿2 error as a function of the viscosity coefficient for the ”Dam-break” problem for (a) Kansa, 

(b) RBF-FD, and (c) RBF-PUM-QR at T = 6 seconds 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2191 
https://internationalpubls.com 

 

 

 

 

 

 

 

 

 

 

 

(a)                                                                            (b) 

Figure 8. (a) 𝐿2 Errors of the RBF-FD method for solving ”Dam-break” problem with respect to the 

number of stencils using 𝜇 = 2 × 10−4. (b) 𝐿2 error in terms of the shape parameter 𝜀. 

 

The solutions presented in this section are more intricate than those in the previous section 4.1-4.2 due 

to the variability of the topography near the boundaries. Consequently, they provide a more rigorous 

validation of the boundary conditions. When 𝑆𝑓 =  0  (indicating bottom friction), the following 

solutions can verify if the friction terms are implemented correctly to maintain steady states. These 

solutions are derived using the procedure introduced by I. MacDonald [52, 53]: given the water depth 

profile and the discharge, the corresponding topography is subsequently computed. The three selected 

solutions illustrate the efficacy of the scheme in addressing stationary states induced by the topography 

and by friction under a wide range of flow conditions. 

At steady states, we have 𝜕𝑡ℎ =  𝜕𝑡𝑢 =  𝜕𝑡𝑞 =  0, thus the mass-conservation equation gives 𝑞 = 

constant and we get the equation 

 

𝜕𝑥𝑏 = (
𝑞2

𝑔ℎ3 − 1) 𝜕𝑥ℎ + 𝑆𝑓.                                           (56) 

 

Where 𝑆𝑓 depends on the friction law chosen. From this relation, one can make as many solutions as 

needed. In this section, we present a few of these solutions obtained for specific value of the domain 

length 𝐿 and fixed parameters, such as the friction law and its coefficient. The water height profile and 

the discharge are provided, and the corresponding topographies are computed by solving the equation 

(56) with a high order iterative method [51]. 

In the following computations, the initial conditions set to ℎ(𝑥) + 𝑏(𝑥)  =  0, 𝑞(𝑥)  =  0  and all 

results are displayed at 𝑇 =  1500𝑠 using 𝑁 =  801, 𝑛𝑙𝑜𝑐  =  100 and 𝑛𝑝  =  400 . The analytical 

solutions shown in [51] are also plotted within the obtained numerical results. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2192 
https://internationalpubls.com 

4.3.1 Subcritical flow 

We consider a 𝐿 =  1000 𝑚 long channel with a discharge of q = 2𝑚2/𝑠. The flow is subcritical at 

inflow and is subcritical at outflow with depth ℎ𝑎𝑛𝑎(1000)  =  0.748409 𝑚 where ℎ𝑎𝑛𝑎 is analytical 

solution. The Manning coefficient for the channel is 𝑛 =  0.033𝑚−1/3𝑠. The boundary conditions set 

to 

 

{
𝑞(𝑥) = 2                                  𝑥 = 0            

ℎ(𝑥) = ℎ𝑎𝑛𝑎(1000)               𝑥 = 𝐿             
                    (57) 

 

The water depth and discharge are shown in Figure 9 for Kansa, RBF-FD with an optimal shape 

parameter ε = 0.25 and RBF-PUM-QR methods using 𝜇 =  10−1 in all cases. The numerical solutions 

is very close to the analytical solution, and we obtained very satisfying results in the approximation of 

the steady discharge and water height as illustrated in Table 4. Figure 10 shows the computed 𝐿2 error 

for different values of the viscosity coefficient 𝜇. From Figure 10, we remark that for Kansa method 

the minimum error is for obtained around 𝜇 =  10−1, for RBF-FD method 𝜇 =  10−1 and for RBF-

PUM-QR method 𝜇 =  10−1. We make an error analysis for the subcritical flow problem using RBF-

FD method (see Figure 11 (a)). Four different stencil sizes are used, 𝑛𝑙𝑜𝑐  =  20, 𝑛𝑙𝑜𝑐  =  100 and 

𝑛𝑙𝑜𝑐  =  150. The result in seen in Figure 11 (a) plotted in logarithmic scales, where the 𝐿2-norm of 

the numerical error is plotted against the grid size, we display in Figure 11 (b) the error of the numerical 

solution at 𝑇 =  1500  seconds in terms of the shape parameter in order to show that stable 

computations were achieved for 𝜀 =  0.25. From the results and comparisons in this example, the 

merits of the proposed meshless method are verified. 

 

4.3.2 Supercritical flow 

 

Considering a 1000 𝑚  long computational domain with a discharge of q = 2𝑚2/𝑠 . The flow is 

supercritical at inflow with depth ℎ𝑎𝑛𝑎(0) =  0.741599 𝑚  and is supercritical at outflow. The 

Manning coefficient for the channel is 𝑛 =  0.0218𝑚−1/3𝑠. The boundary conditions set to 

 

 

{
𝑞(𝑥) = 2,        ℎ(𝑥) = ℎ𝑎𝑛𝑎(0)                    𝑥 = 0,            
 𝑓𝑟𝑒𝑒                                                                   𝑥 = 𝐿.             

                    (58) 

 

The initial free surface is set to zero on the whole domain. We show on Figure 12 the steady state water 

free surface and discharge at 𝑇 =  1500𝑠 and we observe very close agreement between numerical 

results and the reference solution. The numerical examples are performed for Kansa, RBF-FD with an 

optimal shape parameter 𝜀 =  0.1 and RBF-PUM-QR methods using 𝜇 =  10−1 for all cases. The 

error of converged results by different schemes are listed in Table 4. 

 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2193 
https://internationalpubls.com 

4.3.3 Supercritical to subcritical flow 

 

We consider a discharge of q = 2𝑚2/𝑠. Again, the initial free surface and discharge are set to zero on 

the whole domain. The flow is supercritical at inflow with depth ℎ𝑎𝑛𝑎(0) =  0.543853𝑚 and is 

 

Table 4. Error comparison of water level and discharge for different flow conditions (subcritical, 

supercritical and supercritical to subcritical) with friction using Kansa, RBF-FD and RBF-PUM-QR 

methods.  

 

                                                                                    Kansa                    RBF-FD                 RBF-PUM-QR 

                                              Water height ℎ[𝑚]     2.66 × 10−4            2.29 × 10−4              2.55 × 10−4 

         Subcritical 

                                              Discharge 𝑞[𝑚2/𝑠]    1.85 × 10−4            1.39 × 10−4              1.82 × 10−4 

                                              Water height ℎ[𝑚]     8.88 × 10−5            1.07 × 10−5              8.65 × 10−5 

        Supercritical 

                                              Discharge 𝑞[𝑚2/𝑠]    9.55 × 10−5            7.83 × 10−5              6.04 × 10−6 

                                              Water height ℎ[𝑚]     6.45 × 10−4            8.31 × 10−4              5.20 × 10−4 

Supercritical to subcritical 

                                              Discharge 𝑞[𝑚2/𝑠]    2.83 × 10−4            3.21 × 10−4              2.50 × 10−4 

 

subcritical at outflow with depth ℎ𝑎𝑛𝑎(1000)  =  1.334899m . The Manning coefficient for the 

channel is 𝑛 =  0.0218𝑚−1/3𝑠. The boundary conditions set to 

 

{
𝑞(𝑥) = 2,        ℎ(𝑥) = ℎ𝑎𝑛𝑎(0)                            𝑥 = 0,            

  ℎ(𝑥) = ℎ𝑎𝑛𝑎(1000)                                              𝑥 = 𝐿.             
                    (59) 

 

The water depth and discharge are shown in Figure 13 for Kansa, RBF-FD 𝜀 =  0.1 and RBF-PUM-

QR methods using 𝜇 =  10−1. We can clearly see that the flow accurately converges towards the 

steady state and that the shock location is accurately computed. The comparison of the relative error 

of the results by Kansa, RBF-FD and RBF-PUM-QR methods is shown in Table 4. 

 

5. Conclusion 

In this paper, the shallow water equations is presented to simulate the one-dimensional flow of water 

in channels. The resulting system is solved numerically using a globally, localized and stable meshless 

methods based on radial basis functions. In order to stabilize these methods and minimize the non-

physical numerical oscillations near the discontinuities which can solve the strong discontinuity 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2194 
https://internationalpubls.com 

problem in SWEs, an artificial viscosity (Hyperviscosity) technique is implemented. In Kansa’s 

method the system to be solved is dense and therefore it was necessary to refer to localized meshless 

methods such as RBF-FD and RBF-PUM in order to present sparsity and locality. The problem of 

choosing any value of shape parameter ε from the flat limit region which affect the accuracy of the 

solution can be handled by combining the RBF-PUM with RBF-QR approach. 

The proposed meshless methods with hyperviscosity has been tested on systems of shallow water 

equations at different flow regimes. The obtained results indicate good shock resolution with high  

 

 
                             (a) Kansa method                                                      (b) RBF-FD method 

 

(c) RBF-PUM-QR method 

Figure 9. Subcritical flow using (a) Kansa, (b) RBF-FD, and (c) RBF-PUM-QR methods at T = 

1500 seconds with Manning friction 𝑛 =  0.033𝑚−1/3𝑠. 

 

accuracy in smooth regions, and the convergence to the correct steady state solution has been clearly 

verified in flow over a non-flat bottom, which confirm the numerical stability and computational 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2195 
https://internationalpubls.com 

accuracy. Additionally, the main challenge encountered was the selection of hyperviscosity 

parameters, particularly the constant μ and its scaling relations, to ensure adequate stabilization across 

a wide range of flow regimes and discretization densities. 

While our numerical computations have been limited to one-dimensional shallow water problems, the 

current meshless-based RBF methods can be readily extended to two-dimensional shallow water 

problems incorporating source terms. These extensions, along with other related issues, will be the 

focus of future investigations. 

 

 

                      (a) Kansa method                                                            (b) RBF-FD method 

 

                                                               (c) RBF-PUM-QR method 

Figure 10. 𝐿2 error as a function of the viscosity coefficient for the ”Subcritical flow” problem for (a) 

Kansa, (b) RBF-FD, and (c) RBF-PUM-QR at T = 1500 seconds 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2196 
https://internationalpubls.com 

 

(a)                                                  (b) 

Figure 11. (a) 𝐿2 Errors of the RBF-FD method for solving ”Subcritical flow” problem with respect 

to the number of stencils using 𝜇 = 10−1. (b) 𝐿2 error in terms of the shape parameter 𝜀. 

 

 

             (a) Kansa method                                                    (b) RBF-FD method 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2197 
https://internationalpubls.com 

 

(c) RBF-PUM-QR method 

Figure 12. Supercritical flow using (a) Kansa, (b) RBF-FD, and (c) RBF-PUM-QR methods at T = 

1500 seconds with Manning friction 𝑛 =  0.0218𝑚−1/3𝑠. 

 

                            (a) Kansa method                                                      (b) RBF-FD method 

 
(c) RBF-PUM-QR method 

Figure 13. Supercritical to subcritical flow using (a) Kansa, (b) RBF-FD, and (c) RBF-PUM-QR 

methods at T = 1500 seconds with Manning friction 𝑛 =  0.0218𝑚−1/3𝑠. 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2198 
https://internationalpubls.com 

Refrences 

[1] A.B. de Saint-Venant, Théorie du mouvement non permanent des eaux, avec application aux crues 

des rivières et a l'introduction de marées dans leurs lits,  Comptes rendus de l'Académie des 

Sciences de Paris 73, 147--154 and 237-240 (1871). 

[2] M.H. Chaudhry, Open-channel flow, Boston, MA: Springer US (2008). 

[3] O. Delestre, C. Lucas, P.A. Ksinant, F. Darboux, C. Laguerre, et al., SWASHES: a compilation of 

Shallow Water Analytic Solutions for Hydraulic and Environmental Studies, International 

Journal for Numerical Methods in Fluids 72(3), 269-300 (2013). 

[4] Y. Xing, C.W. Shu, Solution of shallow-water equations using least-squares finite-element method, 

Journal of Computational Physics 208(1), 206-227 (2005). 

[5] S.J. Liang, J.H. Tang, M.S. Wu, High order finite difference WENO schemes with the exact 

conservation property for the shallow water equations, Acta Mechanica Sinica 24, 523-532 (2008). 

[6] S.N. Kuiry, K. Pramanik, D. Sen, Finite Volume Model for Shallow Water Equations with 

Improved Treatment of Source Terms, Journal of Hydraulic Engineering 134(2), 231-242 (2008). 

[7] G.H. Duenas, A. Beljadid, A central-upwind scheme with artificial viscosity for shallow-water 

flows in channels, Advances in Water Resources 96, 323–338 (2016). 

[8] Y.C. Hon, K.F. Cheung, X.Z. Mao, E.J. Kansa, Multiquadric Solution for Shallow Water 

Equations, Journal of Hydraulic Engineering 125, 524-533 (1999). 

[9] S.M. Wong, Y.C. Hon, M.A. Golberg, Compactly supported radial basis functions for shallow 

water equations. Applied Mathematics and Computation 127, 79-101 (2002).          

[10] S.M. Wong, Y.C. Hon, M.A. Golberg, Compactly supported radial basis functions for shallow 

water equations. Applied Mathematics and Computation 127, 79-101 (2002). 

[11] E.J. Kansa, Multiquadrics-{A} scattered data approximation scheme with applications to 

computational fluid-dynamics-{II} Solutions to parabolic, hyperbolic and elliptic partial 

differential equations,  Computers and Mathematics with Applications 19(8), 147-161 (1990). 

[12] A. Tolstykh, D. Shirobokov, On using radial basis functions in a finite difference mode with 

applications to elasticity problems, Computational Mechanics 33, 68-79 (2003). 

[13]  B. Fornberg, N. Flyer, A primer on radial basis functions with applications to the geosciences}, 

CBMS-NSF Regional Conference Series in Applied Mathematics, Society for Industrial and 

Applied Mathematics (SIAM), (2015). 

[14] E. Ben-Ahmed, M. Sadik and M. Wakrim, Radial basis function partition of unity method for 

modelling water flow in porous media, Computers and Mathematics with Applications 8(75), 

2925-2941 (2018). 

[15] V. Shcherbakov and E. Larsson, Radial basis function partition of unity methods for pricing 

vanilla basket options, Computers and Mathematics with Applications 71(1), 185-200 (2016). 

[16] E. Ben-Ahmed, M. Sadik and M. Wakrim, RBFPUM with QR Factorization for Solving Water 

Flow Problem in Multilayered Soil, International Journal of Nonlinear Sciences and Numerical 

Simulation, (2018). 

[17] E. Ben-Ahmed, M. Sadik and M. Wakrim, A Stable Radial Basis Function Partition of Unity 

Method with d-Rectangular Patches for Modelling Water Flow in Porous Media, Journal of 

Scientific Computing 84, (2020). 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2199 
https://internationalpubls.com 

[18] B. Fornberg, E. Lehto, Stabilization of rbf-generated finite difference methods for convective 

pdes, Journal of Computational Physics 230, 2270-2285 (2011). 

[19] D. Stevens, H. Power, The Radial Basis Function Finite Collocation Approach for Capturing 

Sharp Fronts in Time Dependent Advection Problems, Journal of Computational Physics 298, 

423-445 (2015). 

[20] H. Ma, Chebyshev-Legendre super spectral viscosity method for nonlinear conservation laws, 

SIAM Journal of Numerical Analysis 35, 893-908 (1998). 

[21] N. Flyer, E. Lehto, S. Blaise, G.B. Wright, A. St-Cyr, A guide to RBF-generated finite differences 

for nonlinear transport: Shallow water simulations on a sphere, Journal of Computational Physics 

231, 4078-4095 (2012). 

[22] N. Flyer, G.A. Barnett, L.J. Wicker, Enhancing finite differences with radial basis functions: 

Experiments on the Navier–Stokes equations, Journal of Computational Physics 316, 39-62 

(2016). 

[23] V. Shankar, G.B. Wright, A. Narayan, A Robust Hyperviscosity Formulation for Stable RBF-FD 

Discretizations of Advection-Diffusion-Reaction Equations on Manifolds, SIAM Journal on 

Scientific Computing 42, A2371-A2401 (2020). 

[24] I. Tominec, M. Nazarov, Residual viscosity stabilized RBF-FD methods for solving nonlinear 

conservation laws, Journal of Scientific Computing 94, (2022). 

[25] M. Al Nuwairan, E. Chaabelasri, Balanced Meshless Method for Numerical Simulation of 

Pollutant Transport by Shallow Water Flow over Irregular Bed: Application in the Strait of 

Gibraltar, applied sciences 12, 6849 (2022). 

[26] Y. Liu, Y. Qiao, X. Feng, A stable radial basis function partition of unity method for solving 

convection-diffusion equations on surfaces, Engineering Analysis with Boundary Elements 155, 

148-159 (2023). 

[27] M. Dehghan, M. Abbaszadeh, The use of proper orthogonal decomposition (POD) meshless RBF-

FD technique to simulate the shallow water equations, Journal of Computational Physics 351, 

478-510 (2017). 

[28] J. Borggaard, T. Iliescu, Z. Wang, Artificial viscosity proper orthogonal decomposition, 

Mathematical and Computer Modelling 53, 269-279 (2011). 

[29] I. Babuska and J.M.~Melenk, The partition of unity method, International Journal for Numerical 

Methods in Engineering 40(4), 727-758 (1998). 

[30] V. Shcherbakov, E. Larsson, Radial basis function partition of unity methods for pricing vanilla 

basket options, Computers and Mathematics with Applications 71, 185-200 (2016). 

[31] D. Shepard, A two-dimensional interpolation function for irregularly-spaced data, in: 

Proceedings of the 23rd ACM National Conference, 517-524 (1968). 

[32] H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions 

of minimal degree, Advances in Computational Mathematics 4(1), 389-396 (1995). 

[33] M.D. Buhmann,  Radial Basis Functions: Theory and Implementation, Cambridge University 

Press, Cambridge (2003). 

[34] H. Wendland, Scattered Data Approximation, Cambridge University Press, Cambridge (2005). 



Communications on Applied Nonlinear Analysis 

ISSN: 1074-133X 

Vol 32 No. 9s (2025) 

  

2200 
https://internationalpubls.com 

[35] C.A. Micchelli, Interpolation of scattered data: Distance matrices and conditionally positive 

definite functions, Constructive Approximation 2, 11-22 (1986). 

[36] E. Larsson, V. Shcherbakov, A. Heryudono, A Least Squares Radial Basis Function Partition of 

Unity Method for Solving PDEs, SIAM Journal on Scientific Computing 39, A2538-A2563 (2017). 

[37] A. Heryudono, E. Larsson, A. Ramage, L. von Sydow, Preconditioning for radial basis function 

partition of unity methods, Journal of Scientific Computing 67, 1089-1109 (2016). 

[38] B. Forenberg, E. Larsson, N. Flyer, Stable computations with Gaussian radial basis functions, 

SIAM Journal on Scientific Computing 33, 869-892 (2011). 

[39] B. Forenberg, C. Piret, A stable algorithm for flat radial basis functions on a sphere, SIAM Journal 

on Scientific Computing 33, 60-80 (2008). 

[40] B. Fornberg, G. Wright, Stable computation of multiquadric interpolants for all values of the 

shape parameter. Computers and Mathematics with Applications 48, 853-867 (2004) 

[41] Q. Liang, F. Marche, Numerical resolution of well-balanced shallow water equations with 

complex source terms, Advances in Water Resources 32, 873-884 (2009). 

[42] T. Zhang, C.X. Zhan, B. Cai, C. Lin, X.M. Guo, An improved meshless artificial viscosity 

technology combined with local radial point interpolation method for 2D shallow water equations, 

Engineering Analysis with Boundary Element 133, 303-318 (2021). 

[43] D. Satyaprasada, S.N. Kuiry, S. Sundar, A shock-capturing meshless method for solving the one-

dimensional Saint-Venant equations on a highly variable topography, Journal of 

Hydroinformatics 25, 17-18 (2023). 

[44] S.  Sarra, A local radial basis function method for advection-diffusion-reaction equations on 

complexly shaped domains. Applied Mathematics and Computation 218, 9853-9865 (2012). 

[45] A. Bermudez, M.E. Vazquez, Upwind methods for hyperbolic conservation laws with source 

terms, Computers and Fluids 23(8), 1049-1071 (1994). 

[46] N. Goutal, F. Maurel, Electriciteé de France Service Applications de l'Electriciteé et 

Environnement & Workshop on Dam-Break Wave Simulation, In: Proceedings of the 2nd 

Workshop on Dam-Break Wave Simulation, Direction des études et recherches, Electricité de 

France, (1997). 

[47] E.F. Toro, Shock-Capturing Methods for Free-Surface Shallow Flows, Wiley and Sons (2001). 

[48] Benkhaldoun. F, Seaid. M, A simple finite volume method for the shallow water equations, 

Journal of Computational and Applied Mathematics 234, 58-72 (2010). 

[49] O. Delestre, Simulation du ruisselement d'eau de pluie sur des surfaces agricoles, PhD thesis 

(2010). 

[50] I. MacDonald, Analysis and computation of steady open channel flow, PhD thesis, University of 

Reading - Department of Mathematics (1996). 

[51] I. MacDonald, M.J. Baines, N.K. Nichols, P.G. Samuels. Analytic benchmark solutions for open-

channel flows. Journal of Hydraulic Engineering 123(11), 1041-1045 (1997).