Autonomous monitoring was carried out by an automated underway seawater and low atmospheric CO₂ molar fraction (xCO₂) measurement system installed inside the engine room of the RENATE P container ship (Curbelo-Hernández et al., 2021). The xCO₂ measurement system, developed by Craig Neill and commercialized by General Oceanics^TM, combines an air and seawater equilibrator with a non-dispersive infrared analyzer for gas detection (Pierrot et al., 2009).

The system was checked automatically on departure and arrival at each port, and periodically every 3 h during the ship’s journey. Four standard gasses (in the order of 0, 250, 400, and 540 ppm, with a ±0.02 ppm accuracy) provided by the National Ocean and Atmospheric Administration (NOAA) and traceable to the World Meteorological Organization (WMO) were used for system check and xCO₂ measurement corrections. Periodic calibrations were automatically performed every 12 h using standard gasses 1 and 4 to adjust the, respectively, zero and span of the infrared analyzer.

The system measures xCO₂ (ppm) in seawater with a frequency of 2–3 min and normally collect 130–150 data in the Strait of Gibraltar. The system also performs three xCO₂ measurements in low atmosphere after each calibration. In both cases, a non-dispersive infrared analyzer built by LICOR^® (initially the 6262 model and after October 2019, a 7000 model) were used for the detection of xCO₂.

The sea surface temperature (SST, in °C) was measured using an SBE38 thermometer with an estimated error of 0.01°C placed at the main seawater intake. In addition, due to the high sensitivity of xCO₂ to temperature changes, an SBE45 thermosalinograph and a Hart Scientific HT1523 Handheld Thermometer, with an estimated error of ±0.01°C each, were used to monitor the temperature near the xCO₂ system and inside the equilibrator, respectively.

Sea surface salinity (SSS) was measured with a manually calibrated SBE45 thermosalinograph, with an estimated error of ±0.005. Atmospheric pressure was measured by the transducer on the deck (these pressure records being different in the order of millibars with the pressure near the xCO₂ system, inside the ship).

Due to technical problems with the measurement equipment aboard the SOOP CanOA-VOS, no measurements were obtained during March 2019 (problems with water intake), during the second half of August 2019 and March–April 2020 (problems with LICOR flow that could not be resolved due to COVID-19 limitations) and on a limited number of return journeys from the port of Barcelona.

This study followed the data collection methodology, quality control, and calculation procedures as published in the updated version of the DOE method manual for ocean CO₂ analysis (Dickson et al., 2007). The fugacity of CO₂ was calculated in seawater (fCO_2,sw) and in the lower atmosphere (fCO_2,atm) from measured and corrected xCO₂ values following Pierrot et al. (2009).

The thermal and non-thermal contribution on fCO_2,sw and variations was studied using the equations presented by Takahashi et al. (2002), with the temperature effects on fCO_2,sw for isochemical seawater of 0.0423°C^–1 determined experimentally by Takahashi et al. (1993). The non-thermal effect (fCO_{2,non–thermal}) was obtained from the seasonal amplitude of the normalized fCO_2,sw to the average temperature (Eqs 1, 2). The thermal effect (fCO_2,thermal) was calculated by applying the observed temperature effect to the average value of fCO_2,sw (Eq. 3) and determining its seasonal amplitude (Eq. 4). The relative importance of thermal and non-thermal effects was expressed by the T/B ratio (ΔfCO_2,thermal/ΔfCO_{2,non–thermal}), with values greater than 1 indicating that the temperature effect dominated over non-physics effects.

The CO₂ fluxes (FCO₂) were determined using Eq. 5 (Broecker and Peng, 1982) with a conversion factor of 0.24 mmol m^–2 d^–1. The solubility (S, mol L^–1 atm^–1) and the difference between seawater and low atmosphere fCO₂ (ΔfCO₂ = fCO_2,sw − fCO_2,atm) were considered. Negative fluxes indicate that the ocean acts as an atmospheric CO₂ sink, while the positive ones indicate that it behaves as a source.

The updated Wanninkhof (2014) parameterization was used in this study due to it considers the most recently advances in the quantification of the input parameters and improves the wind speed products, with k being the gas transfer rate expressed in Eq. 6:

where w is the wind speed (m s^–1) and Sc is Schmidt number (cinematic viscosity of seawater, divided by the gas diffusion coefficient). Daily averages of satellite ocean surface wind speed derived from the Advanced Scatterometer (ASCAT) and obtained from the Satellite Research and Exploitation Center (CERSAT) at IFREMER (Plouzane, France) were considered together with daily measured wind speed averages provided by the Agencia Estatal de Meteorología (AEMET) from the meteorological station located at Tarifa (Cádiz) and corrected to 10 m height (Allen et al., 1998). These wind speed data were interpolated based on the latitude, longitude, and time of the underway xCO₂ measurements. The average ocean surface wind speed for the entire study period calculated from satellite and measured data was 7.77 ± 2.66 and 8.18 ± 2.61 m s^–1, respectively. Due to these wind speed values were consistent, the FCO₂ was calculated in this study by considering the wind speed data measured at the Tarifa meteorological station as was done previously by de la Paz et al. (2011) to assess the air-sea CO₂ exchange in the Strait of Gibraltar.

Autonomous monitoring was combined with 25 discrete surface seawater samples taken manually from the intake seawater line for total alkalinity (A_T, μmol kg^–1) and total dissolved inorganic carbon (C_T, μmol kg^–1) in February 2020 and March 2021 with in situ measurements of SST and SSS. Samples were taken with a temporal frequency of 1–2 h in borosilicate glass bottles, overfilled, preserved with 100 μl of saturated HgCl₂, kept in darkness and analyzed just after arriving at port, in a period less than 2 weeks. The A_T was determined by potentiometrically titration in a VINDTA 3C system following Mintrop et al. (2000) and using CRMs (provided by A. Dickson at Scripps Institution of Oceanography) to test the performance of the titration system and correct the determined values. The CRMs titration gave values with and accuracy of ±1.5 μmol kg^–1.

