Mathematical Problems in Engineering
Volume 2019 (2019), Article ID 8485031, 24 pages
https://doi.org/10.1155/2019/8485031

Research Article

Numerical Modelling of Heavy Metal Dynamics in a River-Lagoon System

F. Torres-Bejarano,¹ C. Couder-Castañeda,² H. Ramírez-León,³ J. J. Hernández-Gómez,² C. Rodríguez-Cuevas,⁴ I. E. Herrera-Díaz,⁵ and H. Barrios-Piña⁶

¹Departamento de Ingeniería Ambiental, Universidad de Córdoba, Montería, Colombia
²Centro de Desarrollo Aeroespacial, Instituto Politécnico Nacional, Mexico
³PIMAS Proyectos de Ingeniería y Medio Ambiente S.C., Ciudad de México, Mexico
⁴Facultad de Ingeniería, Universidad Autónoma de San Luis Potosí, San Luis Potosí, Mexico
⁵Departmento de Ingeniería Agroindustrial, Universidad de Guanajuato, Campus Celaya-Salvatierra, Celaya, Guanajuato, Mexico
⁶Tecnologico de Monterrey, Campus Guadalajara, 45138 Zapopan, Jalisco, Mexico

Correspondence should be addressed to C. Couder-Castañeda; ccouder@ipn.mx

Received 10 January 2019; Revised 1 April 2019; Accepted 18 April 2019; Published 6 May 2019

Academic Editor: Luis Cea

Copyright © 2019 F. Torres-Bejarano et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

Figure 1: Location of study zone and control points (C.P.) in the coastal ecosystem. The centre of lagoon is located at W 94° 01’ 12” N 18° 11’ 6”.

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.

Figure 2: Simulation process flowchart for hydrodynamic module.

Figure 3: Simulation process flowchart for water quality 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 are

where 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 as

where is the horizontal mixing length (m), is a dimensionless constant, and (m) is defined as

where 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 with

where 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 .

Figure 4: Discretization cell. Top view.

Figure 5: Vertical mesh arrangement.

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 as

where 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.

Figure 6: Depiction of reaction metal model parameters for a completely mixed lake. 1 stands for water column and 2 for sediment layer.

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 equation

where 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 by

where 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 by

where 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], as

where 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).

Figure 7: Wind rose used in the hydrodynamic module for dry season.

Figure 8: Wind rose used in the hydrodynamic module for wet season.

(a)

(b)

(c)

(d)

Figure 9: Forcing boundary conditions used in the study domain for dry season: (a) tide condition; (b) hydrological flow, and for wet season; (c) tide condition; (d) hydrological flow. The gauging stations where the hydrological flow and the tide condition were taken are shown in Figure 10(a).

(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

Figure 10: Bathymetry and mesh of the domain.

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].

Table 1: Criteria for the qualitative evaluation of goodness of fit for RMSE.

RMSE
Fit

> 0.70
Not satisfactory

0.60-0.70
Satisfactory

0.30-0.60
Good

0.00-0.30
Very good

Table 2: Considered time spans for wet and dry seasons.

Season
Initial date
Duration (days)

Dry
June, 1st
50

Wet
September, 1st
30

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.

Table 3: WQM parameters for the considered heavy metals that resulted from the validation process.

Variable
Value

Cd
Cr
Ni
Pb

(mg/L)
Initialised from measured data in C.Ps.

()
0.01

(g)
0.2

(m/d)
1.85
1.90
1.90
1.85

(m)
From hydrodynamic module

(m/d)
0.0074
0.0124
0.0114
0.0049

(g)
0.1

(mg/L)
By Control Point

(m/d)
1.5

(g)
0.8

(g)
0.9

(m/d)
0.00035
0.00015
0.00015
0.00035

0.0042

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.

Table 4: Calibrated parameters of the hydrodynamic module [44, 45, 48, 60].

Variable
Value

(m/s)
9.81

(unitless)
0.09

(unitless)
0.41

(unitless)
0.6

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.

Figure 11: (a) Dry season, day 50 snapshot. (b) Wet season, day 30 snapshot.

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.

Table 5: Root Mean Square Error (RMSE) and Nash-Sutcliffe efficiency coefficient (see (25)) obtained as quality control of the simulation for Cadmium transport in dry season. Values are obtained from measured data versus WQM results.

