<figure><img src="image_1.png" /></figure>

HAL
open science

# A Numerical Framework for Pollutant Transport in Shallow-Water Flows: Application to the Niger River in Bamako

Abdoulaye Samaké, Mahamadou Alassane, Amadou Mahamane

► To cite this version:

Abdoulaye Samaké, Mahamadou Alassane, Amadou Mahamane. A Numerical Framework for Pollutant Transport in Shallow-Water Flows: Application to the Niger River in Bamako. Journal of Applied Mathematics, 2025, 2025, 10.1155/jama/6702536. hal-05109537

HAL Id: hal-05109537

https://hal.science/hal-05109537v1

Submitted on 12 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Research Article

# A Numerical Framework for Pollutant Transport in Shallow-Water Flows: Application to the Niger River in Bamako

**Abdoulaye Samaké ID, Mahamadou Alassane ID, and Amadou Mahamane ID**

Faculty of Sciences and Techniques, University of Sciences, Techniques and Technologies of Bamako, Bamako, Mali

Correspondence should be addressed to Abdoulaye Samaké; abdoulaye.samake@usttb.edu.ml

Received 20 September 2024; Accepted 25 April 2025

Academic Editor: Kamal Kumar

Copyright © 2025 Abdoulaye Samaké et al. Journal of Applied Mathematics published by John Wiley & Sons Ltd. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

We propose a unified numerical framework for the transport of passive pollutants by shallow-water flows. The mathematical model we consider for describing this phenomenon results in the coupling of the hydrodynamic shallow-water equations with a two-dimensional advection-diffusion equation governing the pollutant transport. The numerical implementation of this hyperbolic model, based on the finite element method, is achieved using a multiphysics modeling and simulation toolbox featured by Feel++, a versatile C++ library applying the Galerkin methods for solving partial differential equations. Numerical experiments targeted on arsenic, cadmium, and lead, heavy metals among the most harmful to human health, are presented as part of a practical application on the Niger River in Bamako. The framework is validated on the basis of RMSE and MAE metrics, some of the most commonly used error measures in linear regression, using observational data. These indicators, estimated below 5% of the observed mean value, support the reliability and accuracy of the numerical model in capturing pollutant dynamics under flow conditions. The simulation results highlight the predictive effectiveness of this framework and provide better insight into pollution patterns in the scrutinized river section.

**Keywords:** finite element; numerical simulation; shallow-water equations; transport-diffusion; water pollution

## 1. Introduction

Water is a natural resource essential for all life on Earth. Protecting water and the ecosystems that depend on it is therefore of utmost importance for humans and nature. Widespread contamination of water resources observed worldwide in recent years makes the water quality situation worrying.

Water pollution [1] is nowadays a major global challenge that humanity is facing, as it represents an existential threat to human health and environmental balance and hinders economic development [2]. It refers to the release of harmful substances, including chemicals and microorganisms, into water bodies, making them unsafe for human use and related activities and disrupting ecosystem functions [3]. Water pollution, including surface water pollution [4] and groundwater pollution [5], is both natural and anthropogenic in origin. However, it turns out that the most common sources of water pollution are anthropogenic, such as sewage

discharge, agricultural runoff transporting agrochemical (pesticides and fertilizers) residues, urban runoff, and industrial activities. Water pollutants, such as pathogens, organic matter, nutrients, sediments, thermal substances, and radioactive materials, originate from point sources or nonpoint (or dispersed) sources. Point source pollutants are discharged into water bodies from discrete identifiable conveyances, such as factory facilities, sewage treatment plants, or storm drains. In contrast, nonpoint sources of pollution release large quantities of pollutants into water bodies from a broad unconfined area. They are more diffused and are primarily caused by runoff.

Contaminated water and poor sanitation are associated with the spread of several bacterial, viral, and parasitic diseases, such as cholera, diarrhoea, dysentery, Hepatitis A, encephalitis, typhoid, polio, amoebiasis, and giardiasis [6]. It also fuels malnutrition and, in particular, childhood stunting [7]. In 2019, about 1.4 million premature deaths worldwide stemmed from water pollution [8], making it the
second leading cause of premature death in the world after air pollution. It has been reported that approximately 829,000 people worldwide die every year from diseases directly attributable to unsafe drinking water, poor sanitation, and hygiene practices [9]. This number includes about 300,000 children under 5 years of age, representing 5.3% of all deaths in this age group. Developing countries are the most affected by water pollution, since up to 80% of illnesses in this part of the world are attributed to contaminated water and other water-related issues [10]. The African continent, in particular, bears the brunt of the severe consequences of this kind of pollution, overwhelming healthcare infrastructure and hampering efforts to improve human capital and economic growth [11]. As a result, it has become essential to track down and effectively survey the slightest discharge into surface water and groundwater.

Over the decades, major decisions have been taken to address the health and environmental hazards caused by water pollution. They consisted in drawing up and applying legislation and regulations [12] governing water use, implementing pollution control strategies, and developing advanced technical tools for monitoring groundwater and surface water quality [13]. These efforts combine traditional and modern methods to deeply detect pollutants and better identify their sources in order to support authorities in making decisions for preserving resources. Conventional approaches [14] include analysis in laboratory, portable field instrumentation and testing, and visual field inspection [15]. They are far too time-consuming and costly, due to the high maintenance of laboratory facilities and the use of chemical materials, and are particularly inefficient for many practical applications. Modern methods involve advanced and technologically innovative approaches, such as remote sensing [16], Earth observation satellite imagery [17], cyberphysical systems (CPSs) [18], Internet of things (IoT)-based smart sensors [19], and other sophisticated optical techniques [20]. They provide reliable, continuous, and near-real-time data for accurate and high spatial and temporal resolution water quality monitoring. Despite the fact that these smart solutions are more cost-effective, robust, and easier to maintain, they are challenging to implement in developing countries [21], particularly in sub-Saharan Africa, due to the paucity of financial resources to acquire them on the one hand and the lack of technical capabilities to implement them effectively on the other. This clearly demonstrates that these approaches are not concretely viable for these countries, and therefore, it is imperative to investigate more suitable alternative methods.

Relevant alternatives to these conventional approaches include mathematical modeling and numerical simulation [22]. They are fundamental and efficient tools widely used for understanding and predicting the dynamics of pollutants in aquatic environments [23]. They allow for assessing the effects of these substances on public health, environment, and ecosystems and to assist public authorities in formulating effective regulations and policies to reduce pollution, preserve resources, and support sustainable development. They deal with the transport and diffusion of pollutants in water bodies and address chemical processes. Rapid growth

in computational power and recent advances in robust numerical methods [24] have contributed significantly to the remarkable success of these tools.

The Niger River [25] is the principal river in West Africa, see Figure 1. It has a main course length of about 4200 km, making it the third-longest river in Africa after the Nile and Congo rivers. It rises on the mountain range of Fouta Djalon (upper Niger), a highland region of west-central Guinea with an average altitude of about 1100 m. It flows northeastward forming in central regions of Mali, during the wet season, a vast floodplain (inner Delta) which absorbs a large part of its hydraulic potential [26]. On the fringes of the Sahara, downstream of the Inner Delta, the Niger turns and flows due southeast (middle Niger) into Nigeria. It then continues its course (lower Niger) to the Atlantic Ocean where it finally runs into the Gulf of Guinea after its confluence with the Benue River, its largest tributary. The Niger River Basin [27] is the fourth largest in Africa after the Congo, Nile, and Lake Chad and covers an area of about 2.2 million km², including 1.5 million km² of an active hydrological basin. It is shared by nine West and Central African countries: Benin, Burkina Faso, Cameroon, Chad, Côte d'Ivoire, Guinea, Mali, Niger, and Nigeria. This portion of land drained by the Niger is an extraordinary asset for these countries, especially for the three riparian countries, which account for more than 80% of its area: Mali, Niger, and Nigeria. Indeed, the Niger Basin abounds with immense development potential [28], covering irrigation for food production, hydroelectric power generation, transport, industrial and mining development, and other related activities, such as trade, investment in communications, and increased mobility of people and goods [29]. The activities related to irrigation include agriculture, fishing, pastoralism, and market gardening. The resources of this basin are therefore vital for sustainable and integrated regional development and in particular for the riparian population, estimated at over 130 million people [30] with an average annual growth rate of around 3%. However, despite these important social and economic opportunities in the Niger Basin, it is seriously threatened by major phenomena such as the continuing decline in its flows, unsustainable use of its resources, desertification, silting of its bed, and upstream pollution due to various activities, such as industry, mining, agriculture, and urban enterprises, resulting from the combined effects of climate change and strong demographic growth [27]. The assaults of these various anthropogenic polluting activities [31] on the Niger Basin over many years not only weaken ecosystems and biodiversity [32] but also undermine its capacity to meet the expectations of the riparian population. This situation contributes to their further impoverishment and seriously imperils their health [33].

In Mali, the Niger River [27] follows a curve of 1700 km and flows through several major cities of the country, including Bamako, Koulikoro, SAgou, Mopti, Timbuktu, and Gao. It is one of the few permanent sources of water in the heart of the arid and semiarid Sahel [34] and constitutes the principal source of drinking water [35] for the city of Bamako, the political and economic capital of the country. The portion of the Niger Basin in Mali, the largest share in any country,
<figure><img src="image_4.png" /><figcaption>FIGURE 1: A map of the Niger River. Source: Wikimedia Commons (CC BY-SA 3.0).</figcaption></figure>

