Groundwater Quality Evolution Model in the Ring of Cenotes, Yucatan, Mexico

Karst aquifers show dissolution/precipitation processes of the minerals present in the carbonate rocks. The Ring of Cenotes (RC) extends along the edge of the Chicxulub crater, in the limestone platform of the Yucatan Peninsula (YP), where groundwater shows preferential flow paths toward the coast near Celestun and Dzilam Bravo towns. This study aimed to describe the regional hydrogeochemical evolution of groundwater of the RC, and its association with the dissolution/precipitation of the minerals present along its path to the ocean. To achieve this aim, we: a) characterized groundwater's hydrogeochemistry; b) determined the calcite, dolomite, and gypsum saturation indexes (reaction phases with the groundwater) in the study area; c) proposed a hydrogeochemical model developed through PHREEQC using an inverse modelling approach. The model predictions confirmed that there are two evolution pathways of the groundwater consistent with the preferential flow paths suggested in a previous regionalization of the RC. On the western path, where groundwater flows towards Celestun, an important marine intrusion influences the hydrogeochemical processes and represents a risk for the prevalence of freshwater. On the eastern path, where groundwater flows toward Dzilam Bravo, the hydrogeochemistry in the sinkholes correlates well with rainfall, suggesting a higher vulnerability during droughts than during rainy periods.


Introduction
Karstic aquifers represent about 20% of the earth's surface, and they are formed by the dissolution of carbonate rocks [1][2][3]. Coastal aquifers are often the main source of freshwater for human coastal settlements. They have a strong interaction with the sea, and their equilibrium depends both on natural hydrological forces (such as the hurricanes) and on anthropic activities (such as groundwater extractions) [4,5].
Carbonate aquifers, such as the one in the Yucatan Peninsula (YP), show different karstification levels depending on the dominant carbonate rocks dissolution process [6]. The YP aquifer shows characteristic karstic manifestations connected with groundwater (sinkholes, locally named cenotes) and does not exhibit surface runoff [7]. Regional groundwater flow paths in the YP are complex and have been proposed to depend on features and interactions at several scales: (a) regional geologic alignments such as the Ring of Cenotes (RC), and the Holbox Fracture Zone; (b) large dissolution conducts at subregional scale, and (c) small-scale fractures caused by weathering and the dissolution of limestone [8].
The RC is located at the northwest of the YP, and it is formed by a high number of sinkholes that are aligned forming a semicircle, along which groundwater flows preferentially [8][9][10]. The RC is the surficial manifestation of the deep crater structure caused by the Chicxulub meteorite impact 66 Ma ago [11][12][13]. The main characteristic of the sinkholes of the RC is that they intercept the phreatic level of the regional aquifer [14]. Sinkholes are of great importance for human settlements because they have been an important source of water in YP since the Mayan era. Nowadays, sinkholes relate to productive activities, such as tourism, on which the income of populations within RC depends [9]. Besides, the discharge of the regional aquifer into the sea is a major component of the hydrologic balance of the coastal lagoons and mangrove ecosystems of Celestun and Dzilam Bravo [15].
The hydrogeochemistry of the RC suggests that groundwater interacts with the carbonate rocks and exhibits a chemical evolution from the central recharge zone to its discharge ends on the coast [10]. Vertically, the coastal aquifer is characterized by a saline interface (halocline) at which the dissolution activity is maximal, and chemical exchange of calcium with other elements such as magnesium coming from intruded saline waters and carbonate rocks like dolomites occurs [16]. Such a phenomenon is more noticeable at the coastal zone of the RC. On the other hand, gypsum dissolves easier than calcite in the groundwater of the RC, especially in the areas close to the Sierrita de Ticul, which is the main physiographic feature of the YP, reaching 350 m a.s.l. [17][18][19].
Earlier studies on the hydrogeochemistry of the water of the RC focused on the presence of pollutants and hazardous substances of anthropic origin [20][21][22][23], as well as on the presence of sulfides, groundwater hardness and water salinization associated to the natural dissolution of calcite and gypsum, and to the interaction with marine water [18,19,21,24]. Other studies are related to the hydrogeochemical evolution of the water quality of the YP at the regional level using an approach based on the chemical balance of major elements concentration and saturation indexes of the lithological phases [18,19,24]. However, there are few studies on the hydrogeochemistry of the RC region that allow the identification of areas susceptible to natural dissolution of these elements and zones with higher saline intrusion by using inverse hydrogeochemical modelling.
The large number of sinkholes distributed along the RC offers the opportunity to study the hydrogeochemical evolution of groundwater along the regional flow paths using hydrogeochemical approaches. This work aims to understand the interaction of groundwater with the main lithological phases and to assess the dissolution and precipitation of minerals along the RC. Through analysing the concentration of major elements in groundwater, building hydrogeochemical diagrams, and using the inverse modelling software PHREEQC. Finally, we also assessed the proportion of seawater present in the groundwater samples of the sinkholes to understand how marine intrusion modifies the dissolution and precipitation of minerals. Understanding of such processes can help to implement management strategies at the RC for a better use of groundwater.