Error
C.P. 1
C.P. 2
C.P. 3
C.P. 4

RSME
0.0011
0.0006
0.0005
0.0019

NSE
0.3699
0.6818
0.4036
0.2820

Table 6: Root Mean Square Error (RMSE) and Nash-Sutcliffe efficiency coefficient (see (25)) obtained as quality control of the simulation for Cadmium transport in wet season. Values are obtained from measured data versus WQM results.

Error
C.P. 1
C.P. 2
C.P. 3
C.P. 4

RSME
0.0011
0.0011
0.0017
0.0032

N-S
0.3699
0.8224
0.8230
0.2718

Figure 12: Simulated distribution of Cadmium at the end of the dry season (day 50th snapshot).

Figure 13: Behaviour of Cadmium concentration at control points during dry season. The continuous lines are predicted by the WQM; the (orange) asterisks are field measured data.

Figure 14: Simulated distribution of Cadmium at the end of wet season (day 30 snapshot).

Figure 15: Behaviour of Cadmium concentration at control points during wet season. The continuous lines are predicted by the WQM; the (orange) asterisks are field measured data.

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.

Table 7: Root Mean Square Error (RMSE) and Nash-Sutcliffe efficiency coefficient (see (25)) obtained as quality control of the simulation for Chromium transport in dry season. Values are obtained from measured data versus WQM results.

Error
C.P. 1
C.P. 2
C.P. 3
C.P. 4

RSME
0.0026
0.0008
0.0023
0.0103

N-S
0.4121
0.7291
0.9647
0.3132

Table 8: Root Mean Square Error (RMSE) and Nash-Sutcliffe efficiency coefficient (see (25)) obtained as quality control of the simulation for Chromium transport in wet season. Values are obtained from measured data versus WQM results.

Error
C.P. 1
C.P. 2
C.P. 3
C.P. 4

RSME
0.0107
0.0116
0.0031
0.0022

N-S
0.6948
0.6378
0.2016
0.2026

Figure 16: Simulated distribution of Chromium at the end of dry season (day 50 snapshot).

Figure 17: Behaviour of Chromium concentration at control points during dry season. The continuous lines are predicted by the WQM; the (orange) asterisks are field measured data.

Figure 18: Simulated distribution of Chromium at the end of wet season (day 30 snapshot).

Figure 19: Behaviour of Chromium concentration at control points during wet season. The continuous lines are predicted by the WQM; the (orange) asterisks are field measured data.

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.

Table 9: Root Mean Square Error (RMSE) and Nash-Sutcliffe efficiency coefficient (see (25)) obtained as quality control of the simulation for Lead transport in dry season. Values are obtained from measured data versus WQM results.

Error
C.P. 1
C.P. 2
C.P. 3
C.P. 4

RSME
0.0058
0.0031
0.0551
0.0335

N-S
0.4649
0.4906
0.3573
0.2003

Table 10: Root Mean Square Error (RMSE) and Nash-Sutcliffe efficiency coefficient (see (25)) obtained as quality control of the simulation for Lead transport in wet season. Values are obtained from measured data versus WQM results.

Error
C.P. 1
C.P. 2
C.P. 3
C.P. 4

RSME
0.0020
0.0035
0.0089
0.0198

N-S
0.9320
0.8155
0.0163
0.0984

Figure 20: Simulated distribution of Lead at the end of dry season (day 50 snapshot).

Figure 21: Behaviour of Lead concentration at control points during dry season. The continuous lines are predicted by the WQM; the (orange) asterisks are field measured data.

Figure 22: Simulated distribution of Lead at the end of wet season (day 30 snapshot).

Figure 23: Behaviour of Lead concentration at control points during wet season. The continuous lines are predicted by the WQM; the (orange) asterisks are field measured data.

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.

Table 11: Root Mean Square Error (RMSE) and Nash-Sutcliffe efficiency coefficient (see (25)) obtained as quality control of the simulation for Nickel transport in dry season. Values are obtained from measured data versus WQM results.

Error
C.P. 1
C.P. 2
C.P. 3
C.P. 4

RSME
0.0024
0.0015
0.0155
0.0156

N-S
0.3289
0.3101
0.4473
0.0619

