This study was conducted from March 2023 to July 2024 in Gairaca Bay, a protected bay located within TNNP, on the Caribbean coast of Colombia (11°15’33’’ N, 73°24’06’’ W) ( Figure 1 ). The bay is situated between the eastern sector of Taganga and the Piedras River basin, and features diverse geomorphology including rhodolith beds, coral reefs, rocky shores, sandy beaches, shallow sandy-bottom areas, and a transitional zone at the bay’s entrance influenced by open ocean processes ( Figure 1 ).

The region experiences pronounced seasonal variability, modulated by the interplay of coastal upwelling, precipitation, and freshwater inflow (Bayraktarov et al., 2014) (Ricaurte-Villota et al., 2025). The local climate follows a bimodal seasonal pattern, with alternating wet and dry periods throughout the year. The rainy season typically extends from August to November, while the dry season occurs between December and April. During the dry season, intensified northeasterly trade winds promote upwelling of SUW with sea surface temperatures between 19–25 °C and salinities exceeding 36.5 (Paramo et al., 2011; Correa-Ramirez et al., 2020). In contrast, rainy season conditions are characterized by warmer surface waters (average 28.7 °C, reaching up to 30.3 °C) and reduced salinity (~34.7) due to increased freshwater input from the Magdalena river and local tributaries draining from the Sierra Nevada de Santa Marta (Arévalo-Martínez and Franco-Herrera, 2008; Bayraktarov et al., 2013; Alvarado-Jiménez et al., 2024) that alters salinity, total alkalinity (TA), and dissolved inorganic carbon (DIC) within the bay (Ricaurte-Villota et al., 2025). These seasonal changes in water mass properties are driven by the alternating dominance of coastal upwelling and fluvial inputs.

To further characterize the upwelled water mass, we used ARGO profiling float data to identify the presence of SUW in the offshore region of the Colombian Caribbean. The gridded distribution of potential temperature (θ) was derived from validated ARGO profiles (WOD code: 12; originator’s flag set to use = 12), collected in April 2024 approximately 200 km off the coast during oceanographic cruises. The figure was generated using the “Two Scatter Windows” tool in Ocean Data View (ODV) with gridded field interpolation enabled. SUW was identified by its characteristic high-salinity core (~37.16 ± 0.18) located at ~120 m depth (Correa-Ramirez et al., 2020). This thermohaline signature is consistent with SUW properties previously described for the western tropical Atlantic, supporting the notion that this water mass is the primary source of upwelled water in the region during the dry season ( Supplementary Figure SM_1 ).

Three stations were strategically selected to represent different habitats and hydrodynamic conditions: A rhodolith bed located at the center of the bay (sampled at 1, 7, and 15 m depths), the outer bay at the entrance (10 m depth), and the inner bay, a shallow sandy-bottom zone near the beach (sampled at 1 and 6 m depths). The selection of the three sampling sites within Gairaca Bay, was guided by the objective of capturing environmental contrasts within a relatively small spatial scale, such as physical and oceanographic features related to depth, exposure to oceanic exchange, and proximity to terrestrial inputs ( Figure 1 ).

At each station and depth, in situ measurements of temperature, salinity, and pH were taken immediately after sample collection. The pH on the total scale (pH_T) was measured using a SI Analytics HandyLab 100 meter equipped with a WTW Sentix 41 pH electrode. The instrument had a resolution setting of 0.001 pH units (3 decimal places) and an accuracy of ≤ 0.005 pH ±1 digit. pH measurements were recorded in millivolts and later converted to pH_T using the values of Tris buffer (Batch T41) and the in situ temperature, following standard operating procedures (Dickson et al., 2007). The electrode was calibrated monthly using NBS buffers (pH 4.01, 7.00, and 10.00), and during each calibration, the slope was checked and consistently exceeded 97%, indicating a near-Nernstian response (Dickson et al., 2007). Dissolved oxygen was measured in situ using a WTW Oxi 3310 meter, with a manufacturer-reported accuracy of ≤ 0.1 K ±1 digit. Salinity and conductivity were measured using an SI Analytics HandyLab 200 meter, with an accuracy of ≤ 0.5% of the measured value ±1 digit. These measurements were taken immediately after water sample collection at each station and depth to ensure reliable representation of ambient conditions. Additionally, 500 mL seawater samples were collected monthly either by SCUBA diving or using a 7-liter Niskin bottle. Samples for TA analysis were preserved with saturated mercuric chloride to prevent chemical alterations. All samples were transported to the Water Quality Laboratory at the University of Magdalena and stored at temperatures below 17 °C until further analysis.

Precipitation data for the study period were obtained from the Institute of Hydrology, Meteorology, and Environmental Studies (IDEAM), using records from the Simón Bolívar Airport station near Gairaca Bay. Wind speed and wind direction data were retrieved from the Copernicus product Global Ocean Hourly Sea Surface Wind and Stress, with a spatial resolution of 0.125° (DOI: 10.48670/moi-00305). Continuous temperature monitoring was performed in situ at the rhodolith bed (15 m depth) using an Onset HOBO UA-002–64 data logger, programmed to record every 2 hours.

All nutrient analyses were performed following the methodology described by Garay et al. (2003). Dissolved inorganic nutrient concentrations were determined using standard colorimetric techniques. Ammonium (NH₄ ⁺) was quantified using the indophenol blue method, which involves its reaction with phenol and sodium nitroprusside in an alkaline medium in the presence of hypochlorite. Nitrite (NO₂ ^-) was measured via diazotization, using sulfanilamide and N-(1-naphthyl) ethylenediamine as reagents. Nitrate (NO₃ ^-) was first reduced to nitrite using a cadmium column activated with copper sulfate and then analyzed using the same procedure as for nitrite. Inorganic phosphate (PO₄ ^3-) was quantified using the molybdenum blue method, employing a reagent mixture of ammonium heptamolybdate, ascorbic acid, sulfuric acid, and potassium antimonyl tartrate. Absorbance measurements were performed using a UV-Vis spectrophotometer, with wavelengths ranging from 543 to 885 nm depending on the nutrient analyzed.