Study Zone
The Ring of Cenotes (RC) is located in Mexico, at the northwest of the YP, between the parallels 88.5-90.5 °W, and latitudes 20.5-21.5 °N ( ). The aquifer in the RC is unconfined, except for a coastal strip where variable hydrogeological confinement occurs [24,25]. The hydraulic gradient is almost flat, and the regional groundwater flow is towards the coast [8,24,25]. The prevailing climate is wettropical, with an annual average temperature of about 26°C, with a minimum in December and January and maximum during May to August [7]. The rainy season occurs between Jun and October, while the period between November and May is drought. The annual rainfall averages 1218 mm for the whole YP, while the potential evapotranspiration index is 1200 -1400 mm. Within the RC, the annual cumulative rainfall varies from 500 mm in the north coast to 1200 mm in the south of the RC [26].
Because of its water reserves and high biodiversity, the RC has been designated a RAMSAR site. Moreover, the RC is part of a state reserve containing 99 sinkholes that constitute a freshwater source for rural human settlements [9]. The RC lacks surficial runoff (rivers and creeks), and the aquifer has a high transmissivity, particularly along the alignment of karstic structures [8,24]. The regional flow paths of the RC discharge into two coastal lagoons: Ria Celestun and Bocas de Dzilam. This discharge helps maintain the health of the mangrove and seagrasses ecosystems that support a high diversity of wildlife of ecological and commercial importance [15].
The surficial geology in the study area (see [18,[27][28][29][30]) shows that the RC extends over four geological formations: (a) the Chichen Itzá formation composed by carbonate rocks from the Eocene, (b) rocks from the Oligocene, (c) the Carrillo Puerto formation composed by limestone from the Miocene-Pliocene, and (d) quaternary Pleistocene-Holocene limestone along the coast ( ). Geochemically, the RC has been regionalized [10] as shown in : (a) Region 1 (R1), located at the west of the RC, where sulfates, sodium and chlorides increase towards the coast, (b) Region 2 (R2), located at the south side of the RC near the Eocene formation, is where the recharge of the RC occurs, and (c) Region 3 (R3) at the east of the RC characterized by an increase of pH and chloride towards the coast.
The surficial geology and the regional hydrogeochemistry [19,24] shows that calcite, gypsum, and dolomite are the most common minerals of the RC, and these solid phases interact with groundwater and with seawater from the marine intrusion [16]. Apparently, groundwater dissolves calcite, while sulfates and chloride increase along its flow to the coast, due to marine influence [19,24]. On the other hand, Velazquez [19] identified a sulfate enrichment pathway in the Eocene aquifer and the influence of dolomites and gypsum. Graniel-Castro et al. [17] and Perry et al. [18] suggest that the Sierrita de Ticul is the main way whereby the dissolved sulfates -originated from gypsum dissolution due to weathering of the Icaiché formation-travel from the Chichancanab lagoon. In the map, the surficial geology described by [28][29][30]. The sinkholes of the study area split into three regions (R1, R2, R3) according to [10].

Field Measurements and Laboratory Analyses
Water samples from 18 sinkholes selected using the regionalization proposed by [10] were collected during the four seasons: May-June 2006 (SP1), October-November 2006 (SP2), February-March 2007 (SP3), and April-May 2008 (SP4). During each field survey, water samples were collected at three depth levels (0.5, 5.5, and 10.5 m) using a Kemereer bottle (WILCO®). However, to avoid bias due to variations in temperature, dissolved oxygen, pH, and suspended solids at the surficial levels, only the samples from the deeper level (10.5 m) were considered. In all cases, this water level was within the freshwater lens, above the halocline. Physical and chemical parameters were measured in situ using a multiparametric probe (Hydrolab G-Quanta®). The concentrations of major elements (Na + , K + , Ca +2 , Mg +2 , Cl -, HCO3 -, SO4 -2 , NO3 -) were determined according to [31]. All water samples had an ionic balance in the range of ±5% to ±10%.

Evolutionary Trajectories of the RC Groundwater
Groundwater ternary hydrogeochemical diagrams were built to analyse the evolution trajectories of the RC [32]. Dispersion diagrams based on the rSO4 -2 /rClratio were plotted as well to identify the spatial patterns of sulfates and salinization in the CR [19,24]. The conceptual model was based on the interactions between groundwater and the lithologic solid phases along the regional groundwater flow path in the RC.

Determination of saturation indexes (SI)
Because of the low velocity of groundwater flow in the water-rock systems, previous studies have shown that groundwater in the Yucatan aquifer is near-equilibrium regarding calcite and dolomite [18,19,24], so that saturation indexes (SI) can be a good proxy of the saturation level of these minerals [33]. Positive SI mean that the solution is oversaturated, and the mineral will show a tendency to precipitate. When SI is negative, the solution is undersaturated, and the mineral will dissolve. If the index is near zero, then the solution is at equilibrium. Because it is uncommon to find solutions exactly at equilibrium, an equilibrium interval was assumed within 0.00 ± 0.05 % deviation from the log of K, the equilibrium constant [32,34]. SI estimation was carried out using the software PHREEQC 3.0, which calculates SI from the molality and the activity coefficient of the dissolved species [35].
Rainfall (H2O) with a high content of carbon dioxide (CO2) forms carbonic acid (H2CO3), which allows the dissolution of the carbonate rocks when it contacts surficial minerals [32], as described by reactions 1-4: Although early studies showed that groundwater in the Yucatan aquifer is near-equilibrium regarding calcite and dolomite [18,19,24], the rate of dissolution/precipitation of these minerals over their trajectory to the coast has not been quantified yet. We describe this approach in the next sections.

Groundwater mineral precipitation/dissolution rate modelling
The estimation of the dissolution and precipitation rates in the sinkholes water was carried out using the inverse modelling module implemented in the PHREEQC 3.0 software [35]. Inverse modelling is based on the initial and final chemical compositions of a solution along a flow line, as well as on the minerals and gases present in the environment, using the mass balance approach. The results obtained are estimates of seawater mixing proportions and potential mass transfers from the lithological solid phases considered (mmol kg -1 H2O) [35].
Using the RC regionalization proposed by [10], two trajectories were analysed using inverse modelling: (a) toward Celestun and (b) toward Dzilam Bravo. The inverse model was implemented in two segments for each trajectory: the first was from the recharge region (R) to the middle point of each trajectory (MC and MD), and the second from the middle point to the respective coastal zone C and D, Fig. 2). 14 were averaged. For MC, concentrations of sinkholes 6, 7 and 9 were averaged. For MD, concentrations from sinkholes 15-16 were used. In the case of the coastal areas, the average of concentrations from sinkholes 1-5 was used for C, and the average of concentrations from sinkholes 19, 20, 21 for D. Calcite, dolomite, and gypsum were considered as the lithological solid phases [18,19,24]. Additionally, marine water solution was included to simulate the marine influence within the aquifer towards the coast [36]. Finally, O2 and CO2 phases were also included in the model as suggested in Eq. 1 ( Fig. 1 y Fig. 2).
For convergence, the model considered a global uncertainty of 0.10. All possible water-rock interactions were considered. The best solution was selected based on the smaller molal balance residual, previously estimated by the simplex method implemented in PHREEQC [35].

