## 3 Numerical Solution by Meshless Method  
### Formulation of the Local RBF Differentiation  
We employ the radial basis functions method, which is a meshless method originally proposed by [PERSON] [19], for the numerical solution of pollutant transport by the shallow water model. In the following, we outline the main features of this method. Radial basis function interpolation approximates the solution function by an expansion. According to the [PERSON] method, the solution is approximated on a set of \(N\) collocation nodes by a linear combination of local radial basis functions, as follows [20]:  
\[F(\mathbf{x},t)=\sum_{j\in s_{i,n}}\lambda_{F,j}\phi(\|\mathbf{x}-x_{j}\|, \epsilon) \tag{6}\]  
where \(\lambda_{F,j}\) represents the unknown coefficients to be calculated, \(\phi(x)\) is a radial basis function, \(\|\mathbf{x}-x_{j}\|\) is the distance between the nodes \(\mathbf{x}\) and \(x_{j}\), and \(s_{i,n}\) is a stencil that contains the index \(i\) and indices of the \(n\) nearest neighbor of the node \(x_{i}\). The stencil consists of the indices of node points falling within three times the distance to the nearest node. Thus, the number of nodes indexed by the stencil is variable, as it depends on the distribution of nodes in the computational domain. Many radial basis functions have been proposed in the literature. In the current study, we use the infinitely smooth multiquadric radial basis function defined as \(\phi(\|\mathbf{x}-x_{j}\|)=\sqrt{1-\epsilon\|\mathbf{x}-x_{j}\|^{2}}\), where \(\epsilon\
eq 0\) is a shape parameter that plays a substantial role in the resulting accuracy. Generally, the determination of the range of shape values can be obtained by using numerical tests. In this study, we also use the procedure proposed by [PERSON] in [21]. Therefore, an optimal shape parameter has been identified. In all simulations, \(\epsilon=\frac{\epsilon_{0}\sqrt{n}}{d_{min}}\), in which \(d_{min}\) is the minimum distance between two nodes and \(10^{-4}\leq\epsilon_{0}\leq 2\)\(\times 10^{-1}\). The numerical solution of an equation involving partial derivatives can be approximated by using the linear approximation of the partial derivatives at the \(\mathbf{x}\) nodes, based on local radial basis functions as follows:  
\[\mathcal{L}F(\mathbf{x},t)=\sum_{j\in s_{i,n}}\lambda_{F,j}\mathcal{L}\phi(\| \mathbf{x}-x_{j}\|,\epsilon) \tag{7}\]  
This can be rewritten in a more concise form as  
\[\mathcal{L}F^{[i]}=\Phi_{\mathcal{L}}^{[i]}\Lambda_{F}^{[i]} \tag{8}\]  
where \(\Delta_{F}^{[i]}=[\lambda_{F,1},\lambda_{F,2},\ldots,\lambda_{F,n}]^{T}\) is an \(n\times 1\) vector of unknown coefficients to be calculated and \(\Phi^{\mathcal{L}}\) is an \(n\times n\) matrix containing linear differentiation of the radial basis function \(\mathcal{L}\phi(\|x_{i}-x_{j}\|,\epsilon)\). Using Equation (6), the coefficient \(\Lambda_{F}\) can be expressed in terms of the function \(F\) as follows:  
\[\Lambda_{F}^{[i]}=(\Phi^{[i]})^{-1}F^{[i]} \tag{9}\]  
where \((\Phi^{[i]})^{-1}\) is the inverse of the matrix \(\Phi^{[i]}\). By removing the dependence on the expansion coefficients, the equation in (8) becomes:  
\[\mathcal{L}F^{[i]}=\Phi_{\mathcal{L}}^{[i]}(\Phi^{[i]})^{-1}F^{[i]} \tag{10}\]  
Note that in the governing equations the first and second derivatives are calculated for each node by using several matrix operations on \(n\times n\) matrices, which slows down the computations. In order to overcome the problem of the successive inversion of many \(n\times n\) matrices, [PERSON] et al. [20] reformulated the local interpolation and differentiation matrices in terms of a global sparse matrix by appending zeros to the proper entries. This reformulation was also used by [22] for the space-time derivatives. The entries of the global differentiation matrix \(\mathcal{D}_{\mathcal{L}}\) are then defined as follows:  
\[\mathcal{D}_{\mathcal{L}}(i,j)=\left\{ \begin{array}{ll}\Phi_{\mathcal{L}}^{[i]} (\Phi^{[i]})^{-1}&\text{if}j\in s_{i,nt}\\ 0&\text{otherwise}\end{array} \right.