TA was determined via open-cell potentiometric titration with 0.1 mol L^-1 hydrochloric acid buffered in 0.6 mol L^-1 NaCl (Certified Reference Material - CRM, Dickson Laboratory, Batch 205), following the Remarco protocol (Bernal et al., 2021). The titration’s analytical accuracy was monitored using the CRM with an uncertainty of ±10 µmol kg^-1.

Using the measured TA and pH_T values, additional carbonate system parameters, dissolved inorganic carbon (DIC), partial pressure of carbon dioxide (pCO₂), bicarbonate (HCO₃ ^-) and carbonate ions (CO₃ ^2-), and saturation state for calcite (Ω_cal) and aragonite (Ω_arag), were calculated using the Excel-CO2SYS software version 2.5 (Pierrot et al., 2006). Dissociation constants K1, K2 were taken from (Mehrbach et al., 1973) refit by (Dickson and Millero, 1987), KHSO₄ by (Dickson, 1990) and boron concentration following (Lee et al., 2010).

Trends in sea surface temperature, wind velocity, and precipitation were analyzed through time series plots to identify fluctuations associated with different climatic phases (upwelling, non-upwelling, transition pre-upwelling, and transition post-upwelling). Climatic classifications were based on temporal patterns of these environmental drivers.

Descriptive statistics compiling the mean ± standard deviation values of carbonate system variables such as TA, DIC, pH_T, Ω_arag, and CO₃ ^2- were done per site, climatic season, and year (2023–2024). To illustrate the relative temporal deviations in the carbonate system (increases or decreases of the carbonate variables over the study time), minimum and maximum seasonal deltas (Δ - deviations from the mean) were calculated for TA, DIC, Salinity, and Ω_arag across the rhodolith bed bottom, outer bay, and inner bay sites for each climatic season in 2023 and 2024.

Deltas (Δ_i) were computed as:

Where X i is the observed value and X ¯ is the overall mean of each variable, calculated using all available observations across sites, seasons, and years. The mean was chosen as the reference value to center the data around a consistent baseline, facilitating the comparison of variability across sites and time frames. This method allowed for the identification of relative increases or decreases in carbonate system parameters with respect to their average behavior.

Sectional plots of each parameter were then generated in Ocean Data View (ODV) using the Gridded Field option with DIVA gridding, with sampling date (month–year) on the X-axis, depth (m) on the Y-axis, and parameter concentration on the Z-axis. To enhance interpretability, the X-axis scale was set to 50 and the Y-axis scale increased to 350 ‰, a vertical extrapolation of ± 0.5 m was applied, and contour lines were added to highlight gradients throughout the water column.

PERMANOVA was performed to detect differences in the carbonate system parameters (TA, DIC, pCO₂, Ω_cal, Ω_arag, pH_T, salinity) among sites, depths, and seasons. A Euclidean distance matrix was computed using standardized data with the vegdist function (vegan package, R version 4.3.2). PERMANOVA was performed using the adonis2 function, with 9,999 permutations and type II sums of squares. The homogeneity of multivariate dispersion was tested using the betadisper function and confirmed by ANOVA (p > 0.05).

To assess the predictive capacity of temperature on carbonate system variables, we fitted both simple and multiple linear regression models to estimate dissolved inorganic carbon (DIC, µmol kg^-1) and partial pressure of CO₂ (pCO₂, µatm). The analysis used discrete in situ data collected from three representative sites within Gairaca Bay: rhodolith bed bottom, inner bay and outer bay. The following models were applied:

Simple regression model for DIC:

Simple regression model for pCO₂:

Multiple regression model for pCO₂:

Model fitting was conducted in R using the lm () function. Model performance was assessed using the coefficient of determination (R²) and root mean square error (RMSE).

An additional PERMANOVA was performed between seasons at the rhodolith bed bottom site. SIMPER analyses were carried out to identify the contribution of individual variables to the observed dissimilarities between seasons, with 9,999 permutations. Spearman’s rank correlation was used to assess relationships between environment and carbonate system variables (T°C, pCO₂, CO₃ ^2-, salinity, pH_T, DIC, Ω_arag) and nutrient concentrations (ammonium - NH₄ ⁺, nitrate - NO₃ ^-, nitrite - NO₂ ^-, and phosphate - PO₄ ^3-) at the rhodolith bed bottom. All statistical analyses were performed in RStudio (version 2024.12.1 Build 563; Posit Software, PBC, 2025), using the packages vegan (Oksanen et al., 2025), ggplot2 (Wickham, 2016), car (Fox et al., 2024) (Fox and Weisberg, 2019), dplyr (Wickham et al., 2023a), and tidyr (Wickham et al., 2023b).

Seasonal and spatial variability in discrete in situ measurements of temperature, salinity, and dissolved oxygen was evident across the main sampling sites and depths in Gairaca Bay between March 2023 and July 2024. At the rhodolith bed bottom (15 m depth), temperatures ranged from 26.2 °C during upwelling (March 2024) to 30.5 °C in the post-upwelling transition (June 2024), with salinity between 33.5 and 35.9 and oxygen concentrations inversely correlated with temperature (3.50–7.28 mg·L^-1). At the inner bay (6 m depth), temperatures varied from 27.0 to 30.8 °C, salinity ranged from 33.0 to 35.9, and oxygen concentrations declined from a peak of 7.93 mg·L^-1 (March 2023) to 4.03 mg·L^-1 (October 2023). The outer bay (10 m depth) showed intermediate conditions, with temperatures between 26.4 and 30.5 °C, salinity from 32.7 to 35.8, and dissolved oxygen following the seasonal temperature pattern (3.51–7.67 mg·L^-1) ( Table 1 ). Additional in situ data from 1 m and 7 m depths at the rhodolith bed, and from 1 m depth at the inner bay, are presented in Supplementary Material ( Supplementary Table SM_1 ).