Hydrogeochemical Evolution of the Groundwater of the RC
Major ion concentrations in the RC showed two different evolution trajectories toward the coastal areas of Celestun and Dzilam Bravo (Figs. 3 and 4).
In the Celestun trajectory, an increase of the concentration of sulfates segment was observed on R-MC for all the surveys (Fig. 3). This behaviour is coincident with previous findings for the Eocene aquifer near Ticul [19,24]. In the MC-C segment, an increase of chloride concentration was also observed, revealing marine influence and groundwater salinization near Celestun [19].
Figs. 3a-3c show an increasing chlorides concentration towards the coast in the trajectory to Dzilam Bravo (samplings SP1 to SP3), as well as decreasing calcium concentration. Particularly, Fig. 3d  (SId = 0.64), and sample MD from survey SP3 which was oversaturated for both calcite (SIc = 0.67) and dolomite (SId = 1.45). Point C showed chemical equilibrium with calcite (-0.424 < SIc < 0.424), except for the sample from SP2 (SIc = 1.05), while dolomite showed oversaturation in SP2 and SP3 (SId = 2.29 and 0.69, respectively). Point D showed oversaturation with calcite in SP2 and SP4 (SIc = 1.39 and 0.88, respectively) and with dolomite in SP2, SP3, and SP4 (0.56 < SId < 2.87). All samples showed the potential to dissolve gypsum, but it was more notable in the R sample (SIg < -0.229, Fig. 4c). The recharge zone showed chemical balance with calcite and dolomite but towards the coast dolomite may precipitate, while gypsum may be dissolved at any point of the trajectories described.
Marine influence was stronger in Celestun as compared with Dzilam Bravo. Point C showed a seawater proportion between 1.0 and 1.9 % (SP2 and SP3 respectively), while point D showed a percentage range of 0.2 to 0.5 % (SP1 and SP2 respectively Fig. 4). Pérez-Ceballos et al. [10] described such behaviour in their regionalization, and it was also supported by the results observed in SP4.  Finally, the rSO4 -2 /rClratio showed that groundwater hydrogeochemical evolution was consistent with the regional trajectories of the YP's aquifer [19,24]. The trajectory towards Celestun showed sulfate enrichment, likely a product of the gypsum dissolution, while the trajectory towards Dzilam Bravo showed only seawater enrichment (Fig. 5).