makes up around 30% of the total basin area [27] and is home to approximately 75% of the Malian population [36], many of whom live in Bamako on the banks of the Niger. It represents a vital artery for Mali, making an invaluable contribution to the social and economic development, since a large majority of the population depends directly or indirectly on its resources. Despite its immense potential for sustainable development of Mali, the Niger River is severely strained along its course through the country, particularly in the major urban centers, especially Bamako, by several damaging factors. These include overuse of river banks and shallows and pollution generated by urban wastewater. These matter are highly dense in micro-plastics and include domestic sewage, stormwater runoff (urban and agricultural), and wastewater from industrial units and from certain informal economic activities, such as artisanal gold mining, artisanal dyeing and tanning, and sand and gravel extraction [37]. Some of the resulting pollutants, at a certain concentration threshold, can seriously damage human and animal health, harm the aquatic environment, and contaminate consumer products from market gardening and irrigated farming on its banks. They include chemicals, such as phosphate (PO₄³⁻), nitrate (NO₃⁻), ammonium (NH₄⁺), and sulfate (SO₄²⁻); microbials, such as fecal coliforms, especially the major species in this group referred to as *Escherichia coli* (*E. coli*); and heavy metals, such as chromium (Cr), copper (Cu), zinc (Zn), arsenic (As), mercury (Hg), cyanide (CN), lead (Pb), cadmium (Cd), manganese (Mn), and nickel (Ni). As the self-purification capacity [38] of Niger River has been severely impaired by the nature and the extent of the pollution, it is essential to identify, track down, and effectively monitor the multiple sources of pollutant discharge into this river in order to gain in-depth knowledge of their dynamics and thus support and accommodate the authorities' decision-making process for preserving the river and its resources [39].

For this purpose, we introduced in this paper a unified numerical framework for the modeling and simulation of pollutants in the Niger River in Bamako in order to achieve

a better insight into the pollution patterns. In terms of specific contribution to the current literature, this paper provides the scientific community interested in this topic and water pollution issues with a nontrivial application based on new robust algorithms and advanced numerical tools. The model we consider consists of a coupling of the shallow-water equations, which govern river hydrodynamics, with an advection-diffusion equation for the pollutant transport by water flows.

This paper is organized as follows. Section 2 describes a coupled mathematical model of shallow-water equations with pollutant transport. Section 3 deals with the numerical implementation of the model. In Section 4, we present the numerical experiments, discuss simulation results, and analyze their validity. The main conclusions are recapped in Section 5.

## 2. Mathematical Model

Mathematical modeling and computer simulations are crucial tools for predicting and assessing water pollution in natural streams. These powerful numerical tools provide an accurate and reliable estimation of the spatial and temporal distribution of pollutants, essential for an effective prediction of water pollution and the development of better control strategies and risk management. The transport of pollutants in river systems is a complex process [40] that depends on both the flow characteristics of the fluid and the physical and chemical properties of the pollutants. It can be described by a mathematical model with two different but dependent variables, namely, a hydrodynamic variable representing the fluid flow dynamics and a scalar transport variable defining the pollutant concentration.

In the current work, the water flow in the Niger River is governed by the two-dimensional (2D) shallow-water equations [41], often known as the Saint-Venant equations. This hydrodynamic model, fairly well known in the literature, is
derived from the integration of the three-dimensional incompressible Navier-Stokes equations over the fluid depth under the assumptions of hydrostatic pressure distribution and the uniform velocity distribution along the vertical direction. This system of equations [42] is a suited model widely used for describing the flows in environments where the fluid depth is much smaller than the horizontal scale of motion, which absolutely fits the context of this study. On the other hand, for modeling the dynamics of pollutants in water flows, we consider a 2D advection-diffusion equation [43]. Thus, the overall model for fully describing the transport of pollutants by shallow-water flows, regarding concurrently flow and transport phenomena, results in the coupling [44] of the shallow-water equations with the advection-diffusion equation. The resulting system of partial differential equations (PDEs) is written in the so-called conservative form as follows:

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \frac{\partial h}{\partial t} + \operatorname{div}(\mathbf{q}) = 0, \\[1.5ex]
\displaystyle \frac{\partial \mathbf{q}}{\partial t} + \operatorname{div}\left(\frac{1}{h}\mathbf{q} \otimes \mathbf{q}\right) - v\Delta\mathbf{q} + gh\nabla(h + h_b) + (\mathbf{k} \times f_c\mathbf{q}) + \frac{1}{\rho_w}\tau_b - \frac{1}{\rho_w}\tau_s = \mathbf{f}, \\[1.5ex]
\displaystyle \frac{\partial \phi}{\partial t} + \operatorname{div}\left(\frac{1}{h}\mathbf{q}\phi\right) - \operatorname{div}\left(D\nabla\phi - \frac{1}{h}D\nabla h\phi\right) = S_c,
\end{array}
\right. 
\quad (1) $$

where

i. $t$: time

ii. $h$: water depth

iii. $\mathbf{q} = hu$: flow discharge, where $\mathbf{u}$ is the depth-averaged flow velocity

iv. $\mathbf{k}$: unit vector in vertical direction

v. $\tau_b$: bottom shear stress

vi. $\tau_s$: surface stress

vii. $\mathbf{f}$: energy per unit mass due to additional body forces

viii. $h_b$: bottom bathymetry

ix. $v$: kinematic viscosity

x. $f_c$: Coriolis parameter

xi. $g$: gravitational acceleration

xii. $\phi = hc$: pollutant concentration, where $c$ is the depth-averaged concentration of the pollutant

xiii. $\mathbf{D}$: diffusion tensor of the pollutant

xiv. $S_c$: depth-averaged pollutant source term

xv. $\otimes$: dyadic product

All the fields and parameters of Model (1), including those to be introduced later, and their units based on the international system (SI) are listed in Table 1. The first equation of System (1) above represents the mass conservation of

the water. It is also referred to as the continuity equation. The second equation refers to the momentum conservation, while the third equation represents the transport of pollutants by the water flows. Model (1) can easily be extended to several pollutants by replicating the last equation as many times as the number of pollutants considered while updating the field and physical parameters specific to each of them. The pollutants are supposed to be passive (nonreactive). This assumption is based on the fact that the pollutants of interest have low reactivity in water under typical environmental conditions. Indeed, they are not subject to significant chemical reactions or interactions that could alter their behavior as they are transported by the water flow over the timescales of interest. This justifies neglecting reaction processes, much slower than the transport by advection and diffusion. We assume that they do not induce any feedback to the water flow. This assumption is reasonable for this study, since the pollutant concentrations are typically much lower than the total volume of flowing water. Consequently, pollutants do not significantly alter the density or viscosity of the water, primary factors influencing flow dynamics.

The bottom shear stress $\tau_b$ is computed using the Manning law [45] as follows:

$$ \tau_b = \frac{\rho_w}{h^2} c_b \|\mathbf{q}\|\mathbf{q} \quad \text{with} \quad c_b = \frac{gn_b^2}{h^{1/3}}, $$

where $\rho_w$ is the water density, $c_b$ is the bed roughness coefficient, and $\|\mathbf{q}\|$ is the Euclidean norm of the flow discharge. The parameter $n_b$ refers to the Manning roughness coefficient [46]. It defines the resistance, shear, and drag forces included, to flow by the channel. It is essential for the model, since it can significantly influence the numerical results. The estimation of this coefficient is particularly challenging due to the fact that it depends on several factors, such as the surface bed roughness, vegetation density, sinuosity, flow discharge, and potential obstacles [47]. This has been widely investigated and the subject of numerous results, generally based on empirical formulations and field measurements. The setting of this parameter for this application will be discussed in Section 4.2.

The surface stress $\tau_s$ is expressed as a quadratic function of the wind velocity $\mathbf{w}$ as follows [48]:

$$ \tau_s = \rho_w c_s \|\mathbf{w}\|\mathbf{w}, $$

where $c_s$ is the friction coefficient of the wind. It is usually defined by

$$ c_s = 10^{-3} \rho_a (0.75 + 0.067 \|\mathbf{w}\|), $$

where $\rho_a$ is the air density. The additional body forces, causing the specific energy $\mathbf{f}$, include atmospheric pressure gradient, density gradients, radiation stress, and tidal stresses [49]. They will not be considered in this work. The Coriolis parameter $f_c$ is defined by [50]

$$ f_c = 2\omega \sin(\varphi), $$
<table><caption>TABLE 1: Fields and physical parameters of the model.</caption><thead><tr><th>Symbol</th><th>Unit</th><th>Meaning</th></tr></thead><tbody><tr><td>h</td><td>m</td><td>Water depth</td></tr><tr><td>q</td><td>m<sup>2</sup>/s</td><td>Unit discharge</td></tr><tr><td>u</td><td>m/s</td><td>Depth-averaged velocity</td></tr><tr><td>w</td><td>m/s</td><td>Wind velocity</td></tr><tr><td>f</td><td>m<sup>2</sup>/s<sup>2</sup></td><td>Energy per unit mass due to body forces</td></tr><tr><td>τ<sub>b</sub></td><td>kg/m/s<sup>2</sup></td><td>Bottom shear stress</td></tr><tr><td>τ<sub>s</sub></td><td>kg/m/s<sup>2</sup></td><td>Surface stress</td></tr><tr><td>φ</td><td>kg/m<sup>2</sup></td><td>Concentration of pollutant</td></tr><tr><td>c</td><td>kg/m<sup>3</sup></td><td>Depth-averaged concentration of pollutant</td></tr><tr><td>D</td><td>m<sup>2</sup>/s</td><td>Pollutant diffusion tensor</td></tr><tr><td>S<sub>c</sub></td><td>kg/m<sup>2</sup>/s</td><td>Pollutant sources</td></tr><tr><td>k</td><td>1</td><td>Unit vector in the vertical direction</td></tr><tr><td>ν</td><td>m<sup>2</sup>/s</td><td>Kinematic viscosity</td></tr><tr><td>f<sub>c</sub></td><td>1/s</td><td>Coriolis parameter</td></tr><tr><td>g</td><td>m/s<sup>2</sup></td><td>Gravitational acceleration</td></tr><tr><td>h<sub>b</sub></td><td>m</td><td>Bottom bathymetry</td></tr><tr><td>n<sub>b</sub></td><td>s/m<sup>1/3</sup></td><td>Manning roughness coefficient</td></tr><tr><td>c<sub>s</sub></td><td>m<sup>1/2</sup>/s</td><td>Wind friction coefficient</td></tr><tr><td>ω</td><td>rad/s</td><td>Angular speed of Earth rotation</td></tr><tr><td>φ</td><td>rad</td><td>Latitude of the site</td></tr><tr><td>ρ<sub>a</sub></td><td>kg/m<sup>3</sup></td><td>Air density</td></tr><tr><td>ρ<sub>w</sub></td><td>kg/m<sup>3</sup></td><td>Water density</td></tr></tbody></table>

where $\omega$ is the angular speed of the Earth rotation and $\phi$ is the geographic latitude of the site. The diffusion tensor $\mathbf{D}$ is a $2 \times 2$ matrix defined by