Mean wind speeds and water temperatures (the latter continuously measured with a HOBO sensor at the rhodolith bed bottom) exhibited a clear seasonal pattern influenced by climatic and oceanographic processes. Unless otherwise noted, all mean values are presented with their corresponding standard deviations (mean ± SD). In 2023, water temperature and wind speed exhibited marked seasonal variability associated with the upwelling season. The lowest mean water temperature was observed in March (24.73 ± 0.66 °C), corresponding to the peak of the upwelling season, while the highest was recorded in June (28.67 ± 0.98 °C), followed by a slight decrease in July (28.14 ± 0.60 °C), reflecting the transition out of the upwelling phase. Wind speeds during this period peaked in July, reaching a monthly average of 7.80 ± 1.62 m·s^-1, with a maximum of 10.28 m·s^-1. The lowest wind speeds occurred in May and June, with daily minima of 2.04 m·s^-1 and 1.69 m·s^-1, respectively, indicating a relaxation phase before re-intensification ( Figure 2 ).

In 2024, the seasonal pattern remained consistent, but wind speeds were generally stronger and temperatures slightly higher during the post-upwelling transition. The highest mean wind speed was observed in January (10.06 ± 1.16 m·s^-1), and high wind activity persisted through February to April (means between 8.30 and 8.53 m·s^-1), with peak gusts exceeding 13 m·s^-1 in February ( Figure 2 ). As wind strength declined in May (5.83 ± 1.97 m·s^-1), water temperatures rose to a maximum monthly mean of 27.95 ± 0.44 °C, marking the post-upwelling transition. The coldest daily temperature in 2024 was 23.27 °C, recorded in January, slightly earlier than in 2023.

Wind direction also followed a seasonal pattern, with prevailing northeasterly and easterly winds during upwelling months (January to March and December), reinforcing the wind-driven upwelling dynamics. In contrast, during transitional and non-upwelling periods, wind direction became more variable, with increased frequencies from the southeast, south, and southwest, indicating a weakening of the trade wind system and reduced upwelling potential ( Supplementary Figure SM_2 ).

Precipitation in 2023 showed a pronounced seasonal pattern influenced by the upwelling season and associated climatic transitions. During the upwelling season (March to April) only 0.1 mm of rain was recorded on a single day in March, and a moderate onset of rainfall in April, for a total of 81.9 mm over 3 rainy days, with a peak of 54.1 mm. During the post-upwelling transition (May to July), rainfall initially decreased to 12.8 mm in May, but increased considerable in June, reaching 175.6 mm, with an extreme event of 102.4 mm on a single day. Rainfall decreased again in July to 51.1 mm, although distributed over 13 rainy days. In 2024, the first quarter of the year was extremely dry, with only 0.4 mm of rain in February. A slight increase occurred in April (5.9 mm) and May (6.8 mm), with isolated light rain events. June, marked the onset of the rainy season, with 126.9 mm, including 7 days with rainfall exceeding 5 mm, and a maximum daily value of 43.0 mm. In July, rainfall was 72.2 mm, although concentrated in fewer events, with one day of heavy rainfall peaking at 63.1 mm ( Figure 3 ).

A summary of carbonate system variables and in situ environmental parameters across the three main sites, rhodolith bed bottom (15 m), inner bay (6 m), and outer bay (10 m), is shown in Table 1 . A full dataset including additional depths and sites (e.g., rhodolith bed surface (1 m depth) and midwater (7 m depth), inner bay surface (1 m depth)) is provided in Supplementary Material ( Supplementary Table SM_1 ).

The seasonal averages and standard deviations reported below were calculated from the raw data presented in Table 1 and Supplementary Table SM_1 . For clarity and synthesis, values were grouped by site and climatic period (upwelling, non-upwelling, and transitional seasons), to better describe temporal variability across the study period.

Temporal variations were observed in the main physicochemical parameters and carbonate system variables in Gairaca Bay from April 2023 to July 2024. Total alkalinity (TA) exhibited higher values during the upwelling season in both years, reaching 2427.95 ± 29.4 μmol kg^-1 in 2023 at the rhodolith bed bottom ( Table 1 ) ( Figure 4 ). In 2024, the highest TA was recorded at the outer bay, with 2377.46 ± 35.2 μmol kg^-1 ( Supplementary Figure SM_7A ). In contrast, the lowest TA values occurred during the non-upwelling season at the inner bay (1 m depth) in 2023 (2316.3 ± 37.0 μmol kg^-1), reflecting the influence of less alkaline surface waters over sandy bottoms ( Supplementary Figure SM_5A ).

The total pH (pH_T) showed a seasonal pattern consistent with TA fluctuations. The highest values were recorded during the non-upwelling season in 2023 at the rhodolith bed midwater (7 m depth), averaging 7.99 ± 0.03 ( Supplementary Table SM_1 ; Supplementary Figure SM_4 ). A decrease in pH_T was observed during the upwelling season, particularly in 2024, when the lowest value (7.94 ± 0.05) was recorded at the inner bay at 6 m depth ( Table 1 ; Supplementary Figure SM_6 ). However, overall variations in this variable across sites and seasons were relatively small, indicating a general stability of the carbonate system’s pH throughout the study period.

Dissolved inorganic carbon (DIC) showed the highest concentrations during the pre-upwelling transition in 2023 at the outer bay (2134.05 ± 46.5 μmol kg^-1), and during upwelling at the rhodolith bed bottom in both years (2122.60 ± 40.7 μmol kg^-1 in 2023 and 2109.86 ± 40.9 μmol kg^-1 in 2024) ( Table 1 ).

Salinity exhibited both spatial and interannual variability. In 2023, the lowest salinity values were recorded at the inner bay (1 m depth) during the non-upwelling season (33.36 ± 0.29) ( Supplementary Table SM_1 ; Supplementary Figure SM_5C ). In contrast, in 2024, salinity increased significantly at the rhodolith bed midwater, reaching 35.17 ± 0.97 during the post-upwelling transition season ( Supplementary Table SM_1 ; Supplementary Figure SM_4C ).

The aragonite saturation state (Ω_arag) exhibited a seasonal pattern similar to pH_T. In 2023, the highest Ω_arag values were recorded at the rhodolith bed surface during the non-upwelling season (3.65 ± 0.17) ( Supplementary Table SM_1 ; Supplementary Figure SM_3E ). However, a general decline in Ω_arag was observed across all sites in 2024, with the most pronounced decrease at the inner bay (6 m depth) during the upwelling season (average 3.13 ± 0.31), representing the lowest value recorded during the study ( Table 1 ; Supplementary Figure SM_6E ).

