Telemac 2D Modeling of Pollutant Transport in the Medjerda River: Impact of Tamarix Vegetation, Bousselem, Tunisia

The aim of this study is to evaluate the impact of pollutant discharges from Bousselem WWTP in the Medjerda River (Tunisia). We chose to model the transport of a passive pollutant from the Bousselem WWTP using Telemac 2D model. The hydrodynamic model calibration was based on identifying the roughness coefficient derived from existing vegetation, mainly from the Tamarix type. Three scenarios of hydrocarbon dispersion representing potential accidents, with concentrations of 0.25, 10, and 50 kg/m³, were simulated. Three successive downstream sections were chosen to track concentration variability in space and over time. The final concentrations recorded (6 km from discharge point) are 500 μg/l, 12 μg/l, and 2.4 μg/l for scenarios of 50 kg/m3, 1 kg/m³, and 0.25 kg/m³, respectively. These values significantly exceed the international norm for drinking water, set at 1 μg/l, and so causing problems for drinking use. The obtained results serve to ensure environmental safety for hydraulic resources.

Introduction

The expansion of human activities has played a major role in the emergence of river water pollution problems [1]. The population explosion, industrial activities, and the rapid development and use of new chemicals and products pose a global environmental threat. This pollution contributes to water quality deterioration, resulting in a scarcity of water for human use, especially in semi-arid and arid regions vulnerable to anthropogenic pollution from urban, industrial, and agricultural sources [2]. Fluctuations in the physico-chemical quality of surface water in rivers are crucial due to the presence of point-source injections of various pollutants, impacting environmental biodiversity [3].

The Medjerda River is the sole perennial river in Tunisia [4]. originates near Souk-Ahras in Algeria, flowing eastward for 485 km before discharging into the Mediterranean Sea. The Medjerda collects the majority of surface water from northern Tunisia, transporting around one million cubic meters of water annually, representing half of the exploitable surface water resources in the country. However, the Medjerda faces substantial environmental challenges, notably pollution. In the Medjerda catchment area, it is estimated that untreated municipalities and rural sectors are responsible for discharging 13.5 million m³/year of untreated household wastewater into the natural environment. The Oued Medjerda receives urban wastewater discharges estimated to 1.27 million m³/year and treated water discharges from the 19 ONAS treatment plants estimated to 12 million m³/year [5]. Non-compliance of some of these discharges increases the risk of water resource contamination.

Given these concerning data, it becomes imperative to establish alert systems to monitor pollutant propagation. Modeling proves crucial in understanding the dispersion of contaminants in water bodies and identifying risk zones where concentrations exceed established norms. The increasingly frequent occurrence of point-source pollution in rivers has prompted the development of numerical simulation tools that enable to predict accurately the space and time variations of velocity fields and concentrations [6]. These tools provide more complete information on the spatio-temporal development of pollutant load concentrations [7], allowing precise forecasts that can aid in decision-making and developing action plans for extreme cases.

To assess the impact of pollutant discharges on water quality in terms of quantitative and physico-chemical characteristics, we chose to model the transport of a passive pollutant along a section of the Medjerda, starting from the Bousselem wastewater treatment plant using the Telemac 2D model. Although the evolution of pollutants in an aquatic environment is intrinsically linked to the characteristics of the transporting flow [8], and so constructing the model will allow us to understand and anticipate these dynamics in the flow. This approach also enables us to evaluate the extent of environmental impacts, providing an essential basis for developing management and mitigation strategies. By simulating various emission scenarios, modeling becomes a valuable tool for anticipating and assessing the consequences of these discharges, contributing to decision-making in water resource preservation [9].

On the other hand, the Medjerda is recognized for the excessive development of Tamarix vegetation, which influences the dispersion and retention of pollutants, acting as a natural filter by retaining contaminants through adsorption and particle trapping. But these Tamarix that grow fastly, amplify the risk of flooding by narrowing the cross-sections, and they have often been cleaned up by the government. In context, we will consider a scenario assuming the total absence of vegetation. By comparing this result to the existing state, we can deduce the impact of vegetation presence on the distribution of pollutants at the section level.

Methodology

