Regional Hydrogeochemical Evolution of Groundwater in the Ring of Cenotes, Yucat á n (Mexico): An Inverse Modelling Approach

: 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 two preferential ﬂow paths toward the coast near Celestun and Dzilam Bravo towns. The objectives of this study were to describe the regional hydrogeochemical evolution of the groundwater in the RC, and its association with the dissolution/precipitation of the minerals present along its pathway to the ocean. These objectives results were obtained by: (a) characterizing groundwater hydrogeochemistry; (b) calculating calcite, dolomite, and gypsum saturation indexes in the study area; and (c) developing a hydrogeochemical model using PHREEQC (U. S. Geological Survey) inverse modelling approach. The model predictions conﬁrmed that there are two evolution pathways of the groundwater consistent with the preferential ﬂow paths suggested in a previous regionalization of the RC. On the western path, where groundwater ﬂows towards Celestun, marine intrusion inﬂuences the hydrogeochemical processes and represents a risk for the freshwater. On the eastern path, where groundwater ﬂows toward Dzilam Bravo, rainfall has an important effect on the hydrogeochemical processes, evidenced by a higher concentration in sulfates during droughts than during rainy periods. Then, monitoring of marine intrusion and phases dissolution in the RC is highly recommended.


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 on both natural hydrological forces (such as the hurricanes) and 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]. In the YP, regional groundwater flow pathways are complex and depend on karstic manifestations at several scales: (a) regional scale due to geologic alignments such as the Ring of Cenotes (RC), and the Holbox Fracture Zone; (b) large dissolution conduits at subregional scale such as the Aktun-Chen cave in Quintana Roo, and (c) local-scale heterogeneities caused by erosion, tectonics, and the dissolution of limestone [8]. Therefore, RC is an important karstic zone which influences the regional groundwater flow path in YP.
The RC is located at the northwest of the YP (Figure 1), and it is formed by several sinkholes that are aligned forming a semicircle, along which groundwater preferentially flows [8][9][10]. The RC is the surficial manifestation of the deep crater structure caused by the Chicxulub meteorite impact 65 Ma (million years ago) [11][12][13]. Mechanisms such as erosion, tectonics and dissolution of limestones created a morphologic aspect in the Miocene-Pliocene zone that shows a more karstified and inclined setting from Ticul to the RC. In contrast, in the interior of the Yucatan's plain (from RC to the coast), the topography is smoother and shows a less karstified landscape than the former [12]. 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 are a source of water in the YP since the ancient Mayan civilization. Nowadays, sinkholes also relate to productive activities, such as tourism, on which the income of some RC's inhabitants depends [9]. Additionally, 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 quality evolves from the central recharge zone to its final destinations on the coast [10]. Vertically, the study zone is characterized by a saline interface (halocline) located between 50 to 60 m in the recharge zone of the RC, and between 18 to 26 m in the coastal zone ( Figure 2). Halocline shows the maximal dissolution activity and chemical exchange of calcium with other elements such as magnesium coming from intruded saline waters and carbonate rocks [16]. Such a phenomenon is more noticeable at the coastal zone of the RC. On the other hand, gypsum dissolves more easily than calcite in the groundwater, especially in the areas close to Sierrita de Ticul, which is the main physiographic feature of the YP, reaching 350 m above sea level (m. a.s.l.) [17][18][19].
Earlier studies on the hydrogeochemistry 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 with 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 at the regional scale using a chemical-based approach on the chemical balance of major elements' concentration and the saturation indexes of the lithological phases [18,19,24]. However, few studies on the hydrogeochemistry of the RC region identify areas susceptible to natural dissolution of these minerals (calcite, gypsum and dolomite), and zones with higher saline intrusion by using inverse hydrogeochemical modeling.
The large number of sinkholes distributed along the RC offers the opportunity to study the 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 by means of analyzing the concentration of major elements in groundwater, building hydrogeochemical diagrams, and using the inverse modeling 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/precipitation of minerals and interacts with the rainfall. Understanding such processes can help to implement management strategies at the RC for a better use of groundwater.