The carbonate ion concentration (CO₃ ^2-) showed pronounced seasonal and between year variability. In 2023, the highest average CO₃ ^2- concentration was recorded at the inner Bay (6 m depth) during upwelling (274.54 ± 22.1 µmol kg^-1) ( Supplementary Figure SM_6F ). Conversely, in 2024 at the same site and season, the lowest values were observed, averaging 138.59 ± 19.5 µmol kg^-1. Overall, CO₃ ^2- concentrations tended to be higher during non-upwelling and transitional periods, and lower during upwelling events, particularly in 2024 across most sites ( Table 1 ; Supplementary Table SM_1 ).

Seasonal deltas, calculated as deviations from the mean, revealed interannual and spatial variability in carbonate system variables across the three sites: rhodolith bed bottom (15 m depth) ( Figure 5 ), inner bay (6 m depth) ( Figure 6 ), and outer bay (10 m depth) ( Figure 7 ). In 2023, the strongest deviations were observed at rhodolith bed bottom, ΔTA reaching a maximum delta of +92.5 µmol kg^-1 during upwelling and a minimum of –89.5 µmol kg^-1 during non-upwelling ( Figure 5A ). Similar but less pronounced patterns were detected in 2024, with ΔTA ranging from +58.2 to –20.9 µmol kg^-1. At inner bay, ΔTA in 2023 ranged approximately from –27.7 to +88.7 µmol kg^-1, while in 2024 the spread was narrower ( Figure 6A ). Outer bay showed more moderate fluctuations across both years, with most ΔTA values remaining within approximately ±45 µmol kg^-1 (range: −71 to 85 µmol kg^-1) ( Figure 7A ).

ΔDIC at rhodolith bed bottom showed the widest range in 2023, peaking at +88.3 µmol kg^-1 during the transition pre-upwelling and dropping to –130.8 µmol kg^-1 during non-upwelling. In 2024, ΔDIC remained high, with values between +67.4 and –74.1 µmol kg^-1 ( Figure 5B ). At inner bay, ΔDIC in 2023 followed a similar pattern but did not reach the same extremes, with the most negative value (−115.89 µmol kg^-1) recorded in October during the non-upwelling season. In 2024, the variation decreased overall, although a maximum positive ΔDIC of 105.62 µmol kg^-1 was observed during the upwelling season ( Figure 6B ).

The general behavior of salinity across the three sites shows a decreasing trend during the non-upwelling season, especially at outer bay and rhodolith bed bottom. In 2023, rhodolith bed bottom showed the strongest salinity deviations, reaching –1.08 during the wet season and +0.42 during the transition post and pre-upwelling ( Figure 5C ). These patterns persisted in 2024 with similar magnitude. At inner bay, the salinity deltas show a contrasting pattern between 2023 and 2024. In 2023, there was a slight decreasing trend in salinity (mean of -0.25), with a higher proportion of negative values. In 2024, however, a clear increasing trend was observed (mean of 0.41) ( Figure 6C ).

At the rhodolith bed bottom, ΔpH_T values were generally negative during upwelling and positive during non-upwelling. In 2023, they ranged from −0.06 to +0.06, with peaks during the non-upwelling season. In 2024, the largest deviations occurred during the transition post-upwelling (−0.07 to +0.03), while changes during upwelling remained moderate (−0.04 to +0.02) ( Figure 5D ). At inner Bay, ΔpH_T values in 2023 were mostly positive across all seasons, reaching up to +0.18 during upwelling and +0.08 during non-upwelling. However, negative values were observed during the transition pre-upwelling (as low as −0.09). In contrast, in 2024, ΔpH_T values at inner bay were predominantly negative, particularly during the upwelling season, reaching as low as −0.19. During the transition post-upwelling, fluctuations were more moderate, ranging from −0.06 to +0.05 ( Figure 6D ). At outer bay, ΔpH_T values in 2023 were mostly positive, with peaks during the transition post-upwelling (+0.06) and non-upwelling (+0.05) periods. In 2024, however, values exhibited greater variability, with a pronounced minimum of −0.09 during the transition post-upwelling and a maximum of +0.05 during upwelling. Overall, both positive and negative ΔpH_T values remained within ±0.1 units across years ( Figure 7D ).

Ω_arag deltas were generally lowest during upwelling and highest during non-upwelling across all sites. At the rhodolith bed bottom, 2023 showed the largest amplitude, with deltas ranging from −0.36 to +0.52, while 2024 displayed a slightly narrower and more negative range (−0.38 to +0.35) ( Figure 5E ). Inner bay exhibited a different pattern, with much greater variability. In 2023, deltas ranged widely from −0.71 to +1.09, and in 2024 the extremes were even more pronounced (−1.08 to +0.28), indicating fluctuations in carbonate saturation ( Figure 6E ). At outer bay, Ω_arag remained relatively stable in both years, with deltas mostly contained between −0.44 and +0.38 in 2023 and between −0.39 and +0.19 in 2024 ( Figure 7E ).

Although no statistically significant differences were found between sites (PERMANOVA: R² = 0.04, F = 0.76, p = 0.66), the outer bay exhibited the greatest changes in DIC (Δ ≈ 229 µmol kg^-1) and salinity (Δ = 3.1), whereas the inner bay showed the highest variability in pH_T (Δ ≈ 0.37) and aragonite saturation state (ΔΩ_arag ≈ 2.18). Meanwhile, the rhodolith bed bottom recorded the largest variation in total alkalinity (ΔTA ≈ 182 µmol kg^-1). Overall, DIC was the most variable parameter across all sites.

PERMANOVA indicated no significant differences in carbonate system composition among sites or depths (p > 0.05). However, seasonal differences were detected between the non-upwelling and upwelling periods. The transition pre-upwelling season closely resembled the upwelling period. Significant seasonal variation was also observed for temperature (F = 248.42, p < 0.05), salinity (F = 49.02, p < 0.05), TA (F = 7.65, p < 0.00), and DIC (F = 2.54, p < 0.00). The SIMPER analysis identified the carbonate system variables contributing most significantly to seasonal dissimilarities. In the upwelling versus non-upwelling contrast, DIC (38.9%) and TA (28.6%) were the major contributors to observed differences, both statistically significant (p ≤ 0.05). Other variables such as salinity, aragonite saturation state (Ω_arag), and calcite saturation state (Ω_cal) contributed less than 0.4% ( Table 2 ).