Given that our hydrodynamic model developed for the Bousselem part of the Medjerda has been calibrated and yielded positive results in our studies, particularly in “Analysis and forecasting flood risk mapping of the Medjerda river at Bousselem town” [10] and “1D/2D coupling model to assess the impact of dredging works on the Medjerda river floods, Tunisia,” [11] we have extended this hydrodynamic model to simulate pollutants in the Medjerda from the wastewater treatment station in Bousselem. This represents an innovation compared to other works on the Medjerda. We selected the area near the wastewater treatment station in Bousselem, where pollutant injection scenarios were simulated to highlight pollutant dispersion along the stretch. This will provide alert data for potential accidents.

In our study, we opted for hydrocarbons as the choice pollutant, in line with the presence of an industrial zone connected to the Bousselem wastewater treatment plant. Hydrocarbons in our simulation are considered to be conservative passive pollutants that do not undergo significant chemical reactions or major transformations in water throughout the simulation period. The pollutant’s movement is therefore entirely governed by river hydrodynamics.

Our numerical model should thus provide an assessment, diagnostic, and analytical tool for understanding the phenomenon of substance migration in the river. Using different initial concentrations in the discharge, we will model tracer mixing in a river section, involving complex phenomena such as lateral mixing and longitudinal dispersion.

We use the Telemac 2D code to model free-surface flows in two dimensions of horizontal space, and thus to simulate the conservative tracer dispersion along a river section [12]. This model offers the possibility of placing pollutant tracer injection sources (with or without discharge) at any point in the domain. Indeed, it is possible to consider the passive tracer’s numbers transport [13], which may or may not be diffused. Most tracers are subject to diffusion and dispersion in the river environment; these factors will determine their paths in the domain and their concentrations during transport. Diffusion is, therefore, a determining factor in the monitoring of passive tracers such as dyes. Dispersion must also be considered since it is involved in the tracer mixing by physical phenomena [14]. The model outputs will be used to estimate the concentration time series as a function of time and space.

In order to better adapt to reality, a vegetation study in the field was focused along the watercourse in the study area, with the aim of setting the roughness coefficient correctly and thus being able to study the pollutants dispersion in the presence of vegetation. Vegetation is an important characteristic of river areas interacting with water flow [15], and river sections, therefore, have variable roughness along the wetted perimeter [16].

Concentration curves will be estimated in several longitudinal sections, which will enable to analyse the dispersion phenomenon evolution with time in the presence of vegetation. Various accidental release concentration scenarios will be considered, which will enable decision-making support plans to be drawn up in the event of an emergency.

The Considered Study Area in the Medjerda

The Medjerda wadi section passing through Bousselem city was chosen to illustrate the above-mentioned impacts, where the Medjerda wadi receives water from the Bousselem wastewater treatment plant (WWTP) as well as discharges from the existing industrial estate (Fig. 1). Bousselem city, situated along the Medjerda banks at 36° 36′ 40′′ N, 8° 58′ 11′′ E, features a wastewater treatment plant receiving effluents from residential, industrial, and marginal areas. Despite its 2000 commissioning, the Bousselem WWTP, with a 3,000 m³/day capacity and 75% saturation, operates below 60%. Approximately 14,862 residents, with a 98% connection rate, contribute to the plant’s flow rate of 2,043 m³/day. Regrettably, discharges from the plant are non-compliant, posing a risk of resource contamination [17]. The treated water is currently released into the Medjerda River.

Fig. 1. Discharge from the Bousselem WWTP.

In this part of the Medjerda, Tamarix is the dominant invasive plant species along the watercourse. Ten observation points were set up in February 2020 in this study area to describe the existing vegetation distributed along the watercourse (Fig. 2).

Fig. 2. Vegetation observation station in the study area.

Table I outlines the distribution of vegetation formations along the watercourse, with riparian vegetation, covering about 102 ha, being notably prominent. Tamarix, particularly Tamarix gallica (L.), dominates the landscape, followed by Eucalyptus and Acacia species. Despite dredging efforts post-floods, especially in proximity to residential and industrial zones, Tamarix trees persist along the entire watercourse and both banks [18]. Recolonization swiftly occurs in dredged areas, manifesting as 0.5 m to 3 m high bushes with a 25% to 50% cover. Additionally, denser stands, aged over 5 years and reaching 10 m in height with trunks of 20 cm to 30 cm in diameter, occupy certain sections.

Table I. The Tamarix Characterization in the Study Area

Statement no.	C1	C2	C3	C4	C5	C6	C7	C8	C9	C10
Max. vegetation height (m)	6	3	3	4	4	5	15	15	15	15
Tamarix cover (%)	50–80	25–50	25–50	25–50	50–80	25	100	25–50	80	50
Max. height Tamarix (m)	3	3	3	2	3	2	10	3	10	5