Study Zone
The Ring of Cenotes is located in Mexico, at the northwest of the YP, between the parallels 88.5-90.5 • W, and meridians 20.5-21.5 • N ( Figure 1). The aquifer is unconfined, except for a coastal strip where variable hydrogeological confinement occurs [24,25]. The hydraulic gradient typically ranges from 1 to 10 cm km −1 (extremely low). The regional groundwater flow is towards the coast and the RC has been reported as a preferential flow path [8,24,25]. The prevailing climate is wet-tropical, with an annual average temperature of about 26 • C, with a minimum in December and January and maximum from May to August [7]. The rainy season occurs from June to October, while the period between November and May is dry. The annual rainfall averages 1218 mm for the whole YP, while the potential evapotranspiration index is 1200-1400 mm [26]. Within the RC, the annual cumulative rainfall varies from 500 mm in the north coast to 1200 mm in the south [26].
Because of its water reserves and high biodiversity, the RC has been designated a RAMSAR site. Moreover, it is part of a state reserve containing 99 sinkholes that constitute a freshwater source for rural human settlements [9]. The RC lacks of surficial runoff (rivers and creeks), and the aquifer has high transmissivity, particularly along the alignment of karstic structures [8,24]. The flow paths discharge into two coastal lagoons: Ria Celestun and Bocas de Dzilam. This discharge helps in maintaining the health of the mangrove and seagrasses ecosystems which support a high diversity wildlife of ecological and commercial concern [15].
Butterline and Bonet described the study area surficial geology as sedimentary carbonicdolomitic sequences and intercalations of evaporites from Paleocene to Miocene [27]. The RC extends over four geological formations [18,[28][29][30][31]: (a) the Chichen Itzá formation composed of massive fossiliferous limestones from the Eocene, (b) Oligocene stratified sandy limestones with the presence of clays, (c) the Carrillo Puerto formation composed by fossiliferous limestones from the Miocene-Pliocene that covers discontinuously to the Oligocene, and (d) Pleistocene-Holocene littoral and lacustrine sediments with the presence of shells and mollusks along the coast ( Figure 1). Geochemically, the RC has been regionalized [10] into three regions: (a) Region 1 (R1), located at the west, where sulfates, sodium and chlorides increase towards the coast, (b) Region 2 (R2), located at the south side, characterized by slightly acid pH values, lower electrical conductivity, and lower ion (Na + , K + , Mg +2 , Cl − ) concentrations. Here, 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] show that calcite, gypsum, and dolomite are the most common minerals of the RC, and these solid phases interact with groundwater and seawater from the marine intrusion [16]. Apparently, groundwater dissolves calcite, while sulfates and chloride increase its concentrations as the groundwater flows 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 Sierrita de Ticul is the main way along which the dissolved sulfates-originated from gypsum dissolution due to weathering of the Paleocene Icaiché formation-travel from the Chichancanab lagoon (located 110 km southeast of the RC). In the map, the surficial geology described by [29][30][31]. The 18 sinkholes of the study area split into three regions (R1, R2, R3) according to [10]. R is the recharge zone, MC-middle Celestun trajectory, C-Celestun, MD-middle Dzilam trajectory and D-Dzilam. RC's regional groundwater flow was described by [8]. . Selection of the sampled sinkholes followed the regionalization proposed by [10]. During each field survey, water samples were collected at three depth levels (0.5, 5.5, and 10.5 m) using a Van Dorn 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 for analyses. All the collected samples came from the freshwater lens because the halocline was at least 18-25 m below the water table, no saline water was collected neither in the coastal zone (C and D) nor the recharge zone (MD, MC, and R). Temperature (T), pH, electrical conductivity (EC) and dissolved oygen (O2) 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 [32]. All water samples had an ionic balance in the range of ± 5% to ± 10%.