In the transition post-upwelling vs. non-upwelling comparison, DIC accounted for 36.9% of the dissimilarity (p = 0.03), while salinity contributed only 0.7% (p ≤ 0.00), despite its low explanatory power. In the non-upwelling versus transition pre-upwelling contrast, DIC again emerged as the dominant factor (39.0%, p ≤ 0.00), followed by pCO₂ (34.6%, p = 0.01) and TA (25.6%, p = 0.00). Although salinity, Ω_cal, Ω_arag, and pH_T contributed to the differences, their individual contributions were relatively small ( Table 2 ). These results highlight the importance of DIC and TA in the seasonal differentiation of the carbonate system. However, the interaction between different variables can further modulate the variability of the carbonate system.

The variability of the carbonate system cannot be explained solely by individual factors. Temperature explained 39% of DIC variability, showing a statistically significant relationship (R² = 0.39, p < 0.001) ( Figure 8A ). In contrast, no significant relationship was found between temperature and pCO₂ (R² = 0.02, p = 0.32), indicating that temperature alone was not a good predictor of pCO₂ variability in the study area ( Table 3 ).

When considering additional hydrochemical variables, the multiple regression model achieved a markedly improved fit for pCO₂ (adjusted R² = 0.96, RMSE = 18.1 µatm) ( Table 4 ; Figure 8B ). In this model, temperature, salinity, and the DIC/TA ratio emerged as significant predictors, whereas total alkalinity (TA) and DIC alone did not significantly contribute to the explained variance.

At the rhodolith bed bottom site, the carbonate system parameters and salinity exhibited relatively stable average values, with notable seasonal variability ( Figure 4 ). TA had a mean value of 2384.66 ± 49.87 µmol kg^-1, ranging from 2287.30 to 2469.30 µmol kg^-1 ( Figure 4A ), DIC averaged 2097.37 ± 60.64 µmol kg^-1, with a range between 1977.78 and 2196.89 µmol kg^-1 ( Figure 4B ) and salinity averaged 34.58 ± 0.73, fluctuating between 33.5 and 35.9 ( Figure 4C ). The pH_T showed limited variability across seasons, with an overall mean of 7.96 ± 0.04. The maximum value (8.02) was recorded during the non-upwelling season in September 2023, while the minimum (7.89) occurred in July 2024 during the post-upwelling transition ( Figure 4D ). The Ω_arag averaged 3.36 ± 0.28 ( Figure 4E ). The CO₃ ^2- had a mean value of 211.38 ± 13.28 µmol kg^-1, the highest concentration (235.66 µmol kg^-1) was observed in September 2023 during the non-upwelling period, coinciding with elevated Ω_arag and lower DIC values, while the lowest value (185.70 µmol kg^-1) was observed in December 2023 during the pre-upwelling transition, when pCO₂ peaked and Ω_arag declined ( Figure 4F ).

PERMANOVA revealed statistically significant differences among seasons (F = 3.0596; p = 0.01). SIMPER analysis further identified the variables contributing most to these seasonal dissimilarities at the rhodolith bed bottom site. The transition post-upwelling vs. transition pre-upwelling comparison exhibited the strongest dissimilarities, with DIC (41.4%, p = 0.04) and TA (92.7%, p = 0.03) showing statistically significant differences. In the upwelling vs. transition post-upwelling comparison, the major contributors were Ω_arag (34.7%) and Ω_calc (30.1%), followed by DIC (18.8%) and TA (6.7%) with no significant differences (p > 0.1). In the upwelling vs. non-upwelling contrast, DIC (19.6%), Ω_arag (19.1%), and Ω_calc (18.3%) contributed most, followed by TA (16.8%), salinity (12.3%), and pCO₂ (12.1%) with no statistically significant differences detected (p > 0.75).

Similarly, in the upwelling vs. transition pre-upwelling comparison, Ω_calc (31.7%), Ω_arag (29.4%), and pCO₂ (27.2%) were the dominant contributors without significant differences. The transition post-upwelling vs. non-upwelling comparison revealed salinity as the main contributor (21.6%), with no significant differences (p > 0.30), likewise, no significant differences were found between non-upwelling and transition pre-upwelling periods, despite salinity accounting for 100% of the observed dissimilarity (p = 0.80).

The analysis of nutrient concentrations revealed clear seasonal and between year variation. Maximum nitrate concentrations were recorded in June 2023, reaching 0.08 mg·L^-1. Regarding nitrite, the highest values were observed in July 2023 and April 2024, while the lowest concentration was recorded in December 2023 (0.03 mg·L^-1). For ammonium, the highest concentration was detected in August 2023 (0.01 mg·L^-1); however, from June to December 2023 and throughout 2024, ammonium levels were consistently below the detection limit. Phosphate concentrations also varied notably between years. In 2023, the highest concentration occurred in September (0.31 mg·L^-1), whereas the lowest was measured in December (0.22 mg·L^-1). In contrast, 2024 exhibited a significant increase, with a peak concentration of 0.39 mg·L^-1 recorded in May; in all other sampled months, phosphate concentrations were below the detection limit ( Table 5 ). Among the nutrient variables, a positive correlation was found between nitrite (NO₂ ^-) and both CO₃ ^2- (r = 0.59, p = 0.02) and Ω_arag (r = 0.54, p = 0.04). Additionally, nitrate (NO₃ ^-) was positively correlated with phosphate (PO₄ ^3-) (r = 0.53, p = 0.04).

