2.1 Meteorological data

Data on monthly accumulated precipitation were retrieved from the historical database of the Sistema de Información Hidrológica (Hydrological Information System; SIH) of Mexico (^{SIH, 2023}). The data were related to four stations, those closest to the sampling sites (Fig. 1), from November 2018 to March 2021. For characterizing seasonal trends, we defined three seasons based on historical data (1951-2010), as shown in Table I; thereby, the following were proposed: rainy (June to October, ≈ 900 mm precipitation), cold-front (November to February, ≈ 250 mm precipitation), and dry (March to May, <200 mm precipitation).

Fig. 1 Sampling sites, altitude (masl), and codes of the nearest meteorological stations. 

2.2 Hydrochemistry

The chemical composition of rainwater is represented by the major ions (i.e., sodium, magnesium, potassium, calcium, nitrite, nitrate, sulfate, phosphate, and chloride) measured with an 882 Compact IC Plus ion chromatograph (Metrohm) with a conductivity detector. The samples were stored at room temperature for up to three months, then filtered and analyzed. The anions were quantified under the following conditions: Metrosep A4 250/4 column with a pre-column Metrosep R2 Guard 3.5, mobile phase 12.6 mM Na₂CO₃/11.9 mM NaHCO₃. For their part, conditions for cations were Metrosep C4 250/4 column with a pre-column Metrosep R2 Guard 3.5, mobile phase nitric acid 2.7 mM, and dipicolinic acid. Each sample was analyzed twice with the detection limit of 0.01 mg L^-1.

The data normality was evaluated with a Shapiro-Wilk test. Samples that did not have a normal distribution were analyzed with a Kruskal-Wallis test. The rainwater results were analyzed using the volume-weighted mean (VWM) in microequivalents (µeq L^-1). The VWM of the various species present in the rainwater samples were calculated using Eq. (1):

where C _i is the concentration of any given ion (µeq L ^-1 ) and P _i is the precipitation for the event (mm). For our purposes, each precipitation event comprises the total monthly precipitation, as all analyses were completed in composed water samples accumulated by month, i.e., all the rain events that occur in one calendar month. It is important to note that for the Bacalar sampling site, the chemistry data were only available for five months (i.e., from July to November 2019).

A Pearson correlation was completed for the measured ions. The ion ratio SO₄ ^2-+NO₃ ^-/ Ca²⁺+Mg²⁺ was used for assessing the acidity source, where values lower than one (< 1) represent the dominance of alkali ions, and values greater than one (>1) suggest acidic dominance. The NO₃ ^-+Cl^-/SO₄ ^2- ratio relates to the contribution of either acidic source (nitric or sulfuric); values <1 suggest sulfuric acid-derived acidity, while >1 indicate nitric acid-derived acidity (^{Salve et al. 2008}). The neutralization factor per cation (NF _x ) was calculated using Eq. (2) following ^{Oduber et al. (2020)}:

where C _x is the concentration of each cation, and CSO ₄ ^2- and CNO ₃ ^- are the sulfate and nitrate concentrations (µeq L^-1), respectively. We also calculated the ratios between the concentration of elements in the sample (Cx/C _{r sample} ) with C _x for all elements calculated and C _r as the concentration of the reference element. To distinguish the source or origin of the elements in precipitation, we used two sets of calculations; the first ones are seawater (Eq. [3]) and crust (Eq. [4]) enrichment factors (EF), as follows:

The enrichment factors evaluate the same elements in crust (C _x /C _{r crust} ) and seawater (C _x /C _{r sea water} ) relative to the reference concentrations in each of these reservoirs. Conventionally, Na is the reference element for seawater, assuming marine origin (Sp/Na; ^{Keene et al., 1986}); whilst Ca is the lithophilic reference element for continental crust (Sp/Ca; ^{Zhang et al., 2007}). The concentrations utilized in the EF calculations above are expressed in μeq L^-1.

The calculations for the contribution of the second source were performed with the sea salt fraction (%SSF, Eq. [5]) and a non-sea salt fraction (%NSSF, Eq. [6]); the former related to the percentage of seawater from sea spray, and the latter to crustal fractions (i.e., terrestrial dust; ^{Zhang et al., 2007}). The calculations were as follows:

