Modelling, Transport and Assessment of Aquatic Toxic Metals in Coastal Ecosystems

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 enable 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 the Chicozapote and Tonalá rivers, that connect the lagoon with the Gulf of Mexico. We present the results of pollutants transport simulation in the coastal ecosystem, focusing on toxic parameters for two hydrodynamic scenarios: wet and dry seasons. As it of interest in the zone, we study the transport of four metals: Cadmium, Chromium, Nickel and Lead. In order to address our objectives we solved numerically a self-posed mathematical problem,which is based on the measured data. The performed simulations show to characterise metal transport within the acceptable range of accuracy and in accordance with the measured data. The performed simulations show 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 as an excellent option when the distribution of solutes with high accuracy is required, easing environmental, economic and social management of coastal ecosystems.


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 in human health, it is necessary to study heavy metal contamination in aquatic ecosystems [4].
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 is an indicator of anthropogenic activities, such as hazardous wastes derived from industrial activities, mining, agriculture, etc. Nowadays, interest in the pollution of rivers by metals has increased, because rivers serve a medium for transport of dissolved and particulate matter from continents to the ocean. For this reason, heavy metal concentrations in waters have been analysed worldwide [5].
In coastal waters, heavy metals are distributed through the water column 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 flocculation process during estuarine mixing. Thus, heavy metals get bound to these elements and precipitate to the bottom. 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 [6][7][8] (usually characterised by muddy bottoms [9]). 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 [10,11]. Moreover, other metals are known for their nutrient-like behaviour. Thus, flocculation process constitutes an arduous task to model [12], which is not the aim of this work.
For the above reasons, strategies and tools to mitigate the pollution of heavy metals are required [13]. A huge number of mathematical models that intend to predict the transport of heavy metals in flows exist, as 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 [14]. 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 allows to use a minimum of experimental measurements [15].
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,16,17]; 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 [18] (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, has been deeply studied on shallow aquatic systems during this century [19,20]. 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 -reaction equation), considering even the highly nonlinear interactions typical of these regions. Thomann [27] and MIKE Ecolab/ABM. Nevertheless, most of them are based on solving mass balance/advective diffusion equation, but oriented to nutrients (as DO, BOD, PH 4 , Phosphorus) or pathogens (like coliforms). Such models do not contemplate metals transport [28]. A review regarding computational models of water quality can be found in [29], where it is pointed out that 3 of 25 WTM are of great importance for environmental management; nevertheless, they still constitute an open research field. 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 [30], which has been used to model the concentrations of metals in rivers since the decade of 1980.
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-development for research purposes, and has been previously used in different research applications, including: the hydrodynamic-hydrological modelling in flood zones [31], the modelling of flows through vegetation [32], the modelling of the thermal discharges [33,34] and the modelling of fresh water plumes in river-sea interaction [35]. On the other side, in order to accurately simulate turbulent viscosity, we implemented the turbulence model discussed in detail in [36].
The water quality module herein developed 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 [37,38]. 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 [39]. According to our previous knowledge, we did not find an ad hoc module for the transport of these metals, so we 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 o 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, the Chicozapote river and the Tonalá river discharge, which has suffered a serious deterioration due to pollution originated by oil industry activities.
This paper begins by a detailed presentation of the coastal ecosystem we considered in Section 2. It will then go to exposition of the methodology herein applied, including the measurement protocol, the design of the hydrodynamic and water quality modules as well as the numerical strategies for their solution (Section 3). Thereafter, in Section 4 we present 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 to assess the accuracy of the WQM herein developed with respect to measured data. Finally, in Section 5 we present the conclusions of this study as well as some final interesting remarks.