$$ \mathbf{D} = \begin{pmatrix} D_{xx} & D_{xy} \\ D_{yx} & D_{yy} \end{pmatrix}, $$

where $D_{xx}$, $D_{xy}$, $D_{yx}$, and $D_{yy}$ correspond to the diffusion coefficients of the pollutant in each combination of x- and y-directions. The diagonal elements $D_{xx}$ and $D_{yy}$ represent the diffusion coefficients along the principal directions, and the off-diagonal terms $D_{xy}$ and $D_{yx}$ reflect the correlation between diffusion coefficients along the respective pair of principal directions. In general, the components of the diffusion tensor $\mathbf{D}$ account for molecular diffusion and turbulent mixing. They depend on the nature of the pollutant under consideration and the flow characteristics, including water depth, flow velocity, bottom roughness, wind, and vertical turbulence. The turbulent mixing is the dominant component of the dispersion in river flows [40]. It is thereby feasible to neglect the molecular diffusion, a small-scale process that occurs much more slowly than turbulent mixing. Thus, in longitudinal (or streamwise) and transverse reference system, the diffusion tensor can be expressed in the form of a diagonal matrix as follows:

$$ \begin{pmatrix} D_L & 0 \\ 0 & D_T \end{pmatrix}, $$

where $D_L$ and $D_T$ are longitudinal (along the flow) and transverse (perpendicular to the flow) dispersion coefficients, respectively. These coefficients are estimated by

$$ D_L = \alpha_L h u^* \quad \text{and} \quad D_T = \alpha_T h u^*, \quad \text{with} \quad u^* = \frac{n_b \sqrt{g}}{h^{7/6}} \|\mathbf{q}\|, $$

where $u^*$ represents the shear velocity and $\alpha_L$ and $\alpha_T$ are dimensionless constants quantifying the magnitude of longitudinal and transverse dispersion, respectively. Typical values of these constants for river flows [51] are $\alpha_L = 5.93$ and $\alpha_T = 0.15$. In order to compute the components of the diffusion tensor $\mathbf{D}$ in the Cartesian coordinate system, the following expressions hold [52]:

$$ \begin{cases} D_{xx} = D_L \cos^2(\theta) + D_T \sin^2(\theta), \\ D_{xy} = D_{yx} = (D_L - D_T) \cos(\theta) \sin(\theta), & \text{with } \cos(\theta) = \frac{q_x}{\|\mathbf{q}\|} \quad \text{and} \quad \sin(\theta) = \frac{q_y}{\|\mathbf{q}\|}, \\ D_{yy} = D_L \sin^2(\theta) + D_T \cos^2(\theta), \end{cases} $$

where $\theta = \arctan(q_y/q_x)$ represents the angle between the flow direction and the x-axis.

Let $\Omega$ be a bounded domain of $\mathbb{R}^2$, representing the studied area. In order to address the governing Model (1), defined in the spatial domain $\Omega$ over the time interval $(0, T)$, from a practical point of view, it must necessarily be supplemented by suitable boundary and initial conditions. In general situations, the boundary conditions depend not only on the physical nature of the problem and the data availability but also on the flow regime, which is determined by the Froude number defined as

$$ Fr = \frac{\|\mathbf{u}\|}{\sqrt{gh}}. $$

Based on this important dimensionless parameter, the flow is characterized as subcritical (or fluvial) if $Fr < 1$, critical when $Fr \approx 1$ or supercritical (or torrential) if $Fr > 1$.

For natural water flows, two types of boundaries are usually considered, namely, the water (or open) boundaries, which need to be treated carefully, and land (or solid) boundaries. Let $\partial\Omega$ be the boundary of $\Omega$ and $\mathbf{n}$ be the unit outward normal to $\partial\Omega$. We consider a partition of $\partial\Omega$ formed by the subsets $\partial\Omega_I$, $\partial\Omega_O$, and $\partial\Omega_L$, representing the inflow (or upstream) boundary, outflow (or downstream) boundary, and land (or lateral) boundaries, respectively. They can be defined as follows:

$$ \partial\Omega_I = \{\mathbf{x} \in \Omega : \mathbf{u} \cdot \mathbf{n} < 0\}, $$

$$ \partial\Omega_O = \{\mathbf{x} \in \Omega : \mathbf{u} \cdot \mathbf{n} > 0\}, $$

$$ \partial\Omega_L = \{\mathbf{x} \in \Omega : \mathbf{u} \cdot \mathbf{n} = 0\}. $$
The flow regime is supposed to be subcritical in this study, since the Froude number of river flows is typically less than 1. So, according to the literature, the flow boundary conditions commonly prescribed for such a flow regime are defined as follows [49]:

On $\partial\Omega_1$, given discharge

$$q \cdot n = \hat{q} \cdot n.$$

On $\partial\Omega_O$, given water depth

$$h = \hat{h}.$$

On $\partial\Omega_L$, no-slip (or adherence) condition

$$q = 0,$$

where $\hat{q}$ and $\hat{h}$ are specified. These boundary conditions are designed to ensure a realistic model behavior and the stability of numerical methods for approximating the solution. The prescribed discharge at the upstream boundary, driven by the availability of these inflow data, reflects the real-world situation in river systems where the inflow rate is commonly controlled or estimated on the basis of observations or predictions. This constraint, often preferred over prescribing water depth directly, makes the model more stable and physically consistent. On the other hand, the prescribed water depth at the downstream boundary is motivated by the fact that these outflow data are available and that their application is consistent with natural conditions and prevents unrealistic reflections. This ensures the stability and physical relevance of the model. Finally, the no-slip condition at lateral boundaries is a realistic assumption, since it reflects the physical behavior of the water in contact with solid boundaries.

The initial conditions for the flow variables are given by

$$q(\mathbf{x}, 0) = q_0(\mathbf{x}) \quad \text{and} \quad h(\mathbf{x}, 0) = h_0(\mathbf{x}) \quad \forall \mathbf{x} \in \Omega,$$

where the variable $\mathbf{x} = (x, y)$ describes the spatial dimensions and $q_0$ and $h_0$ are specified. Added to these are boundary conditions for pollutant flows from sewage collectors, dumps, and landfills discharging into the river. Let $\partial\Omega_{L_p} \subset \partial\Omega_L$ be the inlet boundaries for pollutant flows from sewage collectors, dumps, and landfills discharging into the river. The boundary and initial conditions associated with the advection-diffusion transport model are given by

$$\phi(\mathbf{x}, t) = \hat{\phi}(\mathbf{x}, t) \quad \forall \mathbf{x} \in \partial\Omega_{L_p} \quad \text{and} \quad \phi(\mathbf{x}, 0) = \phi_0(\mathbf{x}) \quad \forall \mathbf{x} \in \Omega,$$

where $\hat{\phi}$ and $\phi_0$ are specified.

### 3. Numerical Implementation

3.1. Background. The equations of the coupled Model (1) form a hyperbolic system of nonlinear conservation laws with source terms [53]. The numerical solution of such a system remains a major challenge, since it is subject to

strong discontinuities (shocks), mainly due to the nonlinear nature of the system and its hyperbolic character. Thus, the investigation of efficient and robust numerical methods based on hyperbolic systems has attracted increasing interest in recent decades. The main purpose consists of introducing mathematical and numerical techniques capable of effectively computing the solutions, including those of high-order accuracy in time and space over irregular bed topography with complex geometry, of this type of system of PDEs. Accordingly, several discretization methods and numerical schemes [49] for solving shallow-water equations have been developed. Among them, the most commonly used and widely discussed in the literature include the finite element method (FEM) and finite volume method (FVM).

Despite the adherence of the FVM [54] to the physical properties of the flows governed by System (1), coupled with the impressive benefits of this method for solving such systems, including its stability and local and global mass conservation properties, this work is focused on the FEM [55]. Indeed, this method has been used extensively in the study of general hyperbolic systems of conservation laws, involving a wide range of complex engineering and mathematical problems, in particular the shallow-water equations [56]. The pervasiveness of this approach is motivated by numerous advantages it offers, including the flexibility to handle complex and irregular geometric shapes, material properties and boundary conditions, natural support for higher-order discretizations on unstructured meshes, and the ability to allow adaptive mesh refinement for accurate solutions in specific regions of interest. In addition, for many decades, this method has been thoroughly studied, verified, and validated, contributing largely to its success and reliability, especially for highly demanding industrial applications. Furthermore, the implementation of the FEM is strongly supported by several versatile software packages based on powerful and modern computing tools. They include Feel++ [57], FEniCS [58], FreeFem++ [59], deal.II [60], GetFEM [61], and MFEM [62].

3.2. Feel++. The numerical implementation of this framework is performed using Feel++ [63], a versatile open-source finite element library designed for solving PDEs using the arbitrary-order Galerkin methods. This software project provides a domain-specific embedded language (DSEL) in C++, allowing a high-level mathematical expressivity for variational formulations. Feel++ supports both structured and unstructured meshes, continuous Galerkin (CG) and discontinuous Galerkin (DG) formulations, as well as a wider range of finite elements, including Lagrange, Hermite, Nedelec, and Raviart-Thomas elements. It integrates interfaces to advanced linear and nonlinear solvers as well as scalable preconditioners from powerful and flexible scientific computing suites. Feel++ supports seamless high-performance computing (HPC) and enables advanced features, including hp-adaptivity, robust boundary condition handling, isoparametric mapping, and large-scale nonlinear multiphysics simulations. It provides several toolboxes tailored for multiphysics modeling and simulation of various scientific and engineering problems, constituting an alternative approach to the
definition of the underlying equations. These toolboxes are based on two complementary files: JSON and CFG. The JSON file is used to define the physical quantities of the model, such as material properties, loads, constraints, sources, and fluxes. The CFG file is used to configure general options, including mesh file, output directory, and solver.

The numerical implementation proposed here relies on a particular toolbox of Feel++ called coefficient form PDEs (CFPDEs), similar to the interactive interface with the same name of COMSOL Multiphysics [64]. It is based on the representation of the PDEs in terms of their coefficient functions. This tool provides a powerful environment for implementing PDEs expressed in the following generic form:

