## 3 Numerical Solution by Meshless Method
### Numerical Procedure  
In this section, we will describe the implementation of the LRBF method for shallow water equations coupled with the advection-diffusion equation of a non-active pollutant.  
First, a set \(N\) of distinct points, \(\{(x_{i},y_{i})\}_{i=1,N}\), are considered in the study area. The system (3) is approximated by the following equivalent system:  
\[\partial_{t}\mathbf{Q}=-\hat{\mathcal{F}}(\mathbf{Q})-\hat{\mathcal{G}}( \mathbf{Q})+\tilde{\mathbf{S}}_{b}(\mathbf{Q})+\mathbf{S}_{s}(\mathbf{Q}) \tag{14}\]  
where \(\hat{\mathcal{F}}\) and \(\hat{\mathcal{G}}\) are the convective flux and \(\hat{\mathcal{S}}_{b}\) is the term due to bottom variation. Using the local radial basis function meshless method to perform numerical evaluation, we obtain:  
\[\hat{F}(\mathbf{Q})\;\;=\;\partial_{x}\left[\begin{array}{c}\mathcal{D}_{x} q_{x}\\ \mathcal{D}_{x}(\frac{q_{x}}{h}+\frac{1}{2}gh^{2})\\ \mathcal{D}_{x}(\frac{q_{x}dy}{h})\\ \mathcal{D}_{x}\frac{q_{x}d_{x}}{h}+hK_{xx}\mathcal{D}_{xx}(\frac{q_{x}}{h}) \end{array}\right],\tilde{G}(\mathbf{Q})\;\;=\;\partial_{y}\left[\begin{array} []{c}\mathcal{D}_{y}q_{y}\\ \mathcal{D}_{y}(\frac{q_{x}dy}{h})\\ \mathcal{D}_{y}(\frac{q_{y}}{h}+\frac{1}{2}gh^{2})\\ \mathcal{D}_{y}\frac{q_{y}d_{x}}{h}+hK_{yy}\mathcal{D}_{yy}(\frac{q_{x}}{h}) \end{array}\right] \tag{15}\]  
and  
\[\tilde{\mathbf{S}}_{b}\left[\begin{array}{c}0\\ -gh\mathcal{D}_{x}z\\ -gh\mathcal{D}_{y}z\\ 0\end{array}\right] \tag{16}\]  
Thus, the numerical solution consists of the numerical evaluation of all the derivatives involved, and this is performed several times during the same time step. So, for reasons of computational efficiency, we first evaluate the derivation matrices \(Dx\), \(Dy\), \(Dxx\) and \(Dyy\), then the numerical fluxes.  
For time progression, the second-order Runge-Kutta scheme was employed and found to be reasonably stable. The formula used for the solution \(\mathbf{Q}^{n+1}\) at the \(n+1\) time level is as follows:  
\[\begin{split}\mathbf{Q}^{n+1}&=\frac{1}{2}(\mathbf{Q}^{n}+ \mathbf{Q}^{*})-\frac{\Delta t}{2}\big{[}\hat{\mathcal{F}}(\mathbf{Q}^{*})+ \hat{\mathcal{G}}(\mathbf{Q}^{*})-\tilde{\mathbf{S}}_{b}(\mathbf{Q}^{*})- \mathbf{S}_{s}(\mathbf{Q}^{*})\big{]}\\ \mathbf{Q}^{*}&=\mathbf{Q}^{n}-\Delta t\big{[}\hat{\mathcal{F}}(\mathbf{Q}^ {n})+\hat{\mathcal{G}}(\mathbf{Q}^{n})-\tilde{\mathbf{S}}_{b}(\mathbf{Q}^{n}) -\mathbf{S}_{s}(\mathbf{Q}^{n})\big{]}\end{split} \tag{17}\]  
In the above formula, \(\Delta t\) represents the time step; this needs to be chosen carefully to guarantee the stability of the scheme. In all our simulations, we chose \(\Delta t\) using the following formula:  
\[\Delta t=CFL\frac{d_{min}}{max(|\frac{q_{x}}{h}|+\sqrt{gh}I_{t}|\frac{q_{x}}{ h}|+\sqrt{gh})} \tag{18}\]  
where \(d_{min}\) is the smallest distance node and \(CFL\) is the current number taken equal 0.5.  
At the beginning of each step of the RK scheme, the friction effect is evaluated and used implicitly by the splitting scheme proposed in [23]. Keeping only the local acceleration and bed friction terms, the momentum equation is approximated as follows:  
\[\frac{d\mathbf{Q}}{dt}=\mathbf{S}_{f}^{n+1} \tag{19}\]  
Since the friction term is only involved in the momentum equation, only the flow rates need to be evaluated. Thus, the following implicit scheme is used [24]:  
\[q_{i}^{n+1}=q_{i}^{n}+\frac{\Delta t\mathbf{S}_{f}}{1-2\Delta t\mathbf{S}_{f} (q_{i}^{n})^{-1}} \tag{20}\]  
This updated flux is used as an initial condition for the operators in Equation (17).