A non-natural or anthropogenic fraction (Eq. [7]) is commonly assumed as the remnant contribution after computing seawater and crust, expressed as the difference from 100%.

Lastly, a comparison between some rainwater sampling stations worldwide was run to provide the present results with a context.

2.3 Isotopes

The δ²H δ¹⁸O isotopes were measured at the University of Waterloo Environmental Isotopes Laboratory using a Los Gatos Research (LGR) liquid water isotope analyzer (LWIA), model T-LWIA-45-EP instrument with published precisions of δ²H = ± 0.3 ‰ and δ¹⁸O = ± 0.1‰.

The d-excess parameter was obtained as a function of δ¹⁸O by finding d in a line of slope = 8 in the d-space (^{Gat, 2005}; Eq. [8])

Negative or close to zero values for d-excess indicate a possible loss of light isotopes or fractionation of rainwater by evaporation; in this way, two rain samples with a negative value of d-excess were eliminated so that the final number of samples was 100 for subsequent analyses. We carried out a Kolmogorov-Smirnov test for normality. For normally distributed data, we used a two-way ANOVA and a Kruskal-Wallis for non-normally distributed data.

2.4 Backward air parcel trajectories

We used the Hybrid Single Particle Lagrangian Integrated Trajectory model (HYSPLIT; ^{Stein et al., 2015}) to obtain the backward trajectories of the air masses to relate moisture sources with the chemical and isotopic composition of the rainwater samples (^{Burnett et al., 2004}). For the present study, and based on our season’s classification, we selected September, January, and May as representative months for each season (i.e., rainy, cold-front, and dry season, respectively) to explain the trajectory models. We used the trajectory ensemble mode with a 144-h time window between the days 10 to 15 of each of the selected months (mid-month) at three altitudes (250, 500, and 1000 m above the ground), assuming that the cloud base level in the Yucatan peninsula is around 1000 m (^{Nuijens et al., 2014}; ^{NOAA, 2024}). The endpoint of the trajectories corresponds to each sampling location. The seasonal differences in rainwater chemistry, stable isotopes, as well as air parcels were explored based on the year division shown in Table I.

3.1 Hydrochemistry

From the chemical analysis, we found that nitrites (NO₂ ^-) and phosphates (PO₄ ^3-) were always below the detection limit (< 0.01 mg L^-1) in all four sites and that magnesium (Mg²⁺) was also below such limit in Cancún. Figure 2 shows the results of ion concentrations by sampling site and climatic season; the complete dataset is presented in Table SI in the supplementary material. The anion deficiency associated with the exclusion of HCO₃ ^- does not seem very relevant (cations vs. anions correlation coefficient, R²= 0.9508) and might not have had a significant contribution (^{Keresztesi et al., 2019}).

Fig. 2 Box plot of ion concentrations (mg L^-1) at each sampling site. Missing values were below the detection limit (< 0.01 mg L^-1). Chemical species appear without charge. The Shapiro-Wilk normality test indicated that none of the ions had a normal distribution (p < 0.001); therefore, the Kruskal-Wallis test showed that data from November 2018 to December 2019 was different by site, not by season, as shown in Table II.  

The calculated VWM values are shown in Table III. According to these concentrations, the anions in Cancún, Tankah, and Chetumal have the following trend: Cl^- > SO₄ ^2- > NO₃ ^-; whereas in Bacalar, chloride, and sulfate switch order (SO₄ ^2- > Cl^- > NO₃ ^-). For cations, Tankah and Chetumal have the trend Ca²⁺ > Na⁺ > Mg²⁺ > K⁺, Cancún Ca²⁺ > Na⁺ > K⁺ > Mg²⁺, and Bacalar Ca²⁺ > Mg²⁺ > Na⁺ > K⁺. In all cases, Ca²⁺ was the predominant cation. Given that soil dust contains calcite minerals from calcium carbonates present in the soil of the region (^{Bravo et al., 2000}), our results reflect the presence of calcium-rich dust in the Yucatan peninsula’s atmosphere (^{Báez et al., 1997}; ^{Galloway and Lickens, 1978}). Since the calculation of the ion ratio SO₄ ^2- + NO₃ ^-/Ca²⁺ + Mg²⁺ was found to be < 1, it points to alkali dominance (from calcium and magnesium). Despite the alkaline dominance, the small acidity present in the rainwater samples showed dominance from nitric acid (NO₃ ^- + Cl^-/SO₄ ^2- > 1); only Bacalar has sulfuric acid-dominant acidity (Table III).