Table 12: Root Mean Square Error (RMSE) and Nash-Sutcliffe efficiency coefficient (see (25)) obtained as quality control of the simulation for Nickel transport in wet season. Values are obtained from measured data versus WQM results.

Error
C.P. 1
C.P. 2
C.P. 3
C.P. 4

RSME
0.0024
0.0015
0.0155
0.0155

N-S
0.3289
0.3101
0.4473
0.0619

Figure 24: Simulated distribution of Nickel at the end of dry season (day 50 snapshot).

Figure 25: Behaviour of Nickel concentration at control points during dry season. The continuous lines are predicted by the WQM; the (red) asterisks are field measured data.

Figure 26: Simulated distribution of Nickel at the end of wet season (day 30 snapshot).

Figure 27: Behaviour of Nickel concentration at control points during wet season. The continuous lines are predicted by the WQM; the (orange) asterisks are field measured data.

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.

References

P. Kavcar, A. Sofuoglu, and S. C. Sofuoglu, “A health risk assessment for exposure to trace metals via drinking water ingestion pathway,” International Journal of Hygiene and Environmental Health, vol. 212, no. 2, pp. 216–227, 2009.

M. K. Mahato, P. K. Singh, A. K. Tiwari, and A. K. Singh, “Risk assessment due to intake of metals in groundwater of east bokaro coalfield, Jharkhand, India,” Exposure and Health, vol. 8, no. 2, pp. 265–275, 2016.

S. Cao, X. Duan, X. Zhao et al., “Health risks of children's cumulative and aggregative exposure to metals and metalloids in a typical urban environment in China,” Chemosphere, vol. 147, pp. 404–411, 2016.

G. M. Gadd and A. J. Griffiths, “Microorganisms and heavy metal toxicity,” Microbial Ecology, vol. 4, no. 4, pp. 303–317, 1977.

G. M. Gadd, “Heavy metal accumulation by bacteria and other microorganisms,” Experientia, vol. 46, no. 8, pp. 834–840, 1990.

L. Zhang, Y.-W. Qin, Y.-Q. Ma, Y.-M. Zhao, and Y. Shi, “Spatial distribution and pollution assessment of heavy metals in the tidal reach and its adjacent sea estuary of daliaohe area, China,” Huanjing Kexue, vol. 35, no. 9, pp. 3336–3345, 2014.

A. R. Jafarabadi, A. R. Bakhtiyari, A. S. Toosi, and C. Jadot, “Spatial distribution, ecological and health risk assessment of heavy metals in marine surface sediments and coastal seawaters of fringing coral reefs of the Persian Gulf, Iran,” Chemosphere, vol. 185, pp. 1090–1111, 2017.

M. Klavinš, A. Briede, V. Rodinov, I. Kokorite, E. Parele, and I. Klavina, “Heavy metals in rivers of Latvia,” Science of the Total Environment, vol. 262, no. 1-2, pp. 175–184, 2000.

A. Tessier, P. G. C. Campbell, and M. Bisson, “Evaluation of the APDC-MIBK extraction method for the atomic absorption analysis of trace metals in river water,” International Journal of Environmental Analytical Chemistry, vol. 7, no. 1, pp. 41–54, 1979.

R. V. Thomann and J. A. Mueller, Principles of Surface Water Quality Modeling and Control, Harper & Row, New York, NY, USA, 1st edition, 1987.

H. Shimazu, E. Ohnishi, N. Ozaki, T. Fukushima, and O. Nakasugi, “A model for predicting sediment-water partition of toxic chemicals in aquatic environments,” Water Science and Technology, vol. 46, no. 11-12, pp. 437–442, 2002.

A. J. Manning, W. J. Langston, and P. J. C. Jonas, “A review of sediment dynamics in the Severn Estuary: Influence of flocculation,” Marine Pollution Bulletin, vol. 61, no. 1, pp. 37–51, 2010.

L. H. Erikson, S. A. Wright, E. Elias, D. M. Hanes, D. H. Schoellhamer, and J. Largier, “The use of modeling and suspended sediment concentration measurements for quantifying net suspended sediment transport through a large tidally dominated inlet,” Marine Geology, vol. 345, pp. 96–112, 2013.