The determined A_T values were consistent with those estimated according to the global relationship of A_T with SSS and SST in global surface waters described by Lee et al. (2006) for the North Atlantic (presented an average difference of 3.08 ± 6.80 μmol kg^–1). Due to the similarity of both values, A_T was calculated for the time, longitude, and latitude of the CanOA-VOS data using the A_T-SSS relationship obtained from the discrete samples (Eq. 7; r² = 0.87 and Pearson correlation coefficient given hereinafter as ρ is equal to 0.93) (temperature was not found to improve the fitting of the experimental data). It was assumed that the change in A_T with SSS was constant throughout the year in tropical and subtropical latitudes (Lee et al., 2006).

The C_T and pH in total scale (pH_T) were computed with the Excel program CO₂sys, using the carbonic acid dissociation constants of Lueker et al. (2000), the HSO₄ dissociation constant of Dickson (1990) and the value of [B]_I determined by Lee et al. (2010). The longitudinal change in the seasonality of the variables measured in situ (SST and SSS), determined (fCO_2,sw and FCO₂) and computed (C_T and pH_T) was studied in the entire region along the northern and southern sections. The C_T was normalized (NC_T) to a constant SSS of 36.4 (the annual average SSS) to remove the effect of evaporation and precipitation. Similarly, the pH_T values were normalized (pH_T,18) to a constant SST of 18.0°C (annual average SST) to remove the temperature dependence on the variation of pH_T (Santana-Casiano et al., 2001; González-Dávila et al., 2003).

The temporal variation of C_T was studied using the equation provided by Sarmiento and Gruber (2006) (Eq. 8). This equation describes the change in C_T over time produced by air-sea CO₂ exchange (EX), horizontal advection (TRSP), vertical mixing (MIX), and net community production (NCP).

The increase and depletion of NC_T in surface waters was studied during a whole year (February–January) using monthly average NC_T values between February 2019 and February 2021 and considering the variation of the mixed layer depth (MLD). The temporal change in NC_T was calculated as the difference between two consecutive months and was expressed in mmol m^–3, considering the seawater density. Negative values indicate that NC_T decreases in surface waters. The relative contribution of CO₂ exchange in the variation of NC_T (EX term) is given by the relationship between FCO₂ and MLD. The monthly MLD values were obtained with a spatial resolution of 0.028° × 0.028° from the operational Iberian–Biscay–Irish (IBI) Ocean Analysis and Forecasting system based on a (eddy-resolving) NEMO model application and are available at E.U. Copernicus Marine Service Information.¹ The annual MLD cycle was calculated using the harmonic fitting equation and parameters shown in Supplementary Table 1.

The relative contribution of horizontal transport (TRSP term) was calculated by considering a surface transport from the Northeast Atlantic to the Mediterranean Sea. The spatial change in NC_T between westernmost and easternmost waters [ΔNC_T = NC_{T, (West)} − NC_{T, (East)}] was considered. Data available at the Surface Ocean CO₂ Atlas (SOCAT v2020²) referenced to 2019 have been used for the calculation of NC_T at a reference station to the west of the Strait of Gibraltar (6.75 ± 0.25°W) in a similar latitudinal range. The spatial change in salinity was removed by calculating the NC_T at all the selected locations in the Strait of Gibraltar at the constant salinity of the reference station. The horizontal transport contribution to NC_T fluctuations was obtained from the product of ΔNC_T and the change in salinity with time at the locations of interest in the Strait of Gibraltar. Positive values indicate that horizontal transport occurred in a west-to-east direction, while negative values indicated that the transport direction was reversed.

The relative contribution of vertical mixing (MIX term) was estimated following Bégovic and Copin-Montégut (2002) by considering the sum of changes in vertical diffusion across the bottom of the mixed layer and vertical entrainment into the surface layer (Eq. 9). This term was determined using the monthly values of total dissolved inorganic carbon in the MLD obtained from IBI MFC high-resolution biogeochemical forecasts and generated from the PISCES biogeochemical model and the NEMO ocean circulation model. Additional monthly values of temperature and salinity obtained from the operational IBI Ocean Analysis and Forecasting system based on a NEMO model were used to estimate the seawater density. These forecast data have a spatial resolution of 0.028° × 0.028° and are available in E.U. Copernicus Marine Service Information.¹

The difference between NC_T concentration below (NC_T,f) and at the base (NC_T) of the mixed layer was considered to calculate the term of the vertical entrainment. The function θ is equal to 0 when the mixed layer decreases (∂⁡MLD/∂⁡t) ≤ 0 due to the fact that only the deepening of the mixed layer induces mixing with the underlying waters (Peng et al., 1987; Fasham et al., 1990; Gruber et al., 1998). The vertical gradient of NC_T at the base of the mixed layer (∂⁡NC_T/∂⁡z) was considered for the calculation of vertical diffusion. The vertical diffusion coefficient (K_z, m^–2 d^–1) at the base of the mixed layer was determined according to Denman and Gargett (1983) using Eq. 10:

where ρ is the seawater density at the base of the mixed layer, g is the acceleration of gravity and ε is the rate of turbulent energy dissipation. In this study we considered the vertical density gradient (∂⁡ρ/∂⁡z) in the pycnocline and a representative value of 2.0 × 10^–8 m² s^–3 for ε (Gruber et al., 1998).