The correlation coefficients shown in Table IV suggest that potassium might originate from both sea spray and anthropogenic emissions such as biomass burning (^{Oduber et al., 2020}; ^{Raga et al., 2021}). Even though potassium feldspars are present in the carbonaceous shale of this area (^{Lefticariu et al., 2006}), we consider this to be the lesser contribution as it is mainly buried below the actual surface level. Calcium is considered to come from two sources: seawater (the least) and crustal inputs such as soil dust (^{Bravo et al., 2000}); however, calcium can also be incorporated into the atmosphere via fertilizers (^{Zeng et al., 2020}), supported by its correlation with NO₃ ^- (r = 0.954) and sulfates (r = 0.928). The correlation coefficient between SO₄ ^2- and NO₃ ^- (r = 0.98) suggests a common anthropogenic source due to SO₂ and NO_X in the atmosphere (^{Zhang et al., 2007}; ^{Oduber et al., 2020}). There are other ion ratios with good correlation coefficients: Ca²⁺ and SO₄ ^2-; Ca²⁺ and NO₃ ^-; Mg²⁺ and Cl^-; all of them have a common origin that exists in rainwater as a result of atmospheric chemical reactions of acids such as H₂SO₄, HNO₃, and HCl with alkaline compounds rich in calcium and magnesium, found in the atmosphere from wind-blown dust (^{Chandra et al., 2005}; ^{Salve et al., 2008}).

According to the neutralization factor (NF) results (Table V), high NF values for the alkali element Ca²⁺ (from 2.5 to 8.5) are associated with the transport of soil dust with a large quantity of calcite and dolomite due to the calcareous soil of the study area (^{Peikam and Jalali, 2021}).

The enrichment factors of the rainwater samples EF_crust <1 suggest that the element is mainly from crustal origin; values >1 are associated with other sources, either seawater or anthropogenic origin. The EF_seawaters for SO₄ ^2-, Ca²⁺, Mg²⁺, and K⁺ were found to be >1, which indicates enrichment or impacts by sources other than marine; for example, soil or anthropogenic sources (Table VI).

The correlation between Ca and Mg (r = 0.88) indicates crustal origin. In the case of Cl^- its correlations were found to be ≤ 1 in Cancún, Tankah, and Chetumal, indicating an important contribution from seawater. This last correlation is consistent with the relationship between Cl^- and Na⁺ ions (r = 0.97), whose values are very close to those found in seawater (1.2; ^{Möller, 1990}). In sum, our results indicate that chloride, sulfate, sodium, and magnesium are mostly from the sea-salt fraction in the coastal stations (Table VII).

Despite being present in other coastal locations, Mg in Cancún was always below the detection limit. This element has been considered a continental ion tracer linked with economic activities (^{Bravo et al., 2000}). Despite having material extraction sites in the study area, the thickness of the exposed geological formations and the elevated content of calcium carbonate (^{Smart et al., 2006}) might overshadow Mg. Accordingly, calcium is mostly of crustal origin. Potassium, nitrates, and sulfates largely have a non-natural origin (Table VIII); that is, anthropogenic activities in the Yucatan peninsula (^{Bravo et al., 2000}; ^{Raga et al., 2021}).

Variable wet and dry deposition has been previously observed in tropical coastal environments (^{Palani et al., 2012}), which can be attributed to meteorological factors such as extreme rain events and air and moisture parcel trajectories (^{Mullaugh et al., 2013}). The coastal sites have chloride and sodium from the sea salt fraction, whereas sulfates, calcium, magnesium, potassium, and nitrates also have a non-sea salt origin. Similar results were reported by ^{Cerón et al. (2002)}, who found potassium, calcium, and sulfate predominant from the non-sea salt fraction; these authors also stated that nitrate is strongly associated with wind direction and elevated nitrates are related to biomass burning (forest fires); however, supported on the published literature, we consider that sulfate and potassium are the ions most related to biomass burning (^{Oduber et al., 2020}; ^{Raga et al., 2021}). Some aerosol particles, such as metals and metal oxides, as well as airborne microorganisms, pose risks to human health (^{Rodríguez-Gómez et al., 2020}) or the environment (^{Rahim et al., 2019}), whereas hurricanes might represent additional entrance of ions, potentially posing dangerous wet deposition into the ecosystem (^{Qiu and Felix, 2021}).