Correlations between carbonate chemistry variables and nutrient concentrations at the rhodolith bed bottom site revealed strong internal coherence within the carbonate system, but weak associations with nutrient dynamics ( Figure 9 ). Total alkalinity (TA) showed a strong positive correlation with dissolved inorganic carbon (DIC) (ρ = 0.91, p < 0.00), while DIC was also positively correlated with pCO₂ (ρ = 0.67, p = 0.00) and negatively with pH_T (ρ = –0.51, p = 0.04) and CO₃ ^2- (ρ = –0.51, p = 0.04). Aragonite saturation state proxies (pH_T and CO₃ ^2-) were highly correlated with each other (ρ = 0.89, p < 0.00), and both were strongly and negatively correlated with pCO₂ (ρ = –0.96 and –0.89, respectively, p < 0.00). Among the nutrients, only nitrite (NO₂ ^-) showed a significant positive correlation with CO₃ ^2- (ρ = 0.59, p = 0.02). All other correlations between nutrients (NH₄ ⁺, NO₃ ^-, PO₄ ^3-) and carbonate system variables were not statistically significant (p > 0.05). Temperature also exhibited significant relationships with carbonate parameters. It was negatively correlated with TA (ρ = –0.68, p = 0.00) and DIC (ρ = –0.74, p = 0.00), and positively correlated with Ω_arag (ρ = 0.54, p = 0.02) ( Figure 9 ).

This study highlights the significant influence of seasonal variability on carbonate chemistry in tropical coastal ecosystems, particularly those affected by upwelling. Our results show a clear connection between the non-upwelling season, characterized by the rainy period and increased runoff from the Magdalena River (Ricaurte-Villota et al., 2025), and significant fluctuations in the carbonate system parameters in Gairaca Bay. These patterns are comparable to those reported for Moorea reef flats (Kleypas et al., 2011) and suggest that continental runoff plays a central role in modulating water chemistry in Gairaca, setting it apart from other tropical coastal systems where oceanic upwelling exerts a more dominant influence (Sánchez-Noguera et al., 2018). Notably, freshwater inputs contribute substantially to the spatial variability of carbonate chemistry, reinforcing the importance of accounting for runoff in coastal biogeochemical assessments.

During upwelling peaks, DIC and TA increase, while pH_T drops below 7.95 and Ω_arag decreases to ~3.0, similar to conditions observed in the Gulf of Papagayo (Sánchez-Noguera et al., 2018). In contrast, the non-upwelling season reflects a diminished influence of cold, CO₂-rich subsurface waters (Ricaurte-Villota et al., 2025), with pH_T rising from 7.93–7.99 (upwelling) to 8.01–8.03 at sites such as the rhodolith bed surface and inner bay ( Table 1 ; Supplementary Table SM_1 ). This increase in pH_T is accompanied by a decrease in DIC, which drops from 2115.77 µmol kg^-1 in April 2023 (upwelling, inner bay) to around 1973.10 µmol kg^-1 in October ( Table 1 ). Similarly, TA shows a slight decline, such as at rhodolith bed surface, where it decreases from 2487.40 µmol kg^-1 in March to 2268.00 µmol kg^-1 in October ( Supplementary Table SM_1 ). These changes are coupled with a decline in DIC, from 2115.77 µmol kg^-1 in April (inner bay) to ~1973.10 µmol kg^-1 in October, and a moderate decrease in TA, for example from 2487.40 to 2268.00 µmol kg^-1 at the rhodolith bed surface. Conversely, Ω_arag and CO₃ concentrations increase during this period, reaching up to 3.88 and 235.66 µmol kg^-1 at the rhodolith bed bottom, respectively.

Salinity also decreases (e.g., from 34.4 in March to 33.2 in October), while temperature rises to 30.5 °C, compared to 27.0–27.5 °C during upwelling. Dissolved oxygen concentrations tend to decline, as observed at the rhodolith bed bottom (from 7.04 to 4.17 mg·L^-1). These warmer, fresher waters exhibit higher pH_T and Ω_arag, which thermodynamically lowers the energy barrier for CaCO₃ precipitation (Mucci, 1983; Cyronak et al., 2016). Carbonate speciation under elevated pH shifts the equilibrium toward CO₃ ^2-, enhancing local carbonate ion availability (Dickson and Millero, 1987).

The freshwater influx into Gairaca Bay during the rainy season results in a distinct carbonate chemistry response, differentiating it from systems like Papagayo, dominated by oceanic upwelling with minimal freshwater influence (Sánchez-Noguera et al., 2018), and Bocas del Toro, which experiences moderate terrestrial runoff (Pedersen et al., 2024). Gairaca is subjected to pronounced seasonal and interannual variability in runoff, particularly from the Magdalena River during the rainy, non-upwelling season. This input amplifies fluctuations in pH_T, TA, DIC, and Ω_arag, especially during periods of intense rainfall and river discharge, when upwelling is absent ( Table 1 , Figure 4 ). Additionally, the Magdalena River transports organic matter and nutrients (Restrepo et al., 2006), whose remineralization can further alter DIC and pH_T dynamics. Maximum runoff typically occurs from September to November, coinciding with weakened trade winds, enhanced coastal countercurrents, and the suppression of upwelling, thereby strongly influencing the carbonate system (Ricaurte-Villota et al., 2025). Although spatial variability in water chemistry exists within Bocas del Toro due to terrestrial runoff and benthic metabolism (Pedersen et al., 2024), Gairaca shows notably higher temporal variability driven by the strong seasonal freshwater inputs during the non-upwelling season. Nonetheless, the three sites evaluated within Gairaca Bay exhibit similar responses to these runoff and rainy conditions.

Spatial differences among the rhodolith bed bottom, inner bay, and outer bay sites were minimal, likely due to a generally well-mixed water column (see Section 4.3). However, the pronounced seasonal shifts highlight the interplay between oceanic and terrestrial influences, which may enable site-specific buffering mechanisms, such as localized photosynthesis or sediment-driven alkalinity release (Savoie et al., 2022; Ricaurte-Villota et al., 2025). These seasonal shifts in carbonate chemistry not only reflect oceanographic-terrestrial interactions but may also have significant implications for the physiological performance of calcifying organisms (Li et al., 2022).