B. S. Caruso, “Simulation of metals total maximum daily loads and remediation in a mining-impacted stream,” Journal of Environmental Engineering, vol. 131, no. 5, pp. 777–789, 2005.

J. C. Winterwerp, “On the flocculation and settling velocity of estuarine mud,” Continental Shelf Research, vol. 22, no. 9, pp. 1339–1360, 2002.

E. R. Sholkovitz, “The flocculation of dissolved Fe, Mn, Al, Cu, Ni, Co and Cd during estuarine mixing,” Earth and Planetary Science Letters, vol. 41, no. 1, pp. 77–86, 1978.

M. A. Helali, W. Oueslati, N. Zaaboub, A. Added, and L. Aleya, “Chemical speciation of Fe, Mn, Pb, Zn, Cd, Cu, Co, Ni and Cr in the suspended particulate matter off the Mejerda River Delta (Gulf of Tunis, Tunisia),” Journal of African Earth Sciences, vol. 118, pp. 35–44, 2016.

N. Tribovillard, T. J. Algeo, T. Lyons, and A. Riboulleau, “Trace metals as paleoredox and paleoproductivity proxies: an update,” Chemical Geology, vol. 232, no. 1-2, pp. 12–32, 2006.

A. R. Karbassi and S. Nadjafpour, “Flocculation of dissolved Pb, Cu, Zn and Mn during estuarine mixing of river water with the Caspian Sea,” Environmental Pollution, vol. 93, no. 3, pp. 257–260, 1996.

I. D. James, “Modelling pollution dispersion, the ecosystem and water quality in coastal waters: a review,” Environmental Modeling and Software, vol. 17, no. 4, pp. 363–385, 2002.

J. Sun, Z. Wang, R. Li, D. Ma, and B. Liu, “Simulation of metal contents through correlated optimal monitoring metals of dagu river sediments,” Procedia Environmental Sciences, vol. 8, pp. 161–165, 2011.

S. Pintilie, L. Branza, C. Betianu, L. V. Pavel, F. Ungureanu, and M. Gavrilescu, “Modelling and simulation of heavy metals transport in water and sediments,” Environmental Engineering and Management Journal (EEMJ), vol. 6, no. 2, 2007.

T. T. H. Nguyen, W. Zhang, Z. Li et al., “Assessment of heavy metal pollution in Red River surface sediments, Vietnam,” Marine Pollution Bulletin, vol. 113, no. 1-2, pp. 513–519, 2016.

Y. Zhang, C. Chu, T. Li, S. Xu, L. Liu, and M. Ju, “A water quality management strategy for regionally protected water through health risk assessment and spatial distribution of heavy metal pollution in 3 marine reserves,” Science of the Total Environment, vol. 599, pp. 721–731, 2017.

N. R. W. Blessing and S. Benedict, “Computing principles in formulating water quality informatics and indexing methods: An ample review,” Journal of Computational and Theoretical Nanoscience, vol. 14, no. 4, pp. 1671–1681, 2017.

K. W. Chau, “Selection and calibration of numerical modeling in flow and water quality,” Environmental Modeling & Assessment, vol. 9, no. 3, pp. 169–178, 2004.

C. D. W. Canham, J. Cole, and W. K. Lauenroth, Models in Ecosystem Science, Princeton University Press, 2003.

R. V. Thomann, “The future "golden age" of predictive models for surface water quality and ecosystem management,” Journal of Environmental Engineering, vol. 124, no. 2, pp. 94–103, 1998.

R. Lun, K. Lee, L. De Marco, C. Nalewajko, and D. Mackay, “A model of the fate of polycyclic aromatic hydrocarbons in the Saguenay Fjord, Canada,” Environmental Toxicology and Chemistry, vol. 17, no. 2, pp. 333–341, 1998.

Z.-G. Ji, J. H. Hamrick, and J. Pagenkopf, “Sediment and metals modeling in shallow river,” Journal of Environmental Engineering, vol. 128, no. 2, pp. 105–119, 2002.

S. P. Bhavsar, M. L. Diamond, L. J. Evans, N. Gandhi, J. Nilsen, and P. Antunes, “Development of a coupled metal speciation-fate model for surface aquatic systems,” Environmental Toxicology and Chemistry, vol. 23, no. 6, pp. 1376–1385, 2004.