Table IX shows the chemical composition of rainwater in some stations located in tropical and subtropical zones. Notice, on one side, the great contribution of chloride and sodium in Tankah due to the marine influence, and on the other, the unusual nitrate in Chetumal, which may come from fertilizers (^{Zunckel et al., 2003}). Hence, we can conclude that rainwater chemistry experienced a rather local influence from sea-salt spray, land-blown dust, and anthropogenic impacts such as agriculture and biomass burning.

3.2 Isotopes

The local meteoric water line (LMWL) for the coast of Quintana Roo is shown in Figure 3. In comparison with the Global Meteoric Water Line (GMWL), its slightly higher slope (13.802) suggests a greater δ²H enrichment due to evaporative losses (^{Clark and Fritz, 1997}; ^{Sun et al., 2021}). Not surprisingly, the LMWL showed an amount effect (^{Dansgaard, 1964}); that is, an inverse correlation with the amount of rainfall (r = -0.48). Note that the altitude effect is not relevant because all sampling stations are at similar altitudes. We also observed that some samples (Bacalar, October 2019 and June 2020; Cancún, October 2020, Chetumal, May 2019 and June 2020; Table SII in the supplementary material) have a depleted isotopic composition, similar to the water precipitated with hurricanes Delta and Gamma, which comes to show that some monthly accumulated precipitation samples are isotopically similar to tropical storms.

Fig. 3 Local meteoric water line for the Caribbean coast, Quintana Roo, Mexico. 

The high d-intercept value indicated different moisture sources and local climatic patterns (^{Modon et al., 2022}). Some of the most depleted values correspond to hurricanes Delta and Zeta (-8‰ in δ¹⁸O), which were reported as category 2 hurricanes (^{Sánchez-Murillo et al., 2019}). The other main source of moisture comes from North America during the cold-front season (see section 3.3). The local climatic patterns are likely related to convective rains and mixing between two or more air parcels. In coastal locations, convective rains are generally formed from vertically moved ocean-evaporated moisture with high δ¹⁸O (relative to its correspondent δ²H) and have been associated with specific microphysical conditions of rain formation (^{Aggarwal et al., 2016}). The mixing between two or more air parcels does not fit the conventional Rayleigh distillation curve; because of this, the resulting isotopic ratios in tropical and subtropical air parcels are not easily distinguished (^{Galewsky and Hurley, 2010}).

The slope value of the LMWL is slightly greater than those reported for the Cholulteca River basin in Honduras (^{García-Santos et al., 2022}) and the Limbang River basin in Borneo (^{Modon et al., 2022}), which indicates that moisture-generating rain in this area of the world is not formed in equilibrium conditions. The moisture originates in the Caribbean Sea and the Gulf of Mexico with no apparent inputs from the Pacific.

The isotopic composition of each sampling location is shown in Table X. The Kolmogorov-Smirnov test for normality showed that the isotopic composition of the analyzed rainwater samples did not have a normal distribution (p < 0.001), whereas d-excess distribution was normal (K - S = 0.067; p = 0.2). The Kruskal-Wallis test for δ²H and δ¹⁸O reported that the isotope composition was not different by site (δ¹⁸O ‰ p = 0.277 and δ²H ‰ p = 0.425) though it was different by season though (δ¹⁸O ‰ p = 0.012 and δ²H ‰ p = 0.0003), a result opposed to that of the water chemistry, which was different by site, but not by season. This difference is because water chemistry and isotope composition respond to different factors and effectors. Water chemistry responds to the local conditions of the air column below the cloud; thus, when precipitation falls, the aerosols and suspended particles reflect the local conditions, with a large influence of marine aerosols and natural and non-natural dust. On the other hand, the isotope composition of moisture has a lesser local influence since its origin (air parcel trajectory) and the physical condition in which water evaporates and condenses are the most relevant factors for variable isotope composition through the seasons.