Hydrogeochemical Evolution of Groundwater in the RC
The RC's evolution analysis included ternary hydrogeochemical diagrams of the groundwater, using milliequivalents per liter for all samples (sinkholes 1-22) [33]. The rSO4 −2 /rCl − ratio was plotted to identify the spatial patterns of sulfates and the salinization in the RC [19,24]. The conceptual model was based on the interactions between groundwater and the lithologic solid phases along the regional groundwater flow paths. In the map, the surficial geology described by [29][30][31]. The 18 sinkholes of the study area split into three regions (R1, R2, R3) according to [10]. R is the recharge zone, MC-middle Celestun trajectory, C-Celestun, MD-middle Dzilam trajectory and D-Dzilam. RC's regional groundwater flow was described by [8]. . Selection of the sampled sinkholes followed the regionalization proposed by [10]. During each field survey, water samples were collected at three depth levels (0.5, 5.5, and 10.5 m) using a Van Dorn 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 for analyses. All the collected samples came from the freshwater lens because the halocline was at least 18-25 m below the water table, no saline water was collected neither in the coastal zone (C and D) nor the recharge zone (MD, MC, and R). Temperature (T), pH, electrical conductivity (EC) and dissolved oygen (O 2 ) were measured in situ using a multiparametric probe (Hydrolab G-Quanta ® ). The concentrations of major elements (Na + , K + , Ca +2 , Mg +2 , Cl − , HCO 3 − , SO 4 −2 , NO 3 − ) were determined according to [32]. All water samples had an ionic balance in the range of ± 5% to ± 10%.

Hydrogeochemical Evolution of Groundwater in the RC
The RC's evolution analysis included ternary hydrogeochemical diagrams of the groundwater, using milliequivalents per liter for all samples (sinkholes 1-22) [33]. The rSO 4 −2 /rCl − ratio was plotted to identify the spatial patterns of sulfates and the salinization in the RC [19,24]. The conceptual model was based on the interactions between groundwater and the lithologic solid phases along the regional groundwater flow paths.

Determination of Saturation Indexes (SI)
Calcite, dolomite and gypsum saturation indexes (SI) can be a good proxy of the saturation level of these minerals in groundwater [34]. Positive SI means that the solution is oversaturated, and the mineral will tend to precipitate. When SI is negative, the solution is undersaturated, and the mineral will dissolve. If SI is close to zero, then the solution is at equilibrium. Because it is uncommon to find solutions exactly at equilibrium (SI = 0), an equilibrium interval within 0 ± 5% from the log of K was assumed [33,35]. 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 [36]. Average ion concentrations from sinkholes (Table 1, excluding NO3 − ) were used for the inverse modelling. For the recharge area R, ion concentrations from sinkholes 8, 10-14 were averaged. For MC, concentrations of sinkholes 6, 7 and 9 were averaged. For MD, concentrations from sinkhole 15 were used. In the case of the coastal areas, the average concentrations from sinkholes 1-5 was used for C, and the average concentrations from sinkholes 19, 20, 22 for D. Calcite, dolomite, and gypsum were provided 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 [39]. Finally, O2 and CO2 phases were also included in the model as suggested in Equation (1).
For convergence, a global uncertainty of 0.10 was set. 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 [36]. Table 1. Average concentrations for the representative samples. NO3 − was excluded from the inverse modelling. Regions in this table were defined as the description in Figure 1.  Rainfall (H 2 O) with a high content of carbon dioxide (CO 2 ) forms carbonic acid (H 2 CO 3 ), which allows the dissolution of carbonate rocks when it comes in contact with surficial minerals [33], as described by Reactions (1)-(4):

SP Region
Early studies showed that groundwater in the Yucatan aquifer is near-equilibrium with respect to calcite and dolomite because of the low velocity of groundwater [18,19,24]. However, the dissolution/precipitation rate of these minerals is not quantified yet.

Groundwater Mineral Precipitation/Dissolution Rate Modeling
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 [36]. 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 the estimates of seawater mixing proportions and potential mass transfers from the considered lithological phases in millimoles per kilogram of water (mmol kg −1 H2O ) [36]. Using the RC regionalization proposed by [10], two groundwater trajectories, from the recharge region (R), were analysed using inverse modelling: (a) toward Celestun and (b) toward Dzilam Bravo.
The trajectory to Celestun flows through the Eocene and Miocene formations from the recharge region (R) to the midpoint MC, and then continues to the coast through the Miocene only. The trajectory towards Dzilam flows from the recharge region (R) to the midpoint MD passing through the Eocene, Oligocene, and Miocene formations and continues to the coast through the Miocene only ( Figure 1).
Average ion concentrations from sinkholes (Table 1, excluding NO 3 − ) were used for the inverse modelling. For the recharge area R, ion concentrations from sinkholes 8, 10-14 were averaged. For MC, concentrations of sinkholes 6, 7 and 9 were averaged. For MD, concentrations from sinkhole 15 were used. In the case of the coastal areas, the average concentrations from sinkholes 1-5 was used for C, and the average concentrations from sinkholes 19, 20, 22 for D. Calcite, dolomite, and gypsum were provided 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 [39]. Finally, O 2 and CO 2 phases were also included in the model as suggested in Equation (1).
For convergence, a global uncertainty of 0.10 was set. 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 [36].