Lastly, the change in NC_T due to net community production processes (NCP term), considering total changes in photosynthesis/respiration processes, was obtained directly from Eq. 8. Positive contributions of NCP indicated that the remineralization processes increase the NC_T in surface waters while negative contributions shows that biological production decrease it surface concentration.

The raw output data was initially filtered removing data affected by the automatic sampler such as samples measured at low water rates (<2.5 L min^–1). The measured xCO₂ values were corrected following Pierrot et al. (2009) by linearly interpolating the certified standard values (section 3.2) with the xCO₂ measurement times. The three fCO_2,atm values obtained after each calibration were averaged and harmonically interpolated with the times of each fCO_2,sw. The ΔfCO₂ was computed from the difference between the fCO_2,sw and the fCO_2,atm values and was used in the FCO₂ calculations (section 2.2). Additional daily average satellite data for Chlorophyll a obtained from the Operational Mercator Ocean biogeochemical global ocean analysis and forecast system (United States Copernicus Marine Service Information) were used to improve the understanding of the biological influence on the CO₂ system variability. Chlorophyll a data was interpolated with latitude, longitude, and time of the CanOA-VOS line data.

The average physical and biogeochemical variables (y) for each selected station were fitted to a harmonic Eq. 11 as a function of time (x) in order to study the seasonal variability and the thermal and non-thermal effect on fCO_2,sw (Takahashi et al., 2002; Lüger et al., 2004). The a–e coefficients are shown in Supplementary Table 2.

The results were compared to surface data collected by other VOS and research vessels in the oceanographic environment of the Strait of Gibraltar and are available at the SOCAT v2020 database. The SOCAT data between 1999 and 2019 were used to study the monthly frequency and the annual cycle of fCO_2,sw in 2019 by considering an interannual rate of increase of 1.8 μatm per elapsed year (Bates et al., 2014).

The longitudinal west-to-east distribution of the variables of interest was studied (Supplementary Figures 1, 2) and showed differences, especially during summer. The longitudinal variation was analyzed separately in the northern and southern sections: the SST, SSS, and CO₂ system variables and fluxes were annually and seasonally averaged every 0.05° longitude (Figure 2) by considering the following periods: January–March (winter), April–June (spring), July–September (summer), and October–December (autumn). The average values show that SST drives the fCO_2,sw and FCO₂ variations. The greatest latitudinal differences were found toward the east of GC. The average SST, fCO_2,sw, FCO₂, and pH_T values follow an opposite west-to-east trend between the northern and southern sections. The seasonality of SST and fCO_2,sw was maximum (5.29 ± 2.97°C and 45.51 ± 28.11 μatm, respectively) in the westernmost part (6.15–5.95°W) and minimum in the eastern area closest to the Mediterranean Sea (5.30–5.10°W) in both the northern (4.47 ± 2.68°C and 19.11 ± 24.86 μatm, respectively) and the southern (3.39 ± 2.11°C and 23.81 ± 14.96 μatm, respectively) sections. The seasonal change in pH_T (northern and southern sections, respectively) was 0.053 ± 0.028 and 0.033 ± 0.024 U in the westernmost part and 0.019 ± 0.025 and 0.025 ± 0.015 U in the easternmost part.

An annual average west-to-east decrease through the northern and southern routes was obtained for SST [1.70°C (r² = 0.64 and ρ = −0.80) and 1.48°C (r² = 0.83 and ρ = −0.91) per degree of longitude, respectively] and fCO_2,sw [14.21 μatm (r² = 0.65 and ρ = −0.81) and 2.71 μatm (r² = 0.22 and ρ = −0.47) per degree of longitude, respectively] (Figures 2A,E). The seasonality of pH_T increased toward the east in both sections (Figure 2G) and correlated with the decrease in fCO_2,sw (−0.001 U μatm^–1, r² > 0.94, and ρ < −0.97). These longitudinal changes were maximum during summer in the northern section: SST and fCO_2,sw decreased by 2.95°C (r² = 0.66 and ρ = −0.81) and 47.47 μatm (r² = 0.87 and ρ = −0.93) per degree of longitude, respectively, while the pH_T increased by 0.047 U (r² = 0.86 and ρ = 0.93) per degree of longitude.

The annual average FCO₂ (Figure 2E) shows that the CO₂ sink increased toward the east in both the northern (−1.89 mmol m^–2 d^–1 per degree of longitude, r² = 0.68, and ρ = 0.82) and the southern (−1.27 mmol m^–2 d^–1 per degree of longitude, r² = 0.57, and ρ = 0.75) sections. The west-to-east increasing sink gradient was most intense during spring [−2.23 (r² = 0.34 and ρ = 0.58) and −1.38 (r² = 0.64 and ρ = 0.80) mmol m^–2 d^–1 per degree of longitude, respectively] and summer [−6.71 (r² = 0.89 and ρ = 0.94) and −2.69 (r² = 0.60 and ρ = 0.77) mmol m^–2 d^–1 per degree of longitude, respectively]. During summer, both the northern and southern sections behaved as a CO₂ source (3.11 ± 1.76 and 1.24 ± 3.84 mmol m^–2 d^–1, respectively) in the westernmost part (6.15–5.5°W) and as a weak CO₂ sink (−0.87 ± 2.62 and −0.57 ± 2.48 mmol m^–2 d^–1) in the easternmost part (5.5–5.1°W). Maximum outgassing was found in both the northern and southern sections during summer around ES (4.49 ± 2.45 and 2.28 ± 3.91 mmol m^–2 d^–1, respectively) and TN (2.99 ± 0.74 and 2.20 ± 2.59 mmol m^–2 d^–1, respectively). In contrast, maximum ingassing was obtained along the northern section between ES and CS during autumn (−6.01 ± 0.05 mmol m^–2 d^–1) and toward the east of TN along the southern section during winter (−6.57 ± 4.11 mmol m^–2 d^–1).