Groundwater-Rock Interaction in RC
The analysis of the chemical evolution along the trajectories and the SI explained the regional spatial distribution of the water quality in the RC. However, this approach only considered the potential for dissolution of each of the minerals without considering coupled interactions. Conversely, inverse modelling (fed with the information of the lithological phases) allowed prediction of the dissolution/precipitation rates of these minerals, considering the known concentration in the water samples ( Table 1). The optimal hydrogeochemical model was selected using the criterion of smaller molal residual. This model shows the calcite-dolomite-gypsum interaction additionally to the interaction of oxygen and carbon dioxide as phases of the balance. There were fewer model solutions to explain the evolutionary trajectories from R to MC and MD compared with the number of solutions found from MC and MD to C y D, respectively.
In dissolution activity (these minerals even precipitated), and calcite had a potential dissolution rate greater than the remaining lithological phases. On the R-MC, the precipitation of calcite was consistently indicated on all the surveys (SP1-SP4), showing activity rates between -0.42 and -2.048 mmol kg -1 H2O (Fig. 6a). The potential dissolution activity rate of dolomite and gypsum was in the range of 0.04-1.05 and 1.06-1.99 mmol kg -1 H2O, respectively (Fig. 6b y 6c). Finally, on the segment MC-C the calcite showed average potential dissolution rates of 0.71 mmol kg -1 H2O for the surveys SP1 and SP4, while for surveys SP2 and SP3, the average precipitation was 0.23 mmol kg -1 H2O (Fig. 6a). The potential dissolution activity rate of dolomite showed a decreasing trend (0.07 y 0.14 mmol kg -1 H2O, surveys SP2 and SP3), smaller than those observed on the segment R-MC, and we even observed potential precipitation of this mineral in surveys SP1 and SP4 (activity rates between -0.12 and -0.14 mmol kg -1 H2O, Fig. 6b). Similarly, gypsum showed a reduced potential dissolution activity on the segment MC-C, contrasting with the segment R-MC. Gypsum potential dissolution rates ranged between 0.38 to 0.71 mmol kg -1 H2O (surveys SP2 to SP4), and in survey SP1 precipitated at the activity rate of -0.6 mmol kg -1 H2O (Fig. 6c).
On the Dzilam Bravo trajectory, the intermediate segment (R-MD) showed small potential activity and precipitation of gypsum and dolomite. The calcite precipitation showed potential activity between -0.8 and -1.5 mmol kg -1 H2O, except during SP2, which showed a potential dissolution rate equal to 0.056 mmol kg -1 H2O (Fig. 6a). The activity rate of dolomite was in the range of 0.2-1.46 mmol kg -1 H2O and less activity of gypsum (precipitation rate between -0.006 and -0.02 mmol kg -1 H2O, Fig. 6b). During SP4, the dissolution rate of gypsum was 1.16 mmol kg -1 H2O (Fig. 6c). Toward the coast (MD-D), enrichment with calcite was revealed, while the dolomite and gypsum showed less activity, except during SP4 where we observed enrichment and precipitation of these minerals. On this segment, the calcite dissolution activity increased (0.04 to 2.33 mmol kg -1 H2O, Fig. 6a), while the dolomite precipitation activity was between -0.21 and -1.97 mmol kg -1 H2O, and little activity of gypsum (dissolution rate = 0.07-0.18 mmol kg -1 H2O, Fig. 6b). The precipitation rate of gypsum during SP2 was 0.06 mmol kg -1 H2O (Fig. 6c).
There was less calcite activity on the segments R-MC and R-MD (for both trajectories: Celestun and Dzilam Bravo) than on the coastal zone (calcite dissolution increased along MC-C y MD-D). The dissolution of dolomite was consistent on the segments R-MC and R-MD and higher than the activity rates and precipitation towards the coast (MC-C and MD-D). The hydrogeochemical inverse model predicted the Eocene aquifer's influence with inputs of dissolved dolomite in the groundwater, as described by [19]. At the coast, the process of dolomite precipitation (dolomitization) was due mainly to chemical exchange with cations such as Mg +2 available in the ground and from mixing of the regional aquifer with marine water. Such a process has been previously described for the coastal area of YP by [16,37,38].
Regarding the trajectory of gypsum towards Celestun, a constant enrichment of this mineral on the segment R-MC was observed, which matches the trajectories of sulfates dissolution described by [19,24]. On the trajectory to Dzilam Bravo, this mineral showed little activity as compared to the other lithological phases. However, an atypical behaviour in the samples from SP4 survey was found, which showed high gypsum concentrations on the segment R-MD, similar to those in segment R-MC (Fig. 6c).

