Abstract
This paper describes the development of a two-dimensional water quality model that solves hydrodynamic equations tied to transport equations with reactions mechanisms inherent in the processes. This enables us to perform an accurate assessment of the pollution in a coastal ecosystem. The model was developed with data drawn from the ecosystem found in Mexico’s southeast state of Tabasco. The coastal ecosystem consists of the interaction of El Yucateco lagoon with Chicozapote and Tonalá rivers that connect the lagoon with the Gulf of Mexico. The results of pollutants transport simulation in the coastal ecosystem are presented, focusing on toxic parameters for two hydrodynamic scenarios: wet and dry seasons. As it is of interest in the zone, the transport of four metals is studied: Cadmium, Chromium, Nickel, and Lead. In order to address these objectives, a self-posed mathematical problem is solved numerically, which is based on the measured data. The performed simulations show how to characterise metals transport with an acceptable accuracy, agreeing well with measured data in total concentrations in four control points along the water body. Although for the accurate implementation of the hydrodynamic-based water quality model herein presented boundary (geometry, tides, wind, etc.) and initial (concentrations measurements) conditions are required, it poses an excellent option when the distribution of solutes with high accuracy is required, easing environmental, economic, and social management of coastal ecosystems. It ought to be remarked that this constitutes a robust differential equation-based water quality model for the transport of heavy metals. Models with these characteristics are not common to be found elsewhere.
1. Introduction
The concern for water environmental pollution by heavy metals has recently increased due to the negative effects it might have in human beings [1, 2]. Some heavy metals as Cadmium (Cd), Chromium (Cr), and Lead (Pb) may transform into persistent metallic compounds with high toxicity [3]. Due to their damaging effects on the ecological environment and on human health [4, 5], it is necessary to study heavy metal contamination in aquatic ecosystems [6].
Metals are naturally present in small concentrations or traces in earth’s crust; many of them are essential for the growth and development of plants, animals, and human beings. The geo-available origin of these metals occurs from the mother rock to the soils after being released by weathering. In contrast, the presence of high concentrations of metals with respect to the ecological norms is an indicator of anthropogenic activities, such as hazardous wastes derived from industrial activities, mining, and agriculture.
As rivers serve as a medium for transport of dissolved and particulate matter from continents to the ocean, nowadays, interest in the pollution of rivers by metals has increased along with the exponential increment of industrialisation, urbanisation, and agriculturisation of coastal areas. This has substantially increased the concern and level of awareness in this problem [7]. For this reason, heavy metal concentrations in waters have been analysed worldwide, particularly by proposing new numerical approaches [8].
In coastal waters, heavy metals are distributed through the water column (particulate and dissolved) and the bottom sediments. This occurs during the mixing of fresh and marine water, which causes flocculation and sedimentation of organic matter, nutrients, and trace elements from rivers. Actually, dissolved metals come into the particulate phase due to processes as flocculation, water pH, sediment mineralogy, and others during estuarine mixing [9]. Thus, heavy metals get bound to these elements and precipitate to the bottom.
In this work, it is assumed that the partition coefficient does not depend on the concentration of the sorbing solids, according to Thomann and Mueller [10], in which the hypothesis is that the partition coefficient of metal in water is different from the partition coefficient of that metal in the bottom sediment and it is assumed that the decay in the sediment is approximately zero. This analysis applies to rivers where solids are not suffering a net resuspension in the water column; thus, this model was used to evaluate the concentrations of Cd, Cr, Pb, and Ni in the water column. On the other hand, according to Shimazu et al. [11], the sediment-water partition of the chemical mainly depends on sorption to sediment organic matter, sediment inorganic matter, and reaction group.
Flocculation plays a key role in the dynamics of estuarine and coastal environments, controlling the transport of fine-grained cohesive sediments and particulate contaminants throughout these systems [12–14] (usually characterised by muddy bottoms [15]). Nevertheless, it should be pointed out, that during natural estuarine mixing, flocculation process may not occur; actually, salinity plays an important role in the process, depending on the reaction mechanism of a particular metal. For instance, flocculation starts at 10% of salinity during estuarine mixing for Cd [16, 17]. Moreover, other metals are known for their nutrient-like behaviour [18]. Thus, flocculation process constitutes an arduous task to model [19], which is not the aim of this work.
For the above reasons, strategies and tools to mitigate the pollution of heavy metals are required [20]. A huge number of mathematical models that intend to predict the transport of heavy metals in flows exist, for example, the statistical models based on exponential functions (analytic models), which allow the achievement of a simulation in a relatively uncomplicated way. This is the case of the use of sigmoid functions to determine metal concentrations in rivers that can yield to average concentrations in a section [21]. Despite the fact that the power of CPU processors is now much better than decades ago, simulations continue to be carried out through analytic models, which allow using a minimum of experimental measurements [22].
Most common methods to evaluate heavy metals pollution in water bodies are based on quality indexes, which generally use correlation or fuzzy methods for their estimation [2, 23, 24]; these works model pollutant distribution in the area under study, through GIS system modules, which apply interpolation methods. Nevertheless, the obtained spatial distribution through a hydrodynamics module along with transport equation and considering reaction mechanisms produces much more accurate results [25] (although it requires a greater effort to be implemented).
Water quality models (WQM) have increased in number and have improved in recent years, focusing on the study of the water quality as well as pollutant transport in shallow water ecosystems. The behaviour and transport of toxic substances, such as metals and/or hydrocarbons, have been deeply studied on shallow aquatic systems during this century [26, 27]. In this case, the complexity of estuarine coastal systems must be understood in order to clearly pose the solution of the equations representing water hydrodynamics (mass conservation and momentum equations) as well as mass transport of pollutants (advection, diffusion, and reaction equation), considering even the highly nonlinear interactions typical of these regions. Thomann and Mueller [10], Thomann [28], Lun et al. [29], Ji et al. [30], Shimazu et al. [11], and Bhavsar et al. [31] show different approaches of the toxic substances behaviour problem in water columns, as well as their interaction with sediments and air, including mathematical models and solution methods.
The most popular numerical WQM are AQUATOX, Branched Lagrangian Transport Model (BLTM), One-Dimensional Riverine Hydrodynamic, Water Quality Model (EPD-RIV1), QUAL2Kw, Water Quality Analysis Simulation Program (WASP), Water Quality for River-Reservoir Systems (WQRRS), ROMS-ICS [32], MIKE Ecolab/ABM, and IberWQ [33]. Nevertheless, most of them are based on solving mass balance/advective diffusion equation but oriented to nutrients (as DO, BOD, , and Phosphorus) or pathogens (such as coliforms, like the Escherichia coli [34]). Such models do not contemplate metals transport [35]. A review regarding computational models of water quality can be found in the article by Wang et al. [36] who pointed out that most models such as MIKE models, EFDC, and Delft 3D model have been applied to simulate water environmental quality in most cases of environmental impact assessment. However, little information is available on the differences in model results from different models and the suitability and parameter sensitivity of these models.
In recent decades, there has been a boom in the development of hydrodynamic models with coupled water quality models; here we can also mention the POM model, MIKE3, COHERENS, ROMS, and MOHID model.
The MOHID three-dimensional model has the ability to simulate complex estuarine and coastal flows in numerous applications dealing with mesomalar coastal lagoons, tidal channels, and estuarine systems [37]. It also solves the transport equation for salinity and temperature, but the free version is not quite user-friendly and neither is its commercial version.
MIKE 3 model is a general three-dimensional hydrodynamic model for flow simulation in estuaries, bays, and coastal and oceanic areas [38]. Although it is a fairly robust and complete model, only the commercial version is available.
ROMS model (Regional Ocean Modeling System) has been used to simulate the water circulation in different regions of the world’s oceans at different scales (local and basin) [39]. Similarly, POM model (Princeton Ocean Model) is a three-dimensional hydrodynamic model based on semi-implicit finite differences [40].
COHERENS model is a multipurpose three-dimensional model based on finite differences. This model allows the coupling with different submodels that simulate physical and biological processes, as well as the transport and transformation of sediments and pollutants [41].
Most of these last models are commonly used today, but most of them do not focus their strength on the solution of toxic substances transport such as heavy metals; they are more focused on hydrodynamic ocean domain than interaction of coastal-lotic-lentic water bodies.
Although there are available WQM as discussed above, it is common to find models developed by particular researchers that include effects or processes specific to their case of study, such as the MINEQL [42], which has been used to model the concentrations of metals in rivers since the 1980s.
Despite the existing WQM, the aim of this work is to develop a water quality module, which can be coupled to a hydrodynamic numerical model, applicable to the study of hydrodynamics and water quality in coastal ecosystems. The hydrodynamic model adopted in this work is self-developed for research purposes and has been previously used in different research applications, including the hydrodynamic-hydrological modelling in flood zones [43], the modelling of flows through vegetation [44], the modelling of the thermal discharges [45, 46], and the modelling of fresh water plumes in river-sea interaction [47]. On the other side, in order to accurately simulate turbulent viscosity, the turbulence model discussed in detail in the paper by Rodriguez-Cuevas et al. [48] is implemented.
The water quality module developed herein takes into account many parameters, including the advection-diffusion-reaction mechanism. The required parameters cannot be found in a specific existing WQM. These parameters were adopted from Ambrose Jr et al. [49] and Kannel et al. [50]. Four heavy metals, Cd, Cr, Pb, and Ni, were selected for the water quality module because they are required by local institutions as well as because they have strong toxicity to humans [51]. According to the previous knowledge of the authors, they did not find an ad hoc module for the transport of these metals, so they developed one that contemplates the greatest possible number of parameters in the reaction mechanisms. This module can be applied as any WQM (ROMS-ICS, WASP or MIKE), with the advantage that the reaction mechanisms can be modified by calibrating each of the parameters through the equations. The case of study is the ecosystem composed by El Yucateco lagoon, Chicozapote river, and Tonalá river discharge, which has suffered a serious deterioration due to pollution originated by oil industry activities [52, 53].
In this area, several studies have been carried out in the last decade, related to the study of the dynamics of the river-lagoon-sea interaction [54, 55]. On the other hand, the Institute of Marine Sciences and Limnology (ICMyL) from Universidad Nacional Autónoma de México (UNAM) has conducted studies related to the monitoring of heavy metals and toxic substances associated with oil activity for eight continuous years [52, 56]. Based on these previous studies and heavy metal monitoring databases, we made the configuration and development of the proposed model in this work.
This paper begins by a detailed presentation of the coastal ecosystem in Section 2. It then goes to exposition of the methodology applied herein, including the measurement protocol and the design of the hydrodynamic and water quality modules as well as the numerical strategies for their solution (Section 3). Thereafter, Section 4 presents the application of this WQM to the case of study, as well as the detailed results of the simulations for the transport of heavy metals, along with two metrics that allow assessing the accuracy of the WQM developed herein with respect to measured data. Finally, in Section 6, the conclusions of this study are presented as well as some final interesting remarks.
2. Case of Study
2.1. Location of the Region under Study
The coastal ecosystem under study is located at the east of Tabasco State, in the southeast of Mexico. The lagoon is located between LAT 18° 10’ and 18° 12’ N and between LONG 94° 02’ and 94° 00’ W (Figure 1). El Yucateco lagoon interacts with Chicozapote and Tonalá rivers before discharging to the Gulf of Mexico. Nowadays, the lagoon is one of the most important water bodies of the region.
2.2. Evolution of the Region
For hundreds of years, the main activities in the zone have been agriculture, cattle farming, and fishing. Nevertheless, oil field Cinco Presidentes was established in 1963 in the vicinity of this region, with El Yucateco lagoon being the closest water body to the oil facility. Later, a network of artificial channels (approx. 33 Km) was built by an oil company during such decade, with an area of about 130 ha. These channels drain a considerable volume of fresh water to the lagoon product of the drainage of the floodplains (marshes) that surround the area.
In the last 20 years, important changes in hydrology and water quality have occurred, changing the productivity and reducing the number of endemic species. These changes have had social, economical, and commercial impacts for native population. This area is currently under industrial development, where two kinds of activities stand out: oil (extraction and production) and livestock farming, which together account for almost 90% of the productive sectors of Tabasco State [57]. More environmental information about this ecosystem is available in [56].
2.3. Climate
The predominant climate in the region is warm and humid, with abundant rainfall through the whole year, particularly during autumn. Its annual thermal regime oscillates between 25.8°C and 27.8°C; the highest average temperature occurs in May with 29.4°C, while the minimum is registered in January with 23.1°C. September/October is the most rainy period with an average precipitation of 400 mm, while June features the minimum precipitation with an average value of 43.3 mm [58].
2.4. Hydrodynamics
For El Yucateco lagoon, the main hydrodynamic driven forces are winds and tides. The tide and Chicozapote river permanently renew and refresh water in the lagoon (see Section 3.5 for further details). For dry season (see Table 2), current flows head predominantly to north during mornings while mainly to northeast in afternoon hours; in both cases, current flows out of the lagoon with velocities greater than 20 cm/s. For rainy (wet) season (Table 2), current is variable with different directions along the system, heading to northeast in the south part of the lagoon and to the north close to the mouth and to the connection with Chicozapote river. Typical flow velocities about 10 cm/s are observed, which favour the formation of vortices within the lagoon.
In this zone of the Gulf of Mexico, tide presents a mixed type with diurnal influence, whose oscillations are not greater than 30-60 cm. The surge is moderate in E-W direction, with a maximum height of 2 m in normal meteorological conditions.
3. Methodology
The WQM developed herein consists essentially of two parts: a hydrodynamic module and a water quality module. The later one is in charge of transporting heavy metals through the water body by using the hydrodynamic module results, which has been calibrated and tested previously [44, 45, 48, 59, 60]. Later, in Sections 3.3 and 3.4, these modules are appropriately defined. In what follows, the protocol regarding measurement and data analysis is presented.
3.1. Water Quality Measurement Protocol
In measurement campaigns, water and sediment samples were collected following the guidelines of the Standard Methods for the Examination of Water & Wastewaters methodology [61].
For heavy metals in water, 500 ml of sample was taken at each point using a van Dorn water sampler. The samples were filtered (0.45 μm) and the pH was adjusted to 2 using , stored in amber glass bottles and refrigerated at 4.0°C for transportation [61]. Filters for the analysis of particulate metals were conserved in polyethylene bottles, with 2M, while those used to quantify the suspended material were kept in petri dishes, where they were dried and weighed later. The filtered particles were analysed in order to obtain concentration of metals in particulate phase, while the water obtained from this process was analysed to obtain the corresponding dissolved metals in water.
For heavy metals in sediments, samples of 100 g were collected from the surface layer of the bottom sediments (max. 5 cm) with an Ekman dredger. These samples were stored in sterile polyethylene bags and kept at 4.0°C for transport. Particulate and dissolved concentrations of metals were also obtained in sediments.
All the samples were collected in duplicate to determine the precision of tests and sample handling. In order to minimise the effect of hydrological input from rivers, the sediment and water sampling were performed during mornings from 8:00 to 12:00 hours in low tide conditions when there is little penetration from the sea to Chicozapote river and El Yucateco lagoon.
The chemical analyses of heavy metals in water and sediments were performed using atomic absorption spectrophotometry, with spectrophotometer (Shimadzu, Mod. AA 6800).
According to previous in situ studies [56], several metals were initially sampled. Those whose parameters exceeded water and/or sediment quality criteria [62, 63] were selected for diagnosis: Cd, Cr, Pb, and Ni.
Measured data and the chemical analyses of these parameters served to specify boundary and initial conditions to the numerical model developed herein. They also aided in validating the model for the field site.
3.2. Numerical Modeling
The numerical model used in this work is self-developed, strongly based on the proposal of Casulli and Cheng [59]. In this case, 2 modules were applied: the hydrodynamic module and the water quality module (WQM). Figure 2 schematizes the simulation process of the hydrodynamic model to compute the distribution of the velocity components. Once the velocity field is obtained, it constitutes the input for the WQM, which is first executed and then validated with measured data. Then, the WQM is executed iteratively for calibrating some of the coefficients that appear in the equations, until calculations match measured data within an expected error (see Figure 3). Thus, the WQM developed herein is of high computational burden, even though the time step used in the WQM is much greater than that used in the hydrodynamic module.
3.3. Hydrodynamic Module
The hydrodynamic module used in this work is based on the two-dimensional shallow water equations, which can be derived from the Reynolds-averaged Navier-Stokes equations [48]. These equations arewhere and are the depth-averaged velocity components in and directions, respectively (m/s), is the acceleration due to gravity (), is the free surface elevation (m), is the horizontal eddy viscosity (/s), is the water depth (m), and are the wind shear stress terms in and directions, respectively, and and are the bottom shear stress terms in and directions, respectively. It is noted that the units of the wind and bottom shear stress terms are /.
The equation to calculate the free surface elevation () is obtained by integrating the continuity equation over the water depth and by applying a kinematic condition at the free surface, yielding
As mentioned above, in order to achieve more realistic hydrodynamics, mechanical dispersion phenomenon is also considered through the introduction of a model of turbulence. Due to the nature of the water body treated herein, the following mixing-length model is used [48], which contributes to the horizontal eddy viscosity aswhere is the horizontal mixing length (m), is a dimensionless constant, and (m) is defined aswhere is the bathymetry (m), , and and are dimensionless constants (the latter so called von Kármán constant). This hydrodynamic module has been calibrated according to Casulli and Cheng [59] and has been tested in the works of León et al. [60], Ramírez-León et al. [45], Barrios-Piña et al. [44], and Rodriguez-Cuevas et al. [48].
The wind shear stress terms in and directions are calculated withwhere is the air density (kg/), and are the wind magnitude components measured m above the ground level in and directions, respectively, and is the wind drag coefficient obtained from the Garrat formulation , where is the magnitude of the wind velocity vector m above the ground level (m/s). The bottom shear stress terms in and directions are calculated as follows:where is the Chezy coefficient.
The numerical solution scheme of the hydrodynamic module is based on a second-order finite difference formulation in both time and space. The solution method is an adaptation of the semi-implicit Eulerian-Lagrangian scheme proposed by Casulli and Cheng [59]. This method treats the advection and diffusion terms differently. The solution of the nonlinear advection terms of (1) and (2) is given by a Lagrangian formulation through the characteristics method, and the solution of the diffusion terms is given by an Eulerian formulation through the Adams-Bashforth scheme.
A 2D mesh is used for the numerical simulations, based on a staggered cell arrangement, as shown in Figures 4 and 5, where dot at the center of the cell represents the location of any scalar value (free surface elevation, water properties or pollutants, ) and circles on faces indicate the location of the velocity components and .
In Figure 5, is the water depth at point , where the is velocity component and is the water depth at point , where the velocity component is evaluated.
The solution of the system of (1), (2), and (3) is given through linear systems as follows:where the operators and join the advective terms and the turbulent diffusion terms as where and are the explicit velocity components, calculated at time using the characteristics method as where and are the Courant numbers in and directions, respectively, from which and are the integer parts and and are the decimal parts.
3.4. Water Quality Module
The WQM consists of two parts, the one in charge of transport and the one in charge of the reaction mechanisms to which the substance is subject. Within the quality module, the Advection-Diffusion-Reaction equation is solved aswhere is the concentration of a substance (mg/L); and are the horizontal dispersion coefficients (/s) [48]; and is the substance reaction term (mg/L) and are the point sources. Discretization of (16) is given in a similar way to that for the velocity of (1) and (2).
The WQM is initialized considering no reaction () and assuming stationary sediment, constant kinetic coefficients, and suspended solid which is uniformly distributed in space over the river reach. Once the initial concentration distribution is calculated through (16), it is reestimated through a model of reaction of toxic substances. Thomann and Salas [64] present a complete model that considers diffusive exchange, decomposition processes, volatilisation, settling, and resuspension (important processes in a coastal aquatics ecosystem), which is governed by the following ordinary differential equation:where and define water column and sediment bed conditions, respectively, is the concentration of the toxic substance (mg/L) (estimated using the transport equation), is the concentration of the toxicant in sediments (mg/L), is the diffusive exchange of dissolved toxicant between the sediment and the water column (m/d), is the river depth (m), is the dissolved fraction (1), is the particulate fraction (unitless), is the degradation rate of the dissolved toxic substance (), is the overall volatilisation transfer rate (m/d), is the vapour phase concentration (mg/L), is Henry’s constant, is the settling velocity (m/d), is the resuspended velocity (m/d), and is the sediment porosity (unitless). It ought to be mentioned that the concentrations and for any metal are the sum of both particulate and dissolved metals’ concentration [30]. For further reference on these terms, see Figure 6.
In (17), the first term on the right hand is the diffusive exchange of dissolved toxicant between sediment and water column. The second term is the decomposition processes of the dissolved form due to microbial decay, photolysis, hydrolysis, and so forth, in the water column (decay of particulate form is assumed to be zero). The next term is the air water exchange of the toxicant due to volatilisation or gaseous input. The next term represents the settling of the particulate toxicant from the water column to the sediment and the last term is resuspension into the water column of the particulate toxicant from the sediment.
A description of the processes underlying the parameters in (17) is detailed in what follows.
For the decay of the dissolved substances, the most important mechanisms in the degradation rate of the dissolved toxic substance () are represented by the equationwhere is the photolysis rate (), is the hydrolysis rate (), and is the microbial degradation rate ().
The overall exchange rate () estimates the importance of losses due to volatilisation. According to Mackay [65], the application of the two film theory yields an overall volatilisation transfer parameter given bywhere is the liquid film coefficient (m/d) and is the gas film coefficient (m/d). As seen from (19), depends on chemical properties as well as on characteristics of the water body such as water velocity (affecting ) and wind velocity over water surface (affecting both and ).
Henry’s () constant, present in (17) and (19), takes into account the water-atmosphere interaction, and it represents partitioning of the toxicant between water and atmospheric phases. Its dimensionless form is given bywhere is the water solubility (mg/L) and is the vapour phase concentration (mg/L).
The sediment diffusion rate (m/d) reflects the fact that gradients can occur between interstitial sediments and the adjacent water column [66]. An effective model for can be written, according to Manheim [67] and O’Connor [68], aswhere is the molecular weight (g/mol) and is the sediment porosity.
To determine both settling and resuspended velocities, a particle characterisation in the region under study is required. In this work, such characterisation was made with sediment samples taken from field campaigns, where organic matter, particle sizes, porosity, humidity, and so forth were quantified as follows: the organic matter was obtained by the wet oxidation technique using exothermic heating and oxidation of organic carbon of the sediment sample [69]. Humidity was determined by drying method [70], and the percentages of sand, silt, and clay were determined by the Bouyucos method [71].
In coastal ecosystems, a good approximation to the settling velocity is given by Hawley [72] who provides the following empirical relation:where is the particle diameter (m).
For the resuspended velocity , Thomann and Di Toro [73] have proposed the next equation:where is the net rate loss of solids (m/d) and is the settling net rate of the water solids (m/d).
In estuarine and coastal environments, usually characterised by muddy bottoms, flocculation plays a key role [15]; nevertheless, the water body under study has low circulation because there is no important water interchange with Chicozapote river. With this being a quasi-static water body with very small velocities, the sediment transport in the muddy environment at the bottom of the lagoon can be actually neglected.
3.5. Implementation of the Model to the Case Study
The hydrodynamic characteristics of water bodies such as coastal lagoons are governed by a slight balance between tidal forces, currents flow, wind stresses, and density, which induce pressure and friction forces at the bottom [74, 75], in addition to other factors such as the geometry and flow, which is predominantly turbulent with a horizontal length scale much greater than the vertical one [76].
The forcing boundary conditions of the sea-Chicozapote river-El Yucateco lagoon system were set by considering the astronomical tide, winds, and river discharges. The input information into the model as boundary conditions is data obtained from measurements for each season and/or generated by official institutions in the study area [56]. Wind forcing was implemented with data measured directly over El Yucateco lagoon during the field measurements campaigns in both dry and wet seasons (see Figures 7 and 8). Wind forcing was considered spatially invariant along the model domain. Tidal boundary condition was imposed on the northwest edge of the study domain at the river discharge zone by means of the measured variation of water height with respect to the mean sea level (see Figures 9(a) and 9(c)). The hydrological flow of Tonalá river was imposed according to the data shown in Figures 9(b) and 9(d). The points of application of boundary conditions are shown by the arrows in Figure 10(a).
(a)
(b)
(c)
(d)
(a) Bathymetry of the study zone, measured relative to the mean sea level, in metres (m). Arrows indicate gauging stations of hydrological flow and tide boundary conditions
(b) Mesh configuration for the study zone
The bathymetry shown in Figure 10(a) was acquired through an echo sounder and GPS. Several planned transects were carried out to obtain latitude, longitude, and depth in order to fully characterise the complete bathymetry of the water body. Subsequently, these data were interpolated to match the bathymetry with the central points of each cell of the numerical mesh.
The area under study was discretised through a structured mesh. A variable-spacing grid was used to make the model more efficient in zones of interest, with m and m, while m and m. The grid has 213 elements in the -direction and 79 elements in the -direction, for a total of 16827 elements of which 3720 are active elements (see Figure 10(b)). For a staggered cell, free surface elevation and metals concentrations are discretised using central differences, while velocity vector components are considered at the faces of the cell. Also, the time step is related to the Courant number and so to the resolution of the mesh. This guarantees stability and convergence of the numerical solution, while minimising the computational burden.
3.6. Validation Criterion of the Numerical Solution
In order to assess the quality of the numerical solution with respect to field measured data, the Root Square Mean Error and the Nash-Sutcliffe efficiency coefficient (NSE) are presented, with the last given by Horritt [77]: where is an observed/measured datum and is the respective datum obtained from numerical simulation (same place and time). The parameter is the measured data mean. Clearly, if is close to 1, it is possible to consider that numerical simulation data is quite approximate to field measured data.
A value of RMSE = 0 indicates a perfect fit. Some suggested values for decision-making related to the data produced by the RMSE are presented in Table 1 [78].
For NSE coefficient, the range is considered as an acceptable performance with 1 being the optimal value, whereas values indicate that the mean observed value is a better predictor than the simulated value, which indicates unacceptable performance [78].
4. Results and Discussion
4.1. Field Measurements
Tabasco State is the region of Mexico where the climate is very rainy, so seasons of the year are not highly differentiated, with abundant precipitations throughout the whole year. Thus, it is important to differentiate between the typical characteristic climate periods in the study area, such as summertime and wintertime, because the dynamics can significantly change [79]; for instance, they have influence on parameters that define the transport of solutes such as salinity [80] or sediments [31, 81, 82]. For the purpose of yielding to results which can be compared with stationary measurements in situ, wet and dry seasons for this study are selected, as detailed in Table 2 [52, 56]. The time spans selected for wet and dry seasons simulations guarantee that weather remains essentially constant through them. In other words, those periods are highly representative of their respective season [83].
In this work, metals concentrations were measured in four control points (C.P.). Two sampling points are located within the lagoon (C.P. 1 and 2) and two in the river (C.P. 3 and 4) (see Figure 1). For field measurement purposes, the measured data was performed in campaigns each 5 days, for the four referred heavy metals, according to the methodology defined in Section 3.1.
4.2. Validation of WQM
The validation process of the WQM consists in finding the values involved in the considered reaction mechanism (see Section 3.4 and particularly (17)) which yield to the most accurate simulated results with respect to field measured data.
In order to start the validation process, initial conditions are required. Such values are introduced to the computational domain on each control point at the beginning of the simulation. The initial parameters were imposed based on EPA recommendations. Then, a numerical simulation carried out forcing with measured concentrations in C.P. 1, C.P. 2, and C.P. 3 of Figure 1. After convergence, calculated concentrations in C.P. 4 were compared against measurements in the same location. The measurements considered for comparisons were taken at the end of the simulation time period. Then, a trial-and-error procedure was performed up till reaching an error lower than 5% by adjusting the involved parameters (in Table 3) in each simulation. In a second step, a new numerical simulation carried out forcing with measured concentrations in C.P. 1, C.P. 2, and C.P. 4. After convergence, calculated concentrations in C.P. 3 were compared against measurements in the same location. The trial-and-error procedure was performed up till reaching an error lower than 5% by adjusting again the referred parameters. The same procedure with control points C.P. 1 and C.P. 2 was effectuated to complete one calibration cycle. This cycle was repeated iteratively until the best correlation between calculations and measurements was obtained. The parameters that resulted from this validation process are shown in Table 3 and are the ones used for the simulations present in Section 4.4.
It ought to be noted that the values of the system coefficients lie within the typical values reported for them in the literature [64, 66, 73].
4.3. Hydrodynamic Simulations
In order to perform the hydrodynamic simulation, the hydrodynamic module has already been calibrated with the parameters shown in Table 4, as mentioned in Section 3.3.
The hydrodynamic simulations were performed using a time step of s, requiring 2,160,000 and 1,296,000 iterations for dry and wet seasons, respectively. Figure 11 shows the results obtained at the end of simulation.
According to results obtained from field measurements, it is possible to observe that, for dry season, there is an intense recirculation in the lagoon due to the influence of sea. Nevertheless, for wet season, river flow is higher, tide penetration almost does not exist, and lagoon’s behaviour presents lack of recirculation.
4.4. Simulation of Metals Transport
4.4.1. Cadmium Transport
Figures 12 and 14 show the final snapshots of the Cadmium transport simulation for dry and wet seasons, respectively. Figures 13 and 15 depict the behaviour of Cadmium at viewers (control points, see Figure 1). Tables 5 and 6 show RSME and Nash-Sutcliffe efficiency coefficient (see (25)) to evaluate the quality of the simulation against the discrete measured data and the fitted measured data, at every control point for dry and wet seasons, respectively.
4.4.2. Chromium Transport
Figures 16 and 18 show the final snapshots of the Chromium transport simulation for dry and wet seasons, respectively. Figures 17 and 19 depict the behaviour of Chromium at the control points. Tables 7 and 8 show RSME and Nash-Sutcliffe efficiency coefficient (see (25)) to evaluate the quality of the simulation against the discrete measured data and the fitted measured data, at every control point for dry and wet seasons, respectively.
4.4.3. Lead Transport
Figures 20 and 22 show the final snapshots of the Lead transport simulation for dry and wet seasons, respectively. Figures 21 and 23 depict the behaviour of the Lead at the control points. Tables 9 and 10 show RSME and Nash-Sutcliffe efficiency coefficient (see (25)) to evaluate the quality of the simulation against the discrete measured data and the fitted measured data, at every control point for dry and wet seasons, respectively.
4.4.4. Nickel Transport
Figures 24 and 26 show the final snapshots of Nickel transport simulation for dry and wet seasons, respectively. Figures 25 and 27 depict the behaviour of Nickel at the control points. Tables 11 and 12 show RSME and Nash-Sutcliffe efficiency coefficient (see (25)) to evaluate the quality of the simulation against the discrete measured data and the fitted measured data, at every control point for dry and wet seasons, respectively.
5. Discussion
The numerical simulations presented above were performed for wet and dry seasons. The hydrodynamic results show that, for dry season, tidal reversing currents occur affecting the whole ecosystem, favouring the formation of vortices within El Yucateco lagoon. However, for wet season, tidal effects decrease and river flows dominate, causing the weakening of the recirculating vortex.
The metals transport simulations show that, for both season scenarios, dry and wet, the concentration of metals maintains a stable behaviour, which is perturbed by oscillations due to hydrodynamics. The oscillations are greater for dry season, where the hydrodynamics are driven by tides. For the wet scenario, metals concentrations show very slight disturbances, maintaining an almost constant behaviour during the simulation period. In general, the concentrations are higher at control points located in El Yucateco lagoon and Tonalá river discharge (C.P.s 1, 2, and 4) than the one located near the confluence of Tonalá and Chicozapote rivers (C.P. 3). On the other hand, for the wet season, the concentrations of the control points located in the lagoon maintain higher values (C.P.s 1 and 2), due to the small currents that the lagoon possess.
On the other hand, as it can be observed from the structure of the model presented herein, it is of great importance to possess enough information to feed the numerical model with the necessary data in order to validate it and to obtain the results that better reproduce the pollutant transport. In other words, the use of numerical modelling does not exempt the in-depth acquaintance of the dynamics of the system under study. This implies several numerical simulations to achieve the best agreement with field measured data.
This in-depth acquaintance means to know, besides the initial concentrations of the toxicants to be simulated, intrinsic facts of the water body as its detailed bathymetry, contours, and boundary conditions as wind forcing, tides, flow discharges to the water body, and so forth in order to accurately determine the transport of toxic substances. In this sense, the required field work to accomplish the study with this model is enormous, but it provides excellent and more accurate results than most of models.
The procedure adopted to validate the water quality module coupled to the hydrodynamic module is based on a trial-and-error procedure with the aim of adjusting the coefficients of the reaction equation terms for each metal. Concentration field measurements were used to validate the model by adjusting its parameters until acceptable simulation was achieved. Although the validation procedure was designed specifically for the present case of study, it can be straightforwardly extended to other cases. Although this constitutes a robust method to validate the model that yields accurate results, it requires high computational burden and execution times.
6. Conclusions
In this paper, a hydrodynamics-based WQM for the specific evaluation of heavy metals is developed and tested. The model was set for a specific ecosystem located at the east of Tabasco State, Mexico. Chicozapote and Tonalá rivers and El Yucateco lagoon are the parts of this ecosystem where measurements of metals concentrations were obtained from field measurements.
A heavy metal water quality module, based on a laterally averaged two-dimensional hydrodynamics and sediment transport model, was developed and applied to the tidal Chicozapote and Tonalá rivers estuary. The model was validated with measured data, including time-series data and the spatial distribution of the suspended-sediment concentration. The calculated distribution of total Cadmium, Chromium, Lead, and Nickel concentrations along the river and lagoon generally agrees with the field-measured data. Thus, this module shows its potential in the estimation of toxic substances for water quality assessment in aquatic coastal systems. Albeit the application of the model for other ecosystems is limited, it ought to be examined carefully before being applied.
In this way, this model constitutes an excellent option when a distribution of the solutes with a high accuracy is required. Nevertheless, the most challenging fact on carrying out metal transport numerical simulations is the lack of data for validation as well as the lack of information about the values of the reaction equation terms coefficients.
Finally, it ought to be remarked that the necessity of developing this model arose from the fact that it was more convenient to solve and program the differential equations and be able to access and adjust all the model parametrisation to simulate in a better way the hydrodynamics and metals transport. Thus, beside the important assessment of metals transport in the considered estuary, the main contribution of this work is to provide a highly accurate water quality model able to deal with the dynamics of these toxicants in the water body.
Data Availability
The data measured to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.