Seven stations were sampled and used to estimate CO₂ fluxes and carbonate system chemistry in Gorgona National Natural Park from the coast to the open ocean. The closest station to the coast was “Río Guapi - RG”, located approximately 25 km from Guapi-Nariño, Colombia, and the furthest was “North Ocean - ON”, approximately 60 km from the coastline. The stations were distributed along two parallel transects separated by five kilometers, with GPS-WGS84 coordinates. The southern transect was formed by “South Guapi - GS”, “South Reef - AS”, and “South Ocean - OS”, while the northern transect was formed by “North Guapi - GN”, “North Reef - AN”, and “North Ocean - ON” (Figure 1). Monitoring began in November 2021 and extended until July 2022 (nine months), with monthly sampling during the first week of each month. In addition, the “Río Guapi” station was added in January 2022 to characterize its influence.

On site, salinity, temperature, pH millivolts (total protons), percent oxygen, and depth were measured using a Hanna multiparameter sounder (previously calibrated), a CTD-Castaway, and a YSI equipment (accurate to ±2% of reading or ±0.01 PSU, whichever is greater; ± 0.15 °C; ± 0.02 pH/± 0.5 mV; and ranging from 0.00 to 30.00 ppm (mg/L) with an accuracy of ±0.25% of full range.

Discrete water samples were collected at seven stations (GS, GN, AS, AN, OS, ON, and RG; Figure 1) using a 5 L Niskin bottle at depths between 1.8 and 3 meters. The water collected on site was not murky, sampling sites were right after the river plume limit (35 km away from the continent). Seawater samples were stored in 250 mL borosilicate bottles and preserved with 50 µL of a saturated mercuric chloride (HgCl₂) solution for later TA analysis in the laboratory, following recommended best practices (Best Practice Guide for Ocean CO₂ Measurements, SOP 2, SOP-specified concentration range of 0.02-0.05%). Additionally, DIC samples were taken in 50 mL dark bottles (silicone seal), leaving no headspace; and wrapping the cap in parafilm paper. All samples were packed in a Styrofoam refrigerator and sent by plane to our new lab in Bogotá (Javeriana University) in less than four hours from the Pacific region. Then, the samples were kept without exposing them to light and maintaining temperature conditions below 15 °C in the laboratory. TA was measured within two weeks, as well as other parameters such as pH. Regarding DIC samples, those were kept in our lab for two months after the sampling was finished and then sent to the Autonomous University of Baja California in Mexico. They were processed in a maximum of one week after arrival in the UABC.

In the laboratory, total alkalinity (TA) measurements were performed on seawater samples following a standardized protocol (SOP 3b of the open-cell titration method; Dickson et al., 2007) using the GOA-ON titration kit. For titrant, Dickson Lab reference HCl solutions were used (Batch A17: 0.100362 ± 0.000009 mol kg^-1; Batch A24: 0.099922 ± 0.000005 mol kg^-1), buffered with 0.6 M NaCl to match seawater ionic strength. A 50 g sample (~50 g or mL, ± 0.0001 mL) was weighed into a jacketed titration cell (maintained at ~25 °C via a circulating water bath) and titrated. The analytical balance used (Ohaus Pioneer AX225D, ± 0.0001 g) yielded a reproducibility of ±0.001 g for 50 g samples, corresponding to a coefficient of variation<0.003%, well within the precision requirements for TA analysis (0.01 g). The titration involved two stages: first, the addition of a small amount of acid (0.2 mL) with a HandyStep^® Touch S, first by adding to drop the pH to ~3.5 (~0.2100 mV) and waiting 5 minutes to allow for degassing, and second, titration resumed in 0.05 mL steps until a pH close to 3.0 was reached (~0.2300 mV), with pH and temperature measurements taken at each step. Data, including sample volume, titrant volumes (to pH 3.5 and 3.0), salinity, final temperature, and voltage per addition was compiled in an Excel spreadsheet program, which employed a non-linear least squares regression to calculate total alkalinity (Dickson et al., 2007). Approximately 50% of surface samples (n = 31) were analyzed in triplicate. Only triplicate sets with a maximal variation of ±5 µmol kg^-1 TA were retained (± 9 μmol kg^-1, Batch #182: 2230.91 ± 0.71 µmol kg^-1; Batch # 202: 2215.13 ± 0.57 µmol kg^-1; all data set available at NOAA; https://doi.org/10.25921/gfan-3e30).

Dissolved Inorganic Carbon (DIC) was measured at UABC (Mexico) using the UIC C-CM5014 instrument, following the coulometric methodology described by Johnson et al. (1987) and Dickson et al. (2007). Measurements were performed with high precision and accuracy, ensuring that the difference in DIC from the reference values ​​of the standards did not exceed 3 μmol kg^-1 (error of 0.1%). The Dickson Certified Reference Material (CRM) was prepared at the Institute of Oceanography, University of California, San Diego (USA). The substandards and analyzed samples were carried out in parallel at the Oceanographic Research Institute of the Autonomous University of Baja California (Mexico). For a more comprehensive DIC measurement methodology, see Supplementary Table 1.

Similarly, quality “flags” (scale 1 to 5) were reviewed and assigned to the TA and DIC data (equipment accuracy and potential errors in the field or laboratory). Outliers (n=6) were removed from the matrix; this was done by reviewing multiple graphical relationships between temperature, salinity, TA, and DIC.

The Liss and Merlivat (1986) equation allowed for the calculation of the flux magnitude for each of the nine months and seasons in Gorgona. The estimated CO₂ flux in Gorgona was compared with the IPACOA, OOI, NOAA, and SOCAT databases.

The CO₂ fluxes (FCO₂ direction, magnitude, and variability) were calculated following the equation proposed by Liss and Merlivat (1986):

Where, the constant k is the gas transfer coefficient of CO₂ as a function of wind (U₁₀) at 10 meters above sea level (Wanninkhof, 2014), expressed as:

K′ is the solubility of CO₂ and is a function of temperature and salinity (Weiss, 1974); Sc is the Schmidt number and depends on temperature; ΔpCO₂ is the difference between the partial pressure of carbon dioxide (pCO₂) in the surface seawater (pCO_2w) and the pCO₂ of the atmosphere (pCO_2A). ΔpCO₂ determines the direction of the flux, with positive FCO₂ indicating CO₂ moving from the sea to the atmosphere, and negative indicating movement from the atmosphere to the sea. The magnitude of the gas exchange is controlled not only by ΔpCO₂ but also by the gas transfer rate between the two media, which in turn varies with wind speed.

The hourly wind speed from the seven sampling sites - U₁₀ (m s^-1) was extracted from: https://earth.nullschool.net/, taking the 24-hour wind average of the sampling day. For months with sampling of all stations on two consecutive days (adverse weather conditions in two months), 48-hour wind values were averaged. The wind average was also calculated and compared using only the exact sampling hours at each station and by averaging the wind over a month, finding no significant changes in the flow results compared to the daily/monthly wind average (which was finally used).

Since the atmosphere is a well-mixed fluid, the pCO_2A values for the flow calculation were extracted from NOAA/GML https://gml.noaa.gov/ccgg/trends/data.html (NOAA, 2021, 2022), which come from the Mauna Loa climate station, Hawaii, file “Mauna Loa CO₂ monthly mean data.” The pCO_2A data for the sampling months (Nov 2021 to Jul 2022) were extracted.

The pCO_2w was calculated from TA, DIC, salinity, and temperature data resolving the carbonate system equation in seawater using the CO2SYS v3.0 software (Pierrot et al., 2021). In the carbonate system equation solution, the following constants were used: (i) Dissociation constants for K1 and K2 from Millero (2010) for waters ranging from 0 to 40, given that the study area the salinities variation are between 27.03 to 31.29, (ii) KHSO4 dissociation constant from Dickson (1990), (iii) KHF from Perez and Fraga (2003), (iv) Total pH scale (mol-kg SW), (v) [B]T value from Lee et al. (2010), and (vi) EOS-80 standard. The pCO_2w was estimated with an error of ±1.79 concerning the determined fluxes (mmol C m² day^-1). Out of 66 flux data points, six were flagged as possible outliers (flags three onward). The CO2SYS was also used with DIC and TA pairs to estimate the Revelle factor, omega aragonite (Ω_AR), HCO₃, and CO₃. Lastly, the pH in total scale was derived through the CO2SYS using the same specific settings and found to be less variable than in situ measurements, displaying a smaller error. We performed a full uncertainty propagation analysis, which yielded the following average errors: pH_T (± 0.022), pCO_2W (± 30.17 µatm), total alkalinity (TA: ± 9 µmol kg^-1), dissolved inorganic carbon (DIC: ± 3 µmol kg^-1), salinity (± 0.20 PSU), and temperature (± 0.08 °C). These values fall within acceptable ranges for carbonate system studies and support the reliability of our calculated FCO₂ values at Gorgona.

Daily images of Chl-a were downloaded from the Ocean Colour Climate Change Initiative (OC-CCI; version 3.1) for the studied period. The OC–CCI is a merged Level 3 product with a spatial resolution of 4 km, available at http://www.oceancolour.org/. The OC-CCI Chl-a dataset is retrieved by combining the observational data from the MERIS (MEdium spectral Resolution Imaging Spectrometer) sensor of the European Space Agency, the SeaWiFS (Sea-viewing Wide-Field-of-view Sensor) and MODIS-Aqua (Moderate-resolution Imaging Spectroradiometer-Aqua) sensors from the National Aeronautics and Space Administration (NASA-USA), and the VIIRS (Visible and Infrared Imaging Radiometer Suite) from the National Oceanic and Atmospheric Administration (NOAA-USA). The remote sensing reflectance data derived from the sensors were merged by band-shifting and bias-correcting the MERIS, MODIS, and VIIRS data to match the SeaWiFS data. Due to the complexity of the study region, version 3.1 of the OC-CCI product was selected as it improves the performance of the ocean color data in coastal Case-2 waters compared to earlier versions that primarily focus on open ocean waters. Finally, total monthly Chl-a values were obtaining by averaging the Chl-a of the seven sampling stations, resulting in a relative total Chl-a concentration for the entire sampling region.

The analysis began with a Shapiro-Wilk test, which indicated that the flux magnitude data does not follow a normal distribution (S-W, P = 0.01877, W = 0.95, n = 60), repeating this for pCO_2w and ΔpCO₂ components, respectively (S-W, p = 0.147 and 0.134, W = 0.97 and 0.96, n = 60 and 60), indicating that the variables associated with flux do not follow a normal distribution either. Hence, the Box-Cox methodology was applied to transform FCO₂ data into a normal distribution. Furthermore, a Welch F test in the case of unequal variances was applied to all the variables of interest (FCO₂, pCO_2w, ΔpCO₂, SST, SSS, TA, DIC, and Chl-a) which showed that these different data groups have equal variance (W-F, F = 2.512E^-4, df = 200.3, p = 1.855E^-304). Moreover, a t-test (parametric) was used to analyze differences between the two flux periods: Upwelling (January to March of 2022) and Post-upwelling (November to December of 2021 and April to July of 2022). We used a One-way ANOVA to compare different flux magnitudes across the nine months, and then an a posteriori t-test. Additionally, a Spearman linear correlation was run to explore relationships between pCO_2w and SST, SSS, TA, DIC, and Chl-a.

Normalizations to observe the pCO_2w sensitivity were conducted as follows: we used the maximum values for salinity and temperature (32.06 units and 28.1 °C, respectively) setting the maximum value for salinity first, as a constant, while leaving the TA, DIC, and temperature values corresponding to each data point to evaluate the effect of salinity over the pCO_2w, repeating the same process but substituting the temperature for the recorded maximum, and leaving salinity as measured (CO2SYS). Then, we recalculated the pCO_2w and plot the difference between the pCO_2w with in situ salinity against the recalculated pCO_2w that uses the maximum recorded salinity. 49 µatm higher, in average, was the pCO_2w difference by using max salinity effect. Regarding a normalization with maximum temperature, the same procedure was followed where the maximum recorded temperature was constant, and the remaining variables used the corresponding values, revealing an average pCO_2w difference of 31 µatm higher than with the observed temperatures (Supplementary Figure 1) in the CO2SYS (Lewis and Wallace, 1998). Thus, it is clear that salinity has a larger influence (relative percentage) than temperature regarding the pCO_2w variability observed throughout the year. To complement these normalizations, we assessed the uncertainty associated with pCO_2w by incorporating the potential error of ±30.17 μatm into our flux calculations, recalculating them and comparing them to the original estimates.

In addition, to define climatic-oceanographic periods, a Principal Component Analysis (PCA) was performed using the latest version of R Studio software. The PCA used the following variables: month, season, SST, SSS, pCO₂w, pH_T, DIC, Chl-a, to review if there are spatio-temporal FCO₂ groups. Furthermore, to clearly separate both seasons we used several atmospheric/climatic variables, such as total monthly precipitation (mm), average sea surface temperature (SST), average pH (total), average surface water partial CO₂ pressure (pCO_2w), and average CO₂ fluxes (mmol m^-2 d^-1).

The Zonal Ekman Transport (ZET) was obtained following the methodology of Bakun and Nelson (1991) by using weekly wind data obtained from the CCMP product (https://www.remss.com/measurements/ccmp/). For the calculation, we follow the equation:

where ρ_ω is the density of seawater, which is assumed constant at 1025 kg m⁻³; f is the Coriolis parameter and τ is the wind stress for the studied area.

The wind stress (τ) was computed from the weekly wind fields as follows:

where C_d = 0.0015 is the drag coefficient, ρ_a = 1.2 kg m⁻³ is the mean air density, and V₁₀ is the wind speed at 10 m above the sea surface. The wind stress curl was estimated using the zonal and meridional components of the weekly wind stress:

This procedure was performed for each grid point of the wind field by applying the centered finite difference algorithm, and the Ekman pumping velocity was estimated as follows:

where τ_x is the zonal wind stress component, β = (2Ωcosφ/R) is the latitudinal variability of the Coriolis parameter (f = 2Ωsenφ), R = 6,371,000 m is the radius of the earth and ρ_w = 1025 kg m⁻³ is the mean water density. This methodology has been previously used on the Northeast Tropical Pacific (Kessler, 2006) and Panama Bight region (Devis-Morales et al., 2008; Corredor-Acosta et al., 2020).

Previous studies have indicated an interplay between the Ekman transport due to alongshore winds and Ekman pumping due to offshore wind stress curl playing an important role due to their relative contributions for upwelling/downwelling dynamics. For instance, in the northern Chilean upwelling system, the vertical transport induced by coastal divergence (ZET) represented the 60% of the annual total upwelling, however, the Ekman pumping also displayed an important contribution of 40%, indicating that the largest differences between these mechanisms occur at spatial scale. That is, coastal ZET predominated in areas with low orography and headlands, whereas WEK was higher in regions with high orography and the presence of embayments (Bravo et al., 2016). Similar findings were observed in three typical upwelling systems of the South China Sea, showing temporal and spatial differences, highlighting regions where WEK had the same intensity as ZET, even doubling the amount of upwelling predicted solely considering the coastal divergence (Wang et al., 2013). Those physical mechanisms and their coupling are important to be explored because the physical transfer of water properties and nutrients from the deep to the surface ocean is strongly related to the vertical supply, which also involves important effects on primary productivity, changes in the water column mixing/stratification and associated water masses, as well as, in the local chemical processes (e.g., Williams and Follows, 2003; Pasquero et al., 2005). According to this, we assess both mechanisms (ZET and WEK) in this study, in order to relate its contributions to the maintenance of biological production (in terms of Chl-a) and key changes of the carbonate system.

To identify potential processes (physical mixing and biological) and explain the variability observed in carbonate system parameters (e.g., pCO_2w, DIC, TA, and O₂), we employed a three-endmember mixing model following established methodologies in oceanographic studies that developed the analysis of sea water from mixing triangles without assumption of isopycnic mixing (Tomczak, 1981; Paulmier et al., 2011; Kahl et al., 2018). This approach identifies water endmembers, in order to describe the thermohaline variability due to mixing, using conservative tracers (salinity and potential temperature) to construct a mixing triangle. The contribution of the water masses considered (M_k,i) to a given sample ‘i’ can be calculated by solving the followings determined system of three linear equations:

Where k is the water mass endmembers (1,2 and 3) and ‘i’ is the sample number (from 1 to 60). S_k and θ_k are the thermohaline characteristics of the k-water mass endmember (Tomczak, 1981). Finally, it computes expected conservative values due to physical mixing processes (DIC_mix, TA_mix, and O_2mix) for the non-conservative variables measured in situ (DIC, TA, and O₂). DIC_mix is defined as the fraction of dissolved inorganic carbon whose variability (R² from the regression DIC_mix vs. DIC) can predominantly be attributed to physical mixing. The complementary variability (1 – R²) is then attributed to biogeochemical processes, isolated through residual analysis (e.g., ΔDIC = observed DIC – DIC by mixing) and quadrant plots of ΔDIC vs. ΔO₂, where quadrants indicate net production (positive ΔDIC, positive ΔO₂), photosynthesis (negative ΔDIC, positive ΔO₂) ammonification (negative ΔDIC, negative ΔO₂) or respiration/remineralization (positive ΔDIC, negative ΔO₂) (Paulmier et al., 2011; Kahl et al., 2018).

Initially, we constructed T-S diagrams from the full spatiotemporal dataset (13 months, 7 stations, 0–80 m depth) to define three endmembers representing distinct sources of variability and mixing in the area. The river endmember (1; river-influenced) exhibited low salinity (26.81), TA (1750 µmol kg^-1), and DIC (1763 µmol kg^-1), with high temperature (27.94 °C) and O₂ (174 µmol kg^-1). The transitional endmember (2) showed intermediate values: salinity (33.06), TA (2138 µmol kg^-1), DIC (1867 µmol kg^-1), temperature (27.28 °C), and O₂ (123 µmol kg^-1). Finally, the oceanic endmember (3; oceanic influenced) featured high salinity (35.38), TA (1999 µmol kg^-1), and DIC (1989 µmol kg^-1), with low temperature (14.26 °C) and O₂ (40 µmol kg^-1; Preciado et al., in Progress). All three water endmembers incorporate signals from multiple upstream sources, reflecting the complex dynamics of this estuarine system. For instance, the oceanic (3) endmember exhibited influence from at least four diverse origins, including Antarctic Intermediate Water signals (Kawabe and Fujio, 2010); the intermediate endmember (2) blends subsurface and riverine inputs; and the river endmember (1) primarily mixes freshwater with marine waters, potentially including remote eddy-transported surface waters and local coastal currents. We could not directly measure riverine carbonate parameters due to differences in instrumentation and chemistry, river itself even showed tidal intrusion of seawater from below during high tides, further complicating pure freshwater end-member isolation. Instead, riverine influence is inferred through its dilution effects on salinity and biogeochemistry, following a similar approach to Huang et al. (2022) and Sun et al. (2023).

Once the mixing fractions are applied to compute expected DIC_mix values under conservative mixing assumptions. Deviations from these expected values then highlighted non-conservative processes affecting DIC, which is a key objective of our biogeochemical analysis. DIC, TA, and O₂ are used as complementary information for the characterization of the endmembers, enabling us to perform the mixing model analysis and determine the sources of variability (physical-mixing or biogeochemical processes) for the non-conservative variables DIC, TA, and O₂. This approach aligns with standard practices in estuarine and oceanic studies, where salinity and TA are often normalized or used to account for dilution and mixing effects before assessing non-conservative variables like DIC (Cai & Wang, 1998; Oliveira et al., 2017; Courtney et al., 2021). For instance, while TA behaves conservatively in many open-ocean settings and correlates linearly with salinity (Millero et al., 1998; Jiang et al., 2014), the T-S framework was chosen due to its effectiveness in capturing thermohaline distinctions in Pacific Colombian dynamic estuarine system, influenced by high rainfall and runoff. In addition, to validate the implication of the biological removal in pCO_2w concentration, we used a complementary approach, such as ΔDIC together with cumulative monthly in situ Chl-a data to explain CO₂ consumption. These results (biological/temperature implications) were also validated by an analysis to calculate if physical or biological factors are implicated in the pCO_2w concentrations found around Gorgona island, following the methodology described by Takahashi et al. (2002), for more detail on the equations used see Supplementary Table 2.

During the moderate 2021–2022 La Niña event, surface waters around Gorgona displayed a near-neutral net CO₂ flux (0.0105 ± 0.0001 mol C m^-2) over the nine sampled months, with a bimodal pattern driven by upwelling and post-upwelling wet season (Figure 4). Throughout the upwelling season, cooler, saltier waters reached the surface (shallow thermocline), coinciding with the observed low pCO_2w values (below atmospheric levels, mean pCO_2A = 461 ± 92.8 μatm and negative ΔpCO₂, see Figures 4A–E), resulting in near-zero to weakly negative FCO₂ (February–March, Figure 4F). This is additionally in line with elevated DIC and TA (Figures 4G, H). In contrast, during the post-upwelling period, the fluxes were slightly positive (0.0093 ± 0.0001 mol C m^-2) with an average pCO_2w of 567 ± 97.5 μatm (Figure 4D). Similarly, monthly fluxes also differed significantly between periods (Welch test, F = 2.5 × 10^-4, p< 1 × 10^-300; n = 9), peaking in June (2.5 mmol C m^-2 month^-1) and reaching a minimum in March (–0.1 mmol C m^-2 month^-1). Therefore, markedly higher fluxes were observed during the post-upwelling compared to the upwelling season when including the seven stations sampled per month (t = 4.91, critical t = 2.00, p = 7.8 × 10^-6; n = 21 vs. 39; Figure 4F).

Comparably, the average pCO_2w trend (Figure 4D) during the upwelling season evidences surface pCO_2w hovering moderately above atmospheric values in January (near 513 µatm Figure 4D; atmospheric reference value ~418 µatm, horizontal black line), but dips towards/just below atmospheric level by March (minimum average 397 µatm). Correspondingly, a post-upwelling spike in pCO_2w increases values up to ~700 µatm in May, around 1.7 times higher than atmospheric values, which drives a large positive ΔpCO₂ and outgassing (Figure 4E). Thereafter, pCO_2w slightly decreases until July, but remains well above atmospheric values (~580 – 590 µatm).

The ΔpCO₂ values in eight of the nine months were above the atmospheric pCO_2A - (between 50 and 260 micro atmospheres). However, during the upwelling season, there is a reduction in pCO_2w, causing it to be less saturated than the pCO₂ of the surrounding atmosphere. Thus, March evidenced less pCO₂ in the water than in the air (negative ΔpCO₂, Figure 4E) which is reflected with near-zero/negative FCO₂ (Figure 4F). Transition to the post-upwelling (higher rain and streamflow season; April-July) was also evident through the progressive increment in CO₂ peaking during May (ΔpCO₂ ≳ +250 µatm; more pCO₂ in the ocean than in the surrounding air) and then stabilizing, as reflected in the average CO₂ flux (Figure 4F).

The coolest temperatures were recorded during upwelling (~26–26.5 °C). Rapid warming of the surface water was observed in April (~28 °C) lasting until July (>27 °C). The shown warming is concurrent with the sharp rise in pCO_2w (Figures 4A, D). Average salinity was< 31 during the nine months, meanwhile, the highest salinity was recorded during upwelling (31.2 in March). From April to July (post-upwelling months) a steady drop in salinity was measured, coinciding with the peak of precipitation/river discharge (Figure 2). Hence, markedly brackish water was observed (low-alkalinity via dilution), related to pCO_2w increase (Figures 4B, D, H).

The highest pH_T was measured in March (8.01), chemically consistent with reduced pCO_2w values and less riverine/precipitation influence (Figure 4I). Then, pH_T began to decline during April and May (recorded a minimum of 7.81 in May) and stabilized after during June and July. Meanwhile, DIC start increasing between November and January (~1,700–1750 µmol kg^-1) reaching its highest values in April (1,815 µmol kg^-1) and May (~1800 µmol kg^-1). After, there was a general decrease from May onward (~1,790–1,650 µmol kg^-1 by July; Figure 4G).

The dominant DIC species, bicarbonate (HCO₃^-), remained relatively stable from November to February (~1,550–1,600 µmol/kg, Figure 4J) then increased sharply from March–April (~1,720 µmol/kg) in line with the highest DIC/TA values. Thereafter, HCO₃^- decreased gradually during the post-upwelling season (May–July) mirroring the reduction in DIC/TA.

Carbonate ion (CO₃^2-) values showed a clearer seasonal signal. They peaked during the upwelling months (January–March: ~170 and 180 µmol/kg, Figure 4K) in agreement with higher pH_T and Ω_AR (Figures 4I, L). Afterwards, CO₃^2- concentrations decreased during the post-upwelling period, sharply decreasing in April to ~150 µmol/kg, then rebounding slightly during May to ~170 µmol/kg, and finally steadily decreasing to its minimum values from June to July (114–118 µmol/kg) coinciding with lower pH_T, elevated Revelle Factor, and reduced Ω_AR.

Monthly accumulated Chl-a (sum of seven stations) was relatively high during the sampling period (~16–19 mg m^-³, Figure 4M), with average Chl-a of 1.09 – 2.77 mg m^-³ per sampled site (see Supplementary Figure 2). A marked increase from February to March was observed (16 – 19.4 mg m^-³). Afterwards, Chl-a values remain high during May (17.2 mg m^-³) until July with ~18 mg m^-³ (Figure 4).

The Revelle Factor (RF, Figure 4N) and Ω-aragonite (Ω_AR) displayed an inverse relationship. The lowest RF values (9.9) occurred in March, coinciding with weakened or even collapsed stratification (Figure 4C), when Ω_AR reached its peak (2.8). In contrast, RF increased during the post-upwelling season, rising from April and reaching ≥12 by July, concurrent with the lowest Ω_AR values (1.9), indicating suboptimal conditions for calcification. As expected, Ω_AR was inversely proportional to pH_T (Figures 4I, L), while RF varied inversely with SST.

Seasonal variability in dissolved inorganic carbon (DIC), its mixing component (DIC_mix), and related drivers highlighted the coupled influence of physical circulation and biological activity (Figure 5). During the upwelling season, lower ΔDIC values coincided with reduced pCO_2w and CO₂ fluxes, whereas in the post-upwelling season both DIC and pCO_2w increased together with intensified Ekman transport and pumping (Figures 5B–D). Regional chlorophyll-a distributions also followed these seasonal shifts, with peaks during February–March and again in May–July (Figure 5C).

Complementary evidence from a temperature–salinity (T–S) diagram (Figure 6A) underscored the role of mixing in shaping carbonate chemistry. Mixing calculations indicated that up to ~79% of the observed DIC could be explained by conservative mixing of three endmembers, while the residual ~19% reflected non-conservative processes (Figure 6B). The correlation between ΔDIC and ΔO₂ (Figure 6C) suggested that aerobic remineralization and nitrification were the dominant biological drivers of this residual variability. Physical processes dominated DIC variability, as indicated by the strong correlation between DICmix and DIC (R² = 0.79, Figure 6B). In contrast, respiration explained a smaller but measurable fraction of the variance (ΔDIC–ΔO₂, R² = 0.19; Figure 6C). Together, these results highlight that both physical and biological controls on DIC ultimately influenced its relationship with pCO_2w, consistent with the significant negative correlation observed between the two variables (Q = –0.71, p< 0.001, n = 60; Figures 4D, G).

The physical mixing of the water mass contributes significantly to explaining the DIC values ​​in the system (Figure 6A) as validated by the relationship between DIC_mix and DIC (Figure 6B; ~78% of variability) and the seasonal variability of WEK and ZET (Figures 5C, D). Therefore, DIC increased from November to May (Figure 5A) and then decreased from June to July. During La Niña 2021-2022, estuarine dynamics were strongly modulated by horizontal and vertical exchanges between marine and freshwater sources. High-frequency tides (up to 4 per day, with amplitudes of 5 m; Cauca tide table 2021–2022) mixed riverine and oceanic waters, while mesoscale eddies and Ekman-Zonal transport introduced deep offshore waters into the system (Corredor-Acosta et al., 2020; Figures 5C, D). Wind speed and direction generate local coastal currents, with north-south surface flow above the thermocline and south-north flow in subsurface waters (in situ observations); and upwelling and subsidence cycles (downwelling; Figures 5C, D) producing three water masses with different densities along the vertical profile (up to 80 m; see Figure 6A). Peaks in DIC between January and March coincided with intensified WEK and ZET which enhanced vertical injection of CO₂ and TA rich subsurface waters, when freshwater dilution and stratification was weaker (lower relative rain and precipitation; Figures 2B, D). Even so, elevation of DIC values could also be due to a gradual increase of bicarbonate and carbonate in the system (Figures 4J, K), when acidic sub-superficial waters (Ω_AR values sometimes below 1.5, Preciado A., unpublished data) could dissolve the reef structure, bringing the increased DIC signal to the surface. From April to July, the relative decrease in DIC reflected dilution by enhanced river discharge and weaker zonal inputs. The residual variability (~19%) was linked to respiration and remineralization, as evidenced by the correlation between ΔDIC and ΔO₂ where DIC increased but O₂ decreased (Figure 6C; O₂ levels reached 60 µmol/kg) and pCO_2w values were elevated. However, this respiration hypothesis requires further studies in the study site, to understand its contribution in DIC variability.

The seasonal variability of FCO₂ in Gorgona during the study period resulted from the system shifting from post-upwelling, dominated by southerly trade winds, higher precipitation and river dominated system to upwelling, dominated by northern winds, less precipitation, and lower river discharge, and then back to post-upwelling season. During post-upwelling, a slight positive increase in FCO₂ is related to the extensive plume of water from coastal rivers (Figure 2) originating from a 44% increase in river flow during the study period (e.g., the Guapi River with 103 m³ s^-1; Figure 2D), this plume extended for more than 60 km from the coast (as observed by floating logs on the field), forming a large estuary with an average salinity of 29.51 units (SD ± 1.56). Thus, following the definitions of Dai et al. (2013) and McKee et al. (2004), our study area can be classified as a system dominated by riverine inputs. During the post-upwelling season, pCO_2W values ​​were relatively high (Figure 4D), ranging from 378 to 839 μatm (average 567 μatm), peaking in May, in line with high SST (Figure 4A), and ΔpCO₂ ranging from 47 μatm to 418 μatm. Our results agree with Ricaurte-Villota et al. (2025), where rivers via remineralization could explain the high pCO_2w values ​during moderate La Niña years (see Figure 6C). But also, the increased freshwater discharge could also diminish CO₂ fluxes, simultaneously, under La Niña (Reimer et al., 2013; Ricaurte-Villota et al., 2025), via nutrient excess and phytoplankton bloom. During upwelling season, the average water temperature decreased to 26.53 °C, but the average salinity remained below 30, confirming the rivers influence on surface waters. Likewise, the pCO_2W was relatively lower, ranging from 316 to 626 μatm (average 461 ± 92 μatm), with minimum values ​​in March, coinciding with the lowest surface temperature (26.11 °C) and highest salinity (28.39). During this period, FCO₂ reached negative or near neutral values (between 0.3 and -0.1 mmol m^-2 d^-1), with similarly low ΔpCO₂ values, ranging from –101 μatm to 93 μatm. Despite theory indicating that upwelling during a moderate La Niña should increase the CO₂ flux to the atmosphere (Kim et al., 2017), this was not observed in our study site. Similar dynamics i.e., low or negative net fluxes despite upwelling have been reported in tropical estuarine systems where upwelled DIC is rapidly consumed by phytoplankton blooms which counterbalance the physical delivery of dissolved carbon (e.g., Vargas et al., 2007; Bouillon, 2011; Giraldo et al., 2011). Monthly patterns and flux directions were preserved when calculated using the propagated pCO_2w error, including the neutral/slight sink observed in March, with only December and June showing greater variability. Net flux differences over the nine-month period remained negligible (original net flux = 0.0105 mol C m^-2, positive error net flux = 0.0133 mol C m^-2, negative error net flux = 0.0077 mol C m^-2) indicating a net total difference of ±0.0028 mol C m^-2, confirming the robustness of our results.

The temporal variability of pCO_2W was shaped by the combined effects of salinity, temperature, DIC, and ΔO₂ variables were correlated with pCO_2W values. Salinity showed the strongest and most significant correlation with pCO_2w (Q = –0.96; p = 1.77E^–36, n = 60) reflecting the strong influence of river discharge and freshwater dilution, particularly during post-upwelling months, which differs from many oceanic regions where temperature is typically the main driver of FCO₂ variability (Séférian et al., 2013). Ekman-driven processes contributed to FCO₂ variability by vertical mixing processes (upwelling, downwelling; Figure 5C) that reduce salinity up to 40-60m, modifying the depth of the mixed layer over time. Sea surface temperatures recorded at the study site were relatively high (26.1–28.1 °C), as foreseen at equatorial latitudes, due to maximum solar insolation (Saraswat, 2011). On the other hand, sea surface temperature also correlated positively with pCO_2w (Spearman Q = 0.42313, p = 0.0007, n = 61; Figure 4A), in line with previous findings by Zhai et al. (2005), who reported higher pCO_2w values ​​during warmer seasons. During the post-upwelling, the increase in SST from April to June aligned with increases in pCO_2W, ΔpCO₂, and FCO₂. As a result, the observed trend is consistent with the thermodynamic relationships between temperature and CO₂ solubility (Johnson et al., 2010), even though sea surface temperature alone did not fully explain the FCO₂ variability. Furthermore, we assessed the relative influence of salinity and temperature to pCO_2W by performing a sensitivity test which fixed one of the variables (salinity or temperature) while allowing the other to vary in the CO2SYS, recalculating pCO_2w. This approach allowed us to isolate the individual effect of each parameter (see Supplementary Figure 1). Results showed that salinity explained a larger proportion of the variance – particularly in June and July when riverine CO₂ inputs were strongest, supporting the findings of Ricaurte-Villota et al. (2025) who also observed that river influence on salinity in the Colombian Caribbean exerts a greater control on pCO_2w than temperature. The lower salinity observed in June (post-upwelling) was associated with higher pCO_2w values, consistent with dilution by freshwater input and potential organic and inorganic matter increase from riverine mangrove sources (Palacios Peñaranda et al., 2019). The observed negative relationship between DIC and pCO_2W reflects the seasonal interplay of physical and biological processes at Gorgona. During upwelling, physical mixing and water mass advection introduced CO₂-rich subsurface waters, but biological drawdown suppressed pCO_2w provided by the sea and rivers. On the contrary, during the post-upwelling phase, the DIC increases and then decreases, due to the antagonistic relationship between the mixture (ZET, WEK; Figure 5) that, like respiration processes (oxygen reduction; Figure 6), provide pCO_2w (and DIC), and the photosynthetic activity, which consumes it or the excess of fresh water that dilutes the DIC (Figure 4).

Overall, the patterns in Figure 6 suggest that physical mixing and dilution were the primary controls of DIC variability, consistent with estuarine systems where hydrodynamic forcing often dominates carbon system fluctuations (Clark et al., 2022). Nevertheless, biological processes such as respiration and photosynthesis clearly modulated the signal, producing seasonal configurations of DIC–pCO_2W values that reflect the co-occurrence of both physical and biogeochemical drivers. Similar interactions have been reported in estuarine environments by Ahad et al. (2008) and Quiñones-Rivera et al. (2022), where strong physical control is complemented by meaningful, localized biological contributions. As anticipated, pH’s relationship to pCO_2w was inverse (Peng et al., 2013; Hans-Rolf and Fritz, 2023; Ramaekers et al., 2023), and direct with Ω_AR (See Figure 4I), which is similar to the relationship described by Feely et al. (2009), in our case also in agreement with the highest salinity values. In addition, high pH_T values during March were coupled with consumption of CO₂ by phytoplankton (lowered pCO_2w), and higher buffer capacity in the system (high TA; Figure 4), which is similar to cases reported by Cai et al. (2011) and Macdonald et al. (2009) in coastal and estuarine systems (low SSS), despite relatively high DIC values. This condition favors coral reef or marine organisms’ calcification, although the Ω_AR remained below 3, suggesting thermodynamically possible, yet suboptimal conditions (Ries et al., 2009). On the other hand, pH_T decreased during post-upwelling, to its minimum values (May; 7.81), as did TA (~1800 µmol kg^-1, Refer to Figure 4H), when the pCO_2w was at its maximum, thereby lowering the Ω_AR below 2.

The apparent suppression of CO₂ fluxes under upwelling conditions contradicted expectations from classical upwelling theory (e.g., Kim et al., 2017), which generally predicts enhanced CO₂ outgassing as deeper, carbon-rich waters reach the surface. Instead, rapid phytoplankton uptake offset the physical inputs (rivers, subsurface water, advection, atmosphere, and other sources). Evidence for this mechanism comes from strongly negative ΔDIC values in February and March (Figure 6D), indicating intense CO₂ consumption, and from the inverse relationship between surface pCO_2w and cumulative Chl-a concentrations (Figures 4D, M; Figure 5A). This agrees with the findings of Corredor-Acosta et al. (2020) who reported peak productivity around Gorgona during March, based on multi-annual satellite observations. Both their long-term satellite information and our average and accumulated Chl-a data (Figure 4M) confirm that the coastal waters near Gorgona sustain relatively high phytoplankton biomass throughout the year, with average concentrations exceeding 2.45 mg m^-3 in both seasons, and accumulated values surpassing 16 mg m^-3 (Figure 4M; Supplementary Figure 2).

According to Corredor-Acosta et al. (2020) the sustained Chl-a signal on the Colombian Pacific coast is modulated by rising waters along the Panama Bight during the upwelling season and river discharges (see also Rodríguez-Rubio et al., 2003; Devis-Morales et al., 2008; Giraldo et al., 2008; 2008; 2011). Therefore, the productivity observed in Gorgona is attributed to the constant input of nutrients (Giraldo et al., 2011), coming from both, subsurface waters and the nearby Patía–Sanquianga deltaic complex (~50 km), which contributes total suspended solids (23.2 mg L^-1), nitrates (33.55 µg L^-1), and phosphates (<2 µg L^-1; Giraldo et al., 2011). These inputs, according to INVEMAR (2019–2022) result in high Chl-a concentrations near the coast (~0.9 mg m^-3) compared to areas farther from the island (~0.3 mg m^-3). Additionally, the coastal phytoplankton around Gorgona is dominated by large diatoms and dinoflagellates which support a zooplankton biomass up to seven times greater (89 g/100 m³) than the regional average for the Colombian Pacific (12 g/100 m³; Murcia Riaño and Giraldo López, 2007; Giraldo et al., 2011). Therefore, zooplankton also facilitates carbon export through marine snow, contributing to remineralization in the water column and up approximately 90–100 m depth (Alldredge, 1984; Vargas et al., 2007) which further aids in pCO_2w modulation at the surface (WEK, Figure 5C). The persistence of such high productivity effectively modulates surface pCO_2w by enhancing phytoplanktonic CO₂ uptake, thereby reducing the magnitude of air-sea fluxes. Without this sustained biological activity, surface waters around Gorgona would likely accumulate greater CO₂ excesses and display substantially higher efflux to the atmosphere during La Niña years.

The near-neutral to slightly positive CO₂ fluxes observed around Gorgona can be understood as the outcome of several interacting climatic–oceanographic processes. As shown, the system displays marked seasonal variability, yet overall fluxes remained low to neutral compared with other tropical Pacific or estuarine sites. Three main drivers appear to regulate this behavior: wind-driven exchange, vertical stratification, and the buffering capacity of surface waters.

Wind speeds during the study period were consistently low (mean<5 m s^-1; range: 2.1–6.4 m s^-1; Figure 5), which limited the efficiency of gas transfer. Since air–sea CO₂ fluxes are jointly governed by wind speed and the air–sea CO₂ partial pressure gradient (ΔpCO₂), the absence of high-wind events (>10 m s^-1; Prytherch et al., 2010) meant that surface–atmosphere exchange remained weak throughout most of the year. Even when ΔpCO₂ values were favorable for outgassing, limited wind-driven turbulence constrained fluxes to the atmosphere.

Vertical stratification further reinforced this suppression. The water column at Gorgona remained strongly stratified for most of the year, with the thermocline–halocline–pycnocline complex typically lying between 11 and 40 m depth (Figure 4C). This effectively isolated subsurface CO₂-enriched waters from the surface, preventing their sustained transfer into the mixed layer. Only in March, during the shoaling of the thermocline (to ~11 m), did subsurface waters briefly reach the surface, temporarily breaking stratification. Yet even under these conditions, fluxes remained near-neutral rather than shifting strongly positive.

This apparent contradiction can be explained by the buffering capacity of surface waters. Revelle Factor (RF) values were lowest during upwelling (January–March; Figure 4N), indicating a high resistance of surface waters to changes in pCO_2W despite inputs of dissolved inorganic carbon (DIC) from below. As a result, increases in subsurface DIC did not translate into proportional rises in surface pCO_2W, thus preventing strong efflux. By contrast, RF increased during the post-upwelling period (May–July; Figure 4N), when surface waters became more sensitive to DIC changes. Nevertheless, stratification was then strongest, reducing the vertical supply of CO₂-rich subsurface waters and limiting fluxes.

Taken together, these mechanisms explain why Gorgona exhibited suppressed CO₂ efflux during both upwelling and post-upwelling conditions. Wind limitation constrained the physical transfer of CO₂, stratification acted as a barrier separating sources from the surface, and buffering capacity modulated the extent to which subsurface DIC inputs could alter surface pCO_2W. This combination is consistent with earlier findings in the Colombian Caribbean (Reimer et al., 2013; Ricaurte-Villota et al., 2025), where La Niña–driven stratification similarly reduced CO₂ outgassing. At Gorgona, therefore, the interplay of seasonal stratification and surface carbonate chemistry set the stage for fluxes that remained near-neutral, even under physical conditions that would normally favor enhanced CO₂ release.

To enable comparison of our nine-month flux record with annual estimates reported in literature, we estimated the three missing post-upwelling months (August to October 2020) in order to derive a potential annual flux. This approximation was based on our in situ measurements from six post-upwelling months that correspond to the rainy (wetter) season, a period which also encompasses the missing months. From these six months, we calculated a mean monthly flux of 0.00024 ± 0.0001 mol C m². Using this value, we estimated the flux for the three missing months (0.0007 mol C m²; net positive flux) and added it to the nine-month total (0.01049 ± 0.000148 mol C m²) yielding a potential annual flux of 0.0112 ± 0.00015 mol C m² yr^-1. Interpolation of missing FCO₂ data has been applied previously (e.g., Kahl, 2018). Nevertheless, we emphasize that the estimated three-month flux is a hypothetical approximation intended only for comparison purposes, and should be validated with additional in situ measurements.

This approach revealed a neutral to weakly positive net air-sea CO₂ flux at Gorgona Island (FCO₂ = 0.0112 ± 0.00015 mol C m^-2 yr^-1; ~0.031 mmol C m^-2 d^-1) being an order of magnitude lower than those reported by Dai et al. (2022) for non-ENSO years (0.8–1.0 mol C m^-2 yr^-1), and close to the lower limit of the global range reported by Wong et al. (2022) for surface ocean carbon fluxes (0–1.5 mol C m^-2 yr^-1); highlighting potential regional divergence from global modeling outputs (Table 1).

The Gorgona fluxes were also markedly lower than those observed in other tropical and coastal systems. For instance, at similar equatorial latitudes (~3°N), Jersild et al. (2017) and Landschützer et al. (2022) reported fluxes exceeding 0.3 ± 0.1 mol C m^-2 yr^-1 during moderate La Niña conditions. Pacific Ocean continental shelves exhibit values between 0.2 and 0.35 mol C m^-2 yr^-1 (Jersild et al., 2017; Vaittinada Ayar et al., 2022), while the North Pacific shelf, even under non-ENSO conditions, registers around 0.043 mol C m^-2 yr^-1 (Reimer et al., 2013). These regional differences suggest that local environmental controls strongly limit the air-sea CO₂ exchange in the Eastern Tropical Pacific, particularly near Gorgona. Consistent with this interpretation, Ricaurte-Villota et al. (2025) documented similarly low to near-neutral fluxes in the Colombian Caribbean under La Niña conditions, attributing them to high rainfall, enhanced stratification, and strong biological uptake. These drivers mirror those identified in our study, reinforcing the idea that La Niña exerts a suppressive effect on CO₂ release in this region. Comparable mechanisms have also been described in Arctic continental shelves, where freshwater inputs from rivers and ice melts (together with seasonal ice coverage) enhanced stratification, restricting gas exchange and contributing to low or near-neutral fluxes (Else et al., 2013; Miller et al., 2019; Mu et al., 2020).

Furthermore, comparison with other tropical systems reinforces our findings. Similar low or near-neutral fluxes have been reported in Costa Rica’s Gulf of Nicoya (Pfeil et al., 2013) and other areas of the Eastern Tropical Pacific (Laurelle et al., 2013).

Estuarine systems by contrast, tend to show higher fluxes and yearly variability (Chen et al., 2013). For example, the Matla estuary in India exhibits high annual fluxes of 2.3 mol C m^-2 yr^-1 (Akhand et al., 2016) over 10 times greater than the estimated CO₂ flux around Gorgona, and the Patos Lagoon in Brazil ranges from -13.9 to 19.7 mol C m^-2 yr^-1 attributed to changes in freshwater input and phytoplankton activity (Albuquerque et al., 2022). Meanwhile Asian tropical estuaries average around 8.1 mol C m^-2 yr^-1 (Chen et al., 2013) which are also substantially higher than our estimated yearly fluxes. However, comparable near-neutral fluxes were reported in the Arafura and Red Seas under similar La Niña conditions (Hydes et al., 2012).

Although global models (e.g., Dai et al., 2022) and some regional studies do not explicitly separate ENSO phases, this distinction is essential for accurately characterizing flux dynamics in equatorial regions such as Gorgona. Prior literature indicates that La Niña conditions can reduce FCO₂ by 0.2–0.4 mol C m^-2 yr^-1 in the Pacific region, compared to El Niño or neutral years, aligning with the low values observed in this study (Chavez et al., 1999; Landschützer et al., 2016; McKinley et al., 2020; Park et al., 2010).

Our findings challenge the assumption that upwelling-dominated systems and estuarine environment are high CO₂ emitters, especially under La Niña conditions. Since in our tropical coastal area FCO₂ is sensitive to atmospheric conditions (wind velocity, precipitation, solar radiation), freshwater discharge, physical dynamic (vertical-horizontal water movement), vertical stratification, upwelling events, carbonate chemistry and biological uptake. Future research should incorporate higher-resolution temporal data, uncertainty analyses in ΔpCO₂, and ENSO-phase-specific studies (and neutral/non-ENSO years) to improve carbon flux estimates in this and similar tropical estuarine systems.

Our study did not include measurements of δ¹³C-DIC, which limits our ability to trace the biogeochemical sources and transformations of dissolved inorganic carbon. This isotopic parameter is a powerful tool for distinguishing between marine, terrestrial, and anthropogenic carbon inputs, and for assessing the extent of biological activity and remineralization. Its absence constrains our interpretation of carbon cycling dynamics, particularly under varying ENSO conditions. We suggest that future sampling campaigns incorporate δ¹³C-DIC analysis to enhance source attribution and improve the mechanistic understanding of CO₂ flux variability in coastal upwelling systems. Additionally, we did not directly measure nutrients, which limits our ability to consider them in the CO2SYS calculations, specifically those which are proton acceptors (e.g. phosphates, silicates, ammonia), which can have effects on the computation of certain carbonate system parameters under high nutrient loads. We recommend that future sampling efforts incorporate nutrient measurements for improved flux and carbonate system calculations.

It is important to consider methodological differences when comparing with other studies. Although references such as Feely et al. (2004, 2006) and Hydes et al. (2012) offer context, disparities in sampling protocols, spatial resolution, and ecosystem type limit the scope of direct comparison. Finally, we recommend a long-term monitoring program taking water samples to estimate in the lab the couple DIC/TA, this will allow to compare the CO₂ flux variability between ENSO events and neutral conditions. An inexpensive discrete sampling based on our experience for Gorgona National Park could be to take water samples for 2–4 contrasting months.