The longitudinal change in SSS (Figure 2B) was significant in this transition area between Northeast Atlantic and Mediterranean Sea waters and influenced the spatial distribution of C_T. Big differences in SSS and C_T values were observed between the northern and southern sections. In the northern section, the annual average of SSS slightly decreased from the most Atlantic part to GC (0.12 U per longitude, r² = 0.66, and ρ = −0.81) and increased rapidly toward the east of GC (1.40 U per longitude, r² = 0.95, and ρ = 0.97). These longitudinal changes were also observed in the southern section, but both the decrease (0.08 U per degree of longitude, r² = 0.67, and ρ = −0.82) and the increase (0.50 U per degree of longitude, r² = 0.86, and ρ = 0.93) were lower.

The longitudinal distribution of C_T (Figure 2D) was correlated with that of SSS in the northern section (r² = 0.91 and ρ = 0.95), while this control was not observed in the southern section (r² = 0.13 and ρ = 0.36). The average annual NC_T (Figure 2F) increased similarly toward the east in both the northern (11.42 μmol kg^–1 per degree of longitude, r² = 0.64, and ρ = 0.80) and the southern sections (12.64 μmol kg^–1 per degree of longitude, r² = 0.80, and ρ = 0.89). The average NC_T values were maximum during the cold months (January–June) and minimum during the warm months (July–December) in the entire region, except in the northern section during autumn to the east of CS and during summer around TN. The decrease in fCO_2,sw between CS and GC in summer and the effects of remineralization and vertical mixing processes accounted for these changes. The absolute minimum values of NC_T in summer along each route should be highlighted: the minimum value eastward of GC in the northern section (2090.60 ± 4.99 μmol kg^–1) was followed by a reduction of the oceanic CO₂ sink and high values of fCO_2,sw, while that of the southern route in TN (2090.03 ± 2.95 μmol kg^–1) coincided with maximum values of SST, fCO_2,sw, and FCO₂ and may be related to an increase in biological production.

The seasonality of SST, SSS, and CO₂ system variables and fluxes were studied on both the northern and the southern routes (Supplementary Figures 3, 4) by averaging the values obtained on each trip and fitting them to a harmonic Eq. 11. The average values of fCO_2,sw and pH ranged between ∼360 and 430 μatm and between ∼8.03 and 8.09, respectively, in phase with a seasonal temperature change of ∼7.5°C (between ∼15.5 and 23.0°C). Several differences were observed between the northern and the southern sections: the fCO_2,sw and pH trends (Supplementary Figure 3) show a lower seasonality along the southern route despite the fact that the seasonal change in SST is maximum in this section. The fCO_2,sw values were compared with surface fCO_2,sw data available at the SOCAT v2020 database for the Strait of Gibraltar between 1999 and 2019 (Supplementary Figure 5). The CanOA-VOS lines provided a total of 6440 measurements during the 2 years of study (2019–2021), almost the same number of observations (6063 data) as in the SOCAT database for a period of 20 years (1999–2019). The CanOA-VOS line collected data at times of the year when there was no previous data available (February, May–July, and December). The annual cycle of fCO_2,sw provided from both the CanOA-VOS and the SOCAT databases is described by Eqs 12, 13 (r² = 0.67 and r² = 0.64, respectively) in which SST fluctuations were considered (SSS was not found to improve the fitting). Their standard errors of estimate was ±10.2 and ±13.8 μatm, respectively.

The seasonal change of fCO_2,sw with SST obtained through the northern and southern sections of the CanOA-VOS line is given by Eqs 14, 15 (r² = 0.78 and r² = 0.66, respectively) with a standard error of estimate of ±10.3 and ±8.9 μatm, respectively.

When the full data set (SOCAT and CanOA) was considered (more than 12,000 data), Eq. 16 (r² = 0.62), which present a standard error of estimate of ±12.8 μatm, provides an estimation of the interannual trend (2.33 ± 0.06 μatm year^–1), the seasonal variability and the influence of SST as the variable controlling the physical characteristics of the Strait of Gibraltar.

The FCO₂ seasonal trend (Supplementary Figure 3) ranged between −6 and 3 mmol m^–2 d^–1 in the northern section and between −8 and 3 mmol m^–2 d^–1 in the southern section, with a strong behavior as a sink during the cold months and as a weak source during the warm months. The maximum outgassing (2.99 ± 2.70 and 3.60 ± 1.76 mmol m^–2 d^–1 in the northern and the southern section, respectively) occurred in August 2020 and coincided with high values of SST and fCO_2,sw, while the maximum ingassing (−11.48 ± 1.36 mmol m^–2 d^–1) was observed along the southern section in February 2021 and was related with wind speed values higher than 11 m s^–1. The change in FCO₂ was linked to the change in ΔfCO₂, which is controlled by temperature fluctuations: the average values of FCO₂ increased by 1.22 mmol m^–2 d^–1 °C^–1 (r² = 0.71 and ρ = 0.84) in the northern section and 1.03 mmol m^–2 d^–1 °C^–1 (r² = 0.50 and ρ = 0.71) in the southern section. A smoothing spline with a parameter of 0.98 was applied to the daily distribution of SSS (r² > 0.27), SST (r² > 0.82), fCO_2,atm (r² > 0.94), fCO_2,sw (r² > 0.56), and wind speed (r² > 0.49). Daily values of FCO₂ were calculated from these computed variables and fitted to the observed FCO₂ values to improve the evaluation of its temporal distribution. These computed FCO₂ values differ in average by −0.25 mmol m^–2 d^–1 in the northern section and by 0.15 mmol m^–2 d^–1 in the southern section compared to those calculated from the observed data.