$$d \frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot (-c \nabla \mathbf{u} - \alpha \mathbf{u} + \gamma) + \beta \cdot \nabla \mathbf{u} + a \mathbf{u} = \mathbf{f}, \quad (2)$$

where

i. **u**: unknown field

ii. **d**: damping or mass coefficient

iii. **c**: diffusion coefficient

iv. **α**: conservative flux convection coefficient

v. **γ**: conservative flux source term

vi. **β**: convection coefficient

vii. **a**: absorption or reaction coefficient

viii. **f**: scalar or vectorial

### 3.3. Coupled System Implementation.
In order to implement this framework using the CFPDE toolbox, it is imperative to rewrite each equation of System (1) in the form of Equation (2). This crucial step enables one to identify the various coefficients to be introduced into the toolbox implementation JSON file in accordance with the established guidelines. So the first equation of System (1) can be rewritten as

$$\frac{\partial h}{\partial t} = - \operatorname{div}(\mathbf{q}). \quad (3)$$

From (3), the CFPDE coefficients are identified as follows:

$$d = 1,$$

$$\mathbf{f} = - \operatorname{div}(\mathbf{q}).$$

Following the same approach, the second equation of System (1) can be rewritten as

$$\begin{aligned} & \frac{\partial \mathbf{q}}{\partial t} + \operatorname{div} \left( \frac{1}{h} \mathbf{q} \otimes \mathbf{q} \right) - v \Delta \mathbf{q} = \mathbf{f} - gh \nabla (h + h_b) - (\mathbf{k} \times f_c \mathbf{q}) - \frac{1}{\rho_w} \boldsymbol{\tau}_b + \frac{1}{\rho_w} \boldsymbol{\tau}_s \\ & \Rightarrow \frac{\partial \mathbf{q}}{\partial t} + \operatorname{div} \left( -v \nabla \mathbf{q} + \frac{1}{h} \mathbf{q} \otimes \mathbf{q} \right) \\ & = \mathbf{f} - gh \nabla (h + h_b) - (\mathbf{k} \times f_c \mathbf{q}) \\ & \quad - \frac{1}{\rho_w} \boldsymbol{\tau}_b + \frac{1}{\rho_w} \boldsymbol{\tau}_s. \end{aligned} \quad (4)$$

From (4), the following relationships for the CFPDE coefficients hold

$$d = 1,$$

$$c = v,$$

$$\boldsymbol{\gamma} = \frac{1}{h} \mathbf{q} \otimes \mathbf{q},$$

$$\mathbf{f} = \mathbf{f} - gh\nabla(h + h_b) - (\mathbf{k} \times f_c \mathbf{q}) - \frac{1}{\rho_w} \boldsymbol{\tau}_b + \frac{1}{\rho_w} \boldsymbol{\tau}_s.$$

Finally, let us consider the third equation of System (1). It can be rewritten as

$$\begin{aligned} & \frac{\partial \phi}{\partial t} + \operatorname{div} \left( \frac{1}{h} \mathbf{q} \phi \right) - \operatorname{div} \left( \mathbf{D} \nabla \phi - \frac{1}{h} \mathbf{D} \nabla h \phi \right) = S_c \\ & \Rightarrow \frac{\partial \phi}{\partial t} + \operatorname{div} \left( -\mathbf{D} \nabla \phi + \frac{1}{h} \mathbf{D} \nabla h \phi \right) + \frac{1}{h} \mathbf{q} \cdot \nabla \phi + \operatorname{div} \left( \frac{1}{h} \mathbf{q} \right) \phi = S_c \\ & \Rightarrow \frac{\partial \phi}{\partial t} + \operatorname{div} \left( -\mathbf{D} \nabla \phi + \frac{1}{h} \mathbf{D} \nabla h \phi \right) + \frac{1}{h} \mathbf{q} \cdot \nabla \phi \\ & \quad + \left( \frac{1}{h} \operatorname{div}(\mathbf{q}) - \frac{1}{h^2} \mathbf{q} \cdot \nabla h \right) \phi = S_c. \end{aligned} \quad (5)$$

By identifying (5) with the general form of CFPDEs, it follows

$$d = 1,$$

$$c = \mathbf{D},$$

$$\boldsymbol{\alpha} = - \frac{1}{h} \mathbf{D} \nabla h,$$

$$\boldsymbol{\beta} = \frac{1}{h} \mathbf{q},$$

$$\mathbf{a} = \frac{1}{h} \operatorname{div}(\mathbf{q}) - \frac{1}{h^2} \mathbf{q} \cdot \nabla h,$$

$$\mathbf{f} = S_c.$$

Once the coefficients of the various equations have been identified, they can be implemented in a JSON file and completed with loads, including forces and moments, as well as physical constraints (boundary conditions) and initial conditions. The flowchart of the key processes involved in the solution of Model (1) using the CFPDE toolbox is shown in
<figure><img src="image_5.png" /><figcaption>FIGURE 2: Flowchart of the processes involved in the numerical solution of the model.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 3: Map of Bamako including the Niger River. Source: OpenStreetMap.</figcaption></figure>

Figure 2. The first two steps in this flowchart concern the model implementation. The preprocessing step includes mesh generation, preparation of data related to material properties, boundary and initial conditions, and configuration of model-specific physical parameters and simulation settings. The first task in the simulation step involves importing the mesh, configuration files, and other data previously processed. The post-processing includes visualization, analysis, and validation of the numerical results.

## 4. Numerical Experiments

4.1. Computational Mesh Generation. This crucial operation began with the extraction of the map data of the section of the Niger River course in Bamako, as shown in Figure 3, from OpenStreetMap (OSM) [65] using the mirror database Overpass API (Application Programming Interface) [66]. This portion of the riverbed, extended over 17 km with a surface area of 23.58 km², represents the computational domain. The resulting map includes features and geospatial data and are referenced in the World Geodetic System-1984 (WGS 84). It has been processed and edited using JOSM, a powerful standalone and extensible desktop editor for OSM, to correct and smooth out discontinuities on the domain boundaries that could mar the quality of the ensuing mesh. Afterwards, the improved map has been imported into QGIS [67], a cross-platform feature-rich geographic information system application, from which a geometry input file for the mesh generator Gmsh [68], embedding a built-in computer-aided design (CAD) engine and postprocessor, is created using an interactive external plugin. This geometry file, which integrates a customizable parameter to control the element size of the meshes that will be derived from it, is then exported from QGIS according to the Universal Transverse Mercator (UTM) Coordinate Reference System, zone 29 N. Finally, a 2D computational mesh, made up of triangles, over the domain of interest is generated from the
<figure><img src="image_7.png" /><figcaption>FIGURE 4: Computational domain covering the Niger river in Bamako.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 5: A mesh of the computational domain.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIGURE 6: Computational domain and principal wastewater pipeline discharge points.</figcaption></figure>

<table><caption>TABLE 2: Geographic location and discharge of wastewater pipeline outfalls.</caption><thead><tr><td>Point</td><td>Latitude</td><td>Longitude</td><td>Neighborhood</td><td>Commune</td><td>Discharge (m³/day)</td></tr></thead><tbody><tr><td>C<sub>1</sub></td><td>12.5570651</td><td>-8.0453637</td><td>Kabalabougou</td><td>Commune IV</td><td>45,396</td></tr><tr><td>C<sub>2</sub></td><td>12.5890708</td><td>-8.0516760</td><td>Sébénikoro</td><td>Commune IV</td><td>9018</td></tr><tr><td>C<sub>3</sub></td><td>12.5986180</td><td>-8.0490785</td><td>Sébénikoro</td><td>Commune IV</td><td>665</td></tr><tr><td>C<sub>4</sub></td><td>12.6078214</td><td>-8.0358472</td><td>Sébénikoro</td><td>Commune IV</td><td>120,559</td></tr><tr><td>C<sub>5</sub></td><td>12.6164957</td><td>-8.0180427</td><td>Djicoroni</td><td>Commune IV</td><td>70,160</td></tr><tr><td>C<sub>6</sub></td><td>12.6255755</td><td>-8.0086040</td><td>Quartier du fleuve</td><td>Commune III</td><td>52,015</td></tr><tr><td>C<sub>7</sub></td><td>12.6339927</td><td>-7.9901530</td><td>Bozola</td><td>Commune II</td><td>38,927</td></tr><tr><td>C<sub>8</sub></td><td>12.6398227</td><td>-7.9637495</td><td>Zone Industrielle</td><td>Commune II</td><td>12,895</td></tr><tr><td>C<sub>9</sub></td><td>12.6426781</td><td>-7.9500890</td><td>Bougouba</td><td>Commune II</td><td>38,726</td></tr><tr><td>C<sub>10</sub></td><td>12.6421599</td><td>-7.9420692</td><td>Bougouba</td><td>Commune II</td><td>29,569</td></tr><tr><td>C<sub>11</sub></td><td>12.6461825</td><td>-7.9352805</td><td>Sotuba</td><td>Commune I</td><td>13,033</td></tr><tr><td>C<sub>12</sub></td><td>12.6045678</td><td>-8.0305552</td><td>Baco-Djicoroni</td><td>Commune V</td><td>1283</td></tr><tr><td>C<sub>13</sub></td><td>12.6082623</td><td>-8.0143689</td><td>Torokorobougou</td><td>Commune V</td><td>12,297</td></tr><tr><td>C<sub>14</sub></td><td>12.6140275</td><td>-8.0082070</td><td>Torokorobougou</td><td>Commune V</td><td>6436</td></tr><tr><td>C<sub>15</sub></td><td>12.6231479</td><td>-7.9987402</td><td>Badalabougou</td><td>Commune V</td><td>11,015</td></tr><tr><td>C<sub>16</sub></td><td>12.6247446</td><td>-7.9831995</td><td>Badalabougou</td><td>Commune V</td><td>16,661</td></tr><tr><td>C<sub>17</sub></td><td>12.6182011</td><td>-7.9721233</td><td>Magnambougou</td><td>Commune VI</td><td>61,358</td></tr><tr><td>C<sub>18</sub></td><td>12.6264485</td><td>-7.9575455</td><td>Magnambougou</td><td>Commune VI</td><td>14,256</td></tr><tr><td>C<sub>19</sub></td><td>12.6392653</td><td>-7.9207537</td><td>Missabougou</td><td>Commune VI</td><td>17,805</td></tr></tbody></table>