Finally, it is important to consider the analytical uncertainty associated with TA and DIC measurements in this study, which is estimated at ±10 µmol kg^-1. This level of precision corresponds to the “weather goal” defined by the Global Ocean Acidification Observing Network (GOA-ON), which is considered adequate for characterizing short-term variability and spatial gradients in coastal systems (Kortazar et al., 2020). Although not suitable for detecting long-term anthropogenic trends, this uncertainty is acceptable in highly dynamic environments like Gairaca Bay, where variability in TA and DIC often exceeded 100 µmol kg^-1, ensuring that the observed patterns remain robust and interpretable.

During the study period, rainy days were recorded even during months that are typically dry, resulting in unusual hydrological conditions for that time of year ( Figure 3 ). These unexpected rainfall events likely intensified freshwater inputs, leading to dilution effects, changes in salinity and alkalinity, and potential decoupling among carbonate system parameters. Elevated runoff may also have altered water column structure and enhanced biogeochemical fluxes, contributing to the variability observed in carbonate chemistry (Correa-Ramirez et al., 2020; Norzagaray et al., 2020; Cai et al., 2021; Reithmaier et al., 2023).

During the non-upwelling months, which coincide with the rainy season, Gairaca Bay exhibited reduced salinity, lower TA and DIC concentrations, and elevated pCO₂ levels. These shifts likely resulted from increased remineralization and microbial respiration, stimulated by the influx of terrestrial organic matter during intense rainfall. In coastal ecosystems, remineralization refers to the breakdown of organic matter into inorganic constituents. This process, particularly when driven by microbial activity, generates CO₂ and modifies porewater chemistry (Bayraktarov and Wild, 2014; Quintana et al., 2015; Cohn et al., 2024). However, the observed decrease in TA during the rainy season suggests that the dilution effect of freshwater inputs, which are typically low in TA, may outweigh any potential increase in TA from remineralization processes (Pedersen et al., 2024).

Terrestrial–marine interactions, especially those involving organic matter inputs from river discharge, mangroves, and seagrass meadows, are known to modulate carbonate chemistry in Caribbean coastal systems by altering TA and DIC concentrations (Meléndez et al., 2020; Pedersen et al., 2024). These biogeochemical dynamics are ecologically relevant for calcifying organisms such as corals and coralline algae, which are particularly sensitive to fluctuations in pH and carbonate saturation state (Martin and Hall-Spencer, 2017; Silbiger and Sorte, 2018). Increased freshwater runoff can suppress calcification and alter benthic community composition (Fabricius, 2005). However, rhodolith beds may partially buffer these effects by maintaining relatively stable micro-environmental conditions through photosynthesis, calcification, and microbial mediation (Isah et al., 2022).

Delta patterns revealed that climatic seasonality is a major driver of carbonate system variability across Gairaca Bay. Although statistically significant differences between sites were not detected, the magnitude of variation in parameters such as TA, DIC, Ω_arag, and pH_T, differed notably across seasons ( Figures 5 - 7 ).

During the non-upwelling season, the variability between sampling locations was markedly higher for all measured parameters. TA showed the greatest inter-site difference, with a 90.7% relative difference, followed by pH_T (83.2%), DIC (80.5%), and Ω_arag (68.7%). These pronounced variations suggest that localized processes, such as freshwater input, biological activity, and water column stratification, exert greater influence under low-mixing conditions when upwelling is absent. The rhodolith bed bottom recorded the highest TA variability, while the outer bay exhibited the greatest fluctuations in DIC and pH_T. Conversely, the inner bay consistently showed lower variability across most parameters during this period. These findings highlight the complex interplay of regional and local factors in shaping carbonate system dynamics in tropical coastal environments (Sánchez-Noguera et al., 2018; Norzagaray et al., 2020).

Consistent with these observations, temperature emerged as a key driver of DIC variability, due to its effects on solubility and biologically mediated processes such as calcification and respiration. Warmer temperatures generally reduce DIC solubility (Bakker et al., 1999), while biological processes are also temperature-dependent: enhanced respiration increases DIC concentrations, and decreased calcification reduces DIC uptake (Pedersen et al., 2024). In contrast, pCO₂ was influenced by multiple interacting variables. The strong performance of the multiple linear model (R² = 0.96; RMSE ≈ 18 µatm) indicates that including salinity and the DIC/TA ratio significantly improves prediction accuracy, capturing the combined influence of freshwater input, mixing, and net community metabolism. Such complexity is typical in estuarine and coastal systems, where biogeochemical and physical drivers interact dynamically (Cai et al., 2021).

These findings emphasize that, while DIC dynamics are influenced by temperature, pCO₂ is shaped by a broader suite of environmental factors, including carbonate equilibrium and the balance between DIC and TA, which determines buffering capacity and resistance to pH changes (Khan et al., 2020). Understanding these interactions is crucial for predicting the impacts of environmental variability on coastal carbonate chemistry (Carstensen and Duarte, 2019).

Observed reductions in salinity, TA, and DIC during the non-upwelling season support the central role of freshwater inputs in modulating carbonate conditions (Pérez et al., 2015). Biological processes, particularly photosynthesis, may further contribute to DIC drawdown (Isah et al., 2022), while decoupling between production and respiration can amplify pH variability in stratified coastal waters (Carstensen and Duarte, 2019).

In contrast, during the upwelling season, spatial differences across sampling sites diminished considerably. Relative differences in pH_T (48.9%), TA (42.9%), DIC (20.0%), and Ω_arag (18.9%) were all lower, suggesting a homogenizing effect of upwelling, driven by the intrusion of cold, CO₂-rich subsurface and enhanced vertical mixing across the water column, which together affect the entire bay. This influence was most pronounced at the outer bay, where pH_T and Ω_arag variability peaked, while the rhodolith bed bottom and inner bay showed more moderate changes. The predominance of northeasterly and easterly winds during this season plays a key role in sustaining upwelling-favorable conditions, enhancing the upward advection of DIC-enriched subsurface water masses. This mechanism likely contributes to the observed increase in DIC concentrations during upwelling, compared to the more variable and stratified conditions of the non-upwelling season. The stronger divergence observed during the non-upwelling period may be further intensified by unusually high rainfall and riverine discharge (Restrepo and Kjerfve, 2000; Ricaurte-Villota et al., 2025), which enhance the effects of local processes under stratified conditions and limited vertical mixing (Pedersen et al., 2024). On average, Ω_arag remained above the aragonite saturation threshold (>1) at all sites (see Supplementary Table SM_1 ), indicating generally favorable conditions for marine calcification (McGrath et al., 2019). Slightly higher values were observed at the outer bay, followed by the rhodolith bed bottom and inner bay.

