Spatiotemporal variations in surface marine carbonate system properties across the western Mediterranean Sea using volunteer observing ship data

Surface physical and marine carbonate system (MCS) properties were assessed along the western boundary of the Mediterranean Sea. An unprecedented high-resolution observation-based dataset spanning 5 years (2019–2024) was built through automatic underway monitoring by a volunteer observing ship (VOS). The MCS dynamics were strongly modulated by physical–biological coupling dependent on the upper-layer circulation and mesoscale features. The variations in CO₂ fugacity (fCO_2,sw) were mainly driven by sea surface temperature (SST) changes. On a seasonal scale, SST explained 51 %–71 % of the increase in fCO_2,sw from February to September, while total alkalinity (A_T) and sea surface salinity (SSS) explained < 20 %. The processes controlling total inorganic carbon (C_T) partially offset this increment and explain ∼ 23 %–37 % of the fCO_2,sw seasonal change. On an interannual scale, the SST trends (0.26–0.43 °C yr⁻¹) have accelerated by 78 %–88 % in comparison with previous decades. The ongoing surface warming contributed ∼ 76 %–92 % to increasing fCO_2,sw (4.18 to 5.53 µatm yr⁻¹) and, consequently, decreasing pH (−0.005 to −0.007 units per year) in the surface waters. The seasonal amplitude of SST, which is becoming larger due to progressively warmer summers, was the primary driver of the observed slope of interannual trends. The evaluation of the air–sea CO₂ exchange shows the area across the Alboran Sea (14 000 km²) and the eastern Iberian margin (40 000 km²) acting as an atmospheric CO₂ sink of −1.57 ± 0.49 mol m⁻² yr⁻¹ (−0.97 ± 0.30 Tg CO₂ yr⁻¹) and −0.70 ± 0.54 mol m⁻² yr⁻¹ (−1.22 ± 0.95 Tg CO₂ yr⁻¹), respectively. Considering the spatial variability of CO₂ fluxes across the study area, a reduction of approximately 40 %–80 % in the net annual CO₂ sink has been estimated since 2019, which is attributed to the persistent strengthening of the source status during summer and the weakening of the sink status during spring and autumn.

1 Introduction

Semi-enclosed and marginal seas play a relevant role in the global biogeochemical cycles and are highly vulnerable to climate change (IPCC, 2023). These regions encompass extensive coastal and continental shelf and slope areas occupied by multiple diverse ecosystems under anthropogenic pressure. Although these regions present enhanced biogeochemical activity and intensified air–sea CO₂ exchange rates compared to the open ocean (Borges et al., 2005; Cai et al., 2006; Frankignoulle and Borges, 2001; Shadwick et al., 2010), its poor monitoring and assessment have historically excluded them from global studies and models and underestimated them in the Global Carbon Budget (Friedlingstein et al., 2023).

The Mediterranean Sea is a dynamic semi-enclosed system potentially vulnerable to natural and anthropogenic forcing (e.g. Álvarez et al., 2014; Tanhua et al., 2013). The particular oceanography of the Mediterranean Sea, collectively described in several works (e.g. Nielsen, 1912; Robinson et al., 2001; Millot and Taupier-Letage, 2005; Bergamasco and Malanotte-Rizzoli, 2010; Schroeder et al., 2012), has rendered it a “miniature ocean” considered a “laboratory basin” for evaluating physicochemical perturbations that can be extrapolated to larger scales in the global ocean (e.g. Robinson and Golnaraghi, 1994; Bergamasco and Malanotte-Rizzoli, 2010). These perturbations have accelerated since the second half of the 20th century, with temperature and salinity increasing at unprecedented rates of 0.04 °C and 0.015 per decade, respectively (Borghini et al., 2014), impacting the marine carbonate system (MCS). However, the availability of high-quality observation-based data and research in this basin is low due to spatial and temporal limitations in the monitoring and sampling techniques (Millero et al., 1979; Rivaro et al., 2010).