The west-to-east differences were assessed by selecting five stations equidistantly spaced every 0.2° of longitude along the channel (6.0–5.2°W) (Figure 3 and Supplementary Figures 6, 7). The westernmost station at 6.0 ± 0.05°W (S1) is located near ES and was used to study the seasonality in the most Atlantic part over the main gateway of the Mediterranean outflow. The station at 5.8 ± 0.05°W (S2) was selected near the shallower area of CS. The stations at 5.6 ± 0.05°W (S3) and 5.4 ± 0.05°W (S4) were in the longitudinal range around TN and GC, respectively, over the deepest area of the channel. Lastly, the easternmost station at 5.2 ± 0.05°W (S5) was used to study the most Mediterranean part. To consider the observed latitudinal variation, the northern and southern parts of each station were represented separately at all stations except S1 due to upward and backward routes passing through the same latitudinal interval in the westernmost part. A lower temporal resolution was obtained along the northern part of the channel and at S5 (Supplementary Figure 6) due to technical problems with the measurement equipment and subsequent scarcity of data during the first year of observation. Hence, the harmonic fitting is useful to understand the temporal change during periods of lack of data, but the seasonal trend shows a higher uncertainty in these positions.

The latitudinal differences are greater in the most Mediterranean part of the channel. The seasonality of the SST (Figure 3 and Supplementary Figure 6) decreased toward the east through the northern route (in the range between 16–23°C at S1 and 16–20°C at S4) and remained constant with longitude through the southern route (16–23°C). The fCO_2,sw values (Figure 3 and Supplementary Figure 6) were maximum at S1 throughout the year (ranging between 370 μatm during the cold months and 428 μatm during the warm months), and its seasonality remained constant along the southern section between 370 and 415 μatm in phase with temperature and influenced by non-thermal processes. The west-to-east pattern in fCO_2,sw seasonality along the northern section was difficult to characterize due to the lack of values at certain times of the year. However, a slight increase in the seasonal signal was observed from the Atlantic part at S1 (∼57 μatm, between 371 and 428 μatm) to TN at S3 (61 μatm, between 357 and 418 μatm) and a drastic decrease toward the east in GC at S4 (34 μatm, between 362 and 396 μatm).

The pH_T values (Figure 3 and Supplementary Figure 6) were correlated with fCO_2,sw (r² > 0.98 and ρ = −0.99) throughout the region and decreased by 0.001 U μatm^–1 at all stations. The maximum seasonality was found at S1 (8.028–8.080) and in the northern sections of S3 (8.037–8.095). The greatest differences at S2 were observed during the cold months between January and May and could be related with mixing, biological and other non-thermal processes. In GC at S4, the pH_T showed minimum pH_T seasonality ranging from 8.056 to 8.089 in the northern section and from 8.039 to 8.081 in the southern section.

The FCO₂ values (Figure 3 and Supplementary Figure 6) show a slight maximum outgassing and minimum ingassing in the northern part of each station. The outgassing weakened and the ingassing reinforced toward the east, reaching negative FCO₂ values along both the northern and the southern section during the entire first year of monitoring and in the northern part of S4 throughout the study period. A strong CO₂ sink was obtained at all stations during the cold periods, while the system behaved as a weak source during the warm months.

The influence of Mediterranean Sea waters on the distribution of SSS was only observed to the east of GC at S5 through the northern section (Supplementary Figure 7), where the values were maximum throughout the year (∼36.4–36.7). The west-to-east surface circulation from the Northeast Atlantic is reflected in the SSS between ES (S1) and GC (S4): the observed values and seasonality were minimum in the westernmost part at S1 (∼36.3–36.5) and slightly increased through the southern section between S2 and S4 (∼36.3–36.6).

The NC_T (Supplementary Figure 7) decreases between February and September and increases between October and January. The seasonality of C_T and NC_T in ES (S1) (∼30 μmol kg^–1, between 2090 and 2120 μmol kg^–1) increased through the southern section (∼40 μmol kg^–1, between 2090 and 2130 μmol kg^–1) and decreased through the northern section (∼15 μmol kg^–1, between 2100 and 2115 μmol kg^–1) in CS (S2) and TN (S3). In GC at S4, the seasonality decreased in the southern section (∼30 μmol kg^–1, between 2093 and 2123 μmol kg^–1) and increased in the northern section (∼22 μmol kg^–1, between 2098 and 2120 μmol kg^–1). Considering the temporal distribution of Chlorophyll a (Supplementary Figure 7), the seasonal change was reduced and the latitudinal differences increased toward the interior of the channel.

The current study provides a high spatio-temporal resolution evaluation of the surface CO₂ distribution based on high-frequency underway measurements obtained during 2 entire years of observation (February 2019–March 2021). More than 6000 data obtained during those 2 years is the same amount of data available in the SOCAT database for the last 20 years.

According with the high-resolution measurements carried out in this investigation, the seasonal changes in the CO₂ system and fluxes in the Strait of Gibraltar were controlled by the annual cycle of the Atlantic inflow and the position of the AMI. These processes has been shown to control other variables in the region (e.g., Bryden et al., 1994; Bray et al., 1995; Gómez et al., 2000; Skliris and Beckers, 2009). The shallower depth of the AMI favors the deep, cold and CO₂-rich Mediterranean water injection into the upper layers (de la Paz et al., 2009). The AMI elevates toward the northeast throughout the year and reaches minimum depth during late winter (February–March), coinciding with an observed decrease in SST (Figure 2A and Supplementary Figure 1) and increase in SSS and Chlorophyll a (Figure 2B and Supplementary Figure 2). The position of the AMI is also influenced by the tidal amplitude variation and internal wave generation (Gascard and Richez, 1985; Gómez et al., 2001), which favor intermittent vertical mixing events and increase the surface CO₂ during high tides (Wesson and Gregg, 1994; Macias et al., 2006) and especially during the most energetic spring tides (García-Lafuente et al., 2000).