geometry file using Gmsh by specifying a target element size. In Figure 4 and Figure 5, the computational domain and an illustration related to a coarse mesh are presented, respectively.

4.2. *Numerical Simulation Setup.* The problem studied concerns the transport of pollutants by shallow-water flows in the section of the Niger River in Bamako, as shown in
<figure><img src="image_10.png" /><figcaption>FIGURE 7: Daily average concentration of arsenic (As) over November 2023 at the discharge point C<sub>9</sub>. Data source: DNACPN.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 8: Daily average concentration of cadmium (Cd) over November 2023 at the discharge point C<sub>9</sub>. Data source: DNACPN.</figcaption></figure>

Figure 4. This stretch of interest has an average depth of 2.97 m and a mean width of about 1.63 km, according to the data provided by the National Directorate of Hydraulics (DNH) of Mali for the period 2018–2023. In this region, the Niger River flows at an average discharge of 1091.7 m³/s, based on the Niger Basin Authority (NBA) data for the same period as above, from ∂Ω₁ (inlet boundary) toward ∂Ω₀ (outlet boundary) along the land boundaries (see Figure 4). This corresponds to an average annual flow of about 34.43 km³. The mean flow velocity in this section is estimated to 0.23 m/s.

We focus on the dynamics of some heavy metals that are of particular concern to human health and the environment

because of their high toxicity, even at low exposure levels, namely, As, Cd and Pb. These metals occur naturally in the aquatic environment, in particular the Niger River, with a range of natural background concentrations. However, they can become toxic above a critical load threshold, mainly as a result of anthropogenic releases.

The persistent accumulation of these metallic elements in the human body, resulting from prolonged exposure, can have fatal effects. They are associated with a number of serious illnesses, including impaired renal, hepatic, cognitive, and reproductive functions, teratogenic effects, allergic reactions, skin irritation, and cancer. Indeed, based on the International Agency for Research on Cancer (IARC)
<figure><img src="image_12.png" /><figcaption>FIGURE 9: Daily average concentration of lead (Pb) over November 2023 at the discharge point C<sub>9</sub>. Data source: DNACPN.</figcaption></figure>

monographs on the evaluation of carcinogenic risks to humans (Volume 100C) [69], As and Cd have been classified as carcinogenic to humans (Group 1), while Pb has been classified as possibly carcinogenic to humans (Group 2B). As a result, these metallic elements are subject to particularly close monitoring and are strictly regulated throughout the world by international and regional organizations and national governments [70]. We consider pollutants from the river discharge points of 19 major wastewater pipelines, represented in Figure 6. These pipelines drain up to 90% of the effluent from the wastewater collection system of Bamako [71]. This wastewater includes that from the few wastewater treatment plants of the city, which are connected to the main industrial zones, the main administrative and commercial centers, and the major hospitals. The geographic location and discharge of these pipeline outfalls are reported in Table 2.

In Table 2, the flow rate reported refers to the average daily discharge of wastewater released into the Niger River from the corresponding outfalls during the low-water period. These data were provided by the National Directorate for Sanitation and Pollution and Nuisance Control (DNACPN) of Mali. The communes represent the administrative subdivisions of Bamako, while the neighborhoods constitute the administrative entities of the mentioned communes. Of these outflow points, C<sub>8</sub>, C<sub>9</sub>, and C<sub>10</sub> are major outfall pipes discharging most of the heavy metals of interest in this paper, since they are connected to the principal industrial zone, home to more than 70% of Malian industry. The daily average concentrations (ACRs) of As, Cd, and Pb over November 2023 at the outfall C<sub>9</sub> and the standard limits as per World Health Organization (WHO) guidelines [72] are plotted in Figures 7, 8 and 9, respectively.

The wind shear stress on the water surface is assumed to be negligible, which holds that $\tau_s = 0$. The bed topography of

<table><caption>TABLE 3: Model parameter values for numerical simulations.</caption><thead><tr><td>Symbol</td><td>Value</td></tr></thead><tbody><tr><td>ν</td><td>0.8927e-6 (at 25°C)</td></tr><tr><td>g</td><td>9.81</td></tr><tr><td>n<sub>b</sub></td><td>0.033</td></tr><tr><td>ω</td><td>7.2921e-5</td></tr><tr><td>φ</td><td>0.22078</td></tr><tr><td>ρ<sub>a</sub></td><td>1.1839 (at 1.0 atm and 25°C)</td></tr><tr><td>ρ<sub>w</sub></td><td>997.05 (at 25°C)</td></tr></tbody></table>

the river is assumed to be flat, since it is relatively roughless with small variation in altitude, which implies that the bottom bathymetry $h_b$ is a constant and $\tau_b = 0$. The temperature of the river is supposed to be 25°C, the multiyear averaged temperature of tropical rivers [73]. The atmospheric temperature and pressure are assumed to be 25°C and 1.0 atm (101,325 Pa), respectively. The model parameter values considered for the numerical simulations are reported in Table 3. They are related to the mentioned water and atmospheric conditions. The values of the kinematic viscosity $\nu$ and density $\rho_w$ of water are based on reference data [74]. The air density is calculated according to a dedicated equation discussed in [75]. The Manning roughness coefficient $n_b$ is supposed to be uniformly distributed over the computational domain. Its typical value for natural streams ranges from 0.025 to 0.15, as established in [76]. Here, the value taken for this coefficient, as reported in the following table of parameters, refers to a constant widely used for river channels [77]. The value of the angular velocity $\omega$ of the Earth rotation is set following the numerical standards [78] of the International Earth Rotation and Reference Systems
<table><caption>TABLE 4: Finite element parameter values for numerical simulations.</caption><thead><tr><th>Section</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td rowspan="4">Mesh</td><td>Average resolution</td><td>200 m</td></tr><tr><td>Shape of elements</td><td>Triangles</td></tr><tr><td>Number of nodes</td><td>48,545</td></tr><tr><td>Number of elements</td><td>94,539</td></tr><tr><td rowspan="4">Discretization</td><td>Spatial approximation</td><td>P1 Lagrange</td></tr><tr><td>Temporal approximation</td><td>First-order BDF</td></tr><tr><td>Time step</td><td>300 s</td></tr><tr><td>Stabilization method</td><td>GLS</td></tr><tr><td rowspan="3">Time</td><td>Initial time</td><td>November 15, 2023, at 12:00 a.m.</td></tr><tr><td>Final time</td><td>November 16, 2023, at 12:00 a.m.</td></tr><tr><td>Time step</td><td>300 s</td></tr><tr><td rowspan="2">Constraints</td><td>Initial conditions</td><td>Data from NBA and DNACPN</td></tr><tr><td>Boundary conditions</td><td>Data from NBA and DNACPN</td></tr><tr><td rowspan="2">Solver</td><td>Nonlinear solver</td><td>Newton</td></tr><tr><td>Preconditioner</td><td>GASM</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>FIGURE 10: Water depth on November 15, 2023, at 05:00 a.m. and 11:45 p.m.</figcaption></figure>

Service (IERS). The value of $\phi$ is derived from the conversion of the mean latitude of the studied area, which is 12.6137 (in decimal degrees), into radians.

**4.3. Computational Platform.** The numerical simulations have been performed at the Center for Computational Modeling and Simulation (CCMS) of the Faculty of Science and Technology (FST) of Bamako. CCMS hosts a cluster of 10 heterogeneous compute nodes (servers) interconnected by a high-speed QDR InfiniBand network. It consists of five Dell PowerEdge R450 servers, each enabled with two Intel Xeon Silver 4110 processors with eight cores running at 2.10 GHz and 64 GB of RAM and five HPE ProLiant servers, each equipped with two Intel Xeon E5-2623v4 processors with four cores cadenced at 2.60 GHz and 16 GB of RAM. Of these compute nodes, two Dell PowerEdge servers are

featured with two NVIDIA Tesla V100 GPUs. Powered by Volta architecture, the NVIDIA Tesla V100 GPU provides 5120 CUDA computational cores over 84 streaming multi-processors (SMs) and 32 GB onboard memory. CCMS delivers up to 7 and 14 TFLOPS of peak floating-point performance for single and double precision, respectively.

**4.4. Numerical Results and Discussion.** The numerical results presented here are based on a 1-day (model time) simulation, precisely over November 15, 2023, achieved using the CFPDE toolbox, featured by Feel++. The computational mesh has a spatial resolution estimated to 200 m and contains 48,545 nodes and 94,539 unstructured triangular elements. The spatial integration of the model is based on a piecewise linear and continuous finite element approximation, known as P1. Indeed, this discretization corresponds to
<figure><img src="image_14.png" /><figcaption>FIGURE 11: Flow discharge magnitude on November 15, 2023, at 05:00 a.m. and 11:45 p.m.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 12: Concentration of arsenic (As) on November 15, 2023, at 05:00 a.m. and 11:45 p.m.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 13: Concentration of cadmium (Cd) on November 15, 2023, at 05:00 a.m. and 11:45 p.m.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIGURE 14: Concentration of lead (Pb) on November 15, 2023 at 05:00 a.m. and 11:45 p.m.</figcaption></figure>

the finite element pair P1-P1 for the discharge and water depth, solutions of shallow-water equations. It is among the pairs widely discussed in the literature and leading to promising results for this system [79]. The temporal discretization is carried out using a first-order BDF (backward differentiation formula) scheme [80]. The model time-step is set to 300 s. From the simulation toolbox, a nonlinear Newton-based solver and the generalized additive Schwarz method (GASM), both wrapped from PETSc [81], are used for the nonlinear iterations of the numerical approximation of the model and as preconditioners for the solution of the linear systems involved in these iterations, respectively. Given the hyperbolic nature of the system, the Galerkin least-squares (GLS) [63] support from the CFPDE toolbox is enabled for ensuring the stability of the finite element solution. This stabilization technique relies on the introduction of additional least-squares forms of residuals to the standard Galerkin variational formulations for achieving more stable and accurate numerical solutions, especially for such hyperbolic problems, often involving discontinuous solutions. Initial and boundary conditions for the shallow-water hydrodynamic model, relating to the discharge and water depth, are derived from data from the NBA. Initial and boundary conditions for the advection-diffusion transport model are obtained from DNACPN. The inflow boundary conditions concern the pollutant mass fluxes from the upstream boundary of the river and wastewater discharge points into the river. On the other hand, the outflow boundary conditions refer to the prescribed water depth at the downstream boundary. In practice, these source data, preprocessed and stored in Hierarchical Data Format 5 (HDF5), are interpolated in accordance with the finite element approximation and simulation configuration parameters on the computational mesh. The essential finite element parameters are summarized in Table 4.