Rhodolith bed bottom site exhibited the most stable Ω_arag values across seasons, while outer and inner bay sites showed greater variability. This stability supports the buffering potential of rhodolith habitats, consistent with previous findings that suggest rhodolith beds may act as microhabitat refugia under ocean acidification scenarios (Costa et al., 2023). Their capacity to maintain stable chemical conditions through biological and sedimentary processes may offer protection to vulnerable calcifiers. In contrast, shallow inner bay areas, dominated by sandy bottoms, are more susceptible to fluctuations in carbonate availability. These results highlight the interplay between regional climatic drivers and local environmental features in shaping the chemical mosaic of tropical marine systems (Gómez et al., 2023).

Although direct measurements of metabolic or calcification rates were not conducted, the consistently elevated and stable Ω_arag values at the rhodolith bed suggest the presence of biologically mediated buffering. While average values across sites and depths were relatively similar, SIMPER analysis revealed that variables such as TA and DIC contributed significantly to seasonal dissimilarities at the rhodolith bed bottom, particularly between transition phases. These patterns support the notion that, despite limited spatial contrast in mean values, rhodolith beds may play a role in modulating carbonate chemistry under fluctuating conditions, especially during periods of environmental stress such as upwelling pulses or freshwater inputs (Pedersen et al., 2024; Ricaurte-Villota et al., 2025).

Nutrient dynamics at the rhodolith bed bottom revealed clear seasonal trends driven by upwelling, freshwater inputs, and biological processes. Nitrate concentrations peaked during June and July 2023 (0.07 and 0.06 mg·L^-1 respectively), coinciding with the post-upwelling transition. These peaks likely reflect enhanced remineralization following the intrusion of subsurface waters (Schubert et al., 2019; Ricaurte-Villota et al., 2025). In contrast, ammonium remained largely undetectable throughout 2024, while phosphate concentrations increased notably in May 2024 (0.39 mg·L^-1), coinciding with the onset of the rainy season and heightened river discharge.

Interestingly, Ω_arag showed a positive correlation with NO₂ ^- y PO₄ ^3-. While this relationship is not necessarily causal, it may reflect the combined influence of terrestrial runoff and in situ biogeochemical processes in this tropical coastal environment. Seasonal runoff in Gairaca Bay, particularly during the rainy season, introduces nutrients and organic matter that can stimulate biological activity and indirectly affect carbonate chemistry (Aronson et al., 2014). In coastal waters, nutrient enrichment and remineralization can influence pH and DIC, sometimes resulting in complex or even counterintuitive correlations with Ω_arag (Cai et al., 2021). Therefore, the observed positive correlation may be the result of overlapping hydrological and biological processes rather than a direct mechanistic link.

The carbonate system at the rhodolith bed is shaped by a complex interplay between pH_T, DIC, Ω_arag, and nutrient concentrations. A significant inverse relationship was found between DIC and both pH_T and Ω_arag, consistent with acid-base dynamics in marine systems: as DIC increases, pH_T and Ω_arag decrease. This pattern aligns with observations in other tropical regions, such as the northern South China Sea, where similar relationships were reported by Roberts et al. (2021), highlighting the widespread influence of DIC on carbonate chemistry across diverse marine environments.

The positive correlation observed between temperature and aragonite saturation state (Ω_arag) can be attributed to the seasonal dynamics of the studied tropical coastal system. During the rainy season, surface water temperatures rise while upwelling intensity decreases, reducing the input of cold, CO₂-rich subsurface waters and thereby limiting acidification from vertical mixing. Consequently, pH_T and carbonate ion (CO₃ ^2-) concentrations increase, resulting in elevated Ω_arag values. This pattern, characteristic of shallow tropical environments with seasonal forcing, contrasts with systems dominated by persistent upwelling, where the influx of subsurface waters typically lowers Ω_arag despite cooler temperatures (Mucci, 1983; Zeebe and Wolf-Gladrow, 2001; Manzello, 2010). However, the hysteresis effect described by (McMahon et al., 2013) underscores the importance of considering diel variability and metabolic feedbacks when interpreting correlations between Ω_arag and environmental drivers such as temperature.

Nevertheless, most correlations between pH_T or DIC and nutrient concentrations were not statistically significant, indicating that nutrient variability does not directly control carbonate system dynamics in this setting. Although nutrients are essential for biological productivity, their short-term influence on pH_T and DIC appears less pronounced than other drivers such as temperature, salinity, or air–sea CO₂ exchange (Gattuso et al., 1998; Zeebe and Wolf-Gladrow, 2001). In tropical coastal systems, these physical and chemical factors often dominate carbonate chemistry, especially under stratified or runoff-influenced conditions (Salisbury et al., 2008; Cai et al., 2021).

This nutrient enrichment pattern supports the role of freshwater inputs in modulating carbonate dynamics in Gairaca Bay (Ricaurte-Villota et al., 2025). Riverine waters are typically low in TA and DIC but high in nutrients and temperature, promoting seawater dilution and enhanced biological CO₂ uptake through primary production (Borges and Gypens, 2010). Such correlations likely reflect the contribution of remineralization and biologically driven CO₂ uptake in shaping local carbonate dynamics (Borges and Gypens, 2010). These interactions likely explain some of the shifts in carbonate chemistry observed during transition periods, when freshwater delivery and biological activity are both elevated.

Although freshwater inputs from rainfall and runoff are well-documented drivers of biogeochemical variability in this region (Ricaurte-Villota et al., 2025) our data suggest that internal metabolic activity within rhodolith beds may help maintain conditions favorable to calcification, even under elevated pCO₂ and DIC levels. Future studies should investigate diel metabolic variability and long-term biogeochemical trends to better understand the mechanisms sustaining rhodolith bed resilience under changing oceanographic conditions.