The d-excess results showed an average of 13.0 ± 3.2 (5.5 to 24.1), similar to other values reported in Central America (^{García-Santos et al., 2022}). The two-way ANOVA for the d-excess showed differences only by season (F = 12.27, p < 0.001). The post-hoc test Tukey HSD showed that the cold-front season is different (q = 14.77, p < 0.001). When d-excess values are greater than 10‰, they represent potential re-evaporation of precipitation (Fig. 4); however, such precipitation has also been associated with rain events from cold fronts (^{García-Santos et al., 2022}). It also shows enhanced terrestrial recycled or evaporative moisture sources with similar d-excess. Depleted precipitation due to recycled precipitation has also been reported in Nicaragua (^{García-Santos et al., 2022}).

Fig. 4 Deuterium excess values (d-excess) in ‰. The horizontal line in 10 indicates the mean global d-excess value. 

Differences in d-excess vary in response to the event type (frontal systems or low-pressure events) or vertical mixing in convective systems rather than seasonal temperature changes (^{Guan et al., 2013}). The elevated d-excess measured in the present study suggests that clouds may be formed in evaporative conditions, and there is mixing of advective condensation in convective rain systems with moisture from the ocean and probably other sources. This phenomenon might be especially relevant in Bacalar, where, based on the lake’s water isotope composition, an important amount of moisture may evaporate from Bacalar Lake, condense, and precipitate again (^{Evans et al., 2018}).

3.3 Backward air parcel trajectories

Air parcel trajectories were obtained for the entire sampling period (Fig. S1 in the supplementary material). However, in terms of seasonality, we selected for this study three months (September, January, and May) to represent each of the distinctive seasons in the Yucatan peninsula. Mid-September’s rainy season trajectory frequencies (Fig. 5) showed the dominance of east winds from the Caribbean and the mid-Atlantic, passing over Florida and Cuba. These winds probably convey moisture from the mid-Atlantic and the Caribbean into the Yucatan channel, including tropical storms and hurricanes, frequent in this season. The cold-front season trajectory frequencies are exemplified by mid-January in Figure 6, which shows the influence of both north and east winds. The former results from polar air masses moving from the Arctic through North America into the Gulf of Mexico, introducing moisture from the Midwest USA and the Gulf of Mexico into the Yucatan peninsula. This season is characterized by a slightly greater δ²H (and d-excess) than the other two seasons. Cold-front enriched rains have been previously reported by ^{Sánchez-Murillo et al. (2019)} and ^{García-Santos et al. (2022)} in the Caribbean. Finally, the dry season trajectory frequencies (mid-May, Fig. 7) depict dominance of east winds from the Caribbean again, but from the east-southeast direction, moving south of Cuba, from the Greater Antilles and the Yucatan basin. These backward air mass movement trajectories suggest that the predominant air mass that reaches the Yucatan peninsula comes from the Caribbean with contributions of terrestrial sources from the islands; consequently, the rain chemistry retains sea salt influence, while a large part of the ions in precipitation moves from the Caribbean most of the year. We acknowledge that the El Niño-Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO) play an important role in rainwater precipitation at a regional scale; yet, a long record is necessary to find significant changes in rainfall patterns (^{Díaz and Giambelluca, 2012}), which is out of the scope of this research.

Fig. 5 Backward air parcel trajectories for mid-September 2019. 

Fig. 6 Backward air parcel trajectories for mid-January 2019. 

Fig. 7 Backward air parcel trajectories for mid-May 2019. 

As previously observed in the central Mexican highlands, there is continuous and uniform air parcel movement from the Caribbean Sea (^{González Ríos, 2022}). As there are no topographic barriers, rainfall isotopic pulses are associated with moisture origin (^{García-Santos et al., 2022}) with no isotopic effect due to altitude. This variation in air mass movements is also responsible for variations in rain events, changes in sources of atmospheric moisture, and various evaporation effects, as shown by the isotopic composition.

In addition to the back trajectories of air parcels as tools for exploring moisture sources, atmospheric bulk deposition measurement may be assisted or accompanied by remote sensing (e.g., MODIS) or modeling atmospheric transport; these complementary methods of direct measurement can produce information regarding local or long-distance sources of minerals and contaminants (^{Das et al., 2013}).