The spatial variability of the CO₂ system during summer is also influenced by mixing processes driven by local climatology factors, as the wind-induced upwelling in the Gulf of Cádiz, the Strait of Gibraltar and the Alboran Sea at local (e.g., Peliz et al., 2009; Gómez-Jakobsen et al., 2019; Bolado-Penagos et al., 2020) and regional scale (e.g., Richez and Kergomard, 1990; Folkard et al., 1997; Stanichny et al., 2005). The signal of the wind-induced upwelling during the warm months was observed in minimum values of SST around CS and between TN and GC (Figure 2). The lowest fCO_2,sw values obtained during summer along the southern section in the westernmost part are explained by the enhancement of the biological uptake in this area due to wind-induced upwelling events at this time of the year (Stanichny et al., 2005). The opposite occurs to the east of TN, where the wind-induced upwelling weakened and higher fCO_2,sw were encountered along the southern route during summer (Figure 2C). The influence of the physical variability in the seasonality of fCO_2,sw in the area were described in term of a second degree function of SST in Eqs 12–16. The spatio-temporal variability observed in the CO₂ system and fluxes and their causes is consistent with previous studies based on datasets obtained from oceanographic surveys at certain locations and times of the year (Dafner et al., 2001; Santana-Casiano et al., 2002; de la Paz et al., 2008, 2009, 2011; Flecha et al., 2019).

The seasonality of the fCO_2,sw values obtained in this study agreed with that obtained from SOCAT data (Supplementary Figure 5) and both followed SST fluctuations. A high coincidence was found in the annual cycle, with an average difference in fCO_2,sw values between both databases of 6.94 ± 16.93 μatm. The fitting of the CanOA-VOS data and the standard error of estimate in the northern and southern section (Eqs 14, 15) and in the entire region (Eq. 12) were slightly enhanced with respect to the SOCAT data (Eq. 13) due to the higher temporal resolution in the entire annual cycle. The combination of the CanOA-VOS and SOCAT data provided an Eq. 16 that can be used to estimate fCO_2,sw in the surface waters of the Strait of Gibraltar. The fCO_2,sw is increasing in the Strait of Gibraltar at a rate of 2.35 ± 0.06 μatm year^–1, which is close to the Northeast Atlantic ESTOC site (1.92 ± 0.24 μatm year^–1) and similar to the Irminger Sea site (2.37 ± 0.49 μatm year^–1), where also vertical mixing processes contributed to the increased observed trend (Bates et al., 2014).

The coupling between thermal and non-thermal processes drives the seasonal variation of the CO₂ system throughout the region. These non-thermal processes are mainly biological production (including the organic matter production and remineralization) and physical processes that are strongly coupled affecting the fCO₂ distribution with both positive and negative effects (Takahashi et al., 1993). The T/B ratios (Table 2) show that fCO_2,sw seasonality in the Strait of Gibraltar was mainly controlled by thermal processes in the entire region, coinciding with previous studies in the Strait of Gibraltar (de la Paz et al., 2009) and along the northwest African coast (Curbelo-Hernández et al., 2021). However, the seasonality of fCO_2,thermal and fCO_{2,non–thermal} (Figure 3) and the T/B ratios obtained in the northern and the southern sections of each station show a high spatial variability in thermal and non-thermal coupling. The relevance of non-thermal processes in the temporal variation of fCO_2,sw (Figure 3) was greater along the southern section due to greater influence and intensity of upwelling throughout the African coast (Stanichny et al., 2005). The non-thermal effect was maximum and close to canceling the thermal effect to the south of CS, where the tidal dynamics is enhanced and vertical mixing processes are intense (Echevarría et al., 2002). The non-thermal effect decreased from S2 to S5 through the southern section and increased from S3 to S5 through the northern section, which is related to the high-intensity upwelling in the Northwest Alboran Sea and suppression of vertical mixing to the south (Minas et al., 1991; Echevarría et al., 2002; Gómez-Jakobsen et al., 2019).

The influence of thermal processes on the temporal variation of fCO_2,sw was also analyzed through the fCO_2,sw-SST relationship (Table 2). The increase in fCO_2,sw with temperature through the northern route (9.02 ± 1.99 μatm °C^–1, r² = 0.86, and ρ = 0.93) doubled that obtained through the southern route (4.51 ± 1.66 μatm °C^–1, r² = 0.48, and ρ = 0.69). Thus, the control of thermal processes was even higher in the northern section, while in the southern section the relevance of non-thermodynamic processes increased due to the enhancement of vertical mixing and biological activity (Santana-Casiano et al., 2002). The fCO_2,sw-SST relationship is consistent with the T/B ratios obtained at each station (Table 2). The highest correlation and change of fCO_2,sw with SST was obtained in the westernmost part at S1 and in the northern sections of S2 and S3 (>8.5 μatm °C^–1, r² > 0.67, and ρ > 0.82). In the southern sections of CS and TN, the change and correlation of fCO_2,sw-SST was found to be minimal (r² < 0.28 and ρ < 0.53).

The influence of non-thermal processes on the NC_T change was studied through the relationship NC_T-fCO_{2,non–thermal} (Table 2). The NC_T increased with respect to fCO_{2,non–thermal} by 0.50 ± 0.03 μmol kg^–1 μatm^–1 (r² > 0.92 and ρ > 0.96) throughout the region, but decreased with respect to fCO_2,thermal by 0.22 ± 0.06 (r² = 0.82 and ρ = −0.91) and 0.34 ± 0.06 (r² = 0.83 and ρ = −0.91) μmol kg^–1 μatm^–1 through the northern and southern routes, respectively. The NC_T showed a direct relationship and correlated with fCO_{2,non–thermal} (between 0.46 and 0.52 μmol kg^–1 μatm^–1; r² > 0.85 and ρ > 0.92) at all stations. Hence, the NC_T variation was mainly driven by non-thermal processes (Santana-Casiano et al., 2002; de la Paz et al., 2008, 2009).