C. Kim, H. Lim, and C. Cerco, “Three-dimensional water quality modelling for tidal lake and coastal waters with ROMS-ICM,” Journal of Coastal Research, vol. SI, no. 64, pp. 1068–1072, 2011.

L. Cea, M. Bermúdez, J. Puertas et al., “IberWQ: New simulation tool for 2D water quality modelling in rivers and shallow estuaries,” Journal of Hydroinformatics, vol. 18, no. 5, pp. 816–830, 2016.

L. Cea, M. Bermúdez, and J. Puertas, “Uncertainty and sensitivity analysis of a depth-averaged water quality model for evaluation of Escherichia Coli concentration in shallow estuaries,” Environmental Modelling & Software, vol. 26, no. 12, pp. 1526–1539, 2011.

D. Sharma and A. Kansal, “Assessment of river quality models: a review,” Reviews in Environmental Science and Bio/Technology, vol. 12, no. 3, pp. 285–311, 2013.

Q. Wang, S. Li, P. Jia, C. Qi, and F. Ding, “A review of surface water quality models,” The Scientific World Journal, vol. 2013, Article ID 231768, 7 pages, 2013.

F. Braunschweig, P. C. Leitao, L. Fernandes, P. Pina, and R. J. J. Neves, “The object-oriented design of the integrated water modelling system MOHID,” Developments in Water Science, vol. 55, no. 2, pp. 1079–1090, 2004.

Q. Y. Zhang, “Comparison of two three-dimensional hydrodynamic modeling systems for coastal tidal motion,” Ocean Engineering, vol. 33, no. 2, pp. 137–151, 2006.

M. Grifoll, A. Fontán, L. Ferrer, J. Mader, M. González, and M. Espino, “3D hydrodynamic characterisation of a meso-tidal harbour: The case of Bilbao (northern Spain),” Coastal Engineering Journal, vol. 56, no. 9, pp. 907–918, 2009.

W. Huang, “Hydrodynamic modeling of flushing time in a small estuary of North Bay, Florida, USA,” Estuarine, Coastal and Shelf Science, vol. 74, no. 4, pp. 722–731, 2007.

D. Marinov, A. Norro, and J.-M. Zaldívar, “Application of COHERENS model for hydrodynamic investigation of Sacca di Goro coastal lagoon (Italian Adriatic Sea shore),” Ecological Modelling, vol. 193, no. 1-2, pp. 52–68, 2006.

B. M. Chapman, “Numerical simulation of the transport and speciation of nonconservative chemical reactants in rivers,” Water Resources Research, vol. 18, no. 1, pp. 155–167, 1982.

I. Herrera-Díaz, C. Rodríguez-Cuevas, C. Couder-Castañeda, and J. Gasca-Tirado, “Numerical hydrodynamic-hydrological modeling in flood zones containing infrastructure,” Tecnologia y Ciencias del Agua, vol. 6, pp. 139–152, 2015.

H. Barrios-Piña, H. Ramírez-León, C. Rodríguez-Cuevas, and C. Couder-Castañeda, “Multilayer numerical modeling of flows through vegetation using a mixing-length turbulence model,” Water (Switzerland), vol. 6, no. 7, pp. 2084–2103, 2014.

H. Ramírez-León, C. Couder-Castañeda, I. E. Herrera-Díaz, and H. A. Barrios-Piña, “Numerical modeling of the thermal discharge of the Laguna Verde power station,” Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería, vol. 29, no. 2, pp. 114–121, 2013.

H. Ramírez-León, H. Barrios-Piña, F. Torres-Bejarano, A. Cuevas-Otero, and C. Rodríguez-Cuevas, “Numerical modelling of the laguna verde nuclear power station thermal plume discharge to the sea,” Communications in Computer and Information Science, vol. 595, pp. 495–507, 2016.

H. León, H. A. Piña, C. Cuevas, and C. Castaeda, “Baroclinic mathematical modeling of fresh water plumes in the interaction river-sea,” International Journal of Numerical Analysis & Modeling, vol. 2, pp. 1–14, 2005.

C. Rodriguez-Cuevas, C. Couder-Castañeda, E. Flores-Mendez, I. E. Herrera-Díaz, and R. Cisneros-Almazan, “Modelling shallow water wakes using a hybrid turbulence model,” Journal of Applied Mathematics, vol. 2014, Article ID 714031, 10 pages, 2014.