Model Construction

The part studied describes the downstream part of Bousselem City; it has a fairly regular slope of around 25 cm/km (0.0025%) and is fed by the discharge from the Bousselem WWTP (Fig. 2). The area considered for modelling using Telemac 2D extends over a total surface area of almost 985 ha, with a stream length of around 7.2 km and a width of around 2800 m. One difficulty in constructing meshes is selecting the optimum number of meshes to ensure the best representation of the study area [19]. An unstructured triangular surface mesh was developed by the Blue Kenue model, using the finite element method after processing the set of points imported from a LIDAR survey covering the area under consideration (wadi beds and floodplain) that was collected from the Ministry of Agriculture.

In view of the study objective, which is essentially concerned with flows in the minor bed, with relatively low flows that do not reach the floodplains, we chose a variable, more refined mesh for the minor bed, with a mesh size of approximately 5 m longitudinally in order to correctly represent flows in the main direction. The mesh size for the rest of the domain is of the order of 50 m. The grid shown in Fig. 3 consisted of 23286 nodes and 46274 elements.

Fig. 3. Grid and bathymetry of the Medjerda (downstream section of the town of Bousselem).

Model Calibration

The cross-sectional analysis of Medjerda’s cross-section reveals four distinct zones. The first corresponds to the minor bed, where the flow is constant, representing the smoothest area devoid of any vegetation. The second zone, adjacent to the minor bed, is the middle bed, characterized by the absence of vegetation and a frequent tendency to be submerged by water. The third zone encompasses tamarisk vegetation, which becomes submerged during more or less significant flow rates. Finally, the fourth zone, located at a greater distance, exhibits a variety of vegetation cover and is often utilized for agricultural purposes. The latter zone will not be affected by our model, as it remains non-submerged even under high flow rates.

In order to accurately reproduce flow conditions, TELEMAC-2D requires the identification of specific friction forces as a function of geometry. The coefficient of bottom friction, which reflects the roughness of the terrain, is the most significant calibration parameter. In our context, by analysing a high-quality, high-resolution orthophoto (1 m), we have identified the various flow zones (Fig. 4). These zones will be exported in polygonal form, to which a Manning’s coefficient value will be applied. The calibration phase aims to adjust the roughness of the terrain, represented by Manning’s coefficient in our case, in order to obtain water heights corresponding to the actual records associated with the measured flows (Fig. 5).

Fig. 4. Typical distribution of the roughness coefficient in a cross-section as a function of land occupation in the study area.

Fig. 5. Identification and distribution of the roughness coefficient in a cross-section as a function of land use in the study area.

Following the significant dredging works undertaken by the Ministry of Agriculture from 2015 to 2018, the geometric configuration of the riverbed was altered, rendering obsolete the previous rating curve of the Bousselem station. These endeavors aimed to widen the sections, thereby increasing the transit capacity to enhance the city’s resilience against major flood events. It is noteworthy that there have been no high-flow floods after the recent topographic survey campaign in 2020. The Ministry of Agriculture conducted new gauging measurements to calibrate the new heights corresponding to various discharge rates, with the maximum recorded discharge during this measurement campaign being 345 m³/s (Table II). We utilize these measurements as a reference to determine the roughness coefficient specific to each zone, ensuring that the calculated heights in our model align with these observations.

Table II. Calibration Results

	Flow (m3/s)	Water level (m)
Gauging results	345	124.79
Modelling results	345	124.82

Thus, a well-calibrated bed roughness coefficient provides a model that is very close to reality. For our case the values chosen are respectively 0.033 s/m⅓, 0.10 s/m⅓, 0.04 s/m⅓ and 0.06 s/m⅓ for the open minor bed, the middle bed occupied by tamarix, the cleaned middle bed and the major bed with vegetation.

Results and Discussion

After meticulously calibrating and validating the model, we embark on the reconstruction of the hydrodynamic model using Telemac-2D, with a fixed flow rate of 50 m³/s as the initial condition for our study. This flow rate was chosen strategically to ensure a controlled flow along the banks covered with vegetation, preventing any overflow. This methodical approach positions us to achieve our objectives by establishing a precise model, ready to explore predefined scenarios and thoroughly assess the impact of pollutant discharges on water quality in the studied system. This carefully selected flow rate aims to simulate conditions where the flow spreads along the vegetated banks without causing overflow, contributing to a comprehensive understanding of the environmental consequences of pollutant releases.