Case of study
The coastal ecosystem we studied is located at the east of Tabasco State, in the southeast of Mexico. The lagoon is located between LAT 18 o 10' and 18 o 12' N, and between LONG 94 o 02' and 94 o 00' W ( Fig. 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.
The predominant climate in the region is warm and humid, with abundant rainfall during summer. Its annual thermal regime oscillates between 25.8 o C and 27.8 o C; the highest average temperature occurs in May with 29.4 o C while the minimum is registered in January with 23.1 o C. September/october is the most rainy month with an average precipitation of 400 mm while in June the minimum incidence with an average value of 43.3 mm.
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.
For dry season, 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 season, 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. Main activities in the zone are agriculture, cattle farming, fishing and industry. 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 [40]. More environmental information about this ecosystem is available in [41].
Furthermore, Oil field Cinco Presidentes was established in 1963 in the vicinity of this region, being El Yucateco lagoon the closest water body to the oil facility. Later, a network of artificial channels (approx. 33 Km) was built by PEMEX 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.

Methodology
The WQM herein developed consists essentially of two parts: an hydrodynamic module and a water quality module. The later one being in charge of transporting heavy metals through the water body by using the hydrodynamic module results. Later, in Sections 3.2 and 3.3 these modules are appropriately defined. In what follows, the protocol regarding measurement and data analysis is presented.

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 [42].
For heavy metals in water, 500 ml of sample were taken at each point using a Van Dorn water sampler. The samples were filtered (0.45 µm) and acidified (HNO 3 at pH = 2.0), stored in amber glass bottles and refrigerated at 4.0 C for transportation [42]. For heavy metals in sediments, samples of 100 g were collected from the surface layer (max. 5 cm) with an Ekman dredger. These samples were stored in sterile polyethylene bags and kept at 4.0 C for transport. All the samples were collected by duplicate to determine the precision of tests and sample handling. The chemical analyses of heavy metals in water and sediments was performed using atomic absorption spectrophotometry, with a spectrometer ICE 3500 AAS Thermo Fisher Scientifc.
According to previous in situ studies [41], several metals were initially sampled. Those whose parameters exceeded water and/or sediment quality criteria [43,44] were selected for diagnosis: Cd, Cr, Pb and Ni.
Measured data and their chemical analyses of these parameters served to specify boundary and initial conditions to the numerical model herein developed. They also aided to validate the model for the study region.

Hydrodynamic module
The hydrodynamic module is based on the governing equations for shallow waters with constant density and free surface, which can be derived from the Reynolds-averaged Navier-Stokes equations [36]. These equations are: where U and V are the mean components of velocity field in the x and y directions respectively (m/s); g is the acceleration due to gravity (m/s 2 ); ρ is the water density (kg/m 3 ); ρ 0 is the water reference density (kg/m 3 ); η is the free surface elevation (m); P atm is the atmospheric pressure (Pa); and ν TH is the horizontal eddy viscosity (m 2 /s). The equation to calculate the free surface elevation is obtained by integrating the continuity equation over the water depth, which by applying a kinematic condition at the free surface, is given by where h(x, y) is the water depth (m).
As mentioned above, in order to achieve a 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 herein treated, the ml-05 mixing-length model is introduced [36], which contributes to the horizontal eddy viscosity as: where l h is horizontal mixing length.

Water Quality Module
The water quality module can be deem as consisting of a part which is in charge of transporting the substance plus a second one which takes into account the reaction mechanism to which the toxicant is subject. For the transport phase, the Advection-Diffusion-Reaction equation is solved:

∂C ∂t
Rate of change Rate of change due to diffusion Rate of change due to reaction sources (5) where C is the concentration of a substance (mg/L), Ex and Ey are the horizontal dispersion coefficients (m 2 /s) [36]; and Γc is the substance reaction term (mg/L).
The water quality module is initialised considering no reaction (Γ C (t = 0) = 0), and assuming stationary sediment, constant kinetic coefficients and suspended solid's uniformly distributed in space over the river reach. Once the initial concentration distribution C is calculated through Eq. 5, it is re-estimated through a model of reaction of toxic substances. Thomann and Salas presents a complete model which considers diffusive exchange, decomposition processes, volatilisation, settling and re-suspension (important processes in a coastal aquatical ecosystem), which is governed by the following ordinary differential equation, where 1 and 2 represent the water column and sediment respectively, C 1 is the concentration of the toxic substance (mg/L) (estimated using the transport equation), C 2 is the concentration of the toxicant in sediments (mg/L), K f is the diffusive exchange (m/d), h is the River depth (m), f d is the dissolved fraction (1), f p is the particulate fraction (1), K d1 is the degradation rate of the dissolved toxic substance (d −1 ), K L is the overall volatilisation transfer rate (m/d), c g is the vapour phase concentration (mg/L), He is Henry's constant, v s is the settling velocity (m/d), v u is the re-suspended velocity (m/d) and ϕ 2 is the sediment porosity (1). For further reference on these terms, see Figure 2. In Eq. (6), 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, etc., in the water column (decay of particulate form is assumed 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 re-suspension into the water column of the particulate toxicant from the sediment.
A description of the processes underlying the parameters in Eq. (6) 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 (K d1 ) are represented by the equation: where K p is the photolysis rate (d −1 ), K H is the hydrolysis rate (d −1 ) and K B is the microbial degradation rate (d −1 ). The overall exchange rate (K L ) estimates the importance of losses due to volatilisation. According with [46], the application of the two film theory yields to an overall volatilisation transfer parameter given by, where k l is the liquid film coefficient (m/d) and k g is the gas film coefficient (m/d). As seen from Eq. (8), K L depends on chemical properties as well as on characteristics of the water body such as water velocity (affecting k l ) and wind velocity over water surface (affecting both k l and k g ). Henry's (H e ) constant, present in Eqs. (6) and (8), takes into account the water-atmosphere interaction, and it represents a partitioning of the toxicant between water and atmospheric phases. Its dimensionless form is given by where c w the water solubility (mg/L). The sediment diffusion rate K f reflects the fact that gradients can occur between interstitial sediments and the adjacent water column. Di Toro et al. has summarised some data respect to this coefficient in the following expression, where MW is the Molecular weight (g/mol) and ϕ 2 the sediment porosity.
To determine both settling and re-suspended 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, etc. were quantified. In coastal ecosystems, a good approximation to the settling velocity v s is given by Hawley, who provides the following empirical relation, v s = 0.2571(d) 1.138 (11) where d is the particle diameter (m).
For the re-suspended velocity v u , Thomann and Di Toro have proposed the next equation: where v d is the net rate loss of solids (m/d) and v n 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 [9]; nevertheless, the water body under study have low circulation because there is not important water interchange with Chicozapote River. Being this 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. Figure 3 schematises the simulation process to obtain the distribution of heavy metals. Hydrodynamics is calculated for both seasons (wet and dry) based on the application of appropriate boundary conditions (see Section 3.5). Once obtained the velocity field (U, V), it constitutes the input for the water quality module, which is executed to validate the measured data. Then, the water quality module is executed iteratively readjusting the parameters, until it matches the measured data within an expected error. Thus, the WQM herein developed is of high computational burden, even though the time increment (∆t) used in the WQM is much greater than the used in the hydrodynamic module.

Boundary conditions
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, that induce pressure and friction forces at the bottom [50], 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 [51].
The modelling of the Sea -Chicozapote River -El Yucateco lagoon system was carried out by using the astronomical tide, atmospheric (wind direction and velocity) forces and river discharges as forcing boundary conditions. This methodology allow us to reproduce the dynamics of the hydrosystem and gives us a better understanding of how the transport processes occur within the water body. 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 area of study (see Figures 4, 5 and 6) [41]. 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. Wind forcing was implemented with data measured directly over El Yucateco lagoon during the field measurements campaigns in both seasons. Wind forcing was considered spatially invariant along the model domain. The hydrological flow of Tonalá river was imposed according to the data shown in Figures 5b and 6b. The points of application of boundary conditions are shown in the arrows in Figure  7.    The bathymetry shown in Figure 7 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 (see Section 3.6).

Numerical solution methods
The numerical solution scheme 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. This method treats the advection and diffusion terms differently. The solution of the advection terms 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. This procedure is applied for the solution of the momentum equations for U and V velocity components (Eq. (1) and (2)) and for the solution of metal's transport equation (Eq. (5)).
The area under study was discretised through a structured mesh, composed of 213 × 79 staggered cells (see Figure 8). 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, time step ∆t is selected to guarantee stability and convergence of the numerical solution, while minimising the computational burden (see Section 4.

Quality of numerical solution checkup
In order to assess the quality of the numerical solution with respect to field measured data, we present the root square mean error (RSME) as well as the Nash-Sutcliffle efficiency criteria, the last given by [53], where ϕ obs is an observed/measured datum, ϕ sim is the respective datum obtained from numerical simulation (same place and time). The ϕ parameter is the measured data variance. Clearly, if R is close to 1, it is possible to consider that numerical simulation data is quite approximate to field measured data.

Field measurement campaigns
Tabasco State is the region of Mexico where the climate is much 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 [54]; for instance, they have influence on parameters that define the transport of solutes such as salinity [55] or sediments [26,56,57]. For the purpose of yielding to results which can be compared with stationary measurements in situ, we selected wet and dry seasons for this study, as they are detailed in Table 1. 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 [58]. 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 Fig. 1). The location of control points depicted in Figure 1 is shown in Table 2. 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. In order to compare the measured data with the numerical solution of the water quality module, a fitting of the data for each Control Point, heavy metal and season was best performed by a quartic polynomial model. Table 2. Considered control points (depicted in Figure 1).

Validation of WQM
The validation process of the WQM consists in finding the values involved in the considered reaction mechanism (see Section 3.3 and particularly, Eq. (6)) that 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 was carried out forcing with measured concentrations in Control Points C.P.1, C.P.2, C.P.3 of Figure 1. After convergence, calculated concentrations in Control Point 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 to reach an error lower than 5%, by adjusting the involved parameters in each simulation. In a second step, a new numerical simulation was carried out forcing with measured concentrations in Control Points C.P.1, C.P.2, C.P.4. After convergence, calculated concentrations in Control Point C.P.3 were compared against measurements in the same location. The trial and error procedure was performed up to reach an error lower than 5% by adjusting again the referred parameters. The same procedure with 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.

Hydrodynamic simulations
The hydrodynamic simulations were performed using a time step of ∆t = 2.0 s, requiring 2,160,000 and 1,296,000 iterations for dry and wet seasons respectively. Figure 9a and Figure 9b show 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 re-circulation 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 re-circulation. Figures 10 and 12 show the final snapshots of the Cadmium transport simulation for dry and wet seasons respectively. Figures 11 and 13 depicts the behaviour of Cadmium at viewers (control points, see Figure 1). Table 4 and 5 show RSME and Nash-Sutcliffe efficiency coefficient (Eq. (13)) 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.   Figure 13. Behaviour of Cadmium concentration C at Control Points during wet season. The continuous (blue) line is predicted by the WQM, the (red) asterisks are field measured data and the dotted (orange) line fits the measured data. Table 5. Root Mean Square Error (RMSE) and Nash-Sutcliffe efficiency coefficient (Eq. 13) obtained as quality control of the simulation for Cadmium transport in wet season.      Table  8 and 9 show RSME and Nash-Sutcliffe efficiency coefficient (Eq. (13)) 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 10 and 11 show RSME and Nash-Sutcliffe efficiency coefficient (Eq. (13)) 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.

Conclusions and final remarks
In this paper we develop and test a hydrodynamics-based WQM for the specific evaluation of heavy metals. 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 agree 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.
We ought to remark the importance of possessing enough information to feed the numerical model with the necessary data in order to validate, to calibrate 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 agree with field measured data.
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 calibrate 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.
The numerical simulations were carried out 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 domain, 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 is 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. 1 and 2), due to the small currents that the Lagoon possess.
The numerical simulations herein performed allow to characterise the evolution of metals in a specific water system, which is a complex and expensive task to achieve only by field measurements. The application of Computational Fluid Dynamics techniques to real-world environmental problems has increased due to the computer capacities (power and storage), allowing to simulate flows over complex topographies including many parameters. In this way, numerical modelling constitutes an excellent option when a distribution of the solutes with a high accuracy is required. The most challenging fact on carrying out metal transport numerical simulations is the lack of data to validate as well as the lack of information about the values of the reaction equation terms coefficients.