R. B. Ambrose Jr., T. A. Wool, and T. O. Barnwell Jr., “Development of water quality modeling in the United States,” Environmental Engineering Research, vol. 14, no. 4, pp. 200–210, 2009.

P. R. Kannel, S. R. Kanel, S. Lee, Y.-S. Lee, and T. Y. Gan, “A review of public domain water quality models for simulating dissolved oxygen in rivers and streams,” Environmental Modeling & Assessment, vol. 16, no. 2, pp. 183–204, 2011.

P. G. Campbell and J. Gailer, “Effects of non-essential metal releases on the environment and human health,” in Metal Sustainability: Global Challenges, Consequences, and Prospects, vol. 221, John Wiley & Sons, 2016.

F. Villanueva and A. V. Botello, “Vigilancia y presencia de metales tóxicos en la laguna el yucateco, tabasco, méxico,” Golfo de México Contaminación e Impacto Ambiental: Diagnostico y Tendencias, pp. 407–430, 2005.

F. E. C. Casanova, “Efecto de la contaminacin por metales pesados en los ecosistemas costeros del sureste de mxico,” Kuxulkab, vol. 19, no. 37, pp. 65–68, 2013.

C. Couder-Castaeda, L. H. Ramírez, and D. H. E. Israel, “Numeric optimization of the hydrodynamic model yaxum/3d,” in Numerical modeling of coupled phenomena in science and engineering, M. S. Arriaga, J. Bundschuh, and F. Dominguez-Mota, Eds., chapter 34, Taylor & Francis Group, 1st edition, 2009.

H. A. Barrios-Piña, H. Ramírez-León, and F. Torres-Bejarano, Dinmica de la interaccin ro-laguna-mar. el caso de la laguna el yucateco, tabasco, mxico, chapter 19, Instituto Politcnico Nacional, 2012.

A. V. Botello, “Monitoreo ambiental integral de los impactos de la actividad petrolera en la laguna El Yucateco, Tabasco, Mxico,” Segundo Informe tcnico 12, ICMyL, UNAM, México, 2004.

A. Yáñez-Arancibia and J. W. Day, “The gulf of mexico: towards an integration of coastal management with large marine ecosystem management,” Ocean & Coastal Management, vol. 47, no. 11, pp. 537–563, 2004.

Servicio Meteorolgico Nacional (SMN), Clicom Daily Climate Data, Government Mexican Agency National Metereological Service (SMN), 2017.

V. Casulli and R. T. Cheng, “Semi-implicit finite difference methods for three dimentional shallow water flow,” International Journal for Numerical Methods in Fluids, vol. 15, pp. 629–648, 1992.

H. R. León, C. R. Cuevas, and E. H. Daz, “Multilayer hydrodynamic models and their application to sediment transport in estuaries,” in Current Trends in High Performance Computing and Its Applications, pp. 59–70, Springer, 2005.

A. Eaton, M. Franson, A. P. H. Association, A. W. W. Association, and W. E. Federation, Standard Methods for the Examination of Water and Wastewater, Standard Methods for the Examination of Water and Wastewater, American Public Health Association, 2005, https://books.google.com.mx/books?id=buTn1rmfSI4C.

C. D. H. Congreso de la Unin, “DOF (Diario Oficial de la Federacin),” Federal Rights Law of Mexico, December 2005.

Canadian Council of Ministers of the Environment, “Canadian sediment quality guidelines for the protection of aquatic life,” in Canadian Environmental Quality Guidelines, vol. 1299, Canadian Council of Ministers of the Environment, 2002.

R. Thomann and H. Salas, “Modelos del destino de sustancias txicas,” in Manual de evaluacin y manejo de sustancias txicas en aguas superficiales, CEPIS-OMS, 1988.

D. Mackay, Environmental and Laboratory Rates of Volatilization of Toxic Chemicals from Water, vol. 1, Hazard Assessment of Chemicals: Current Developments, 1981.

D. M. Di Toro, D. J. O'Connor, R. V. Thomann, and J. P. St John, “Analysis of fate of chemicals in receiving waters, phase 1,” in Analysis of fate of chemicals in receiving waters, p. 260, HydroQual, Inc, Mahwah, NJ, USA, 1st edition, 1981.