To be able to analyse the space-temporal development of concentrations, it is important to understand the velocity fields distribution in the river. In pollution transport, the velocity evolution is important and necessary to obtain a very precise mathematical model with all the physical processes involved [20].

The results in Fig. 6 of the model’s hydrodynamic simulation show that water velocity varies along the length of the river course. This variation can be explained firstly by the variation in the slope and typology of the wadi, as there is a considerable change in velocities in the meanders, where there is additional resistance caused by changes in the flow direction. The sinusoidal shape with several meanders increases the total length of the river, which reduces the water’s overall flow velocity.

Fig. 6. Water velocity distribution and heights in the Bousselem studied area: (a) Speed distribution and (b) Water height.

Transversally, the water velocity is higher at the center of the river and decreases towards the banks with values ranging from 0.05 m/s on the banks to 1.2 m/s at the center of the flow (Fig. 6). The velocity decreases on the wadi banks as the water depth is less important, which increases the roughness effect of the wadi bottom. Another reason for the reduction in velocity on the banks is the presence of vegetation, which increases roughness.

Given that the plant is connected to an industrial zone, it was decided to reconstruct three hydrocarbon dispersion scenarios that could be reached accidentally. Three different discharges with different initial concentrations (Ci) of the order of 0.25, 10, and 50 kg/m³ were carried out over a period of 4 hours. The discharge flow rate is 50 l/s, and the wadi flow rate is 50 m³/s.

Fig. 7 shows the pollutant propagation as a function of time for a continuous discharge with a flow rate equal to 50 l/s and an initial concentration of 50 kg/m³ in Medjerda. The flow rate is considered fixed, equal to 50 m³/s. The simulation results show that the average pollutant concentration decreases as a function of time and distance downstream in a river. As the simulation models the pollutant propagation over the field adjacent to the injection point, we observe a transverse variability in concentrations, which can be explained by the significant two-dimensional nature of dispersion in this upstream zone. We also observe that the influence of this factor diminishes with distance from the discharge point because in the far field, where the dispersion is one-dimensional, the stain occupies the entire cross-section, we are no longer interested in the dispersion details, but in the overall decrease in the average concentration as a function of time and distance downstream. The final concentrations recorded at T = 4 h in our modelled area (6 km from the initial injection point) are respectively 500 µg/l, 12 µg/l and 2.4 µg/l for the 50 kg/m³, 1 kg/m³, and 0.25 kg/m³ scenarios. Bearing in mind that the international standard for drinking water is 1 µg/l, these values are much higher than the standards, which puts any operator of this water at risk.

Fig. 7. The hydrocarbon concentration evolution as a function of time in the wadi (Q = 50 m³/s, Qi = 50 l/s, Ci = 50 kg/m³): (a) T = 30 m, (b) T = 60 m, (c) T = 180 m, (d) T = 240 m.

To understand this cross-sectional variability in concentration, three successive downstream sections were chosen to plot the concentration variability in space and as a function of time. These sections are respectively 0.6 km apart for the first section S1, 1.3 km apart for the second section S2, and 2.2 km apart for the third section S3.

The simulation results show significant spatial variations in pollutant concentrations. The curves show that the maximum concentration is in the deepest part of the wadi and decreases on the banks. This corresponds to the transverse variation in flow speeds, which reaches a maximum in the flow central part and decreases on the banks. It can also be seen that maximum concentrations decrease as a function of the distance from the discharge point and as a function of the initial concentration of the pollutant injected since the pollutant’s passage time is greater in a section located downstream than in a section located upstream of the flow [21].

Fig. 8 shows that the maximum concentrations are reached in the zones where the depths are greatest. It should be noted that the maximum depths are sometimes eccentric due to the meander’s presence, which explains why the maximum concentrations are also eccentric according to the wadi bed topography.

Fig. 8. The hydrocarbon concentration evolution on cross-sections in the wadi: (a) Ci = 0.25 kg/m³, (b) Ci = 1 kg/m³, and (c) Ci = 50 kg/m³.

The Fig. 9 shows the variation in concentration on S2 for different times corresponding to the pollutant propagation. Since our discharge is continuous, the average concentration tends to increase until it reaches a stable maximum value.