The water depth, flow discharge magnitude, and the concentrations of As, Cd, and Pb obtained from simulations are plotted in Figures 10, 11, 12, 13, and 14, respectively, for the model time of 05:00 a.m. and 11:45 p.m. of the stated date.

<table><caption>TABLE 5: Computational costs of simulations.</caption><thead><tr><td>Number of CPU cores</td><td>CPU time (in hours)</td></tr></thead><tbody><tr><td>16</td><td>35.81</td></tr><tr><td>32</td><td>21.53</td></tr><tr><td>64</td><td>13.47</td></tr></tbody></table>

They provide insight into the river hydrodynamics and the dynamics of heavy metals of interest between the two reported time steps. They clearly demonstrate that pollutants are carried downstream by water flow. These substances spread as quickly as the flow rate increases, illustrating the clear dependence of this transport process on the hydrodynamic characteristics of the river flow, as expected by the coupled model proposed here. It can be seen that in high-flow areas, pollutants are more rapidly diluted and disperse less, while in low-flow areas, they tend to stabilize, resulting in higher concentrations and greater dispersion. It can also be observed that the highest concentrations are located close to the points of discharge and decrease progressively as the flow continues and the pollutants move away from these sources. These behaviors are consistent with physical observations and are justified by well-established physical principles of advection and diffusion, essential components of the passive transport model considered in this framework and river hydrodynamics.

In terms of real-world applications, this framework is very useful for predicting the pollution patterns of the Niger River resulting from an accidental spill or industrial releases and providing a precise estimate of when and where the pollution level will reach a critical threshold. This enables authorities and communities to make informed decisions and to better plan and optimize intervention strategies for mitigating the impacts of such disasters on public health, human activities, and ecosystems. It can also be used to assess the effectiveness of these interventions and identify measures for possible readjustment.
<figure><img src="image_18.png" /><figcaption>FIGURE 15: Average concentration of arsenic (As) over the computational domain on November 15, 2023.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIGURE 16: Average concentration of cadmium (Cd) over the computational domain on November 15, 2023.</figcaption></figure>

In terms of computational costs, the simulation presented above has been achieved in 13.75 h (CPU time) using 64 processor cores of the CCMS cluster. For further illustrating the scalability and efficiency of the implementation, simulations based on the same configuration parameters as above have been performed using different numbers of processor cores. Estimates of CPU time consumed as a function of the number of CPU cores employed are reported in Table 5. These clearly show that this framework scales well as the number of computing resources used for simulations increases.

According to the same simulation setups, the ACRs of As, Cd, and Pb over the computational domain on November 15,

2023, as a function of time are illustrated in Figures 15, 16, and 17, respectively. They display ample evidence of the extent of the pollution of the Niger River in Bamako, since the pollutant loads indicated far exceed WHO standard limits. This confirms our intuitive hunch about the situation of this river. However, in order to establish the validity of these numerical results in concrete terms, a comparative analysis with observational data is carried out in Section 4.5.

It is well established in the literature that mesh convergence is a fundamental concept of finite element analysis [55]. It refers to the behavior of the numerical solution, in terms of accuracy, as the computational mesh becomes finer. This process provides a metric to ensure that, as the mesh is
<figure><img src="image_20.png" /><figcaption>FIGURE 17: Average concentration of lead (Pb) over the computational domain on November 15, 2023.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 18: Mesh convergence study based on relative L2-norm differences between concentrations from successive mesh refinements.</figcaption></figure>

progressively refined, the numerical solution achieves a better approximation that no longer depends on the mesh resolution. Clearly, effective mesh refinement results in more computational effort. In order to analyze the convergence of a numerical model and eventually determine an optimal level of refinement that guarantees a suitable balance between numerical accuracy and computational cost, a mesh convergence study is usually conducted.

In this context, to ensure that the numerical results presented here satisfy the convergence properties as expected by

the FEM, we perform numerical simulations based on the same configuration parameters as above using a series of three computational meshes resulting from successive and consistent refinements, where the resolution is halved at each iteration. It follows that halving the mesh resolution does not necessarily mean doubling the number of elements. These meshes, from the coarsest to the finest, have resolutions of 400, 200, and 100 m, respectively. Given the absence of an exact solution, which cannot be expressed in this framework, and of sufficient observation data throughout the computational domain, we refer to the relative L2-norm differences between pairs of concentrations of each of the heavy metals of interest, including As, Cd, and Pb, obtained from simulations based on successive levels of mesh refinement at specific points in the domain. These are plotted as a function of mesh resolution in Figure 18. These plots clearly show that the difference between numerical solutions from successive mesh refinements decreases as the mesh is refined, as expected in finite element analysis. Indeed, it would have been very useful to extend this study to even finer meshes and confirm the clear evidence of mesh convergence that emerges. However, we are constrained by computational resources to be able to achieve such simulations within a reasonable time.

## 4.5. *Numerical Verification and Validation.*
This important step pursues to assess the validity and reliability of the numerical results produced by the simulation model presented above. For this purpose, we perform a comparative analysis of the numerical results with observational data at a gauge point, marked A, targeted downstream of the computational domain (see Figure 6). This approach is motivated by the fact that a sufficient amount of pollution-related data is available at this point, which hosts a reliable IoT-based water quality monitoring system, consisting of multiple smart sensors, deployed in the water body, a
<figure><img src="image_22.png" /><figcaption>FIGURE 19: Concentrations of arsenic (As) from simulation and observation at the Gauge Point A on November 15, 2023.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>FIGURE 20: Simulated concentration of arsenic (As) as a function of observed concentration and linear regression error metrics.</figcaption></figure>

microcontroller, interfaced with these sensors, and connectivity devices. The sensors measure essential water quality parameters, including common physicochemical and microbiological as well as most hazardous heavy metals. The sensor data are collected by the microcontroller at half-hourly intervals for initial processing and recording. This collection frequency represents the highest configurable for this system. Later, the processed data are transmitted, on a six-hourly basis, from the microcontroller to a central server

through a high-speed Wi-Fi module for storage and further processing. This central controller is hosted at DNACPN, the operator of the system. Typically, such a frequency of data collection, also known as temporal resolution, is not sufficient for critical applications requiring real-time decision-making for water pollution monitoring, since it is unable to capture shorter-term fluctuations in pollution patterns that can be subject to uncertainties affecting the accuracy and reliability of the data provided. However, it is
suited for the purposes of this validation, since the data processed provide sufficient insights for the comparative analysis with numerical results.

Based on these collected data, the concentrations of As from numerical simulation (ASM) and As from observation (ASO) as well as the WHO permissible limit are displayed in Figure 19. These plots clearly show the similarity of trends between simulation results and observed data, supporting the numerical experiments achieved and the proposed framework. To further support these hopeful qualitative results, we have introduced quantitative measures of errors between numerical results and observational data on the basis of RMSE and MAE metrics, some of the most commonly used error measures in linear regression. The scatter (SC) plot of simulated concentration of As as a function of observed concentration, the linear regression (LR) line, and error metrics are displayed in Figure 20. These metrics, both estimated below 5% of the observed mean value, support the reliability and accuracy of this framework in capturing pollutant dynamics under flow conditions and highlight its predictive effectiveness.

## 5. Conclusion

We presented a unified numerical framework for the transport of passive pollutants by shallow-water flows. A mathematical model governing this phenomenon and consisting of the coupling of the hydrodynamic shallow-water equations with an advection–diffusion equation for pollutant transport has been described in detail, including the underlying assumptions, the classification of flow regimes, and the initial and boundary conditions. The numerical implementation of the resulting hyperbolic coupled model, including the time discretization and spatial integration, the computational aspects and strategies, and the essential software tools, has been discussed. Substantial effort has gone into numerical experiments centered on As, Cd, and Pb, some of the most toxic heavy metals and the most scrutinized in water quality monitoring and assessment. They have been completed as part of a nontrivial practical application on the Niger River in Bamako, where this resource is particularly vital for the population and constitutes the principal source of drinking water. Particular emphasis was placed on the design of the computational domain, a complex-shaped geometry based on a map of the Niger River in Bamako formatted from OpenStreetMap, and the generation of the ensuing simulation meshes using powerful APIs and computer software. The numerical results from simulations, achieved using realistic configurations, were illustrated and discussed. The numerical validation of these results with reference to observational data has been assessed on the basis of RMSE and MAE metrics, some of the most commonly used error measures in linear regression. These achievements, including numerical results and their validity, have clearly demonstrated the suitability of the proposed framework to reproduce the pollutant dynamics and to provide a better insight into the pollution patterns of the Niger River in Bamako.

In the short term, this work will be followed by numerical simulations taking into account the bottom bathymetry of the Niger River and the bottom shear and wind shear stresses as soon as the data required for this purpose are available. Future developments will integrate a module for physical and chemical processes for addressing reactive pollutants. Finally, this framework will be extended to advanced artificial intelligence (AI) algorithms for enhancing its predictive capability in the near future.

## Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## Conflicts of Interest

The authors declare no conflicts of interest.

## Author Contributions

**Abdoulaye Samaké**: conceptualization, investigation, methodology, software, validation, visualization, writing – original draft. **Mahamadou Alassane**: Methodology; validation; writing – review and editing. **Amadou Mahamane**: methodology, validation, writing – review and editing.

## Funding

No funding was received for this research.

## Acknowledgments

The authors would like to thank Christophe Prud'homme and Vincent Chabannes for many fruitful discussions. They acknowledge the support in terms of computing resources from the Center for Computational Modeling and Simulation (CCMS) of the Faculty of Science and Technology (FST) of Bamako.

## References

[1] P. Goel, *Water Pollution: Causes, Effects and Control* (New Age International, 2006).