F. T. Manheim, “The diffusion of ions in unconsolidated sediments,” Earth and Planetary Science Letters, vol. 9, no. 4, pp. 307–309, 1970.

D. O'Connor, ““Modeling frameworks, toxic substances notes,” in Manhattan College Summer Institute in Water Pollution Control, Manhattan College, Bronx, NY, USA, 1st edition, 1985.

H. E. Gaudette, W. R. Flight, L. Toner, and D. W. Folger, “An inexpensive titration method for the determination of organic carbon in recent sediments,” Journal of Sedimentary Research, vol. 44, no. 1, pp. 249–253, 1974.

ASTM, “Standard test methods for moisture, ash, and organic matter or peat and other organic soils,” in Annual Book of ASTM Standard, vol. 2974, pp. 31–33, American Society for Testing and Materials, West Conshohocken, PA, USA, 2000.

ASTM, “Standard test method for particle-size analysis of soils,” in Annual Book of ASTM Standard, vol. 42, pp. 1–8, American Society, for Testing and Materials, West Conshohocken, PA, USA, 1998.

N. Hawley, “Settling velocity distribution of natural aggregates,” Journal of Geophysical Research: Oceans, vol. 87, no. 12, pp. 9489–9498, 1982.

R. V. Thomann and D. M. Di Toro, “Physico-chemical model of toxic substances in the great lakes,” Journal of Great Lakes Research, vol. 9, no. 4, pp. 474–496, 1983.

J. Imberger and G. dl Silvio, “Mixing processes in a shallow lagoon,” in Proceedings of the 23rd Coastal Engineering Research Council of the COPRI (Coasts, Oceans, Ports, Rivers Institute) of the American Society of Civil Engineers, Coastal Engineering 1992, pp. 1867-1868, Venice, Italy, 1993.

L. Cea and M. E. Vázquez-Cendón, “Unstructured finite volume discretisation of bed friction and convective flux in solute transport models linked to the shallow water equations,” Journal of Computational Physics, vol. 231, no. 8, pp. 3317–3339, 2012.

I. K. Tsanis, J. Wu, H. Shen, and C. Valeo, “Environmental hydraulics: hydrodynamic and pollutant transport models of lakes and coastal waters,” in Developments in Water Science, vol. 56, Elsevier, 2007.

M. S. Horritt, “Parameterisation, validation and uncertainty analysis of CFD models of fluvial and flood hydraulics in the natural environment,” Computational Fluid Dynamics: Applications in Environmental Hydraulics, pp. 193–213, 2005.

D. N. Moriasi, J. G. Arnold, M. W. Van Liew, R. L. Bingner, R. D. Harmel, and T. L. Veith, “Model evaluation guidelines for systematic quantification of accuracy in watershed simulations,” Transactions of the ASABE, vol. 50, no. 3, pp. 885–900, 2007.

M. Brocchini, M. Postacchini, C. Lorenzoni et al., “Comparison between the wintertime and summertime dynamics of the Misa River estuary,” Marine Geology, vol. 385, pp. 27–40, 2017.

L. Dong, J. Su, L. Wong, Z. Cao, and J.-C. Chen, “Seasonal variation and dynamics of the pearl river plume,” Continental Shelf Research, vol. 24, no. 16, pp. 1761–1777, 2004.

I. Herrera-Díaz, F. Torres-Bejarano, J. Moreno-Martínez, C. Rodriguez-Cuevas, and C. Couder-Castañeda, “Light particle tracking model for simulating bed sediment transport load in river areas,” Mathematical Problems in Engineering, vol. 2017, Article ID 1679257, 15 pages, 2017.

D. S. Van Maren and P. Hoekstra, “Seasonal variation of hydrodynamics and sediment dynamics in a shallow subtropical estuary: The Ba Lat River, Vietnam,” Estuarine, Coastal and Shelf Science, vol. 60, no. 3, pp. 529–540, 2004.

A. Yáñez-Arancibia and J. W. Day, “Environmental sub-regions in the Gulf of Mexico coastal zone: The ecosystem approach as an integrated management tool,” Ocean & Coastal Management, vol. 47, no. 11, pp. 727–757, 2005.