Fig. 9. The hydrocarbon concentration evolution as a function of time in the second section (Q = 50 m³/s, Qi = 50 l/s, Ci = 50 kg/m³).

The figure clearly shows this phenomenon, as the concentration increases at every point in the section until it reaches a maximum value at every point in the section. It can also be seen that the concentration approaches maximum values in the flow central part well before the banks. This can be explained by the roughness factor, which becomes increasingly important on the banks due to the lower water level and the vegetation presence.

To better understand the vegetation presence impact on the concentration’s variability in a section, a fourth scenario was developed, in which the total vegetation absence on the wadi banks was assumed, thus imposing a homogeneous roughness coefficient corresponding to the bottom smooth part.

Fig.10 shows that pollutant concentrations are around 0.001 kg/m³ higher on the banks when vegetation is present. We can deduce from this that the vegetation presence increases pollutant concentrations as they slow down velocities due to their higher roughness coefficient. In the central part, our fourth scenario presents higher concentrations of the order of 0.001 kg/m³ since the flow rate and initial concentration of the pollutants injected are the same for both scenarios. Another reason is that the vegetation presence did not exceed 10% of the total wetted surface area spread over the shallow banks, which did not influence the water level in the wadi.

Fig. 10. Hydrocarbon concentration trends in the second section in the presence and absence of vegetation (Q = 50 m³/s, Qi = 50 l/s, Ci = 50 kg/m³).

Conclusion

In this study, our decision to model the transport of a pollutant along a section of the Medjerda from the Bousselem wastewater treatment plant aims to assess the impact of pollutant discharges on water quality. Although the dynamics of pollutants in aquatic environments are inherently linked to flow characteristics, the model construction provided us with an advanced understanding of this dynamic. This approach has yielded significant results, contributing to informed decisions regarding water resource preservation. we employed Telemac 2D to construct a hydrodynamic model covering a 7 km stretch at the outlet of the Bousselem station, maintaining a constant flow of 50 m³/s. Our findings indicate that water velocity is highest in the river’s central region, gradually diminishing towards the banks, ranging from 0.05 m/s to 1.2 m/s.

To further investigate the dynamics, we conducted three distinct flow discharges, each at 50 l/s, with varying concentrations of 0.25, 10, and 50 kg/m³ over a 4-hour period. We selected three successive downstream sections to monitor concentration variability both spatially and temporally. These sections, denoted as S1, S2, and S3, are positioned at intervals of 0.6 km, 1.3 km, and 2.2 km, respectively. The maximum concentrations observed across the three sections decrease with increasing distance from the discharge point and the initial concentration of the injected pollutant. Notably, the concentration peaks earlier in the central part of the flow than at the banks. This phenomenon is attributed to the heightened roughness factor on the banks, resulting from lower water levels and the presence of vegetation.

At T = 4 h, the final concentrations in our modeled area (6 km from the initial injection point) were recorded at 500 µg/l, 12 µg/l, and 2.4 µg/l for scenarios with concentrations of 50 kg/m³, 1 kg/m³, and 0.25 kg/m³, respectively. It is crucial to note that these values significantly surpass the international standard for drinking water, set at 1 µg/l, and so posing a substantial risk to any operator utilizing this water.

Furthermore, we explored a fourth scenario, specifically investigating the consequences of the total absence of Tamarix vegetation on the wadi banks. In this scenario, we imposed a uniform roughness coefficient aligned with the smooth bottom part of the riverbed. The intention was to isolate the impact of vegetation on pollutant dispersion dynamics. The outcome of this scenario revealed pollutant concentrations approximately 0.001 kg/m³ higher on the banks when Tamarix vegetation was entirely absent.

This intriguing finding suggests that the presence of Tamarix vegetation plays a mitigating role in pollutant concentrations. The absence of vegetation appears to contribute to a slight elevation in pollutant levels along the banks. This deduction implies that Tamarix vegetation, with its unique characteristics, may act as a natural buffer or filter, influencing the transport and dispersion of pollutants within the water system. The heightened roughness and structural characteristics of Tamarix, when present, seem to create conditions that attenuate pollutant concentrations, potentially by impeding the movement or enhancing the settling of pollutants along the riverbanks.

In conclusion, our hydrodynamic model and subsequent analyses provide valuable insights into the flow dynamics and pollutant dispersion, highlighting the importance of considering environmental factors such as vegetation in assessing water quality and potential risks for water operators.