The MCS dynamics have been evaluated in the north-western Mediterranean basin (Bégovic and Copin-Montégut, 2002; Copin-Montégut and Bégovic, 2002; Copin-Montégut et al., 2004; Coppola et al., 2020; Hood and Merlivat, 2001; Mémery et al., 2002; Merlivat et al., 2018; Touratier and Goyet, 2009; Ulses et al., 2023), mainly at the time series DYFAMED (43.42° N, 7.87° E; Marty, 2002) and BOUSSOLE sites (43.37° N, 7.90° E; Antoine et al., 2006, 2008a, b). These investigations have shown that the seasonal cycle of the surface CO₂ is primarily governed by thermal fluctuations and the behaviour of the area as a relatively weak sink for atmospheric CO₂ on an annual scale. Long-term changes estimated by Merlivat et al. (2018) reported increases in the surface CO₂ fugacity (fCO_2,sw) and pH of ∼ 40 µatm and ∼ 0.04 units, respectively, since the 1990s. The interannual trends given for fCO_2,sw (2.3 ± 0.23 µatm yr⁻¹; Merlivat et al., 2018) and pH (0.002–0.003 units per year; Yao et al., 2016) were in agreement with those encountered in the north-eastern Atlantic at the ESTOC site (2.1 ± 0.1 µatm yr⁻¹ and 0.002 ± 0.0001 units per year, respectively; González-Dávila and Santana-Casiano, 2023). Long-term variations in the MCS within the north-western Mediterranean occur at rates exceeding those anticipated from chemical equilibrium with atmospheric CO₂, which has been attributed to the intense deep-convection processes in this area (Copin-Montégut, 1993; D'Ortenzio et al., 2008; Cossarini et al., 2021) and the substantial input of anthropogenic carbon from the North Atlantic (Merlivat et al., 2018; Palmiéri et al., 2015; Schneider et al., 2010; Ulses et al., 2023). Based on a high-resolution regional model, Palmiéri et al. (2015) estimated that ∼ 25 % of the anthropogenic carbon storage in the Mediterranean Sea comes from the Atlantic. The water exchange processes in the Strait of Gibraltar become the western boundary of the Mediterranean Sea in a crucial region for MCS variability, which significantly modulates the basin-wide anthropogenic carbon inventory and ocean acidification trends in the Mediterranean basin and could significantly affect the general circulation and composition of seawaters in the North Atlantic. Additionally, this region is subject to variability related to (1) the intense deep-water convection in the adjacent north-western area of the Mediterranean Sea and (2) the unique circulation patterns shaped to the irregular coastlines and islands, which form quasi-permanent eddies and other (sub-)mesoscale features (Alberola et al., 1995; Bosse et al., 2021, 2016; Bourg and Molcard, 2021).

This research focuses on the surface spatiotemporal variations of the MCS and air–sea CO₂ fluxes on the western boundary of the Mediterranean Sea. High-resolution and reliable data were obtained through autonomous underway monitoring of the ocean surface from February 2019 to February 2024 by a volunteer observing ship (VOS). This systematic strategy represents a powerful tool for analysing the distribution and changes in physical and MCS properties in highly different areas such as coastal transitional zones, where the availability of data has historically been scarce. The cruise track (Fig. 1) followed the southern and eastern geographically rugged coastline of the Iberian Peninsula and allowed characterization of the Alboran Sea (∼ 2–5.1° W) that was separate from the eastern coastal shelf area between the Cape of Gata (Almería) and Barcelona (∼ 36.5–41.3° N). The changes observed in the MCS on seasonal and interannual timescales (even considering the limitations of 5 years of data), the mechanisms controlling their variations, and the changes in the air–sea CO₂ exchange have been addressed in this study.

Figure 1 (a) Map of the western boundary of the Mediterranean Sea with the ES-SOOP-CanOA tracks between February 2019 and February 2024 (red) and the locations of the stations of interest along the southern (S1–S5) and eastern (E1–E6) sections. The main capes and gulf along the geographically rugged Iberian coastline are shown. The schematic diagram summarizes the classical circulation patterns: in the Alboran Sea (blue), the Atlantic Jet (AJ) surrounds the Western Anticyclonic Gyre and Eastern Anticyclonic Gyre (WAG and EAG, respectively) and forms Modified Atlantic Water (MAW), while along the eastern Iberian margin (purple) the Mediterranean Water (MW) is transported from the north-western Mediterranean basin along the path of the Northern Current. The northward spreading of MAW during summer and southward spreading of MW during winter are depicted with the dashed arrows. The thermal front that formed in the axis of the Pyrenees during summer is depicted with a black dashed line. (b) SST maps built with reanalysis monthly data (0.042° × 0.042°) for February and September 2023 from the Med MFC physical multiyear product (Escudier et al., 2020, 2021; Nigam et al., 2021) available at the Copernicus Marine Data Store (https://data.marine.copernicus.eu/products; last access: 15 May 2025).

2 Material and methods

2.1 Study area

The western boundary of the Mediterranean Sea encompasses the Alboran Sea, land-loaded by the southern Iberian Peninsula coast, the northern African coast, and the coastal transitional area along the eastern Iberian margin (Fig. 1a). The classical surface circulation pattern in the Alboran Sea (e.g. Bormans and Garrett, 1989; Peliz et al., 2013; Sánchez-Garrido et al., 2013; Sánchez-Garrido and Nadal, 2022; Speich et al., 1996; Whitehead and Miller, 1979), together with the Atlantic Jet (AJ) following a wavelike path of the quasi-permanent Western Anticyclonic Gyre (WAG) and Eastern Anticyclonic Gyre (EAG) and constituting the Modified Atlantic Water (MAW; Lopez-García et al., 1994; Viúdez et al., 1998), drives west–east variations in physical and biogeochemical terms. The intensity and direction of the AJ, depending primarily on sea level pressure and local wind fluctuations, vary on different timescales and govern the circulation patterns in the Alboran Sea, influencing the biogeochemistry (Sánchez-Garrido and Nadal, 2022; Solé et al., 2016). On a seasonal scale, the AJ oscillates between the two main circulation modes (García Lafuente et al., 2002; Macías et al., 2008, 2016; Vargas-Yáez et al., 2002), which is detectable by reanalysis-data-based SST signals (Fig. 1b): a high-intensity AJ flowing north-eastward during spring and summer and a lower-intensity AJ flowing with a more south-easterly direction during autumn and winter. The stronger AJ during the warm months feeds the classical two-gyre configuration in the Alboran Sea, while the weak AJ only allows the existence of the WAG (Renault et al., 2012). The AJ forms a filament flowing from the Iberian coastal upwelling in the north-western Alboran Sea and surrounding the eastern edge of the WAG, which is most frequently present during summer (Gómez-jakobsen et al., 2019; Millot, 1999). The westernmost part of the Alboran Sea is affected by the shallow position of the Atlantic–Meridional Interface (AMI) layer (Bray et al., 1995; Lacombe and Richez, 1982), which promotes the injection of deep water into the surface (Echevarría et al., 2002; Gómez-jakobsen et al., 2019; Minas et al., 1991).

The eastern Iberian margin is influenced by the path of the Northern Current transporting Mediterranean Water (MW; Pinot et al., 1995), which originates around the Gulf of Lion, where the forcing of the north-easterly winds is frequently strong and flows southward along the eastern coastline of the Iberian Peninsula (Conan and Millot, 1995; Millot, 1999; Sammari et al., 1995). The seasonality of the Northern Current (Millot, 1999) infers meridional variations in the thermal signals between cold and warm months (Fig. 1b). The enhanced wind forcing during winter intensifies the Northern Current, which fits the Iberian continental slope and recirculates offshore of the Cape of Nao (Millot, 1999), while a low-intensity branch progresses southward of the Cape of Nao and reaches the eastern Alboran Sea. During summer, the weakening of the wind forcing forms a surface thermal front in the axis of the Pyrenees that was detectable on the reanalysis-based SST map (Fig. 1b). This front changes the path of the Northern Current further away from the Iberian coast (Lopez-García et al., 1994), which allows the MAW to reach its northernmost spreading. The interaction of the Northern Current with the variety of mesoscale features (mainly meanders and eddies) and the variations in stratification within the annual cycle introduced spatiotemporal differences into the biogeochemical properties (Bosse et al., 2021; Millot, 1999). Additionally, although terrestrial and riverine inputs have a less pronounced impact on biogeochemistry compared to the eastern Mediterranean basin (Cossarini et al., 2015), they can act as a source of local variability. Most significant in this area is the Ebro River runoff, which peaks in March–May due to the combined action of precipitation during winter and snowmelt in the upper river basins during spring (Zambrano-Bigiarini et al., 2010). It feeds the coastal area around the Ebro Delta with fresh and cool waters (see the minimum SST compared to the adjacent waters in February; Fig. 1b).

2.2 Data collection

A high-spatiotemporal-resolution dataset spanning 5 years was constructed based on weekly physicochemical observations of the surface western boundary of the Mediterranean Sea between February 2019 and February 2024. Data were collected automatically by the VOS MV JONA SOPHIE (IMO: 9144718, called RENATE P before November 2021), a container ship managed in Spain by Nisa Maritima, which links the Canary Islands with Barcelona. This VOS line was designed and is maintained by the QUIMA research group at IOCAG-ULPGC, and it has operated within the framework of the Integrated Carbon Observation System (ICOS; https://www.icos-cp.eu/; last assess: 15 May 2025) as a Ship of Opportunity (SOOP) Ocean Station (station ID: ES-SOOP-CanOA) since 2021 (upgraded to an ICOS Class-1 Ocean Station in May 2024). Therefore, the measurement equipment and underway data collection techniques verify the ICOS high-quality requirements and methodological recommendations.

The ES-SOOP-CanOA station allows monitoring of a coastal transitional zone transect across the western Mediterranean Sea (Fig. 1), together with a north-eastern Atlantic subtropical area (Curbelo-Hernández et al., 2021a) and the Strait of Gibraltar (Curbelo-Hernández et al., 2021b). In the Alboran Sea, the vessel advanced eastward and longitudinally crossed the WAG through its northern part and followed the northern path of the EAG. The irregular south-eastern and eastern coastlines of the Iberian Peninsula caused local differences in the oceanographic features and variances in the distance to land of the vessel track.

The system operates fully unattended in underway mode, with biweekly (the time required to complete a round trip) routine maintenance at the port of Las Palmas de Gran Canaria (28.13° N, 15.42° W). Data are automatically transferred to a server when the vessel docks at each of the ports along the usual route (Las Palmas de Gran Canaria, Santa Cruz de Tenerife, Arrecife, Sagunto, and Barcelona). A total of 92 routes were completed in the Mediterranean Sea (Fig. 1).

2.3 Monitoring routines

The autonomous underway monitoring of CO₂ in the surface ocean (water intake placed at 5 m depth) and low atmosphere (air intake placed at 8 m above sea level) and the data collection routines followed the recommendations described by Pierrot et al. (2009) to ensure comparable and high-quality datasets. An automated underway CO₂ molar fraction (xCO₂, ppm) measurement system, developed by Craig Nail and commercialized by General Oceanics, was installed in the engine room of the vessel and described by Curbelo-Hernández et al. (2021a, b).

The xCO₂ measurement system combines an air and seawater equilibrator placed inside the wet box, with a non-dispersive infrared analyser for gas detection placed inside the dry box. The analyser used for xCO₂ detection was built by LICOR^® (initially the LI-6262 model and, after October 2019, the LI-7000 model). The nominal accuracy of the LICOR infrared gas analyser given by the manufacturer is 1 % for CO₂ concentrations within the range 0–3000 ppm. The system performs, in-loop at 3 min intervals, 5 measurements of atmospheric xCO₂ (xCO_2,atm) and 80 measurements of surface seawater xCO₂ (xCO_2,sw). The xCO_2,atm data were consistent with daily xCO_2,atm records from the Izaña Atmospheric Research Center (IZO) site located in Tenerife, Canary Islands, Spain (28.3090° N, 16.499° W, located 2372.9 m above sea level; https://gml.noaa.gov/dv/site/site.php?code=IZO, last access: 14 May 2025), which is operated by the Spanish Meteorological Agency (AEMET) and forms part of several major international atmospheric monitoring networks (Fig. S1). Daily xCO_2,atm data from IZO are available through the National Ocean and Atmospheric Administration (NOAA) Global Monitoring Laboratory (GML) dataset (https://gml.noaa.gov/data/dataset.php?item=izo-co2-flask; last access: 14 May 2025). From 2019 to 2024, xCO_2,atm measurements from the ES-SOOP-CanOA station were, on average, 1.14 ppm higher than those recorded at IZO (Fig. S1), which may be due to the fact that air sampling at IZO is conducted at approximately 2400 m above sea level in a remote location far from major urban or industrial areas and above the atmospheric inversion layer, which shields the station from surface level pollution. In contrast, the ES-SOOP-CanOA measurements are conducted in the lower atmosphere, near the sea surface, and closer to greenhouse gas emission sources (particularly when the vessel operates near the coast in the Mediterranean basin).

The LICOR^® analyser is automatically calibrated on departure and arrival at each port and periodically every 3 h using four standard gases. Additionally, the system is zeroed and spanned (with standard gases 1 and 4, respectively) every 12 h to properly interpolate the standard values and correct for instrument drift. The four standard gases, with an accuracy of ± 0.02 ppm, were provided by the NOAA and are traceable to the World Meteorological Organization (WMO). They were in the order 0, 250, 400, and 550 ppm until January 2021, when the gas bottles for standards 2 to 4 were exchanged for a new set with concentrations in the order 300, 500, and 800 ppm as provided by the ICOS central analytical laboratories.

The sea surface temperature (SST, °C) was monitored using a SBE38 thermometer placed at the primary seawater intake in the engine room, with a reported instrument error of ± 0.01 °C. The high sensitivity of xCO₂ to temperature fluctuations required monitoring of temperature at different locations across the system. A SBE45 thermosalinograph and a Hart Scientific HT1523 handheld thermometer, with reported instrument errors of ± 0.01 °C, were used to monitor the seawater temperature at the entrance of the wet box and inside the equilibrator, respectively. The measured SST was analysed in conjunction with SST monthly reanalysis data (0.042° × 0.042° spanning 24 years from 1 January 2000 to 1 March 2024) from the Med MFC physical multiyear product (Escudier et al., 2020, 2021; Nigam et al., 2021) available at the Copernicus Marine Data Store (https://data.marine.copernicus.eu/products; last access: 15 May 2025). The SST reanalysis data were interpolated to the coordinates of the ES-SOOP-CanOA data to perform direct comparison of their dynamics.

The sea surface salinity (SSS) was measured by the SBE45 thermosalinograph, whose instrumental error is of the order of ± 0.005. Lastly, pressure is measured within ± 0.0002 atm at the deck box transducer close to the air intake (used as atmospheric pressure), in the wet box inside the equilibrator at the time of equilibration, and in the dry box to be used by the LICOR analyser to correct the analogue signal for any pressure effects.

Discrete surface seawater samples were collected manually with in situ records of SST and SSS during three round trips in February 2020, March 2021, and October 2023 (a total of 102 were collected in the Mediterranean Sea). The discrete sampling was performed along the vessel track from the seawater supply line every 1–2 h in borosilicate glass bottles, overfilled, and preserved with 100 µL of saturated HgCl₂. Samples were kept in the dark and analysed just after arriving at port, in a period of less than 2 weeks, for total alkalinity (A_T, µmol kg⁻¹).

The ongoing observational dataset exhibits a gap of 1 year between September 2021 and September 2022 due to the temporary cessation of the measurement system for vessel maintenance in dry dock. During this period, the measurement system was sent for calibration and maintenance to General Oceanics, Miami, USA. There are also several gaps of less than 1 month related to different technical issues with the measurement equipment that were addressed during the routine maintenance visits to the vessel (i.e. problems with the pump and seawater intake and with the LICOR analyser, depletion of gas bottle supplies, and electrical issues in the engine room). Resolution of certain technical issues encountered during 2020 was delayed due to the constraints imposed by COVID-19.

2.4 Calculation procedures

2.4.1 CO₂ system variables

Data processing for fCO₂ calculations

The present investigation 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 correction of the measured xCO₂ and calculation of the fugacity of CO₂ (fCO₂) in surface seawater (fCO_2,sw) and the atmosphere (fCO_2,atm) followed the procedure described by Pierrot et al. (2009). This procedure avoids significant uncertainties in the determination of fCO₂ arising from differences in pressure and temperature conditions between sampling (atmospheric pressure and SST) and equilibration (pressure and seawater temperature inside the equilibrator once equilibration is reached). By calibrating the instrument with standard gases ranging from 0 to 800 ppm (which encompasses the measurement range of 300 to 600 ppm) and actively minimizing temperature and pressure drift through continuous monitoring (see Sect. 2.3 for the standard gas, temperature, and pressure accuracies), the system achieved the target accuracies of ± 0.2 µatm for fCO_2,atm and ± 2 µatm for fCO_2,sw (Pierrot et al., 2009). The full set of standard gases was linearly interpolated to the time of the observations to generate the calibration curve used for xCO₂ correction before calculating fCO₂.

The raw output data were initially filtered by removing data affected by the automatic sampler, such as samples measured at low water rates (< 2.0 L min⁻¹) and/or samples in which the difference in temperature between the seawater intake and the equilibrator was higher than 1.5 °C. The outliers, assumed to be elements of more than 3 local standard deviations from the local mean over a window length of 50 elements, were also removed from the dataset. The xCO₂-measured values in the low atmosphere after each calibration were averaged and interpolated at the time of each xCO₂ observation in seawater by applying a piece-wise polynomial-based smoothing spline.

A_T determination and reconstruction

Discrete seawater samples were analysed for A_T using a VINDTA 3C and following the procedure detailed by Mintrop et al. (2000). The VINDTA 3C was calibrated through the titration of certified reference materials (CRMs, provided by Andrew G. Dickson at Scripps Institution of Oceanography), giving values with an accuracy of ± 1.5 µmol kg⁻¹.

An empirical salinity-based relationship was developed to reconstruct A_T specifically for the monitored transect. The A_T–SSS linear relationship (Eq. 1), derived from 46 discrete samples, is statistically significant at the 99 % level of confidence (p value < 0.01) and presents a high degree of correlation (r²=0.99) and a root mean square error (RMSE) of ± 5.6 µmol kg⁻¹.

$$\begin{array}{}\text{(1)}& A_{\mathrm{T}}=\mathrm{100.5}\left(\pm \mathrm{2.9}\right)\mathrm{SSS}-\mathrm{1271}\left(\pm \mathrm{108}\right)+\mathit{\epsilon }\end{array}$$

Although the reconstruction of A_T from its linear relationship with SSS does not account for biological processes that cannot be traced with salinity (Wolf-Gladrow et al., 2007) or the input of dissolved carbonate minerals and bicarbonate–carbonate species from river runoff, sediments, and water mixing, it has been widely used. It provides useful general approximations in regions with stable conditions and limited influence from these processes. The strong correlation of the A_T–SSS linear relationship and its consistency with those proposed for various regions of the Mediterranean Sea (Schneider et al., 2007; Copin-Montégut and Bégovic, 2002; Jiang et al., 2014; Cossarini et al., 2015) indicate a dominant control of salinity-driven processes over A_T variability. In the western Mediterranean, these processes primarily include the conservative effects of evaporation–precipitation balance and the inflow of cooler, less saline Atlantic waters (Jiang et al., 2014). River runoff also causes local changes in SSS, but its role as a salinity-independent source of A_T variability is minimal compared to the eastern Mediterranean basin (Jiang et al., 2014; Cossarini et al., 2015). With respect to biogeochemical processes, they cause negligible changes in salinity and induce small direct changes in A_T.

To account for minor variations in A_T that are independent of salinity, a non-conservative term (ε) was included in the linear model (Eq. 1), representing the residuals (the difference between measured A_T values and those predicted from SSS) and capturing variability not explained by salinity alone. A_T was calculated at the times of the SSS observations (Curbelo-Hernández et al., 2021a, b, 2023) using Eq. (1). A Monte Carlo simulation was employed to propagate the uncertainties from both the linear model parameters (slope and intercept) and the ε term: 10 000 simulations were conducted for each SSS observation by generating random realizations of the regression parameters from normal distributions centred at their best-fit values with standard deviations equal to their respective standard errors. Similarly, random realizations of ε were performed from a normal distribution characterized by the mean and standard deviation of ε. This yielded a full probability distribution of estimated A_T values for each observed SSS, from which the ensemble mean, standard deviation, and 95 % confidence intervals were computed. The propagated uncertainty in A_T estimates, considering the errors in A_T determination and SSS measurements (Sect. 2.3) and the linear model uncertainty, was approximately ± 5.7 µmol kg⁻¹. This error in A_T estimation falls within the accepted uncertainty range of ± 10 µmol kg⁻¹ for A_T when used as an input variable alongside fCO_2,sw (when its uncertainty is up to ± 2 µatm) for the calculation of other MCS variables aligning with the criteria for the “weather goal” level of measurement quality (Steinhoff and Skjelvan, 2020). This new A_T–SSS relationship is applicable for estimating A_T in surface seawaters within the western Mediterranean Sea that are predominantly influenced by conservative processes and where non-salinity factors such as biological activity and fluvial inputs are limited. This includes waters with salinities in the range of 36 to 38.5.

pH and C_T calculation

The pH and C_T were calculated at the times of the underway observations using the CO_2SYS program developed by Lewis and Wallace (1998) and run with the MATLAB software (van Heuven et al., 2011; Orr et al., 2018; Sharp et al., 2023). The fCO_2,sw and A_T were used as input CO₂ system variables. The set of constants used for the computations includes the carbonic acid (H₂CO₃) dissociation constants of Lueker et al. (2000), the sulfuric acid (HSO₄) dissociation constant of Dickson (1990), the hydrofluoric acid (HF) dissociation constant of Pérez and Fraga (1987), and the total boron concentration ([B]_T) determined by Lee et al. (2010). The effect of temperature on pH was removed by computation at a constant temperature of 19 °C, which is the mean temperature within the observational period (referred to as pH₁₉). Further data adjustments and statistical procedures are detailed in Appendix A.

2.4.2 Thermal and non-thermal fCO_2,sw

The relative influence of the thermal and non-thermal processes on the variation of fCO_2,sw has been addressed. The non-thermal processes mainly include the biological and carbonate pumps, circulation patterns, and air–sea gas exchange (De Carlo et al., 2013). The collectively known methodology presented by Takahashi et al. (2002) with the experimentally determined temperature effects on pCO₂ for isochemical seawater of 0.0423 °C⁻¹ (Takahashi et al., 1993) was used. This procedure was previously applied to ES-SOOP-CanOA data and detailed by Curbelo-Hernández et al. (2021a, b). An alternative procedure recently introduced by Fassbender et al. (2022), detailed by Rodgers et al. (2023), and modified from the Takahashi et al. (2002, 1993) framework was also applied in this investigation. This updated method addresses the slight variations in the thermal sensitivity of fCO_2,sw due to background chemistry (Wanninkhof et al., 1999, 2022), which introduce a slight difference between the observed seasonal cycle of fCO_2,sw and the calculated seasonal cycle through the sum of its thermal and non-thermal components. The Takahashi et al. (2002) and Fassbender et al. (2022) procedures are referred to hereinafter as T'02 and F'22, respectively.

The new approach in F'22 for the thermal component of fCO_2,sw (fCO_{2, T FASS}) was computed from the annual means (denoted with the subscript “AM”) of SSS, A_T, and C_T at in situ temperature (Eq. 2) using the CO_2SYS program (Lewis and Wallace, 1998) for MATLAB (van Heuven et al., 2011; Orr et al., 2018; Sharp et al., 2023). Then, the thermally driven change in fCO_2,sw (fCO_{2, T anom}) can be calculated as the difference between the thermal component of fCO_2,sw (fCO_{2, T FASS}) and the annual mean of fCO_2,sw (Eq. 3).

$$\begin{array}{}\text{(2)}& {\displaystyle }f{\mathrm{CO}}_{\mathrm{2},\mathrm{T}\mathrm{FASS}}={\mathrm{CO}}_{\mathrm{2},\mathrm{SYS}}\left(C_{\mathrm{T},\mathrm{AM}}{A}_{\mathrm{T},\mathrm{AM}},{\mathrm{SSS}}_{\mathrm{AM}},\mathrm{SST}\right)\text{(3)}& {\displaystyle }f{\mathrm{CO}}_{\mathrm{2},\mathrm{T}\mathrm{anom}}=f{\mathrm{CO}}_{\mathrm{2},\mathrm{T}\mathrm{FASS}}-f{\mathrm{CO}}_{\mathrm{2},\mathrm{AM}}\end{array}$$

The new approach in F'22 for the non-thermal component (fCO_{2, NT FASS}) is given by the difference between fCO_2,sw at the times of the observations and fCO_{2, T anom} (Eq. 4). The difference between fCO_{2, NT FASS} and the annual mean of fCO_2,sw provides the change in fCO_2,sw that is explained by non-thermal processes (fCO_{2, NT anom}) (Eq. 5).

$$\begin{array}{}\text{(4)}& {\displaystyle }f{\mathrm{CO}}_{\mathrm{2},\mathrm{NT}\mathrm{FASS}}=f{\mathrm{CO}}_{\mathrm{2},\mathrm{sw}}-f{\mathrm{CO}}_{\mathrm{2},\mathrm{T}\mathrm{anom}}\text{(5)}& {\displaystyle }f{\mathrm{CO}}_{\mathrm{2},\mathrm{NT}\mathrm{anom}}=f{\mathrm{CO}}_{\mathrm{2},\mathrm{NT}\mathrm{FASS}}-f{\mathrm{CO}}_{\mathrm{2},\mathrm{AM}}\end{array}$$

Considering the seasonal amplitudes of fCO_2,T and fCO_2,NT (dfCO_2,T and dfCO_2,NT), the relative importance of thermal and non-thermal processes was expressed by the $T/B$ ratio (dfCO_2,T $/$ dfCO_2,NT), with values greater than 1 indicating that the temperature effect governs the fCO_2,sw variations.

2.4.3 Factors controlling the seasonal amplitude of fCO_2,sw

The changes in the surface fCO_2,sw result from the combined variation in the physical and biochemical seawater properties. The seasonal variability of the surface fCO_2,sw was addressed by attending to the partial contributions of SST, SSS, C_T, and A_T (e.g. Takahashi et al., 2014). The influence of each driver was quantified by assuming linearity and employing a first-order Taylor series deconvolution (Sarmiento and Gruber, 2006) given in Eq. (6) and previously used for pCO₂ (Doney et al., 2009; Lovenduski et al., 2007; Takahashi et al., 1993; Turi et al., 2014) and pH (Fröb et al., 2019; García-Ibáñez et al., 2016; Pérez et al., 2021; Takahashi et al., 1993; Curbelo-Hernández et al., 2024a).

Seasonal changes in each driver (SST, SSS, C_T, and A_T) in Eq. (7) $\left(\frac{\mathrm{d}X}{\mathrm{d}t}\right)$ were assumed as their differences per month elapsed between the times of the year in which fCO_2,sw was at its minimum and maximum (seasonal amplitudes). Seasonal amplitudes were calculated between monthly means (based on observations and computed data) for February and September (when the minimum and maximum fCO_2,sw values were observed). An error propagation based on standard deviations for February and September was performed to calculate the error of the seasonal change.

Due to the high relevance of the evaporation and precipitation processes in the Mediterranean Sea and in order to avoid the influence of freshwater fluxes, the most recent equation (Eq. 7) given by Pérez et al. (2021) with salinity-normalized C_T and A_T (NC_T and NA_T) was used. The normalization was performed to a reference salinity (SSS₀) of 37.4 (NC_T= SSS₀ × $C_{\mathrm{T}}/$ SSS; NA_T= ((A_T−728.3) $/$ SSS) × SSS₀+728.3, following Friis et al., 2003), which is the average SSS for the entire monitored area.

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}f{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{d}t}}& ={\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial \mathrm{SST}}}{\displaystyle \frac{\mathrm{dSST}}{\mathrm{d}t}}+{\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial \mathrm{SSS}}}{\displaystyle \frac{\mathrm{dSSS}}{\mathrm{d}t}}\\ & +{\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial C_{\mathrm{T}}}}{\displaystyle \frac{\mathrm{d}{C}_{\mathrm{T}}}{\mathrm{d}t}}+{\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial A_{\mathrm{T}}}}{\displaystyle \frac{\mathrm{d}{A}_{\mathrm{T}}}{\mathrm{d}t}}\end{array}\text{(7)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}f{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{d}t}}& ={\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial \mathrm{SST}}}{\displaystyle \frac{\mathrm{dSST}}{\mathrm{d}t}}+({\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial \mathrm{SSS}}}+{\displaystyle \frac{{\mathrm{NC}}_{\mathrm{T}}}{{\mathrm{SSS}}_{\mathrm{0}}}}{\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial C_{\mathrm{T}}}}\\ & +{\displaystyle \frac{{\mathrm{NA}}_{\mathrm{T}}}{{\mathrm{SSS}}_{\mathrm{0}}}}{\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial A_{\mathrm{T}}}}){\displaystyle \frac{\mathrm{dSSS}}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{SSS}}{{\mathrm{SSS}}_{\mathrm{0}}}}{\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial C_{\mathrm{T}}}}{\displaystyle \frac{{\mathrm{dNC}}_{\mathrm{T}}}{\mathrm{d}t}}\\ & +{\displaystyle \frac{\mathrm{SSS}}{{\mathrm{SSS}}_{\mathrm{0}}}}{\displaystyle \frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial A_{\mathrm{T}}}}{\displaystyle \frac{{\mathrm{dNA}}_{\mathrm{T}}}{\mathrm{d}t}}\end{array}\end{array}$$