Discussion
The RC showed two different hydrogeochemical evolutions along the trajectories from the recharge zone towards the coastal zones of Celestun (R-MC-C) and Dzilam Bravo (R-MD-D). These hydrogeochemical groundwater evolutions were suggested by the ternary diagrams and the ratio rSO4 -2 /rCl -, and confirmed through inverse modelling.
Celestun trajectory showed an increase of sulfates associated with the dissolution of the gypsum's Eocene aquifer and towards the coast also likely related to seawater mixing caused by marine intrusion. The dissolution of gypsum in R-MC was indicated in all of our surveys, which is consistent with the results of Perry et al. [18,24,27], who posed that sulfates originate from gypsum dissolution at the Chichancanab lagoon and move through the Ticul fault. The segment R-MC showed chemical equilibrium for calcite and dolomite, and undersaturation for gypsum (Fig. 4). Additionally, gypsum and dolomite dissolution as well as calcite precipitation were observed. This combination of processes is also known as dedolomitization [39], and is closely related to coupled variations of Mg +2 and Ca +2 concentrations in the groundwater of the YP's aquifer [23]. In the coastal segment of this trajectory (MC-C), gypsum and dolomite showed a decreasing activity (and even precipitation during SP1, Fig. 6b y 6c); but also calcite dissolution and marine water influence ( Fig. 4d y 6a). Conceptual models [16,37] also predict dolomite precipitation (dolomitization) in the presence of brackish water through lattice replacement with Mg +2 from seawater [40]. Overall, the inverse modelling implemented in this work is in agreement with earlier conceptual models proposed by [16,18,19,23,24,37,39] and supports them.
Dzilam Bravo trajectory also showed salinization toward the coast, although less intense than in the Celestun trajectory (Fig. 4d). Calcite was the lithological phase showing higher activity along the Dzilam trajectory, contrasting with gypsum and dolomite. While calcite apparently precipitated on the R-MD segment, later it dissolved on the MD-D segment (Fig. 6a), which was likely favoured by the higher influence of seawater towards the coast [16]. On the other hand, gypsum activity, which was overall lower in this trajectory, increased to levels as high as those found along the Celestun trajectory ( Fig. 4 y 6) during SP4, when rainfall was lowest at Dzilam Bravo (Fig 7). This showed a pattern -consistent with [18,19,24,27], who suggested that gypsum dissolved by weathering of the older formations (Eocene and Paleocene) travels toward the RC-which is apparently blurred by the stronger rainfall along this trajectory. The importance of rainfall for Dzilam Bravo is further outlined by the significant and strong inverse correlation found for seawater mixing with rainfall in this trajectory (Fig.  8b), showing the effect of rainfall variations on the marine intrusion there. In contrast, although seawater mixing was higher at Celestun, it did not correlate significantly with rainfall. Similarly, the SI of gypsum did not correlate significantly for Celestun trajectory, but it did with rainfall for the Dzilam trajectory (Fig. 8, SIg determination coefficient R 2 > 0.5). Calcite dissolution/precipitation correlated directly with rainfall, while inverse correlations were obtained for dolomite and gypsum (Fig 9, R 2 >0.5). This suggests that the aquifer's recharge associated with the rainfall it is an important effect on the hydrogeochemical model for the Dzilam Bravo trajectory, similarly to the findings of [23,21]. However, Celestun showed a higher potential transference rate from the lithological phases as compared with Dzilam Bravo due to higher influence of marine water. Thus, the hydrogeochemistry of Celestun is more closely related to seawater mixing, as opposed to Dzilam Bravo, where rainfall is more important. This higher seawater mixing towards Celestun relative to Dzilam Bravo could be related to coastal differences among these two sites. Along the Yucatán's coast, the confinement effect influences the interaction between marine water and the aquifer, affecting the coastal aquifer's salinity [25,41,42]. Because the coastal confinement effect extends about 20 km towards Celestun but about 3 km toward Dzilam Bravo [24,43], it likely allows a higher interaction between groundwater and seawater in Celestun, which results in the higher precipitation rate of gypsum and dolomite, and the higher seawater mixing as compared to Dzilam Bravo. Therefore, saline intrusion on the Celestun trajectory could have a significant effect on the RC water quality, adding up to the impacts caused by increasing human population, tourism and industrial development, as well as sea rise level driven by climate change [44,45]. Because of this higher vulnerability of the RC groundwater on the Celestun trajectory, we here highlight the importance of establishing a permanent salinization monitoring program on the Celestun trajectory, with higher priority during the winter months, when atypical meteorological tides are caused by cold fronts [10,42]. In this regard, chlorides should be monitored on segment MC-C (about 70 km from the shoreline), while sulfates are of particular concern in segment MC-R (cenotes 6,7,9). We also propose monitoring sulfates during the dry season for Dzilam Bravo because of the increase of sulfates activity during lower rainfall periods. We recommend monitoring along the MD-D segment (cenotes 17 and 19) located at 20 km from the shoreline, and in R-MD (cenotes 10 to 16) located 30 to 62 km from the Dzilam Bravo's shoreline.
The high concentrations of sulfate we found in the recharge zone (R) of the RC have a natural origin in the dissolution of gypsum from the Eocene formations. However, earlier studies reported high concentrations in wells drilled near Ticul, and that agricultural practices in the area likely also contribute to the increasing sulfates concentration [17,23], which could in short time become unacceptable for human populations [46]. For example, it is known that sulfates may have a laxative effect when combined with Ca +2 and Mg +2 [46,47], and both ions are also present in the RC. Therefore, the cenotes at the recharge zone must also be monitored, especially during the dry season, when gypsum interaction is higher (Fig. 8 y 9).
Finally, we would like to briefly outline the advantages of the inverse model here used over the SI analyses derived from its capability to analyse the lithological phases in a coupled way, producing more detailed predictions. For example, the SIg suggested gypsum dissolution all along the Celestun trajectory (Fig. 4) but the inverse model predicted a decrease of gypsum's activity and even its precipitation at the MC-C segment. Additionally, the inverse model allowed identifying the occurrence of dolomitization at R-MC. Altogether, our results sustain the usefulness of this approach for assessing the evolution of groundwater´s hydrogeochemistry and quality.