[2] R. P. Schwarzenbach, T. Egli, T. B. Hofstetter, U. Von Gunten, and B. Wehrli, "Global Water Pollution and Human Health, Annual Review of Environment and Resources 35, no. 1 (2010): 109–136, https://doi.org/10.1146/annurev-environ-100809-125342.

[3] L. Schweitzer and J. Noblet, "Water Contamination and Pollution," in *Green Chemistry* (Elsevier, 2018), 261–290, https://doi.org/10.1016/B978-0-12-809270-5.00011-X.

[4] D. Walker, D. J. Baumgartner, C. P. Gerba, and K. Fitzsimmons, "Surface Water Pollution," in *Environmental and Pollution Science* (Elsevier, 2019), 261–292, https://doi.org/10.1016/B978-0-12-814719-1.00016-1.

[5] S. Kurwadkar, S. R. Kanel, and A. Nakarmi, "Groundwater Pollution: Occurrence, Detection, and Remediation of Organic and Inorganic Pollutants," *Water Environment Research* 92, no. 10 (2020): 1659–1668, https://doi.org/10.1002/wer.1415.
[6] A. Prüss-Ustün, J. Wolf, J. Bartram, et al., "Burden of Disease From Inadequate Water, Sanitation and Hygiene for Selected Adverse Health Outcomes: An Updated Analysis With a Focus on Low-and Middle-Income Countries," *International Journal of Hygiene and Environmental Health* 222, no. 5 (2019): 765-777, https://doi.org/10.1016/j.ijheh.2019.05.004.

[7] D. Vilcins, P. D. Sly, and P. Jagals, "Environmental Risk Factors Associated With Child Stunting: A Systematic Review of the Literature," *Annals of Global Health* 84, no. 4 (2018): 551, https://doi.org/10.29024/aogh.2361.

[8] R. Fuller, P. J. Landrigan, K. Balakrishnan, et al., "Pollution and Health: A Progress Update," *Lancet Planetary Health* 6, no. 6 (2022): e535-e547, https://doi.org/10.1016/S2542-5196(22)00090-0.

[9] UNESCO, *The United Nations World Water Development Report 2021: Valuing Water* (United Nations, 2021).

[10] M. A. Abedin, A. E. Collins, U. Habiba, and R. Shaw, "Climate Change, Water Scarcity, and Health Adaptation in Southwestern Coastal Bangladesh," *International Journal of Disaster Risk Science* 10, no. 1 (2019): 28-42, https://doi.org/10.1007/s13753-018-0211-8.

[11] A. O. Fayiga, M. O. Ipinmoroti, and T. Chirenje, "Environmental Pollution in Africa," *Environment, Development and Sustainability* 20, no. 1 (2018): 41-73, https://doi.org/10.1007/s10668-016-9894-4.

[12] I. Adam, T. R. Walker, J. C. Bezerra, and A. Clayton, "Policies to Reduce Single-Use Plastic Marine Pollution in West Africa," *Marine Policy* 116 (2020): 103928, https://doi.org/10.1016/j.marpol.2020.103928.

[13] A. M. Freeman, "Water Pollution Policy," in *Public Policies for Environmental Protection* (Routledge, 2016), 97-149.

[14] U. Ahmed, R. Mumtaz, H. Anwar, S. Mumtaz, and A. M. Qamar, "Water Quality Monitoring: From Conventional to Emerging Technologies," *Water Supply* 20, no. 1 (2020): 28-45, https://doi.org/10.2166/ws.2019.144.

[15] K. Acharya, A. Blackburn, J. Mohammed, A. T. Haile, A. M. Hiruy, and D. Werner, "Metagenomic Water Quality Monitoring With a Portable Laboratory," *Water Research* 184 (2020): 116112, https://doi.org/10.1016/j.watres.2020.116112.

[16] X. Wang and W. Yang, "Water Quality Monitoring and Evaluation Using Remote Sensing Techniques in China: A Systematic Review," *Ecosystem Health and Sustainability* 5, no. 1 (2019): 47-56, https://doi.org/10.1080/20964129.2019.1571443.

[17] L. Agnoli, E. Urquhart, N. Georgantzis, et al., "Perspectives on User Engagement of Satellite Earth Observation for Water Quality Management," *Technological Forecasting and Social Change* 189 (2023): 122357, https://doi.org/10.1016/j.techfore.2023.122357.

[18] C. Alexandra, K. A. Daniell, J. Guillaume, C. Saraswat, and H. R. Feldman, "Cyber-Physical Systems in Water Management and Governance," *Current Opinion in Environmental Sustainability* 62 (2023): 101290, https://doi.org/10.1016/j.cosust.2023.101290.

[19] J. O. Ighalo, A. G. Adeniyi, and G. Marques, "Internet of Things for Water Quality Monitoring and Assessment: A Comprehensive Review," in *Artificial Intelligence for Sustainable Development: Theory, Practice and Future Applications* (Springer, 2021), 245-259, https://doi.org/10.1007/978-3-030-51920-9_13.

[20] K. Murphy, B. Heery, T. Sullivan, et al., "A Low-Cost Autonomous Optical Sensor for Water Quality Monitoring," *Talanta* 132 (2015): 520-527, https://doi.org/10.1016/j.talanta.2014.09.045.

[21] H. Wang, T. Wang, B. Zhang, et al., "Water and Wastewater Treatment in Africa-Current Practices and Challenges," *CLEAN-Soil, Air, Water* 42, no. 8 (2014): 1029-1035, https://doi.org/10.1002/clen.201300208.

[22] J. A. Sokolowski and C. M. Banks, *Modeling and Simulation Fundamentals: Theoretical Underpinnings and Practical Domains* (John Wiley & Sons, 2010), https://doi.org/10.1002/9780470590621.

[23] K. Karroum, Y. Lin, Y.-Y. Chiang, et al., "A Review of Air Quality Modeling," *Mapan* 35, no. 2 (2020): 287-300, https://doi.org/10.1007/s12647-020-00371-8.

[24] V. K. Singh, D. Gao, and A. Fischer, *Advances in Mathematical Methods and High Performance Computing* (Springer, 2019), https://doi.org/10.1007/978-3-030-02487-1.

[25] R. Welcomme and H. L. Dumont, *The Niger River System, the Ecology of River Systems* (Springer, 1986), https://doi.org/10.1007/978-94-017-3290-1_2.

[26] T. Abrate, P. Hubert, and D. Sighomnou, "A Study on Hydrological Series of the Niger River," *Hydrological Sciences Journal* 58, no. 2 (2013): 271-279, https://doi.org/10.1080/02626667.2012.752575.

[27] I. Andersen and K. G. Golitzen, *The Niger River Basin: A Vision for Sustainable Management* (World Bank Publications, 2005).

[28] J. Yang, Y. C. E. Yang, H. F. Khan, et al., "Quantifying the Sustainability of Water Availability for the Water-Food-Ecosystem Nexus in the Niger River Basin," *Earth's Future* 6, no. 9 (2018): 1292-1310, https://doi.org/10.1029/2018EF000923.

[29] A. O. Abdoulaye, H. Lu, Y. Zhu, and Y. A. Hamoud, "Future Irrigation Water Requirements of the Main Crops Cultivated in the Niger River Basin," *Atmosphere* 12, no. 4 (2021): 439, https://doi.org/10.3390/atmos12040439.

[30] O. Seidou, C. Ringler, S. Kalcic, L. Ferrini, T. A. Ramani, and A. Guero, "A Semi-Qualitative Approach to the Operationalization of the Food-Environment-Energy-Water (FE²W) Nexus Concept for Infrastructure Planning: A Case Study of the Niger Basin," *Water International* 46, no. 5 (2021): 744-770, https://doi.org/10.1080/02508060.2021.1956231.

[31] S. Ba, E. U. Onyeabor, and A. N. Moneke, "The Current Legal Framework for Pollution Control in the Niger River Basin Relative to SDG 6.3," *Water International* 47, no. 8 (2022): 1217-1234, https://doi.org/10.1080/02508060.2022.2073756.

[32] C. Béné, L. Evans, D. Mills, et al., "Testing Resilience Thinking in a Poverty Context: Experience From the Niger River Basin," *Global Environmental Change* 21, no. 4 (2011): 1173-1184, https://doi.org/10.1016/j.gloenvcha.2011.07.002.

[33] A. Ogilvie, G. Mahé, J. Ward, et al., "Water, Agriculture and Poverty in the Niger River Basin," *Water International* 35, no. 5 (2010): 594-622, https://doi.org/10.1080/02508060.2010.515545.

[34] R. van der Wijngaart, J. Helming, C. Jacobs, P. G. Delvaux, S. Hoek, and S. Gomez y Paloma, *Irrigation and Irrigated Agriculture Potential in the Sahel: The Case of the Niger River Basin* (Publications Office of the European Union, 2019).

[35] L. O. Sangaré, S. Ba, O. Diallo, D. Sanogo, and T. Zheng, "Assessment of Potential Health Risks From Heavy Metal Pollution of Surface Water for Drinking in a Multi-Industry Area in Mali Using a Multi-Indices Approach," *Environmental*
*Monitoring and Assessment* 195, no. 6 (2023): 700, https://doi.org/10.1007/s10661-023-11258-x.

[36] A. du Bassin du Fleuve Niger, *Rapport sur L'état du Fleuve Niger au Mali* (Rapport Final, 2018).

[37] O. Abayomi, T. E. Olayemi, and T. Ogungbade, "Environmen-
tal Pollution and Its Ecological Consequences on the Niger
Delta: A Review of the Literature," *African Journal of Environ-
ment and Natural Science Research* 4 (2021): 27-42.

[38] D. Orange, A. Palangié, Y. Xu, and B. Usher, "Assessment of Water Pollution and Risks to Surface and Groundwater Resources in Bamako, Mali," in *Groundwater pollution in Africa* (CRC Press, 2006), 8, https://www.taylorfrancis.com/chapters/edit/10.1201/9780203963548-24/assessment-water-pollution-risks-surface-groundwater-resources-bamako-mali-orange-palangi%C3%A9.

[39] J. Grijsen, C. Brown, A. Tarhule, Y. Ghile, and M. Taner, *Climate Risk Assessment for Water Resources Development in the Niger River Basin* 12, Water Partnership Program (WPP)/TWIWA, HydroPredict, 2013).