The annual cycle of C_T was studied at each of the selected stations by considering its temporal change between consecutive months. The C_T was normalized to a constant salinity of 36.4 to remove the effect of precipitation and evaporation and study the other non-conservative processes involved in the surface variation of C_T (i.e., primary production, oxidation of organic matter, precipitation/dissolution of CaCO₃, air-sea CO₂ fluxes, and physical processes) (e.g., Chen and Pytkowicz, 1979; Takahashi et al., 1993; Wanninkhof and Feely, 1998; Lee et al., 2000).

The results allow to establish the relative contribution of biological activity, air-sea exchange, horizontal transport, and vertical mixing (Eq. 8) on the increase and depletion of NC_T in surface waters (Figure 4 and Table 3). The greatest seasonal changes of NC_T were found at S1 and throughout the southern section, especially at S2 and S3 (Table 3). A high spatial variability was found in the seasonality of NC_T. The depletion of NC_T in the southwestern part occurred from February–March to September, while it was difficult to establish a seasonal pattern through the northern section and in the easternmost part due to the high variability of the surface inflow in the tidal time-scale (de la Paz et al., 2009).

The seasonal change of NC_T was controlled by NCP processes at most of the stations except in the northern part at S2 (0.4–10.6%; Figure 4B), where the air-sea exchange and vertical mixing dominated (62.7–68.4 and 29.7–36.7%, respectively) and horizontal transport influenced (8.0–15.6%). It is related to the generation of internal bores at CS (Boyce, 1975; Armi and Farmer, 1988; Bruno et al., 2002) and the minimum intensity of the wind-induced upwelling toward the northwest of the channel (Peliz et al., 2009). The biological effect was especially significant along the African coastal upwelling (Stanichny et al., 2005) in the southern part of S2 (81.8–83.9%) and S3 (94.2–98.4%). The upwelling effect in the Northwest Alboran Sea (Minas et al., 1991; Echevarría et al., 2002; Gómez-Jakobsen et al., 2019) was detected at the easternmost stations (Figures 4D,E). The change in NC_T explained by NCP exceeded the total change of NC_T by 24.4–30.0 and 21.9–69.3% in the northern part at S4 and S5, respectively. The biological effect in these locations was counteracted by the contribution of air-sea exchange (12.3–33.6% in the northern part at S4 and S5 and 69.04–69.2% in the southern part at S5) and the influence of vertical mixing and horizontal transport (<12 and <10%, respectively). The contribution of NCP is strongly linked to the vertical mixing processes in this region, which inject deeper and remineralized CO₂-rich water in surface layers. Positive contributions of vertical mixing indicated that the NC_T increase in the surface layer due to the deep-water injection, while negative values indicate that the concentration of NC_T on the surface is higher than in depths layers due to a minimal biological production on the surface. In some periods, the amount of injected CO₂ to the surface is not counterbalanced by primary production and positive values contributed to the observed monthly variability. Accordingly, the seasonal and spatial variation of the AMI and the deep water injection processes in the surface layers drove by the wind (wind-induced upwelling) through the southernmost part and by tides (tidal-induced upwelling) through the northernmost part could account for the observed positive and negative contributions of both NCP and vertical mixing. The contribution of vertical mixing was especially significant in certain locations due to the advection of C_T in upper layers toward the east of CS during tidal-induced upwelling events (de la Paz et al., 2008). The diffusion processes of C_T through the pycnocline increased its surface concentration by an average of 3.55 ± 2.69 mmol m^–3 during the entire annual cycle throughout the region. The vertical entrainment increased the surface C_T by an average of 0.58 ± 0.52 mmol m^–3 between August and January and was canceled between February and July due to the elevation of the mixed layer and subsequent increase in stratification. About 5.7–5.8% of the seasonal change in NC_T at S1 (Figure 4C) was explained by vertical mixing processes. The vertical mixing inversely affected the change of NC_T at S3 and S5. During the periods of increase in NC_T, between 0.02 and 9.10% of the NC_T was removed from the surface by vertical mixing processes, while during periods of depletion, between 2.90 and 8.20% of the surface NC_T was injected from deep waters. The contribution of vertical transport was maximum in the northern part of S2 (∼29.7–36.7%) and at S4 (∼10.3–11.9%) (Figures 4B,D) due to the effect of internal waves in CS (Bruno et al., 2002; Macias et al., 2006) and the elevation of the pycnocline in the most Mediterranean part (Echevarría et al., 2002). The influence of the seasonality of the wind-induced upwelling on the change of NC_T by vertical mixing was detected in S2, where a difference of 1.34 mmol m^–3 between increase and depletion periods was found.

The horizontal transport in the west-to-east direction decreased the surface NC_T (0.8–15.6%) in the entire annual cycle due to the lower A_T and C_T content of the Atlantic inflow compared to the Mediterranean outflow (de la Paz et al., 2008). An increase in NC_T by horizontal transport was only observed at S1 (0.4%) and in the northern part of S4 (2.6%) and S5 (9.6%).