Rainfall Analysis
To analyse the rainfall influence in the study area during the sampling period (SP1-SP4), we used diurnal and monthly accumulated precipitation data from ten meteorological stations from the National Commission of Water (CONAGUA): four stations for the recharge zone (Cuzama, Cantamayec, Sotuta, and Kantuninkil), two for Celestun (Kinchil and Celestún), and four for Dzilam Bravo (Tunkas, Temax, Buctzotz, and Dzilam Bravo). We averaged the accumulated precipitation measured during each sampling period (two months).
We analysed the rainfall's influence on the seawater mixing, chemical activities, and saturation indexes through a simple linear regression model. The regression considered the rainfall as the explanatory variable and seawater mixing, chemical activities, and saturation indexes as response variables. Pearson's correlation was used for most of the data that met normality, and Spearman's method was used for those with non-normal distribution. All tests were performed applying a significance level of 0.05. Table 1. Average concentrations for the representative samples. NO 3 − was excluded from the inverse modelling. Regions in this table were defined as the description in Figure 1.

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 (Figures 3 and 4). Marine influence was stronger in Celestun as compared to 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 Figure 4). Pérez-Ceballos et al. [10] described such behaviour in their regionalization, and this was also supported by the results observed in SP4.  Finally, the rSO4 −2 /rCl − ratio in contrast to rCl 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 consequence of the gypsum dissolution, while the trajectory towards Dzilam Bravo showed only seawater enrichment ( Figure 5). During the period SP4, the sulfates and chloride enrichment increased towards the Dzilam trajectory, likely due to reduced rainfall. Celestun's trajectory showed an increase of sulfates associated with gypsum's dissolution from the Eocene aquifer on the middle trajectory (R-MC). This behaviour is consistent with previous studies for the Eocene aquifer near Ticul [19,24]. In the MC-C segment, we found an increase in chloride and a decrease in sulfate concentration (Figure 3), revealing marine influence due to the reduction in the freshwater lens's thickness above the halocline that caused groundwater salinization near Celestun [19]. The freshwater lens reduction also caused decreased sulfates because of seawater mixing.
The Dzilam Bravo trajectory showed an increasing chloride concentration toward the coast (samplings SP1 to SP4) and a decreasing calcium concentration (Figure 3). Figure 3d shows an increase in R-MD sulfate concentration, similar to the Celestun R-MC trajectory only during SP4, probably due to meteorological reasons.
The saturation indexes of calcite, dolomite, and gypsum (SI c , SI d and SI g , respectively) are shown in Figure 4a-c. SI showed equilibrium for calcite (−0.424 < SI c < 0.424) and dolomite (−0.416 < SI d < 0.416) at all samples from R. On the middle points (MC and MD), most of the samples were also in equilibrium with calcite (−0.424 < SI c < 0.424) and dolomite (−0.416 < SI d < 0.416); the only exceptions were the MC average sample, which was undersaturated for dolomite (SI d = −0.74) during survey SP2 and oversaturated during SP3 (SI d = 0.64), and MD samples from survey SP3, which were oversaturated for both calcite (SI c = 0.67) and dolomite (SI d = 1.45). Point C showed chemical equilibrium with calcite (−0.424 < SI c < 0.424), except for the sample from SP2 (SI c = 1.05), while dolomite showed oversaturation in SP2 and SP3 (SI d = 2.29 and 0.69, respectively). Point D showed oversaturation with calcite in SP2 and SP4 (SI c = 1.39 and 0.88, respectively) and with dolomite in SP2, SP3, and SP4 (0.56 < SI d < 2.87). All samples showed the potential to dissolve gypsum, but it was more evident in the R sample (SI g < −0.229, Figure 4c). The recharge zone showed chemical balance with calcite and dolomite but towards the coast dolomite may precipitate, while gypsum could be dissolved at any point of the trajectories described.
Marine influence was stronger in Celestun as compared to 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 Figure 4). Pérez-Ceballos et al. [10] described such behaviour in their regionalization, and this was also supported by the results observed in SP4.
Finally, the rSO 4 −2 /rCl − ratio in contrast to rCl 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 consequence of the gypsum dissolution, while the trajectory towards Dzilam Bravo showed only seawater enrichment ( Figure 5). During the period SP4, the sulfates and chloride enrichment increased towards the Dzilam trajectory, likely due to reduced rainfall.

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 one of the minerals without considering coupled interactions. Conversely, inverse modelling (fed with the information of the lith-

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 one 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, MD to C, D, respectively.
In the trajectory towards Celestun, our results indicate calcite precipitation and gypsum and dolomite dissolution on the segment (R-MC). In contrast, on the second segment (MC-C), the gypsum and dolomites showed a decrease 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 (Figure 6a). The potential dissolution activity rate of dolomite and gypsum was in the range of 0.33-1.05 and 1.06-1.99 mmol kg −1 H2O , respectively (Figure 6b,c). Finally, on 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 (Figure 6a). The potential dissolution activity rate of dolomite showed a decreasing trend (0.07 and 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 survey SP1 (activity rates of −0.14 mmol kg −1 H2O , Figure 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 (Figure 6c).
range of 0.33-1.05 and 1.06-1.99 mmol kg −1 H2O, respectively (Figure 6b, 6c). Finally, on 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 (Figure 6a). The potential dissolution activity rate of dolomite showed a decreasing trend (0.07 and 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 survey SP1 (activity rates of −0.14 mmol kg −1 H2O, Figure 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 (Figure 6c). On the Dzilam Bravo trajectory, the first segment (R-MD) the calcite precipitation showed potential activity between −0.2 and −1.5 mmol kg −1 H2O, except during SP2, which showed a potential dissolution rate equal to 0.056 mmol kg −1 H2O (Figure 6a). The dolomite activity rate was in the range of 0.2-1.46 mmol kg −1 H2O (Figure 6c). The least activity of gypsum was observed in SP1-SP3, with precipitation rates between −0.006 and −0.02 mmol kg −1 H2O (Figure 6c). During SP4, the dissolution rate of gypsum was 1.16 mmol kg− 1 H2O (Figure 6c). On the second segment (MD-D), enrichment with calcite was revealed, while the precipitation dolomite and gypsum showed less activity. However, during SP1, we observed a further enrichment of calcite and dolomite precipitation due to higher rainfall.
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 and 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 was mainly due to chemical exchange with cations such as Mg +2 available in the aquifer and from mixing with marine water. Such a process has previously been described for the coastal area of YP by [16,40,41].
Regarding the trajectory towards Celestun, a constant enrichment of gypsum on the segment R-MC was observed, which matches the evolution 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 was found, in the samples from SP4 survey which showed high gypsum concentrations on the segment R-MD, similar to those in segment R-MC (Figure 6c).  On the Dzilam Bravo trajectory, the first segment (R-MD) the calcite precipitation showed potential activity between −0.2 and −1.5 mmol kg −1 H2O , except during SP2, which showed a potential dissolution rate equal to 0.056 mmol kg −1 H2O (Figure 6a). The dolomite activity rate was in the range of 0.2-1.46 mmol kg −1 H2O (Figure 6c) . The least activity of gypsum was observed in SP1-SP3, with precipitation rates between −0.006 and −0.02 mmol kg −1 H2O (Figure 6c). During SP4, the dissolution rate of gypsum was 1.16 mmol kg −1 H2O (Figure 6c). On the second segment (MD-D), enrichment with calcite was revealed, while the precipitation dolomite and gypsum showed less activity. However, during SP1, we observed a further enrichment of calcite and dolomite precipitation due to higher rainfall.
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 and 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 was mainly due to chemical exchange with cations such as Mg +2 available in the aquifer and from mixing with marine water. Such a process has previously been described for the coastal area of YP by [16,40,41].
Regarding the trajectory towards Celestun, a constant enrichment of gypsum on the segment R-MC was observed, which matches the evolution 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 was found, in the samples from SP4 survey which showed high gypsum concentrations on the segment R-MD, similar to those in segment R-MC (Figure 6c).

Discussion
The groundwater 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), as suggested by the ternary diagrams and the ratio rSO 4 −2 /rCl − . The inverse modelling approach confirmed such hydrogeochemical tendencies.
From the recharge area to Celestun, the increase of sulfates is likely associated with the dissolution of the gypsum's Eocene aquifer. The dissolution of gypsum in R-MC was present in all surveys, which is consistent with the results of Perry et al. [18,24,28], who argued 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 ( Figure 4). Additionally, gypsum and dolomite dissolution, as well as calcite precipitation, were observed. This combination of processes is also known as dedolomitization [42], and is closely related to coupled variations of Mg +2 and Ca +2 concentrations in the groundwater [23].
In the coastal segment of this trajectory (MC-C), gypsum and dolomite showed a decreasing activity (and even precipitation during SP1, Figure 6b,c); but calcite dissolution and marine water influence were also evident ( Figure 4d or Figure 6a). Conceptual models [16,40] predict dolomite precipitation through lattice replacement with Mg +2 from seawater [43]. The observed decrease of sulfates ( Figure 3) was due to seawater mixing as a consequence of the reduction in the thickness of the freshwater lens above the halocline. Overall, the inverse modelling implemented in this work is in agreement with earlier conceptual models proposed by [16,18,19,23,24,[40][41][42] and supports them.
The Dzilam Bravo trajectory also showed salinization toward the coast, although less intense than in the Celestun trajectory (Figure 4d). The enrichment of chloride and sulfate was more noticeable during the dry season of 2008 (SP4, Figure 3d or Figure 5d), likely being a consequence of the scarce rainfall. 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 (Figure 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 (Figures 4 and 6) during SP4 when rainfall was lowest at Dzilam Bravo (Figure 7). The importance of rainfall for Dzilam Bravo is further outlined by the strong inverse relation found for seawater mixing with rainfall in this trajectory (R 2 = 0.94; p = 0.03, Figure  8), showing the effect of rainfall variations on the marine intrusion there. In contrast, although seawater mixing was higher at Celestun, it did not relate significantly with rainfall (R 2 = 0.26; p = 0.49, Figure 8).
Although the regression analyses were non-significant for most of the analysed variables (Figure 9), our findings suggest a groundwater input, non-dependent on the rainfall coming from the south towards Dzilam, as suggested by Long et al. [23]. Earlier studies suggest that gypsum dissolves by weathering downstream in the older formations (Eocene and Paleocene), traveling toward the RC [18,19,24,28] and, according to our findings, later to the coast. This suggests that gypsum concentration depends on the previously dissolved sulfates in the Eocene and Icaiché formations and not necessarily on the local recharge by rainfall. The importance of rainfall for Dzilam Bravo is further outlined by the strong inverse relation found for seawater mixing with rainfall in this trajectory (R 2 = 0.94; p = 0.03, Figure 8), showing the effect of rainfall variations on the marine intrusion there. In contrast, although seawater mixing was higher at Celestun, it did not relate significantly with rainfall (R 2 = 0.26; p = 0.49, Figure 8). For the Celestun trajectory, gypsum was inversely related with the rainfall, although such regression was statistically non-significant (Figure 9a). Similarly, no statistically significant regression of the dissolution/precipitation rates for dolomite and gypsum regarding the rainfall were present for Dzilam (Figure 9b). These results suggest, according to the hydrogeochemical model, that the rainfall causes an important effect on the aquifer's recharge only for the Dzilam Bravo trajectory, similarly to the findings of [21,23]. 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 hydro- Although the regression analyses were non-significant for most of the analysed variables (Figure 9), our findings suggest a groundwater input, non-dependent on the rainfall coming from the south towards Dzilam, as suggested by Long et al. [23]. Earlier studies suggest that gypsum dissolves by weathering downstream in the older formations (Eocene and Paleocene), traveling toward the RC [18,19,24,28] and, according to our findings, later to the coast. This suggests that gypsum concentration depends on the previously dissolved sulfates in the Eocene and Icaiché formations and not necessarily on the local recharge by rainfall.
For the Celestun trajectory, gypsum was inversely related with the rainfall, although such regression was statistically non-significant (Figure 9a). Similarly, no statistically significant regression of the dissolution/precipitation rates for dolomite and gypsum regarding the rainfall were present for Dzilam (Figure 9b). These results suggest, according to the hydrogeochemical model, that the rainfall causes an important effect on the aquifer's recharge only for the Dzilam Bravo trajectory, similarly to the findings of [21,23]. 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 the rainfall effect is more important. This higher seawater mixing towards Celestun relative to Dzilam Bravo could be related to hydrogeological coastal differences among these two sites. Along the Yucatán's coast, the confinement influences the interaction between seawater and the aquifer, affecting the coastal aquifer's salinity [25,44,45]. The coastal confinement effect extends about 20 km towards Celestun, but only about 3 km toward Dzilam Bravo [24,46]. This means For the Celestun trajectory, gypsum was inversely related with the rainfall, although such regression was statistically non-significant (Figure 9a). Similarly, no statistically significant regression of the dissolution/precipitation rates for dolomite and gypsum regarding the rainfall were present for Dzilam (Figure 9b). These results suggest, according to the hydrogeochemical model, that the rainfall causes an important effect on the aquifer's recharge only for the Dzilam Bravo trajectory, similarly to the findings of [21,23]. 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 the rainfall effect is more important.
This higher seawater mixing towards Celestun relative to Dzilam Bravo could be related to hydrogeological coastal differences among these two sites. Along the Yucatán's coast, the confinement influences the interaction between seawater and the aquifer, affecting the coastal aquifer's salinity [25,44,45]. The coastal confinement effect extends about 20 km towards Celestun, but only about 3 km toward Dzilam Bravo [24,46]. This means that Celestun is more confined towards the continent than Dzilam Bravo, and hypothetically, Celestun is more influenced by the seawater and its physical processes, 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 [47,48]. 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 with higher priority during the winter months, when atypical meteorological tides are caused by cold fronts [10,45]. 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. 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 seasonal monitoring along the MD-D segment (sinkholes 17 and 19) and in R-MD (sinkhole 10 to 16) located a 20 km and 30 to 62 km, respectively, 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 of sulfate in wells drilled near Ticul, likely due to extensive use of fertilizers in local agriculture that eventually reach the aquifer contributing to the increase of sulfate concentration [17,23], which could in the short term become unacceptable for human populations [49]. For example, it is known that sulfates may have a laxative effect when combined with Ca +2 and Mg +2 [49,50], 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 (Figures 8 and 9).

Conclusions
Groundwater in the RC flows along two different trajectories with different evolution. Both initiate 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 influences the dissolution and precipitation processes of the lithological phases showing dedolomization in the segment R-MC of the trajectory to Celestun and dolomite precipitation in the coastal zone (MC-C).
In the trajectory towards Dzilam Bravo, calcite precipitation/dissolution processes were more active than gypsum and dolomite. Groundwater with higher gypsum concentration due to the weathering on the Eocene formations was observable only during the lower rainfall sampling (SP4).
The marine intrusion (up to 2% of seawater) in the Celestun trajectory was higher than in the Dzilam Bravo trajectory. Marine intrusion in the Celestun trajectory is not affected by the rainfall, while in the coastal zone of Dzilam Bravo, it strongly varies inversely with rainfall.
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 SI g suggested gypsum dissolution all along the Celestun trajectory (Figure 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 to identify the occurrence of dolomite precipitation at R-MC. Altogether, our results sustain the usefulness of this approach for assessing the evolution of groundwater hydrogeochemistry and quality.
Finally, both regional flow trajectories of the RC are vulnerable to climate change effects, such as the sea level rise, because the hydrogeochemical processes are affected by the interaction with seawater. However, Celestun is more likely susceptible to the rising sea level. In contrast, Dzilam Bravo can be affected by changes in the rainfall patterns because the hydrogeochemical processes taken place, depend more on the aquifer recharge. Monitoring of marine intrusion and phases dissolution in the RC is highly recommended to assess the evolution of the system and could be done simply through chlorides and sulfates measurements.

Data Availability Statement:
The chemical data is presented in Table 1, data for Figure 2 can be found on [37,38].