[40] Z.-G. Ji, *Hydrodynamics and Water Quality: Modeling Rivers, Lakes, and Estuaries* (John Wiley & Sons, 2017), https://doi.org/10.1002/9781119371946.

[41] D. Bresch, "Chapter 1 - Shallow-Water Equations and Related Topics," in *Handbook of Differential Equations: Evolutionary Equations* 5, Elsevier, 2009), 1-104, https://doi.org/10.1016/S1874-5717(08)00208-9.

[42] P. García-Navarro, J. Murillo, J. Fernández-Pato, I. Echeverribar, and M. Morales-Hernández, "The Shallow Water Equations and Their Application to Realistic Cases," *Environmental Fluid Mechanics* 19, no. 5 (2019): 1235-1252, https://doi.org/10.1007/s10652-018-09657-7.

[43] J. F. C. Meyer and G. L. Diniz, "Pollutant Dispersion in Wetland Systems: Mathematical Modelling and Numerical Simulation," *Ecological Modelling* 200, no. 3-4 (2007): 360-370, https://doi.org/10.1016/j.ecolmodel.2006.08.005.

[44] S. Li and C. J. Duffy, "Fully-Coupled Modeling of Shallow Water Flow and Pollutant Transport on Unstructured Grids," *Procedia Environmental Sciences* 13 (2012): 2098-2121, https://doi.org/10.1016/j.proenv.2012.01.200.

[45] D. Vanzo, A. Siviglia, and E. F. Toro, "Pollutant Transport by Shallow Water Equations on Unstructured Meshes: Hyperbolization of the Model and Numerical Solution via a Novel Flux Splitting Scheme," *Journal of Computational Physics* 321 (2016): 1-20, https://doi.org/10.1016/j.jcp.2016.05.023.

[46] G. J. Arcement and V. R. Schneider, *Guide for Selecting Man-
ning's Roughness Coefficients for Natural Channels and Flood
Plains* (Water Supply Paper 2339, US Geological Survey,
1989).

[47] Y. M. Choo, G. S. Yun, and T. H. Choo, "A Research on the Estimation of Coefficient Roughness in Open Channel Applying Entropy Concept," *Environmental Earth Sciences* 77, no. 17 (2018): 624, https://doi.org/10.1007/s12665-018-7809-4.

[48] F. Benkhaldoun, S. Daoudi, I. Elmahi, and M. Seaid, "Numerical Modelling of Sediment Transport in the Nador Lagoon (Morocco)," *Applied Numerical Mathematics* 62, no. 12 (2012): 1749-1766, https://doi.org/10.1016/j.apnum.2012.05.010.

[49] C. B. Vreugdenhil, *Numerical Methods for Shallow-Water Flow* 13, Springer Science & Business Media, 2013).

[50] P. K. Kundu, I. M. Cohen, and D. R. Dowling, *Fluid Mechanics* (Academic press, 2015).

[51] M. Morales-Hernández, J. Murillo, and P. García-Navarro, "Diffusion-Dispersion Numerical Discretization for Solute Transport in 2D Transient Shallow Flows," *Environmental Fluid Mechanics* 19, no. 5 (2019): 1217-1234, https://doi.org/10.1007/s10652-018-9644-2.

[52] F. Yang and D. Liang, "Random-Walk Simulation of Non-Conservative Pollutant Transport in Shallow Water Flows," *Environmental Modelling & Software* 134 (2020): 104870, https://doi.org/10.1016/j.envsoft.2020.104870.

[53] P. Öffner, *Approximation and Stability Properties of Numerical Methods for Hyperbolic Conservation Laws* (Springer Nature, 2023), https://doi.org/10.1007/978-3-658-42620-0.

[54] A. Kurganov, "Finite-Volume Schemes for Shallow-Water Equations," *Acta Numerica* 27 (2018): 289-351, https://doi.org/10.1017/S0962492918000028.

[55] A. Ern and J. Guermond, “Theory and Practice of Finite Elements,” in *Applied Mathematical Sciences* 159, Springer, 2004), https://doi.org/10.1007/978-1-4757-4355-5.

[56] E. Hanert, D. Y. L. Roux, V. Legat, and E. Deleersnijder, "An Efficient Eulerian Finite Element Method for the Shallow Water Equations," *Ocean Modelling* 10, no. 1-2 (2005): 115-136, https://doi.org/10.1016/j.ocemod.2004.06.006.

[57] C. Prud'homme, V. Chabannes, V. Doyeux, M. Ismail, A. Samake, and G. Pena, "Feel++: A Computational Framework for Galerkin Methods and Advanced Numerical Methods," *ESAIM: Proceedings* 38 (2012): 429-455, https://doi.org/10.1051/proc/201238024.

[58] A. Logg, K.-A. Mardal, and G. Wells, *Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book* 84, Springer Science & Business Media, 2012), https://doi.org/10.1007/978-3-642-23099-8.

[59] F. Hecht, "New Development in FreeFem++," *Journal of Numerical Mathematics* 20, no. 3-4 (2012): 251-266.

[60] D. Arndt, W. Bangerth, B. Blais, et al., "The Deal. ii Library, Version 9.2," *Journal of Numerical Mathematics* 28, no. 3 (2020): 131-146, https://doi.org/10.1515/jnma-2020-0043.

[61] Y. Renard and K. Poulios, "GetFEM: Automated FE Modeling of Multiphysics Problems Based on a Generic Weak Form Language," *ACM Transactions on Mathematical Software (TOMS)* 47, no. 1 (2021): 1-31, https://doi.org/10.1145/3412849.

[62] R. Anderson, J. Andrej, A. Barker, et al., "Mfem: A Modular Finite Element Methods Library," *Computers & Mathematics with Applications* 81 (2021): 42-74, https://doi.org/10.1016/j.camwa.2020.06.009.

[63] C. Prud'homme, V. Chabannes, T. Saigre, et al., *feelpp/feelpp: Feel++ Release V111 Preview.10* (v0.111.0-review.10) (Zenodo, 2024), https://doi.org/10.5281/zenodo.12742155.

[64] C. Multiphysics, *Introduction to Comsol Multiphysics*® (COMSOL Multiphysics, 1998), accessed Feb 9.

[65] P. Mooney, "An outlook for openstreetmap," in *OpenStreetMap in GIScience: Experiences, Research, and Applications*, eds. J. Jokar Arsanjani, A. Zipf, P. Mooney, and M. Helbich (Springer International Publishing, 2015), 319-324, https://doi.org/10.1007/978-3-319-14280-7_16.

[66] R. M. Olbricht, "Data Retrieval for Small Spatial Regions in Openstreetmap," in *OpenStreetMap in GIScience: Experiences, Research, and Applications* (Springer, 2015), 101-122, https://doi.org/10.1007/978-3-319-14280-7_6.

[67] QGIS Development Team, *QGIS Geographic Information System* (QGIS Association, 2022).
[68] C. Geuzaine and J.-F. Remacle, “Gmsh: A 3-d Finite Element Mesh Generator With Built-in Pre-and Post-Processing Facilities,” *International Journal for Numerical Methods in Engineering* 79, no. 11 (2009): 1309–1331, https://doi.org/10.1002/nme.2579.

[69] International Agency for Research on Cancer, “Arsenic, Metals, Fibres, and Dusts. Volume 100c. A Review of Human Carcinogens,” in *IARC Monographs on the Evaluation of Carcinogenic Risks to Humans* (International Agency for Research on Cancer, 2012), 41–94, Arsenic and arsenic compounds.

[70] United Nations Food and Agriculture Organization, *Heavy Metals Regulations Legal Notice no 66/2003* (Food and Agriculture Organization (FAO), 2003).

[71] S. Ba, B. Dicko, and M. Sidibé, *Cartographie du Réseau D'égout de Bamako et Évaluation des Déversements des Eaux Usées de la Ville Dans le Fleuve Niger* (ÉTUDE CART'EAU: Rapport Technique, 2020), https://www.pseau.org/outils/biblio/resume.php?d=9275.

[72] World Health Organization, *Guidelines for Drinking-Water Quality: Incorporating the First and Second Addenda* (World Health Organization, 2022).

[73] S. Liu, Z. Xie, B. Liu, et al., “Global River Water Warming Due to Climate Change and Anthropogenic Heat Emission,” *Global and Planetary Change* 193 (2020): 103289, https://doi.org/10.1016/j.gloplacha.2020.103289.

[74] M. L. Huber, R. A. Perkins, A. Laesecke, et al., “New International Formulation for the Viscosity of H2O,” *Journal of Physical and Chemical Reference Data* 38, no. 2 (2009): 101–125, https://doi.org/10.1063/1.3088050.

[75] S. Rehman and N. M. Al-Abbadi, “Wind Shear Coefficient, Turbulence Intensity and Wind Power Potential Assessment for Dhulom, Saudi Arabia,” *Renewable Energy* 33, no. 12 (2008): 2653–2660, https://doi.org/10.1016/j.renene.2008.02.012.

[76] [V. Te Chow], *Open-Channel Hydraulics*, Civil Engineering Series (McGraw-Hill, 1959), https://books.google.ml/books?id=JG9_PwAACAAJ.

[77] J.-S. Kim, C.-J. Lee, W. Kim, and Y.-J. Kim, “Roughness Coefficient and its Uncertainty in Gravel-Bed River,” *Water Science and Engineering* 3, no. 2 (2010): 217–232, https://doi.org/10.3882/j.issn.1674-2370.2010.02.010.

[78] G. Petit and B. Luzum, *IERS Conventions (2010)* (Germany: Bundesamt für Kartographie und Geodäsie Frankfurt am Main, 2010).

[79] D. Y. Le Roux, A. Staniforth, and C. A. Lin, “Finite Elements for Shallow-Water Equation Ocean Models,” *Monthly Weather Review* 126, no. 7 (1998): 1931–1951, https://doi.org/10.1175/1520-0493(1998)126%3C1931:FEFSWE%3E2.0.CO;2.

[80] U. M. Ascher and L. R. Petzold, *Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations* (SIAM, 1998), https://doi.org/10.1137/1.9781611971392.

[81] S. Balay, S. Abhyankar, M. F. Adams, et al., *Petsc Users Manual* (Argonne National Laboratory, 2019), https://publications.anl.gov/anlpubs/2019/12/155920.pdf.