It is important to note that the changes in NA_T and NC_T are linked to biogeochemical processes which have different influences: the processes involved in the organic carbon pump contribute to strongly changing the NC_T that weakly affects the NA_T, while those involved in the carbonate pump affect the NA_T twice as much as NC_T. The positive values of $\frac{\mathrm{d}f{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{d}t}$ and $\frac{\partial f{\mathrm{CO}}_{\mathrm{2}}}{\partial X}\frac{\mathrm{d}X}{\mathrm{d}t}$ indicate an increase in fCO_2,sw from February to September, while the negative values do the opposite.

2.4.4 Air–sea CO₂ fluxes

The air–sea CO₂ fluxes (FCO₂) were determined using the bulk formula in Eq. (8):

$$\begin{array}{}\text{(8)}& {\mathrm{FCO}}_{\mathrm{2}}=\mathrm{0.24}{K}_{\mathrm{0}}{K}_{\mathrm{660}}\mathrm{\Delta }f{\mathrm{CO}}_{\mathrm{2}},\end{array}$$

where K₀ is the solubility of CO₂ in seawater, K₆₆₀ is the gas transfer velocity, and ΔfCO₂ represents the difference between fCO_2,sw and fCO_2,atm. A conversion factor of 0.24 was used to express FCO₂ values (mmol m⁻² d⁻¹). K₀ was calculated using the equation and coefficients given by Weiss (1974) and measured the SST and SSS which fall within the valid application limits. Considering the fitting error from the original parameterization of K₀ ($\pm \mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ mol L⁻¹ atm⁻¹; Weiss, 1974) and the instrumental errors of the SST and SSS measurements (Sect. 2.3), the uncertainty associated with the solubility estimation had a negligible impact on the calculation of FCO₂. K₆₆₀ was calculated through its quadratic dependency on wind speed (Eq. 9) using the parameterization given by Wanninkhof (2014):

$$\begin{array}{}\text{(9)}& K_{\mathrm{660}}=\mathrm{0.251}\cdot w^{\mathrm{2}}\cdot {\left({\displaystyle \frac{\mathit{Sc}}{\mathrm{660}}}\right)}^{-\mathrm{0.5}},\end{array}$$

where w is the wind speed and Sc is the Schmidt number (kinematic viscosity of seawater, divided by the gas diffusion coefficient). ERA5 hourly wind speed reanalysis data at 10 m above sea level and with a spatial resolution of 0.25° × 0.25° (Hersbach et al., 2023) were used to calculate K₆₆₀. ERA5 for the global climate and weather is available at the Copernicus Climate Data Store (https://cds.climate.copernicus.eu/; last access: 15 May 2025). The uncertainty in K₆₆₀ reported by Wanninkhof (2014) when using wind speeds ranging between 3 and 15 m s⁻¹ is ± 20 %. The error in the determination of fCO_2,sw and fCO_2,atm (Sect. 2.4.1) propagates into the calculation of ΔfCO₂ and constitutes an additional source of uncertainty. The statistical procedure used to quantify the uncertainty in FCO₂ arising from the uncertainty in ΔfCO₂ is described in Appendix B. The mean absolute error in FCO₂ due to the propagated uncertainty of ΔfCO₂ (± 2.01 µatm) was ± 0.14 mmol m⁻² d⁻¹, which in relative terms is ± 0.05 %. Negative FCO₂ values indicate that the ocean acts as an atmospheric CO₂ sink, while the positive values indicate that it acts as a source.

3 Results

A total of 157 984 data for surface ocean xCO₂ were collected during the study period (34 015 data during 2019, 28 590 data during 2020, 33 288 data during 2021, 19 102 data during 2022, 39 738 data during 2023, and 3251 data during January and February 2024). This exceeds the total number of data points available in the historical record for the western Mediterranean (34.8–43.1° N, 5.5° W–4.7° E) since 1999 (146 094 data) that are available in SOCAT v2024 (Bakker et al., 2016, 2024). The total number of data points in this region included in the SOCAT v2024 database since 2019 is 44 520.

Due to differences in the spatial distribution of the observations, two subregions (referred to as sections) were identified along the vessel track (Fig. 1): the longitudinally distributed southern section (hereinafter S section), encompassing the Alboran Sea (∼ 2–5.1° W), and the latitudinally distributed eastern section (hereinafter E section), following the eastern coastline of the Iberian Peninsula (∼ 36.5–41.3° N). The spatiotemporal distribution of fCO_2,sw and the total number of data points available in each dataset for sections S and E are shown in Fig. S2. In the S section, fCO_2,sw values from the ES-SOOP-CanOA station are consistent with those in SOCAT v2024, although the limited number of cruises covering this section in SOCAT v2024 makes a direct comparison difficult and prevents robust characterization of spatial and seasonal variability patterns. In the E section, some differences between the two datasets are observed (i.e. during spring–summer 2021, fCO_2,sw was higher in SOCAT v2024 than in the ES-SOOP-CanOA dataset). These differences are mainly explained by the distinct sampling trajectories in SOCAT v2024, with some routes extending further eastward, including coastal areas around the Balearic Islands.

The spatial distribution of the average values allowed us to identify heterogeneity in the annual cycle of each variable along both sections (Figs. 2 and S3). The standard deviation of the spatially averaged variables is presented in Table S2. A strong west–east increasing gradient in SST was observed in summer through the S section (∼ 5.5 °C), which led to an increment in fCO_2,sw of ∼ 57.5 µatm and a depletion in pH of ∼ 0.040 units eastward across the Alboran Sea. Despite the approximately constant SST through the S section during the rest of the year (less than 1.5 °C difference between the westernmost and easternmost parts), an eastward decrease in fCO_2,sw of less than 18 µatm accompanied by an increase in pH of less than 0.030 units was observed between October and March.

The latitudinal gradient of SST through the E section was weaker throughout the year, keeping the fCO_2,sw and pH spatially stable. The maximum change in SST occurs during winter, in which a northward decrease of less than 2 °C explained the minimum seasonal average temperatures and fCO_2,sw through the cruise track (14–15 °C and 350–360 µatm, respectively). This contrasts with the maximum average temperatures and fCO_2,sw encountered during summer (25.0–26.5 °C and 450–470 µatm, respectively). These results indicated that the maximum amplitude of the seasonal cycle of SST, fCO_2,sw, and pH occurs along the eastern coastline of the Iberian Peninsula and especially over the continental shelf between Valencia and Barcelona (northernmost part of the E section), while the minimum seasonal amplitude occurs near the Strait of Gibraltar (westernmost part of the S section).

The spatial variations in C_T were significant throughout the year along both sections (Fig. 2). C_T increases eastward in the order 20–45 µmol kg⁻¹ along the S section throughout the year. This increment accelerated along the E section from the Cape of Gata to the Cape of Nao and became approximately stable from the Cape of Nao to the port of Barcelona. The spatial distribution of C_T was greatly influenced by the progressive salinification observed along the S section. The SSS increased during the entire annual cycle from 36.3–36.5 around the eastern part of the Strait of Gibraltar to 37.7–38.1 around the Cape of Nao (Fig. S3). Removing the effect of salinity, NC_T (Fig. S3) presents a weaker spatial variation through the vessel track that is mainly led by biological and mixing processes.

Figure 2 Spatial distribution of the average SST, fCO_2,sw, pH, and C_T calculated on seasonal and annual bases every 0.1° in longitude along the S section (left panels) and every 0.25° in latitude along the E section (right panels). The 3-month periods January–March, April–June, July–September, and October–December were considered to be winter, spring, summer, and autumn, respectively. Note the different scales used for C_T due to significant variations in the S and E sections. The standard deviations are provided in Table S1 and indicate the range of variability in the study period.

The surface physicochemical properties show heterogeneities during some seasons of the year in several key locations along the sections (Figs. 2 and S3). The heterogeneities in the temporal evolution of the SST, SSS, and CO₂ system variables were assessed by the strategic selection of five stations along the S section (stations S1–S5) and six stations along the E section (stations E1–E6), which are depicted geographically in Fig. 1. S1 (4.95 ± 0.05° W) occupied the easternmost part of the Strait of Gibraltar, S2–S4 (4.35 ± 0.05° W, 3.85 ± 0.05° W, and 2.95 ± 0.05° W) were located in the central Alboran Sea, and S5 (2.45 ± 0.05° W) was located south of the Cape of Gata. The stations along the E section include E1 (37.1 ± 0.2° N) in the Gulf of Mazarron, E2 (37.6 ± 0.2° N) to the east of the Cape of Palos, E3 (38.2 ± 0.2° N) in the Gulf of Alicante, E4 (38.7 ± 0.2° N) to the east of the Cape of Nao, E5 (39.3 ± 0.2° N) in the Gulf of Valencia over the continental slope, and E6 (40.2 ± 0.2° N) near the Ebro estuary over the continental shelf.

The temporal variations of each variable at S1–S5 and E1–E6 are depicted in Figs. 3, 4, S4, S5, and S6. The seasonal amplitudes and interannual trends are summarized in Table 1. The seasonal amplitude of SST (minimum values in February–March around 14–17 °C and maximum values in August–September around 20–26 °C) increased eastward through the S section, with however a local decrease at S2 (Figs. 3 and S4, Table 1). The seasonal changes were larger through the E section (∼ 14 to ∼ 28 °C) and show weaker spatial variations (Figs. 4 and S5, Table 1). The SSS (Fig. S6) does not exhibit a seasonal cycle that is correlated well with the harmonic function in Eq. (A1) (r² < 0.5; Table S2). The lower and spatially more stable SSS values were observed along the S section during the entire period (around 36.0–37.5), with an increase with latitude through the E section (around 36.7–38.1).

Table 1 Seasonal amplitudes and interannual trends of SST, SSS, fCO_2,sw, pH, pH₁₉, C_T, and NC_T. The seasonal changes were calculated as the amplitude of Eq. (A1) fitted to the weekly average data at each station. The error of the seasonal amplitudes was assumed to be the result of the standard error of the estimate given by the harmonic function multiplied by 2. The interannual changes were based on linear regressions and are given for each station and for the entire S and E sections (considering the total amount of average data at S1–S5 and E1–E6, respectively) during the cold and warm seasons. The interannual trends of SST during 2000–2019 (based on monthly reanalysis data from the Med MFC physical multiyear product (Escudier et al., 2020; 2021; Nigam et al., 2021) detailed in Sect. 4.2) were included for comparison. The trends were obtained from the linear regressions of the seasonally detrended weekly average data and include their standard error of the estimate. ${}^{***}$ denotes that the trends are statistically significant at the 99 % level of confidence, ${}^{**}$ at the 95 % level of confidence, and ^* at the 90 % level of confidence.

The seasonal amplitudes of fCO_2,sw (from ∼ 340 to ∼ 460 µatm in the S section and from ∼ 340 to ∼ 470 µatm in the E section) and pH (from ∼ 8.00 to ∼ 8.12 units in the S section and from ∼ 8.00 to ∼ 7.98 to ∼ 8.13 units in the E section) were strongly linked to those of SST. They exhibit a west–east increment through the S section, with the exception of S2 (Figs. 3 and S4, Table 1), and remain approximately constant through the E section (Figs. 4 and S5, Table 1). These spatial heterogeneities in the seasonal cycles were found to be caused by the different rise in SST during late summer along each section, as minimal spatial differences were observed during the rest of the year.

C_T (Fig. S6) decreased seasonally from January–February to September–October (from ∼ 2180 to ∼ 2085 µmol kg⁻¹ in the S section and from ∼ 2260 to ∼ 2105 µmol kg⁻¹ in the E section), in phase with the enhancement of biological production. The seasonal amplitude of C_T increased eastward through the S section and northward through the E section, following the salinification gradient (Fig. S6, Table 1). Once the effect of salinity is removed, the seasonal cycle of NC_T shows minimal differences in the S section between the western and easternmost parts, while in the E section the NC_T and its seasonal amplitude increased northward (Fig. S6, Table 1). The enhanced adjustment (correlation) of NC_T with Eq. (A1) (0.47 < r² < 0.61 in the S section and 0.70 < r² < 0.88 in the E section) compared to C_T (0.28 < r² < 0.56 in the S section and 0.45 < r² < 0.73 in the E section) emphasizes the relevance of the processes varying the salinity. The lower correlations encountered through the S section show the higher impact of occasional processes (i.e. changes in the evaporation, precipitation, river runoff, and mesoscale features), locally modifying the surface carbon system in this area and introducing spatial heterogeneities into their seasonal cycles.

4 Discussion

4.1 Spatial characterization of the CO₂ system and its seasonality

4.1.1 The Alboran Sea

The seasonal variability of the AJ (García Lafuente et al., 2002; Macías et al., 2008, 2016; Vargas-Yáez et al., 2002) modified the SST signature in the S section, thus influencing the fCO_2,sw and pH. The maximum intensity of the AJ during summer (Peliz et al., 2013; Renault et al., 2012) caused a more intense warming and salinification of MAW while advancing into the Mediterranean Sea and mixing with the fraction of MW which surrounds the Cape of Gata and recirculates westward (Millot, 1999; Sánchez-Garrido et al., 2013). It explained the eastward increase in fCO_2,sw and the decrease in pH at this time of the year (Fig. 2; Sect. 3.1).

The relatively low SST and fCO_2,sw around S1 (20.68 ± 2.20 °C and 401.68 ± 27.13 µatm) and S3 (21.15 ± 2.11 °C and 407.30 ± 26.20 µatm) were mainly due to the highest intensity of the wind-induced upwelling along the northern coast of the western Alboran Sea during the warm season. This cooled the surface and enhanced the biological drawdown (e.g. Bolado-Penagos et al., 2020; Folkard et al., 1997; Gómez-Jakobsen et al., 2019; Peliz et al., 2009; Richez and Kergomard, 1990; Stanichny et al., 2005) while favouring the formation of the cold and nutrient-rich filament separating the WAG and EAG (Gómez-Jakobsen et al., 2019; Millot, 1999). Differences in the influence and strength of this filament may have contributed to the observed heterogeneities in SST, fCO_2,sw, and pH at S1 during the warm season (Fig. 3), which in turn accounts for the reduction in the model's fitting performance. Additionally, the shallowest position of the AMI during late winter (de la Paz et al., 2009; Echevarría et al., 2002; Gómez-Jakobsen et al., 2019; Minas et al., 1991) fed the surface with CO₂-rich waters from deeper areas of the Mediterranean Sea (de la Paz et al., 2009; Echevarría et al., 2002; Gómez-Jakobsen et al., 2019; Minas et al., 1991), elevating fCO_2,sw around S1. The increase in C_T and NC_T during summer around S3 (Figs. 2 and S3), which contributed to reducing their seasonal amplitudes in this area (Fig. S6, Table 1), suggests that the upwelled waters transported by the filament at this time of the year were not remineralized enough to compensate for the SST-driven decrease in fCO_2,sw. As a consequence, the western and eastern edges of the WAG presented the shortest seasonal amplitudes along the S section for SST, fCO_2,sw, and pH (Fig. 3, Table 1).

Figure 3 Time series of SST, fCO_2,sw, and pH at S1, S3, and S5 along the eastern Iberian margin within the 5 years of observations. The weekly average data were fitted to the harmonic Eq. (A1). The thermal and non-thermal terms of the average fCO_2,sw calculated by following the procedures of Takahashi et al. (2002) (T'02), Fassbender et al. (2022) (F'22), and pH₁₉ are depicted. The coefficients a–f, standard errors of the estimate, and r² as given by Eq. (A1) are presented in Table S1.

Conversely, the increase in SST and fCO_2,sw during summer around S2 (22.63 ± 2.05 °C and 429.98 ± 24.86 µatm), S4 (23.89 ± 2.03 °C and 438.25 ± 25.22 µatm), and S5 (24.05 ± 1.61 °C and 441.67 ± 16.22 µatm) contributed to extending their seasonal amplitudes in these zones (Figs. 2 and 3; Table 1). This suggests that, during the warm season, the increase in fCO_2,sw led by the surface warming near the core of the gyres was not compensated by the biological drawdown occurring at this time of the year (which caused weak decrements in C_T and NC_T at S2; Figs. 2 and S3).

Figure 4 Time series of SST, fCO_2,sw, and pH at E1, E4, and E5 in the Alboran Sea within the 5 years of observations. The weekly average data were fitted to the harmonic Eq. (A1). The thermal and non-thermal terms of the average fCO_2,sw calculated by following the procedures of Takahashi et al. (2002) (T'02), Fassbender et al. (2022) (F'22), and pH₁₉ are depicted. The coefficients a–f, standard errors of the estimate, and r² as given by Eq. (A1) are presented in Table S1.

4.1.2 The eastern Iberian margin

The eastern coastal transitional area of the Iberian Peninsula was subject to variability related to changes in the intensity, morphology, and path of the Northern Current (Fig. 1b). The SST decreased in the northernmost part of the E section from Sagunto to Barcelona throughout the year (north of S5; Fig. 2). The cooling of this area intensified during the cold season due to the mixing of warm waters in the wind-shielded area north of the Cape of Nao, with cool and salty MW transported by the Northern Current. However, it weakened during the warm season due to the northward spreading of MAW favoured by the formation of the thermal front in the axis of the Pyrenees, changing the path of the Northern Current (López-García et al., 1994). In the southernmost part of the section, the enhanced northward spreading of MAW and less wind stress during summer drive the warming observed from the Cape of Gata (at S5) to the Cape of Nao (at E4), while a low-intensity branch of the Northern Current transporting MW and progressing southward of the Cape of Nao weakly cool the area during winter (López-García et al., 1994; López-Jurado et al., 1995).

The local decrease in SST and fCO_2,sw observed during the warm seasons at E4 traced the offshore recirculation of the Northern Current at the Cape of Nao (Millot, 1999) and separated the E section into its northernmost and southernmost areas. This division was also evident based on the C_T and NC_T signatures (Figs. 2 and S6): the northernmost part of the section received remineralized MW transported by the Northern Current, which elevates C_T and NC_T, while the southernmost part was supplied with recent MAW with relatively low C_T and NC_T. Additionally, Ulses et al. (2023) recently suggested that the convective area in the Gulf of Lion acts as a source of natural and anthropogenic carbon for the intermediate waters of the western Mediterranean, which can enter the surface through vertical mixing and account for the observed high amount of C_T and NC_T.

Although the spatial heterogeneities and the northward cooling during the cold season (Fig. 1) increased the seasonal changes in SST, the seasonal amplitudes of fCO_2,sw were approximately constant within E1–E6 (Figs. 3 and S5; Table 1). The location of station E5, away from the influence of the Northern Current during the warm months, explained its locally lower seasonal amplitudes compared to the adjacent waters in the northernmost part of the E section. Nevertheless, these heterogeneities were minimal and did not cause differences in the seasonal amplitude of pH (Table 1).

In the cases of C_T and NC_T (Fig. S6, Table 1), the enhancement in the mixing of MAW with MW during winter increased the seasonality northward from E1 to E4. In the northernmost part, the seasonal variations in C_T and NC_T became shorter due to their increment during the cold season. This was caused by the combined action of the enhanced arrival of remineralized MW at this time of the year and the mesoscale structures locally favouring injections of CO₂-rich deeper waters into the surface (Bosse et al., 2021; Millot, 1999). The Ebro River runoff peaks in late winter and spring (Zambrano-Bigiarini et al., 2010) can also act as a source of variability around E5–E6.

4.2 Warming and interannual trends of MCS variables

The monitoring of the surface western Mediterranean basin allowed the identification of interannual trends for physical and MCS properties (Tables 1 and 2). The SST increased at rates of 0.38 ± 0.05 °C yr⁻¹ in the S section and 0.30 ± 0.04 °C yr⁻¹ in the E section. The rate of increase in SST that intensified locally at S2 (0.50 ± 0.09 °C yr⁻¹) may be due to the transport and accumulation of surface waters towards the core of the WAG. Its variability, migration, and progressive collapse can also account for the rapid warming of the area (Sánchez-Garrido et al., 2013; Viúdez et al., 1998; Vélez-Belchí et al., 2005).

The SST trends based on ES-SOOP-CanOA data were of the same order of magnitude as those derived from reanalysis data for the period 2019–2024 but were 1 order of magnitude higher than the reanalysis-based trends for 2000–2019, indicating a reinforcement of sea surface warming by approximately 80 %–90 % (Table 1). The ES-SOOP-CanOA data-based interannual SST trends were found to be reinforced during summer by 55.2 % in the S section and by 32.4 % in the E section compared to winter. The Northern Current cooling the northernmost part of the E section decelerated the warming in comparison to the S section. The ES-SOOP-CanOA data-based trends reported cumulative increases in SST from 2019 to 2024 of 1.91 ± 0.26 °C in the Alboran Sea (S section) and 1.52 ± 0.22 °C along the eastern Iberian margin (E section). These cumulative increments were 48.3 % and 34.94 % higher than those estimated for the global surface ocean from 1850–1900 to 2001–2020 (0.99 ± 0.12 °C; IPCC, 2023). This aligns with projections from climate models for both terrestrial and marine environments at the middle latitudes, particularly within the Mediterranean region, which is a consequence of human-induced global warming detailed by Hoegh-Guldberg et al. (2018) in the AR6 Synthesis Report (IPCC, 2023).

The warming contributes to modifying the MCS dynamics, mainly accelerating the increase in fCO_2,sw and the acidification. The interannual trends of fCO_2,sw and pH (Table 1) were more than twice (except for the trends at S1) those reported for the north-western Mediterranean at the DYFAMED site based on the difference between the average observation-based data for the periods 1995–1997 and 2013-2015 (2.30 ± 0.23 µatm yr⁻¹ and −0.0022 ± 0.0002 units per year; Merlivat et al., 2018) and for the north-eastern Atlantic at the ESTOC site based on in situ measurements since 1995 (2.1 ± 0.1 µatm yr⁻¹ and 0.002 ± 0.0001 units per year, respectively; González-Dávila and Santana-Casiano, 2023). The interannual rates accelerated eastward along the S section and northward along the E section (Table 1). The stronger trends at S3 compared to adjacent waters (S2 and S4) may be due to the transport of CO₂-rich waters from the southern Iberian coast through the filament. The trends in the S section were conducted by the larger rates of change encountered during the warm season compared to the cold season. The opposite occurred in the E section, where an intense increase in fCO_2,sw accompanied by a drawdown in pH occurred during winter, and the trends were reversed during summer (Table 1).

These spatial differences between the cold and warm seasons were mainly linked to variations in the biological production, remineralization, and mixing and were independent of the surface ocean warming. Hence, they were required to be assessed together with the NC_T trends for better understanding. The NC_T decreases interannually throughout the region (Table 2). The rapid depletion in the S section during winter in comparison to summer could be due to, first, an interannual weakening in remineralization processes and/or inputs of CO₂-rich water to the area during the cold months and, second, an interannual strengthening of the biological uptake during the warm months. However, these variations were insufficient to compensate for the increase in fCO_2,sw and subsequent decrease in pH induced by warming during the cold months, and even more so during the warm months. Conversely, in the E section, the variations in lateral and vertical advection as well as primary driven variations in the (sub-)mesoscale structures (Alberola et al., 1995; Bosse et al., 2021; 2016; Bourg and Molcard, 2021) were of high relevance and introduced differences into the annual cycle of NC_T. The interannual variations during winter were minimal (Table 1, Fig. S6), likely due to non-significant changes in remineralization and the dissolved CO₂ concentration of waters transported into the area. The decrease in NC_T intensified during summer (Table 1, Fig. S6), which was likely caused by the enhancement in biological production together with the weakening lateral advection (this may be related to a reinforcement in the front formed in the axis of the Pyrenees due to the increasingly higher SST of the MAW).

Table 2 Means, seasonal amplitudes, and interannual rates of change of the thermal and non-thermal components of fCO_2,sw (fCO_2,T and fCO_2,NT, respectively) calculated by following Takahashi et al. (2002) and Fassbender et al. (2022) (T'02 and F'22, respectively). The seasonal changes were calculated as the amplitude of Eq. (A1) fitted to the weekly average data at each station. The error of the seasonal amplitudes was assumed to be twice the standard error of the estimate given by the harmonic function. The trends were obtained by the linear regressions of the seasonally detrended weekly average data and include their standard error of the estimate. ${}^{***}$ indicates that the trends are statistically significant at the 99 % level of confidence. ${}^{**}$ is at the 95 % level of confidence and ^* at the 90 % level of confidence.

Once the effects of temperature were removed, the interannual pH₁₉ trends became negligible and were not statistically significant in the S section (< −0.001 units per year and p values > 0.1). This suggests that warming is directly driving the acidification (and indirectly the rising fCO_2,sw), while the progressively enhanced biological productivity partially compensates for the expected fall in pH. In the E section, pH₁₉ was reduced by 63 % (−0.002 ± 0.001 units per year and p values < 0.01) in comparison to the pH trends, which explains why the increase in SST is contributing more than half of the acidification due to the atmospheric fCO₂ increase only. The negative pH₁₉ trends were reinforced in the E section by 47 % during the cold season due to the enhancement in remineralization. The pH₁₉ trends reverted to being positive during the warm season due to the important role of biological production in actively decreasing fCO_2,sw and increasing pH at this time of the year.

However, despite the high statistical confidence in the trends and the consistency found with the reanalysis products, the acceleration in the surface warming and the consequent observed changes in fCO_2,sw and pH may be linked to isolated extreme events such as marine heat waves and are not necessarily indicative of prolonged behaviours over time. The globally increased frequency and magnitude in marine heat waves in phase with warming (Oliver et al., 2018; Hoegh-Guldberg et al., 2018; Frölicher et al., 2018; Smale et al., 2019) could provide feedback and hence continue to expedite the surface ocean warming. The influence of these extreme events is especially relevant in semi-enclosed seas as the Mediterranean, recognized as one of the most affected marine areas in the yearly Copernicus Ocean State Reports (OSR; EU Copernicus Marine Service; https://marine.copernicus.eu/access-data/ocean-state-report; last access: 15 May 2025) since 2016 (OSR1–OSR7).

4.3 The relative contributions of thermal and non-thermal processes to the surface fCO_2,sw

The temporal evolution of fCO_2,sw due to thermal and non-thermal effects (fCO_2,T and fCO_2,NT, respectively) showed a high degree of agreement between the T'02 and F'22 methodologies (Figs. 3 and 4). The average fCO_2,T and fCO_2,NT values differed by less than 5 µatm between the two methodologies. The consistency with the widely employed T'02 engenders confidence in the validity and reliability of the most updated F'22 method.

The seasonal variations in fCO_2,sw were close to 2 times as high in the E section as in the S section (Table 1). The thermally driven seasonal changes (dfCO_2,T) were found to be approximately double those independent of temperature (dfCO_2,NT) throughout the region (Table 2). The $T/B$ ratios demonstrated the control of thermal processes on the seasonality of fCO_2,sw throughout the region (Table 2). The $T/B$ ratios in the westernmost part of the S section (ranging between 1 and 2) were consistent with previous studies on the Strait of Gibraltar (Curbelo-Hernández et al., 2021b; de la Paz et al., 2009). The $T/B$ ratios increased eastward as the AJ advanced in the Alboran Sea and were caused by the intense increase in dfCO_2,T compared to dfCO_2,NT. They exceeded 2 in S4–S5 and E1–E6, which demonstrated the larger control of SST on fCO_2,sw in areas less influenced by the input of surface Atlantic Water.

The interannual trends show the control of thermal processes on the increase in fCO_2,sw during 2019–2024 (Figs. 3 and 4; Table 2). The strong and statistically significant interannual fCO_2,T trends show the important role of warming in elevating fCO_2,sw. The weak and non-significant fCO_2,NT trends suggest that spatiotemporal variations in the biological processes, circulation patterns, and air–sea gas exchange introduced local differences into the distribution of fCO_2,sw. It is difficult to assess the impact of the non-thermal processes on an interannual scale at each of the stations. The interannual trends of fCO_2,T and fCO_2,NT for the entire S and E sections (Table 2) were statistically significant at more than the 95 % level of confidence, and their coupling described, with a difference of less than 0.3 µatm yr⁻¹ (< 1 %), the interannual rates of fCO_2,sw during 2019–2024 (Table 1; Sect. 4.2).

The thermal processes govern the changes in fCO_2,sw on an interannual scale, with contributions ranging between ∼ 76 %–92 % in the S section and ∼ 73 %–83 % in the E section. The contributions for fCO_2,NT were between ∼ 8 %–25 % and ∼ 17 %–27 %, respectively. The decrease in fCO_2,NT compensated by ∼ 6 %–30 % the increase in fCO_2,sw at S1–S5 and E1–E2, while its increase contributed ∼ 24 %–53 % to raising fCO_2,sw at E3–E6. The negative fCO_2,NT trends in the S section were related to progressive enhancement in the biological uptake (mainly during spring/summer) not compensated by remineralization and/or vertical–lateral advection of remineralized waters (mainly during autumn and winter) in areas influenced by recent MAW. Conversely, the interannual increase in fCO_2,NT in the E section suggests that the supply of cool and remineralized MW along the path of the high-intensity Northern Current surpasses the biological drawdown of surface CO₂ and accelerates the increase in fCO_2,sw on an interannual scale.

4.4 Mechanism controlling the seasonal cycle of fCO_2,sw

To infer the causes of variations in the seasonal cycle of fCO_2,sw in the study period, the seasonal rates of change in fCO_2,sw ($\frac{\mathrm{d}f{\mathrm{CO}}_{\mathrm{2},\mathrm{sw}}}{\mathrm{d}t}$, hereinafter dfCO₂) were decomposed into their individual components ($\frac{\partial f{\mathrm{CO}}_{\mathrm{2},\mathrm{sw}}}{\partial X}\frac{\partial X}{\mathrm{d}t}$, hereinafter dfCO${}_{\mathrm{2}}^X$) as described in Sect. 2.4.3 (Eqs. 6 and 7). The results of the solved Eq. (7) for each year at S1–S5 and E1–E6 are depicted in Fig. 5. The uncertainty associated with the difference between the monthly means for each term and year was obtained through error propagation considering their individual standard errors and is presented in Table S3. The dfCO₂ resulted from the cumulative sum of the individual terms in Eq. (7) (indicated with the subscript “sum”) matching the dfCO₂ calculated directly from observations between both seasons (indicated with the subscript “obs”), which instills confidence in the methodology (Fig. 5).

Figure 5 Temporal evolution of the seasonal rates of fCO_2,sw explained by each of their drivers within the 5 years of observations. The differences between the monthly average data for February and September (where the minimum and maximum SST and fCO_2,sw were encountered) were considered to compute the seasonal trends. The standard deviation of the monthly average data was considered in the calculation of the seasonal changes and inferred errors in the computation of fCO_2,sw that are summarized in Table S3. The cumulative fCO_2,sw changes ($\frac{\mathrm{d}f{\mathrm{CO}}_{\mathrm{2},\mathrm{sw}}}{\mathrm{d}t}\left(\mathrm{sum}\right))$ resulting from the distinct drivers were consistent with the observed seasonal fCO_2,sw trends ($\frac{\mathrm{d}f{\mathrm{CO}}_{\mathrm{2},\mathrm{sw}}}{\mathrm{d}t}\left(\mathrm{obs}\right)$), thereby instilling confidence in the methodology.

The SST was identified as the main driver of dfCO₂ throughout the stations, accounting on average for 51 %–71 % of its values over the study period. At some stations and in specific years, this contribution occasionally decreased to ∼ 45 % or increased to ∼ 83 %. In the S section (Fig. 5a), dfCO${}_{\mathrm{2}}^{\mathrm{SST}}$ increased westward as MAW was warmed in the Alboran Sea, while the incursion of the filament cooled the surface locally and decreased dfCO${}_{\mathrm{2}}^{\mathrm{SST}}$ at S3. In the E section (Fig. 5b), dfCO${}_{\mathrm{2}}^{\mathrm{SST}}$ increased northward and reached its maximum north of the Cape of Nao (at E4–E6), particularly during 2021–2022 (32.0–32.5 µatm month⁻¹), due to the greater influence of warmed MW.

A_T described on average < 18 % of the dfCO₂ in the entire region, occasionally increasing up to 22 %. As fCO_2,sw changes inversely with A_T, the weakly negative dfCO${}_{\mathrm{2}}^{\mathrm{AT}}$ found for some years along the S section shows fluctuations in the periods of increment and decrement for A_T that are likely related to changes in the mixing processes. The A_T contribution becomes negligible at E6 (< 2 % throughout the study period) due to the minimal seasonal amplitude of A_T and NA_T (Fig. S6). The approximately constant A_T and NA_T levels throughout the year may be due to the bicarbonate and carbonate content from the Ebro River runoff being neutralized by that in MW and MAW, which spread into the area during winter and summer, respectively. dfCO${}_{\mathrm{2}}^{\mathrm{AT}}$ has tended to decrease since 2020–2021 in S1–S3, S5, and E1 due to the progressive weakening of the NA_T depletion from February to September. The opposite occurred north of the Cape of Palos, where the seasonal cycle of NA_T reaches its maximum amplitude (20–28 µmol kg⁻¹ at E3 and E4). The interannual dealkalinization in the S and E sections (Table 1) acts as a source of heterogeneities: the interannual negative NA_T trends during the cold months (p values < 0.01) were stronger than during the warm months (p values > 0.1) and were consistent in both sections. The spatial differences in the summer trends (weaker in the S section compared to the E section) account for an enhanced reduction in the seasonal amplitude of NA_T in the S section.

The dfCO${}_{\mathrm{2}}^{\mathrm{SSS}}$ values were minimal in both the S and E sections (< 0.6 and < 1.9 µatm per month, respectively) and show the weak impact of SSS on dfCO₂ (< 3.5 %). The entrance of MAW and its mixing with the saltier MW in the Alboran Sea do not allow us to identify a seasonal pattern in SSS (Fig. S6), thus explaining the negligible contribution of SSS in the S section (∼ 2.0  % at S1, which fell to < 0.7 % at S2–S5). The larger seasonal amplitudes of SSS at E2–E5 (Fig. S6) led to a relatively major influence of SSS (∼ 1.0 %–2.3 % during most of the years). The low seasonal amplitude of SSS and A_T at E6, likely related to an approximately constant influence of the Northern Current at this location throughout the annual cycle, caused a negligible variation in dfCO₂ (< 0.4 %).

The depletion in C_T, mainly driven by the increased biological production from February to September, had a significant inverse impact on dfCO₂ (23 %–37 %). In 2019, at stations S1–S3 (and in 2020 only at S1), the inverse contribution of C_T reached 39 %–47 %. This suggests that the influence of the C_T cycle on dfCO₂ in the westernmost Alboran Sea increasingly resembles that observed in the rest of the Mediterranean. These findings indicate that the seasonal drawdown of C_T offsets between one-fourth and one-half of the expected increase in dfCO₂ driven by SST and is slightly prompted by A_T. In the S section (Fig. 5a), the lower increase observed from 2019 to 2023 in dfCO${}_{\mathrm{2}}^{\mathrm{CT}}$ (4–6 µatm per month) compared with dfCO${}_{\mathrm{2}}^{\mathrm{SST}}$ (6–9 µatm per month) demonstrated that fluctuations in C_T were increasingly insufficient for counterbalancing the warming-driven increase in dfCO₂, even at S2–S4, where the biological production was enhanced and the dfCO${}_{\mathrm{2}}^{\mathrm{CT}}$ hence reinforced after 2020. In the westernmost part of the S section, the influence of C_T offsetting dfCO₂ was at its maximum during 2019–2020 at S1 (> 84 %), S2 (67.3 %), and S3 (86.1 %) and diminished towards 2023 (37.1 %, 38.3 %, and 45.1 %, respectively). In the easternmost part, this compensation was around 33 %–44 % at S4–S5 throughout the period (as at S2 and S3 since 2020), except for 2023 at S5, in which dfCO${}_{\mathrm{2}}^{\mathrm{CT}}$ weakened and offset only 22.8 %. In the E section (Fig. 5b), the progressive strength in the processes depleting C_T throughout the period at E1–E4 and since 2020 at E5–E6 compensated the dfCO${}_{\mathrm{2}}^{\mathrm{SST}}$ by 33 %–46 %, which changes inversely to dfCO${}_{\mathrm{2}}^{CT}$. The lowest compensations found in 2019 at E5 (28.8 %) and E6 (18.4 %) were likely related to possible injections of remineralized waters along the Northern Current path that offset the biological uptake of C_T and the elevated dfCO${}_{\mathrm{2}}^{\mathrm{CT}}$.

4.5 Air–sea CO₂ exchange across the western boundary of the Mediterranean Sea

The continuous observation of MCS variables enabled the calculation of FCO₂ at an unprecedentedly high spatiotemporal resolution in the western Mediterranean Sea. The FCO₂ was found to be governed by fluctuations in ΔfCO₂ (Fig. 6) that were mainly controlled by the broader variability of fCO_2,sw (325–500 µatm) compared to fCO_2,atm (390–425 µatm). The SST fluctuations play a relevant role by primarily controlling fCO_2,sw (Sect. 4.3) and modulating the solubility of CO₂ at the air–sea interface. The entire monitored area was undersaturated for CO₂ with respect to the low atmosphere between late October and June (ΔfCO${}_{\mathrm{2}}=-\mathrm{35.30}$ ± 8.97 µatm), acting as an atmospheric CO₂ sink (−2.56 ± 0.55 mmol m⁻² d⁻¹) which peaks in winter (−4.53 ± 0.44 and −3.29 ± 0.31 mmol m⁻² d⁻¹ in the S and E sections, respectively). During summer, the area was supersaturated for CO₂ (ΔfCO₂= 36.43 ± 0.35 µatm) and acted as a source that was about 3 times more intense along the E section (1.70 ± 0.43 mmol m⁻² d⁻¹) compared to the S section (0.57 ± 0.35 mmol m⁻² d⁻¹).

Figure 6 Temporal variations of FCO₂ (blue; left axis), ΔfCO₂ (orange; right axis), and wind speed (gray; left axis) at (a) S1–S5 and (b) E1–E6. A piece-wise polynomial-based smoothing spline was applied to the weekly average data (represented by dots). Gaps were covered by the harmonic fitting (Eq. A1; dashed line). The black lines represent the interannual increase in FCO₂. The seasonally detrended interannual rates of change in FCO₂ and ΔfCO₂ are shown in each panel. ${}^{***}$ indicates that the trends are statistically significant at the 99 % level of confidence, ${}^{**}$ at the 95 % level of confidence, and ^* at the 90 % level of confidence. The wind speed does not show statistically significant interannual trends (p values > 0.1).

The spatial differences in SST during the warm months introduced heterogeneities into the seasonal outgassing in both sections: the higher SST during summer in the E section reduced the solubility and contributed to a greater increase in fCO_2,sw with respect to fCO_2,atm (ΔfCO₂=49.83 ± 0.32 µatm) compared to the cooler S section (ΔfCO₂=16.35 ± 0.14 µatm). The seasonality in the formation of the CO₂ sink and source in the Alboran Sea was consistent with previous studies of the Strait of Gibraltar (Curbelo-Hernández et al., 2021b; de la Paz et al., 2011, 2009) and the north-western African coastal transitional area in the north-eastern Atlantic (Curbelo-Hernández et al., 2021a; Padin et al., 2010) and agreed with the seasonal pattern that is characteristic of tropical and subtropical regions (Bates et al., 2014; Takahashi et al., 2002). The warming during summer at S1 was insufficient to lead to supersaturated conditions (ΔfCO${}_{\mathrm{2}}=-\mathrm{5.56}$ ± 0.26 µatm) and thus acted as a CO₂ sink throughout the year (−2.83 ± 1.77 mmol m⁻² d⁻¹ during the cold months and −0.52 ± 0.02 mmol m⁻² d⁻¹ during the warm months), which coincided with the behaviour observed in the Strait of Gibraltar during 2019 (Curbelo-Hernández et al., 2021b). The sink and source status during the cold and warm months encountered in the eastern Iberian margin agrees with FCO₂ evaluations based on observations in the Mediterranean basin through its north-western (Wimart-Rousseau et al., 2023, 2021, 2020) and eastern parts (Sisma-Ventura et al., 2017) and confirms previous estimations based on satellite data and models (D'Ortenzio et al., 2008; Taillandier et al., 2012).

The variations in FCO₂ during the period of study were addressed by averaging the data across seasons and years at each of the selected stations (Fig. 7). The same procedure was applied to ΔfCO₂ and wind speed (Figs. S7 and S8). The evolution of the seasonal in-gassing and outgassing was evaluated by computing interannual trends for the average FCO₂ and ΔfCO₂ (Fig. 7). The interannual FCO₂ trends showed the progressive strength of the summer source in the S section, which was accelerated at S2 in response to the enhanced warming around the WAG (detailed in Sect. 4.2) and at S4–E1 due to their exposure to increasing wind forcing (Figs. S7 and S8). This was caused by the increase in fCO_2,sw during the warm months that was not offset by biological drawdown, which elevated ΔfCO₂. By contrast, the localization of E2–E6 over the eastern Iberian continental shelf and slope allowed the relevant biological uptake at this time of the year to compensate for the influx of CO₂-rich water. It introduced heterogeneities in ΔfCO₂ between years which do not allow us to identify statistically significant trends.

Figure 7 Temporal evolution of the average FCO₂ calculated on seasonal and annual bases for each year (2019–2023) at S1–S5 and E1–E6. The same representation for ΔfCO₂ and wind speed is available in Figs. S5 and S6. The 3-month periods January–March, April–June, July–September, and October–December were considered to be winter, spring, summer, and autumn, respectively. The legend includes the interannual trends for FCO₂ (mol m⁻² yr⁻¹ yr⁻¹) based on linear regression of the seasonal and annual means. ${}^{***}$ denotes that the trends are statistically significant at the 99 % level of confidence, ${}^{**}$ at the 95 % level of confidence, and ^* at the 90 % level of confidence. The standard deviations are presented in Table S4.

During spring and autumn, the increase in ΔfCO₂, mainly driven by warming and accompanied by the decreasing wind stress (Figs. S7 and S8), led to the positive interannual FCO₂ trends at S2–S5 and E1–E6 (Fig. 7). They show the weakening of the in-gassing during autumn and the achievement of a near-equilibrium state with the atmosphere during spring by the end of the study period. The FCO₂ reverted to being weakly positive during the spring of 2023 in the E section, which prolonged the seasonal source period's relevant impact on the net annual FCO₂. During winter, the increasing wind forcing compensated for the reduction in the in-gassing expected by the rise in ΔfCO₂ (Figs. S7 and S8). However, the variability in the wind speed and other processes involved in the non-thermal change in fCO_2,sw between the years does not allow identification of statistically significant rates of change in the CO₂ sink status. In particular, the relatively high wind speed during the winter of 2021 may have contributed to accelerated horizontal transports, increasing fCO_2,sw and hence ΔfCO₂ (Figs. S7 and S8).

The predominantly negative FCO₂ trends during most of the year led to net annual CO₂ sink behaviour. The positive FCO₂ trends during summer, spring, and autumn have forced the annual average CO₂ invasion to decrease by 44 %–65 % at S2–S5 (ranging from −0.66 ± 0.06 and −0.84 ± 0.04 mol m⁻² during 2019 to −0.27 ± 0.09 and −0.47 ± 0.09 mol m⁻² during 2023) and by 60 %–80 % at E1–E6 (ranging from −0.32 ± 0.09 and −0.53 ± 0.09 mol m⁻² during 2019 to −0.11 ± 0.10 and −0.13 ± 0.09 mol m⁻² during 2023). The unique hydrodynamics of the Strait of Gibraltar strongly influenced the air–sea CO₂ exchange at S1: the in-gassing during summer partially compensated for the reduction in the annual influx and resulted in a lower increase in FCO₂ (23 %) from 2019 (−0.77 ± 0.02 mol m⁻² yr⁻¹) to 2023 (−0.60 ± 0.06 mol m⁻² yr⁻¹).

Considering the annual average FCO₂ for the S and E sections, the net in-gassing has decreased at rates of 0.11 ± 0.02 mol m⁻² yr⁻¹ yr⁻¹ (p value < 0.01) in the Alboran Sea and 0.08 ± 0.02 mol m⁻² yr⁻¹ yr⁻¹ (p value < 0.01) in the eastern Iberian margin. This contrasts with the strength of the CO₂ sink across the western Mediterranean basin recently reported by Zarghamipour et al. (2024) for 1984–2019 based on a combination of observational data and model simulations (0.007 ± 0.001 mol m⁻² yr⁻¹ yr⁻¹). Additionally, Zarghamipour et al. (2024) noted the reduction in the annual net CO₂ source behaviour of the central Mediterranean basin at an estimated rate of 0.003 ± 0.001 mol m⁻² yr⁻¹ yr⁻¹. The findings suggest that the accelerated increase in fCO_2,sw induced by the rapid warming, together with the progressive reduction in solubility, is reversing the interannual FCO₂ trends compared to previous decades, which may be causing the study area to resemble the central and eastern Mediterranean basins in terms of air–sea CO₂ exchange. The reduction in the net annual invasion was consistent with previous estimations in such coastal and shelf environments across the eastern tropical and subtropical South Atlantic during 2002–2018 (between 0.03 ± 0.01 and 0.09 ± 0.02 mol m⁻² yr⁻¹ yr⁻¹; Ford et al., 2022) and toward mid-latitudes over the Scotian Shelf (with the average FCO₂ ranging from −1.7 mol m⁻² yr⁻¹ yr⁻¹ in 2002 to −0.02 mol m⁻² yr⁻¹ yr⁻¹ in 2006; Sisma-Ventura et al., 2017). The continuation of this decreasing rate for net annual in-gassing would imply the reversion of the study area to net annual CO₂ source behaviour before 2030.

The net CO₂ invasion was calculated by integrating the annual cycle of FCO₂ during 2019–2023. The net FCO₂ in the Alboran Sea was −1.57 ± 0.49 mol m⁻² yr⁻¹, which represented a strength in the CO₂ sink in comparison with the adjacent surface areas across the Strait of Gibraltar (between −0.82 and −1.01 mol m⁻² yr⁻¹ during 2019–2021; Curbelo-Hernández et al., 2021b) and the eastern Iberian upwelling (−1.33 mol m⁻² yr⁻¹; Chen et al., 2013). The net FCO₂ along the eastern Iberian margin was −0.70 ± 0.54 mol m⁻² yr⁻¹, which falls within the range of that modelled for the deep-convection area around the Bay of Marseille (north-western Mediterranean basin) during 2012–2013 (−0.5 mol m⁻² yr⁻¹; Ulses et al., 2023) and estimated based on observations during 2017–2018 (between −0.26 and −0.81 mol m⁻² yr⁻¹; Wimart-Rousseau et al., 2020). However, this was contrary to the net outgassing across the eastern Mediterranean basin (0.85 ± 0.27 mol m⁻² yr⁻¹ during 2009–2015; Sisma-Ventura et al., 2017). The net CO₂ sinks for the monitored area across the Alboran Sea (14 000 km²) and the eastern Iberian margin (40 000 km²) were −0.97 ± 0.30 Tg CO₂ yr⁻¹ (−0.26 ± 0.08 Tg C yr⁻¹) and −1.22 ± 0.95 Tg CO₂ yr⁻¹ (−0.33 ± 0.25 Tg C yr⁻¹). These findings greatly contributed to the assessment of the air–sea CO₂ exchange in the Mediterranean basin and the global coastal and shelf areas.

5 Conclusions

The 5 years of automatically underway observations at the ES-SOOP-CanOA ocean station provided a high-spatiotemporal-resolution dataset which includes the surface physical and MCS properties across the western margin of the Mediterranean Sea. It allowed characterization, with an improved degree of certainty, of mechanisms involved in the MCS dynamics in the Alboran Sea and eastern Iberian coastal transitional area on seasonal and interannual timescales.

The variations in fCO_2,sw were found to be strongly controlled by temperature fluctuations. On a seasonal scale, the thermally driven variations intensified as the AJ advanced eastward in the Alboran Sea and MAW was formed, moving northward along the eastern Iberian margin and mixing with MW. In the Alboran Sea, the high intensity of the AJ during summer warms the surface layer towards the core of the WAG and EAG, driving larger seasonal changes in SST, fCO_2,sw, and pH that increased during the study period. The eastern Iberian margin was meridionally separated at the Cape of Nao by the path of the Northern Current: the northernmost part, fed with cool, salty, and remineralized MW during the cold season and influenced by the northward spreading of MAW during the warm season, shows the largest seasonal amplitudes for SST, fCO_2,sw, and pH compared to the southernmost part, which is supplied with recent MAW during most of the year and by a weak and relatively warm branch of the Northern Current during winter. Driver analysis has shown that 51 %–71 % of the increase in fCO_2,sw from February to September within the entire monitored area was explained by SST and < 20 % by A_T and SSS, while the processes controlling C_T offset this increment, contributing 23 %–37 %.

The changes in the seasonal cycles were driven primarily by the increasing contribution of temperature (due to the larger seasonal amplitude of SST), and secondarily by the decreasing contribution of C_T (due to the decreasing remineralization/production ratio). On an interannual scale, the SST increased at rates ranging between 0.26 and 0.43 °C yr⁻¹ and drove a rapid increase in fCO_2,sw within 4.18 and 5.53 µatm yr⁻¹ and a decrease in pH within −0.0049 and −0.0065 units per year. The ∼ 76 %–92 % of the interannual increase in fCO_2,sw was described by warming. In the Alboran Sea, and extending northwards to the Cape of Palos, non-thermal processes, which are primarily biological drawdowns during spring blooms, compensated for up to one-third of the expected increase in fCO_2,sw due to warming. The opposite occurred north of the Cape of Palos, where non-thermal processes, mainly the inflow of CO₂-rich MW during the cold season, accounted for the increase in fCO_2,sw.

The assessment of the air–sea CO₂ exchange shows the western boundary of the Mediterranean basin to be undersaturated and acting as a significant sink of atmospheric CO₂ during most of the year, while the presented supersaturated conditions led to a CO₂ source status during the warm months. The entire monitored area acted as a net annual CO₂ sink that is weakening at statistically significant rates ranging between 0.06 and 0.13 mol m⁻² yr⁻¹ yr⁻¹ (40 %–80 % since 2019). These trends would cause the area to shift towards becoming a net annual CO₂ source before 2030 if the current climate conditions persist. The weakening in the net annual CO₂ sink was driven by the ongoing strength of the summer outgassing (mainly in the Alboran Sea) and the weakening in the autumn and spring in-gassing (throughout the region). Integrating the annual cycle of FCO₂ during the entire study period, the net CO₂ in-gassing values calculated for the Alboran Sea and eastern Iberian margin were −1.57 ± 0.49 and −0.70 ± 0.54 mol m⁻² yr⁻¹.

This study highlights the need for systematic observation strategies to characterize the physicochemical properties of seawater in the Mediterranean, an effort that has been required by the scientific community for the last few decades. It demonstrates the effectiveness of SOOP/VOS in monitoring surface physical and biogeochemical variables, especially in highly variable and anthropogenically pressured areas such as coastal and semi-enclosed seas. The findings enhance our understanding of MSC dynamics in a key coastal transitional area of the western Mediterranean, which has high environmental and socioeconomic importance and implications for the regional climate. Likewise, they contribute to a more accurate understanding of the role of coastal areas in the context of global change at both the basin and global scales. Despite the relatively short study period, this research captured shifts likely driven by isolated events caused by climate change, offering insights into future ocean conditions.

Appendix A: Data adjustments and statistical procedures

The temporal evolution of the physicochemical data was analysed by weekly averaging (the time required by the vessel to complete a trip) at different locations along the vessel track. The average values (y) were fitted to Eq. (A1) as a function of time (yearly fraction). This equation updates the one used by Curbelo-Hernández et al. (2021a, b) to study seasonal cycles through the addition of the b(year-2019) term, which provides the interannual rate of change of each seasonally detrended variable between 2019 and 2024. The coefficients a–f and the standard errors of the estimate given by Eq. (A1) for the variables considered are available in Table S1.

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}y& =a+b\left(\text{year-}\mathrm{2019}\right)+c\cdot \cos \left(\mathrm{2}\mathit{\pi }\mathrm{year}\right)\\ & +d\cdot \sin \left(\mathrm{2}\mathit{\pi }\mathrm{year}\right)+e\cdot \cos \left(\mathrm{4}\mathit{\pi }\mathrm{year}\right)\\ & +f\cdot \sin \left(\mathrm{4}\mathit{\pi }\mathrm{year}\right)\end{array}\end{array}$$

The errors in the weekly averages were determined by dividing the standard deviation by the square root of the number of data points used to calculate the means $(\mathrm{Standard}\mathrm{Deviation}/\sqrt{n})$. The coefficient b in Eq. (A1) represented the interannual variation rates for each variable, which coincided with the slope derived from linear regressions of the detrended average values over time. The standard errors of these slopes were calculated by propagating the errors from the annual mean values.

The strength and direction of the linear regressions and the significance of the interannual trends were evaluated through a Pearson correlation test. This test yielded correlation coefficients (r²) and corresponding p values to determine statistical significance. Trends with p values ≤ 0.01 were statistically significant at the 99 % confidence level, those with p values ≤ 0.05 were significant at the 95 % confidence level, and those with p values ≤ 0.1 were significant at the 90 % confidence level. Trends with p values > 0.1 were not statistically significant but still provided an estimate of the temporal evolution of the variables within their respective layers.

Appendix B: Uncertainty in FCO₂ explained by the propagated error in ΔfCO₂

The uncertainty in ΔfCO₂ was calculated by applying standard error propagation rules for the difference of two independent measurements with associated uncertainties (Eq. B1):

$$\begin{array}{}\text{(B1)}& {\mathit{\sigma }}_{\mathrm{\Delta }f{\mathrm{CO}}_{\mathrm{2}}}=\sqrt{{\mathit{\sigma }}_{f{\mathrm{CO}}_{\mathrm{2},\mathrm{sw}}}^{\mathrm{2}}+{\mathit{\sigma }}_{f{\mathrm{CO}}_{\mathrm{2},\mathrm{atm}}}^{\mathrm{2}}},\end{array}$$

where ${\mathit{\sigma }}_{f{\mathrm{CO}}_{\mathrm{2},\mathrm{sw}}}$ and ${\mathit{\sigma }}_{f{\mathrm{CO}}_{\mathrm{2},\mathrm{sw}}}$ are the uncertainties for fCO_2,sw and fCO_2,atm, respectively (see Sect. 2.4.1). The absolute error in FCO₂ (${\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}}}$; mmol m⁻² d⁻¹) associated solely with uncertainty in ΔfCO₂ was estimated for each data point using Eq. (B2):

$$\begin{array}{}\text{(B2)}& {\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}}}=K_{\mathrm{660}}{K}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{\Delta }f{\mathrm{CO}}_{\mathrm{2}}}.\end{array}$$

To represent the average magnitude of uncertainty in the estimated FCO₂ over the entire dataset (with n being the total number of data), the mean absolute FCO₂ error was calculated using Eq. (B3) and the mean relative FCO₂ was estimated with Eq. (B4):

$$\begin{array}{}\text{(B3)}& {\displaystyle }\stackrel{\mathrm{‾}}{{\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}}}}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n{\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}},i},\text{(B4)}& {\displaystyle }\overline{{\displaystyle \frac{{\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}}}}{{\mathrm{FCO}}_{\mathrm{2}}}}}={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n\left|{\displaystyle \frac{{\mathit{\sigma }}_{{\mathrm{FCO}}_{\mathrm{2}},i}}{{\mathrm{FCO}}_{\mathrm{2}}}}\right|\cdot \mathrm{100}.\end{array}$$