The average seasonal change of NC_T for all selected locations across the Strait of Gibraltar was 29.0 ± 9.96 mmol m^–3. Biological processes controlled the seasonal change of surface NC_T (93.5–95.6%), increasing its concentration by 27.72 mmol m^–3 in periods of higher remineralization influence and decreasing it by 27.11 mmol m^–3 when production was the controlling process. The contribution of air-sea CO₂ exchange (0.9–3.3%) was higher during the increasing periods due to the strong behavior as a sink of the channel (0.94 and 0.27 mmol m^–3 during the increase and depletion periods, respectively). Vertical mixing increased the NC_T by 1.32 mmol m^–3 and decreased it by 1.02 mmol m^–3, with a total contribution of 4.1%. Finally, the effect of horizontal transport was minimal (<2.8%) and decreased the NC_T by 1.08 and 0.52 mmol m^–3 in both total increase and depletion periods, respectively.

The temporal distribution of FCO₂ calculated from observed and computed data allowed the characterization of the air-sea CO₂ exchange in the Strait of Gibraltar throughout the study period. The monthly average FCO₂ values (Supplementary Figure 8) proved the strong behavior of the region as a CO₂ sink from November to July and as a weak source between August and October. This seasonal pattern is consistent with previous studies in the Strait of Gibraltar (de la Paz et al., 2009, 2011) and along the coastal to open-ocean transitional area of the Northeast Atlantic (Curbelo-Hernández et al., 2021). In the northern section, the CO₂ sink behavior remained throughout the first year of observation (February 2019–February 2020) even in the warm months, while a late-summer CO₂ source behavior was found in the second year of observation (March 2020–March 2021). In contrast, the southern section behaved as a late-summer CO₂ source in both entire years considered.

The FCO₂ values show a direct relationship with ΔfCO₂ in both the northern and the southern sections [0.11 (r² = 0.87 and ρ = 0.93) and 0.14 (r² = 0.72 and ρ = 0.84) mmol m^–2 d^–1 μatm^–1, respectively]. Despite the important role of the wind speed and direction controlling the seasonal and spatial distribution and intensity of upwelling events across coastal areas in the Strait of Gibraltar (Stanichny et al., 2005), its short-term variability is not significant on an annual time scale in this region (de la Paz et al., 2011) and specially along the northern section. The difference between the averages FCO₂ calculated with satellite and measured wind speed data was 0.14 mol C m^–2 year^–1. The changes in wind speed influence the intensity of CO₂ exchange between low atmosphere and surface ocean and could modified the physical properties of the upper layers, but do not determine the formation of CO₂ source/sink. Northeasterly winds bring on the deep-water injection in upper layers and subsequent cooling of the surface seawater in the southwest area, favoring the solubility of atmospheric CO₂. Accordingly, a more intense annual CO₂ sink was found in this area compared to the northwest area (Figure 5). The assessment of the spatio-temporal changes of FCO₂ in this study show that the formation of the seasonal CO₂ sink and source is controlled by the ΔfCO₂ variations influenced by changes in the physical properties of the surface ocean, which were driven by the local climatology and the unique hydrology of the channel (Supplementary Figure 8). Therefore, Eq. 16 which provides fCO_2,sw data for the area using SST as unknown variable (available from satellite data) together with satellite winds, can be used to estimate long term CO₂ fluxes in the Strait of Gibraltar.

The annual negative values of FCO₂ obtained along both the northern and the southern section show that the Strait of Gibraltar behaved as a net CO₂ sink (Figure 5). The annual average FCO₂, ΔfCO₂, and wind speed were calculated at all the selected stations, along the northern and southern sections and in the entire region (Table 4). The average ingassing was found to be maximum in the eastern area (Figure 5). The total annual FCO₂ in most of the Northeast Atlantic and Mediterranean Sea parts of the channel was −1.02 and −1.30 mol C m^–2 year^–1, respectively, which is related with lower SST values observed throughout the year and the enhancement of the biological uptake toward the northeastern part of the channel (Echevarría et al., 2002) that increased the ΔfCO₂.

The net CO₂ sink obtained in this study between 2019 and 2021 for the entire region was −0.98 mol C m^–2 year^–1, which falls in the range of the first annual FCO₂ provided by Santana-Casiano et al. (2002) (−2.5 mol C m^–2 year^–1) and the annual FCO₂ reported in most recent studies of de la Paz et al. (2009) (−0.28 mol C m^–2 year^–1) and de la Paz et al. (2011) (−0.02 ± 0.13 mol C m^–2 year^–1). The annual FCO₂ was calculated in these studies from limited datasets and using the parameterization of Wanninkhof (1992) (W92). When applying this parameterization to the data collected, an annual value of FCO₂ of −1.21 mol C m^–2 year^–1 was obtained. This average ingassing is more than two times higher than the values reported with W92 for the Gulf of Cadiz (−0.4 mol C m^–2 year^–1; Borges et al., 2006) and for the Northwest African Atlantic coast (−0.59 mol C m^–2 year^–1, using W92; Curbelo-Hernández et al., 2021). However, it falls within the annual FCO₂ estimated by Chen et al. (2013) for the Northeast Atlantic (−1.04 mol C m^–2 year^–1) and for the Iberian upwelling (−1.33 mol C m^–2 year^–1). The net CO₂ sink reported in this study is consistent with that obtained from SOCAT data referenced to 2019 (−0.97 mol C m^–2 year^–1; using W92).

The net FCO₂ for the entire area (1650 km²) was −7.12 Gg CO₂ year^–1 (1 Gg = 10⁹ g) (−1.94 Gg C year^–1). This result is in agreement with previous studies of the behavior of the coastal regions and continental shelves as a net CO₂ sink (e.g., Borges et al., 2006, 2005; Chen and Borges, 2009; Laruelle et al., 2010, 2013, 2018; Cai et al., 2006; Cai, 2011; Chen et al., 2013) and can be integrated into the estimated net CO₂ sink for the European coastal regions (−68.1 × 10³ Gg C year^–1; Borges et al., 2006) and global continental shelves (−0.4 × 10⁶ Gg C year^–1; Chen et al., 2013).