Conclusions
Groundwater in the RC flows along two different evolution trajectories. Both begin at the RC's recharge zone and travel towards the coastal towns of Celestun and Dzilam Bravo.
In the Celestun trajectory, gypsum inputs from the Eocene aquifers were detected, which then precipitate due to mixing with seawater towards the coast. The marine intrusion there influences the dissolution and precipitation processes of the lithological phases showing dedolomization in the middle segment of the trajectory (R-MC) and dolomization in the coastal one (MC-C).
In the trajectory towards Dzilam Bravo, calcite precipitation/dissolution processes were more active than for gypsum and dolomite, and correlated directly with the rainfall. Groundwater with higher gypsum concentration due to the weathering on the Eocene formations was observable only during the lower rainfall sampling.
The higher intensity of marine intrusion (up to 2% of seawater) in the Celestun trajectory than in the Dzilam Bravo trajectory, in possibly due to a differential confinement effect along the coast. Marine intrusion in Celestun trajectory barely changes with rainfall, while in the coastal zone of Dzilam Bravo it clearly variates inversely with rainfall.
Both regional flow trajectories of the RC are vulnerable to climate change effects, such as sea level rise, because hydrogeochemical processes there are affected by the interaction with seawater, but Celestun is likely more susceptible to the rising sea level. In contrast, Dzilam Bravo can be affected by changes in the rainfall patterns because hydrogeochemical processes there depend more on the aquifer's recharge. Monitoring of marine intrusion and phases dissolution in the RC is highly recommended to assess the evolution of the system, and could simply be done through chlorides and sulfates measurements.