This study focuses on three PERLE cruises: PERLE0, PERLE1, and PERLE2 (Figure 1). These cruises were carried out in the EMed between 2018 and 2019. At all stations, a CTD-Rosette was deployed (1) to acquire data with sensors (Conductivity Temperature and Depth–CTD and associated parameters) along vertical profiles and (2) to collect discrete seawater samples from Niskin bottles for chemical analysis. Over the 11, 31, and 125 casts performed during the PERLE0, PERLE1, and PERLE2 cruises, seawater was sampled from 1, 12, and 17 casts, respectively, for carbonate parameter analysis (see Supplementary Table 1 and Supplementary Figure 1). Details for the cruises and parameters measured during each PERLE cruise are summarised in Table 1.

Systematic primary quality control of the measured data was performed on each PERLE dataset. No significant problems have been detected for Winkler oxygen and pH measurements. During PERLE1, for a few casts, a CTD pump dysfunction significantly altered the quality of the CTD oxygen: oxygen measurements from these casts were disregarded. A systematic quality control procedure for A_T and C_T was conducted based on internal consistency tests between A_T, C_T and pH_T (see details in the Supplementary Material). Following these steps, only 15 PERLE2 casts were validated, leading to the loss of ca. 60% of the PERLE2 A_T/C_T dataset. All the A_T/C_T PERLE1 dataset was lost. A comparison of the quality controlled PERLE dataset with previously collected data does not reveal systematic biases for A_T, C_T, or pH_T²⁵ (Figure 2A–C).

Relationships between years and carbonate parameters (A_T, C_T, and pH_T²⁵) and between A_T and salinity were computed using a linear regression model. Linear regression statistics, including the standard error of the slope (i.e., the error of the estimated trend), the coefficient of determination (r²) and the significance of the trend (p-value) were calculated using the R software. Linear relationships have been tested using the Pearson coefficient for parametric test (Sokal and Rohlf, 1969) with a significance level of 95%.

Parameters derived from the A_T-S linear relationship were tested against previously published A_T-S relationships in the area using a Student’s t-test for the slope and intercept. The null hypothesis, H₀, was that our observations were not significantly different from these linear models.

Absolute salinity (S_A), conservative temperature (Θ) and potential density (σ_θ) were derived from practical salinity, temperature and pressure and the geographic position based on the TEOS-10 (The International Thermodynamic Equation of Seawater-2010). In this study, following the recommendations of the Intergovernmental Oceanographic Commission (Valladares et al., 2011), S_A and Θ were used to study the hydrological context (Θ−S_A diagrams). Calculations were made with the “oce” R package (Kelley et al., 2017). Note that practical salinity (labelled Salinity) and in situ temperature (labelled Temperature) were used in this study to facilitate comparisons with previous studies in particular, for A_T-S relationships.

Apparent Oxygen Utilisation (AOU–μmol.kg^–1) was calculated from the difference between oxygen solubility concentration (at P = 0 dbar) estimated with the “Benson and Krause coefficients” (Garcia and Gordon, 1992) and in situ [O₂]_mes.

A density threshold of 0.03 kg.m^–3 with a reference depth of 10 dbars was used to compute the Mixed Layer Depth (MLD) (D’Ortenzio et al., 2005).

Salinity data were used to reconstruct an A_T time-series using the sub-surface A_T-S relationship proposed by Hassoun et al. (2015a) (see discussion in section “Total Alkalinity and Salinity Relationships Within the Mixed Layer”). In this study, the PERLE1 and the Argo float A_T datasets were reconstructed following this A_T-S relationship. Considering the standard deviation of the A_T-S relationship proposed by Hassoun et al. (2015a), the accuracy of the calculated A_T values is ± 19 μmol.kg^–1.

Salinity-normalised changes in A_T (NA_T^39.3) and C_T (NC_T^39.3) were calculated dividing by in situ salinity and multiplying by 39.3 (i.e., the mean PERLE salinity above 200 dbars).

Seawater carbonate system parameters were derived from A_T and C_T values. Calculations were made with the software program CO2SYS-MATLAB (van Heuven et al., 2011) using silicate and phosphate concentrations. When nutrient data was not available, silicate and phosphate mean concentrations for each depth were used. As recommended for the MedSea by Álvarez et al. (2014), the carbonic acid dissociation constants K₁ and K₂ from Mehrbach et al. (1973) as refitted by Dickson and Millero (1987) and the dissociation constant for HSO₄^– from Dickson (1990) were used. Uppström (1974) was used to calculate the ratio of total boron to salinity and Dickson and Riley (1979) to calculate the hydrogen fluoride constant K_F.

The buffer factors γA_T (γC_T), βA_T (βC_T) and ωA_T (ωC_T) provide an estimation of the seawater’s ability to buffer changes in the aqueous CO₂ [CO₂], protons [H⁺] and the carbonate saturation state (Ω) when A_T (C_T) changes at constant C_T (A_T) (Egleston et al., 2010). The calculations were performed following the formula proposed by Álvarez et al. (2014).

Net Ecosystem Production (NEP) is defined as the sum of biotic and abiotic carbon fluxes in the ecosystem (Borges et al., 2008). Net Ecosystem Calcification (NEC) is a measure of the balance between CaCO₃ formation (calcification) and dissolution (Smith and Kinsey, 1978). Based on the NA_T^39.3 and NC_T^39.3 plot, the reaction path can take on variable slopes depending on the ratio of different processes, such as photosynthesis/respiration, carbonate dissolution/formation and CO₂ release/invasion (Zeebe, 2012). Temporal changes in NA_T^39.3 (ΔNA_T^39.3) and NC_T^39.3 (ΔNC_T^39.3) between each PERLE cruise can be calculated according to NEP and NEC processes as:

Following equation (2), NEP can be expressed according to NEC as:

Then, by replacing the NEP term in equation (1) by equation (3), NEC can be calculated as:

NEC and NEP are expressed in μmolC.kg^–1.d^–1. Salinity-normalised A_T and C_T values “exclude” the “precipitation-evaporation” influence in the layer where biological activity is at a maximum. It is assumed that the layers considered (MLD-200 dbars) to estimate the NEP and NEC processes are not influenced by air-sea CO₂ fluxes, which were therefore not considered.

CARIMED (CARbon, tracer and ancillary data In the MEDsea) aims to be an internally consistent database containing inorganic carbon data relevant for this basin (Álvarez et al., in preparation). Ancillary (hydrographic, inorganic nutrients and dissolved oxygen), CO₂ (pH, A_T, and C_T) and transient tracer (CFC-11 and 12, Tritium, SF₆, Neon, CCl₄, and ΔHe³) data from several cruises in the MedSea from 1976 until 2018 were assembled. Primary and secondary quality control procedures following the GLODAP (Global Ocean Data Analysis Project) philosophy (Tanhua et al., 2010) are locally adapted to this marginal sea. This work only uses data collected in the Levantine basin (Supplementary Table 2).

All vertical profiles for A_T, C_T and pH_T²⁵ measured during the PERLE cruises are presented in Figures 2A–C, respectively. All the A_T profiles presented maximum values in the surface, minimum values between 500 and 700 dbars and remained almost constant (or slightly decreasing) below 1000 dbars. Most of the C_T vertical profiles presented the lowest values in surface waters, reaching maximum values between 500 and 700 dbars and then remaining relatively invariable below 1000 dbars. pH_T²⁵ presented maximum values at the surface (with values around 8.060 measured during PERLE1 cruise), minimum values close to 700 dbars and nearly constant values under 1000 dbars (Figure 2C). The main water masses are identified in Figures 2D–F and detailed in Supplementary Figure 1.

Intermediate waters (mostly LIW) were located around the 29.0 kg.m^–3 isopycnal layer (Lascaratos and Nittis, 1998; see Supplementary Figure 1) and were characterised by an A_T maximum evolving from 2,600 to 2,640 μmol.kg^–1 (Figures 2A,D). As observed by Álvarez et al. (2014), the LIW was located above the layer of maximum organic matter mineralisation in the EMed and was associated with low C_T concentrations (ca. 2,290 μmol.kg^–1) and high pH_T²⁵ values (ca. 8.000) in contrast to the deepest water masses. It can be observed that slightly colder, more haline and denser Cretan Intermediate Waters (Velaoras et al., 2019) were detected during PERLE2 in the western part of the Cretan Sea with the highest A_T value for PERLE2 cruise (ca. 2,660 μmol.kg^–1, Figure 2A).

In the deep-water layer (i.e., EMDW), both AeDW and AdDW presented similar C_T values (Figure 2E) while slightly higher pH_T²⁵ (Figure 2F) and A_T (Figure 2D) values were measured in the AeDW (see Supplementary Figure 2). On the Cretan shelf, deep waters were comprised of dense EMDW with high A_T (≈ 2,650 μmol.kg^–1) and C_T values (≈ 2,350 μmol.kg^–1). Deep waters of the Cretan Sea were filled with CDW with low pH_T²⁵ (≈ 7.950) values resulting from relatively low A_T and high C_T content (Figures 2D–F).

This description of the carbonate chemistry in the deep and intermediate water masses in the PERLE area is in good agreement with previous studies (Schneider et al., 2010; Álvarez et al., 2014). However, the PERLE strategy based on an intense observation period over a year is not appropriate to describe changes in deep-water masses. For the rest of this study, in order to tackle the seasonal dynamics of the surface waters, only data in the NWLB (Figure 1) where all three PERLE cruises were conducted, will be discussed further.

The highest spatial and temporal variability in carbonate chemistry parameters was encountered in the upper water layer which has been defined to be approximately the first 200 dbars. Discrete pH_T²⁵ values (measured and calculated), taken from the southern part of the PERLE sampling area (the NWLB) illustrate the seasonal variability of the carbonate chemistry in the upper layer (Figure 3A). The pH_T²⁵ was the most measured carbonate parameter in this study and, when normalised to 25°C, can be considered as an indicator of the carbonate chemistry status by including the changes in A_T and C_T. An overview of the upper layer seasonal dynamics is also presented for temperature, salinity, fluorescence, and AOU profiles in Figures 3B–E, respectively.

The lowest pH_T²⁵ values were encountered in March 2019 during the PERLE2 cruise and correspond to the relatively higher C_T values and lower A_T values. During this cruise, a significant range in the MLD was encountered with the deepest values observed. This cruise coincided with the abrupt stratification observed in the EMed after the deepening of the MLD from November to February-March (D’Ortenzio et al., 2005). Increased fluorescence values were observed in shallow waters at the end of the cruise (in the eastern part of the area) in comparison to the beginning of the cruise (in the western part).

Intermediate pH_T²⁵ values were measured in June 2018 during the PERLE0 cruise corresponding to increased surface alkalinity and a moderate depletion in inorganic carbon. The PERLE0 cruise is an early summer cruise characterised by a shallow MLD. The highest fluorescence values were recorded during this cruise well below the MLD (ca. 90 dbars) and light oxygen supersaturation (AOU ≈ −20 μmol.kg^–1) just beneath the MLD.

Finally, high pH_T²⁵ values (>8.000) were measured up to 100 dbars during the PERLE1 cruise, probably in association with a high A_T content due to evaporation. During this late summer cruise, the deepest Deep Chlorophyll Maximum (DCM) with the lowest fluorescence values but also the deepest negative AOU concentrations were encountered. Moreover, during this cruise, the mesoscale Ierapetra Eddy (IE) was crossed (see Supplementary Figure 3 and Ioannou et al., 2019). The core of this warm and salty eddy (Figures 3B,C) was characterised by a deepening of the MLD associated with a deep DCM and negative AOU values. Nonetheless, no clear IE signal was observed on the pH_T²⁵ values (Figure 3A).

In the EMed, spring and autumn seasons need to be considered as short transition periods between the summer and winter, which come later than on the continent (Özsoy et al., 1989). Moreover, in the EMed, summer is characterised by maximum heat in the surface layer that can remain up until November, whereas winter is identified with minimal heat that can occur until April. Considering each cruise as representative of a period within the annual cycle, the PERLE0 cruise (June 2018) associated with intermediate pH_T²⁵ values corresponds to the early summer period with decreasing biological activity associated with the strengthening of stratification. PERLE1 (October 2018) is associated with the highest pH_T²⁵ values and corresponds to the end of the summer period characterised by a warm and stratified water column with deep and low fluorescence maximum. PERLE2 (March 2019), associated with the lowest pH_T²⁵ values and shallow fluorescence maximum, corresponds to the end of the winter period, with the beginning of the seasonal stratification of the water column in the eastern part. These features agree with the analysis of the seasonal patterns of surface chlorophyll a concentration (Chl a) (based on remote sensing). The lowest values of surface Chl a were observed during the summer period, whereas an increase in surface Chl a was observed in winter, concomitantly to the deepening of MLD (Bosc et al., 2004; D’Ortenzio and Ribera d’Alcalà, 2009).

When no A_T values were available (see section “Primary Quality Control of the Measured Data”), A_T can be estimated based on an A_T-S relationship. In the MedSea, several linear relationships between A_T and salinity in the surface waters have been proposed for different sub-basins (e.g., Schneider et al., 2007; Cossarini et al., 2015; Hassoun et al., 2015a; Gonzaìlez-Daìvila et al., 2016).

During the PERLE cruises, in the NWLB, AT was significantly (n = 14, p-value = 0.014, r² = 0.36) influenced by salinity variations within the mixed layer (Figure 4). Figure 4 also displays the A_T-S distribution in the Cretan Sea (grey dots on Figure 4). The mixing of high alkalinity Black Sea waters (values of ca. 2,967 μmol.kg^–1; Hiscock and Millero, 2006) in the Cretan Sea shifts the A_T-S characteristics of surface waters in agreement with Schneider et al. (2007) who demonstrated that freshwater and Black Sea inputs affect the A_T-S relationship. More pronounced deviations from the expected linear A_T-S relationship are observed for stations with deeper mixed layers (Figure 4). This might be the result of the mixing of water masses with different A_T-S relationships during winter mixing. As A_T values were available only for PERLE0 and PERLE2 cruises, the A_T-S relationship derived for the PERLE cruises in the mixed layer (and in the NWLB) have been based on a very limited number of data. The PERLE A_T-S linear relationship was tested against the Hassoun A_T-S linear model (Hassoun et al., 2015a). No significant differences were found on either the slope (t-test = 1.86, n = 14, p < 0.05) or the intercept (t-test = 0.27, n = 14, p < 0.05). Therefore, the annual time-series were reconstructed based on the A_T-S linear relationship measured by Hassoun et al. (2015a) in the surface waters (0–25 m) of the eastern Mediterranean sub-basin, and A_T has been estimated based on this relationship.

During the PERLE cruises, the NWLB exhibited a greater range in A_T than C_T values within the mixed layer (see section “Total alkalinity control on the seasonal air-sea CO₂ exchanges”). A_T ranged between 2,610 and 2,693 μmol.kg^–1 whereas C_T ranged between 2,292 and 2,332 μmol.kg^–1. Over an annual scale, the ratio of the range in A_T variations to the range in C_T variations (ΔA_T/ΔC_T) can be used to infer the sensitivity to A_T and C_T changes in the upper ocean. Over the period studied, in the NWLB, the ratio ΔA_T/ΔC_T is equal to 2.1. In the global ocean, long-term time-series ΔA_T/ΔC_T ratios are lower than 1.0 (Table 2).

The reasons for these apparent and rather unique ranges of A_T and C_T over the year in the NWLB can be attributed to several factors: (1) The main drivers of the C_T gradient in the water column are, primary production transforming the C_T into organic carbon in the photic layer, and respiration transforming the organic carbon into C_T. As the EMed is an area of low productivity (Moutin and Raimbault, 2002), the vertical C_T gradient is lower than in other oceanic areas. Consequently, the C_T range in surface waters, driven by C_T consumption during the stratified period and replenishment via vertical mixing with sub-surface waters enriched in C_T, is greatly reduced. (2) The high levels of evaporation that affect the MAW in the EMed during the summer season increases salinity by nearly 1 g.kg^–1 (Figure 2) between the end of winter (PERLE2) and the end of summer (PERLE1). The A_T and C_T parameters should be equally affected by evaporation in a closed system. However, when reported on a A_T/C_T diagram (with normalised axes—see Figure 5), a higher range of A_T variation compared to C_T is observed. This indicates that when salinity increases in surface waters, a concomitant consumption of C_T must occur to compensate for the C_T increase due to evaporation to maintain an apparent stability in C_T concentrations. The biological consumption of C_T will be discussed in the next section as a possible mechanism to explain this low C_T variability.

To understand the overall impact of biological processes on the seasonal variations in the carbonate system in the NWLB, changes in A_T and C_T need to be considered independently from the changes induced by dilution and evaporation. For this purpose, salinity-normalised changes of A_T and C_T in the upper 200 dbars are plotted in Figure 5. To differentiate waters affected by air-sea exchanges from sub-surface waters, the upper 200 dbars of water column has been divided into two layers: within and below the mixed layer (0 dbars—MLD and MLD—200 dbars). The barycentre of all observational points, defined as the coordinate of the mean A_T and C_T values during each cruise, is reported and considered to be representative of the “season” sampled.

The barycentres are spread along the photosynthesis-respiration line between the three cruises, reflecting the effects of biological processes on the carbonate system over the year. From the early summer period (PERLE0—red dots on Figure 5) to the end of the summer period (PERLE1—blue dots on Figure 5), for both layers, the barycentre shift was a signature for increased photosynthetic processes compared to respiration processes. The deepening of the DCM observed between the PERLE0 and PERLE1 cruises and the negative AOU values recorded during these cruises supported this observation. The deepening of the DCM is a signature to the downward displacement of primary producers related to surface nutrient depletion (Sigman and Hain, 2012), and negative AOU values reflect oxygen production. All these elements indicate that autotrophic processes dominate the upper water column between early and late summer. Based on these assumptions, between the end of the summer period (PERLE1) and the end of the winter period (PERLE2—green dots on Figure 5), the barycentre shift indicates that heterotrophic processes were dominant in the upper water column. Whilst observations cannot be time related, it can be assumed that between the late winter period of PERLE2 and the early summer period of PERLE0, the “theoretical” shift of the barycentre indicates a balance in favor of autotrophic processes during this period. When considered together, these seasonal changes in normalised A_T and C_T confirm that during periods of high evaporation, autotrophic processes are consuming C_T and increasing A_T. This can explain the apparent C_T stability and the important change in A_T over an annual cycle.

Based on the assumption that, below the mixed layer, the PERLE sampling area is a closed system (unimpacted by air-sea CO₂ fluxes), the temporal evolution in NA_T^39.3 and NC_T^39.3 was used to calculate NEP and NEC fluxes. From the end of the bloom period (PERLE0) to the end of the summer period (PERLE1), daily NEP and NEC values of 0.53 and 0.01 μmolC.kg^–1.d^–1, respectively, were estimated whereas from the end of the summer period (PERLE1) to the start of the bloom period (PERLE2), negative daily NEP and NEC values of −1.02 and −0.04 μmolC.kg^–1.d^–1, respectively, were estimated. In the MedSea, the MLD seasonal variability is characterised by a deepening from November to February-March (D’Ortenzio et al., 2005). Therefore, it can be assumed that the water masses below the mixed layer remain isolated from surface CO₂ inputs between the PERLE0 and PERLE1 cruises. However, due to the late winter deepening of the MLD (Figure 3), between the end of the summer period (PERLE1) and the late winter period (PERLE2), NEC and NEP could be biased by air-sea exchanges.

The seasonal NEP values estimated in this study confirm previous estimations based on oxygen concentration changes monitored with short-time incubations during the stratified period. In June 2006, Regaudie-de-Gioux et al. (2009) reported a positive NEP value of 0.22 ± 1.30 mmol O_2.m^–3d^–1 in waters above 100 meters in the EMed and in summer 2008, Christaki et al. (2011) reported positive NEP values of 4 ± 14 mmol O_2.m^–2d^–1. As previously observed by Schneider et al. (2007), the contribution of calcification and dissolution processes to variations in the carbonate system could be assumed to have a minor role in the MedSea. The NEC values calculated in the NWLB confirm this. The spreading of PERLE2 data points along the CaCO₃ formation/dissolution line in Figure 5 (green dots) might be associated to the spatial changes in alkalinity content across the geographical distribution of sampling sites during this cruise rather than to calcification and dissolution processes.

To address the question of the control of A_T and C_T changes on the “source” (pCO₂^SW > pCO₂^ATM) or “sink” (pCO₂^SW < pCO₂^ATM) of CO₂ in the NWLB, PERLE’s A_T and C_T values are reported in Figure 6. The temperature range in the area has been used to draw the red and blue “iso pCO₂^SW-lines” as representative of the pCO₂^SW values encountered during the winter and summer PERLE cruises. Considering a mean atmospheric partial pressure (pCO₂^ATM) value of 403 μatm (recorded at Lampedusa site from October 2018 to December 2019; Dlugokencky et al., 2021), the upper seawaters encountered at the warm end of summer with high alkalinity (PERLE1) were a “source” of CO₂. In contrast, the cold and low alkalinity end of winter (PERLE2) surface waters were a “sink” of CO₂ with pCO₂^SW.

Although the C_T content remained almost stable between the PERLE cruises, the A_T variability was noticeable with the lowest A_T values measured at the end of the winter period (PERLE2) and the highest A_T values estimated during PERLE1, at the end of the summer period. When considering the large pCO₂^SW variations due to the temperature variability represented by the shift between the red and blue isolines, the high alkalinity seawater at the end of summer (PERLE1–blue dots on Figure 6) induces low pCO₂^SW values when seawater starts to cool and therefore highlights the potential for surface waters to absorb atmospheric CO₂. In the NWLB, the variability of the A_T content of the surface waters over an annual cycle impacts the air-sea CO₂ exchanges. The “classical” vision that the pCO₂^SW variability is not driven by temperature change but by the biological control on C_T, must be largely revisited in light of the important effect that variations in A_T have on the pCO₂^SW regulation capability in the EMed.

In order to estimate the effect of the A_T variability on the pCO₂^SW over an annual cycle, alkalinity was derived from salinity data from an Argo float that cycled in the NWLB for over a year. The temperature and total alkalinity (derived from salinity) values recorded by the float in the upper 20 dbars of the water column representative of the surface mixed layer affected by air-sea exchanges are presented in Figure 7. The cruise data within the mixed layer are also reported. In Figure 7, the red “iso pCO₂^SW-line” indicates the pCO₂ equilibrium between the ocean and the atmosphere. This isoline was derived at constant C_T, based on the assumption that the pCO₂^SW is, apart from temperature, controlled by A_T rather than by C_T in the NWLB. The distribution of data above and below this line highlights the “source” or “sink” status of the NWLB for atmospheric CO₂, respectively.

The float derived data agreed with data measured during the PERLE cruises and indicate a penetration of atmospheric CO₂ into the EMed from December to April, and a release of CO₂ into the atmosphere from May to November. It must be noted that these estimates are sensitive to the C_T value used. Indeed, by considering a high C_T content (grey isoline labelled “C_T max” in Figure 7), the period of CO₂ “sink” for the atmosphere will be shorter (from February to April). Conversely, if the lowest C_T mean value is considered (black isoline labelled “C_T min” in Figure 7), the area will act as a “sink” from December to May. The observed “iso pCO₂^SW-lines” shift (grey and black isolines in Figure 7) from the “iso pCO₂^SW-line” at mean C_T (red isoline in Figure 7) due to the C_T variability over a year induces a temporal change in the status of “source” or “sink” of the upper water masses. Moreover, by considering the accuracy of ± 19 μmol.kg^–1 associated to the A_T estimation (according to Hassoun et al., 2015a), the uncertainty of the estimated pCO₂^SW has been calculated (Orr et al., 2018) and ranged between the two “iso pCO₂^SW-lines” deduced from the maximum and minimum C_T values (red area on Figure 7). Although the displacement of the air-sea pCO₂ equilibrium might shift considering the A_T uncertainty, the temporal succession of the “sink” or “source” status for atmospheric CO₂ throughout a year in the NWLB is evidenced. It confirms that the A_T content of the surface waters is a significant driver of the air-sea CO₂ fluxes in the NWLB.

These are, to the best of our knowledge, the first estimates of the succession of the “sink” and “source” status in the NWLB based on in situ data. Previous estimates based on satellite observations of sea surface properties, and on a model characterising the evolution of the mixed layer pCO₂^SW (D’Ortenzio et al., 2008; Taillandier et al., 2012) are confirmed by this study. Moreover, coastal observations in the South eastern Levantine basin close to the Israeli shelf have also reported a CO₂ source for the atmosphere in summer (from May to December) and a sink of atmospheric CO₂ in winter (from January to April) (Sisma-Ventura et al., 2017).

Based on historical observations from the CARIMED dataset and observations from the PERLE cruises, temporal changes in carbonate chemistry between 2001 and 2019 in the surface NWLB have been assessed to study the mechanisms that could explain the carbonate system changes over the last twenty years (Figure 8). The surface layer has been defined to a depth of 50 dbars to include sufficient data. Due to the seasonal changes in surface salinity in the EMed (Grodsky et al., 2019), salinity-normalised A_T (NA_T^39.3) and C_T (NC_T^39.3) were used to facilitate the comparison between the different datasets across space and time. Indeed, due to the strong salinity dependency of alkalinity, by normalising by salinity, a significant part of the seasonal signal for alkalinity is removed.

While being higher (even when salinity-normalised) than the trends observed in the North Western MedSea (i.e., 1.40 ± 0.15 μmol.kg^–1.a^–1; Merlivat et al., 2018), the temporal C_T increase in the NWLB surface waters (Figure 8A) is consistent with other trends measured in the eastern Levantine basin (i.e., 5 ± 2 μmol.kg^–1.a^–1; Hassoun et al., 2019). However, when compared to other time-series over the global ocean, the trends measured in the surface NWLB waters are 3.7–1.5 times higher (if the NC_T^39.3 trend is considered) than the global ocean range which lies between 0.78 μmol.kg^–1.a^–1 (Munida South Pacific time-series) and 1.89 μmol.kg^–1.a^–1 (CARIOCA time-series; Bates et al., 2014). This suggests that distinct mechanisms explaining the increasing C_T trend exist in the NWLB.

While A_T is considered insensitive to atmospheric CO₂ penetration (Zeebe, 2012), positive trends in C_T and negative trends in pH_T²⁵ (Figures 8A,E) can be explained, at least partially, by the increase in atmospheric CO₂. Indeed, between 2006 and 2018, a mean annual increase of 2.2 ± 0.08 ppm.a^–1 in xCO₂^ATM (mole fraction of CO₂) was recorded at the Lampedusa site (equivalent to the trend recorded on a global scale; Dlugokencky et al., 2021). To estimate the sensitivity of the estimated trends to the increase in atmospheric CO₂, the increase in xCO₂^ATM was assumed to be equivalent to a surface ocean increase in pCO₂^SW. Based on the estimated trends in pCO₂^SW, NA_T^39.3, and NC_T^39.3, annual changes in carbonate chemistry pCO₂^SW have been calculated by solving thermodynamic equations (Table 3). The observed annual decrease in pH_T²⁵ (Figure 8E) and increase in C_T (Figure 8A) lies between the values estimated with and without an A_T increase. This suggests that an A_T increase must exist to compensate for the decrease in pH and the increase in C_T or, in other words, that the high observed C_T trend is the consequence of the observed A_T increase. Although a positive A_T trend has been observed elsewhere in a coastal site of the MedSea (Kapsenberg et al., 2017), it remains unexplained. These changes could be related to changes in riverine inputs or changes in Black Sea water inputs (Schneider et al., 2007).

It is worth noting that the CARIMED database, by merging data measured over the past 20 years, has a large over-representation of the spring season (Supplementary Figure 1 and Supplementary Table 2). Moreover, the spatial distribution of the sampled stations was different for each cruise. The scarcity of observations in the NWLB precludes the estimation of the seasonal variability on the observed trends. Due to the observed influence of seasonal conditions on the carbonate parameters during the PERLE cruises, time-series that would include observations of the peculiar conditions observed in the late summer (high surface pH_T²⁵ associated with high A_T values during PERLE1—Figure 2C) or winter could modulate the observed temporal trends. Nonetheless, when data collected during “not spring” cruises are not considered to estimate the trends, despite shifting the temporal trend values, tendencies remain significant for each parameter. Thus, the conclusion that a decadal A_T increase must exist to counterbalance the pH decrease associated to the C_T increase remains coherent and valid.

In the projected warmer MedSea (Nykjaer, 2009), increased stratification but also reduced nutrient inputs from river discharge caused by more frequent drought periods could increase the oligotrophy of the MedSea (e.g., Moon et al., 2016; Pagès et al., 2019, 2020). As this study suggests that the magnitude of the annual C_T variation in surface waters is reduced in the EMed due to the low C_T vertical gradients, all processes that could decrease primary production in the future could reduce the C_T contribution to the air-sea exchanges.

Even if internal thermohaline oscillation needs to be considered to draw solid conclusions about salinity trends, over the past 30 years, a positive long-term trend in salinity for the LSW and LIW has been recorded (Ozer et al., 2017). Because of the salinity impact on alkalinity concentrations (Figure 4) and of the A_T impact on the air-sea CO₂ fluxes (Figure 7), if the PERLE1 conditions are exacerbated in the future with marine heatwaves extending over longer periods of the year, even more alkaline waters can be expected at the end of the summer. An even greater potential pCO₂^ATM sink will result when surface seawaters cool. The gyres (such as the IE), which have a higher A_T content due to their saltier waters, might be even more efficient at catching atmospheric CO₂ when seawater cools. The control of air-sea CO₂ exchange by alkalinity that is suggested in this study could be enhanced in a future warmer and less productive EMed. However, as C_T and A_T are equally affected by evaporation and as, in the future less productive EMed, the C_T biological consumption will be less efficient, the mechanisms leading to stable inorganic carbon content described in this study might be altered.

In an attempt to quantify the sensitivity of the carbonate system to future C_T and A_T changes, estimated buffer factors within the MLD for each PERLE cruise are presented in Table 4. At a comparable period of the year (March–April for PERLE2 cruise), the estimated buffer factors are in good agreement with former estimates (Álvarez et al., 2014) whereas the estimated buffer factors for PERLE0 and PERLE1 cruises during summer are significantly higher. Higher absolute buffer values imply higher buffering capacity and lower changes in [CO₂], pH or Ω for a given change in A_T or C_T. Assuming that the PERLE1 conditions will be exacerbated in the future (Darmaraki et al., 2019), the EMed surface water is moving toward an overall increase in its buffering capacity (relative to changes in A_T and C_T).

It is worth noting that, when atmospheric CO₂ dissolves in seawater, the CO₂ concentration in solution changes due to the carbonate ion buffering effect. The future effects of the decadal trends measured in the NWLB on the buffering capacities of the carbonate ion can be discussed using three different perspectives: (1) By considering the observed decrease in pH_T²⁵, the carbonate ion availability will decrease accordingly, reducing the atmospheric CO₂ uptake by the MedSea. (2) The greater increase in C_T in comparison to the increase in A_T will reduce the carbonate ion availability, but, nevertheless, will compensate for the impact of a pH decrease on the carbonate ion content, so allowing the CO₂ uptake into the atmosphere. (3) The positive trend in A_T, and its impact on the CO₂ atmospheric uptake and on mitigating the decreasing pH trend, may indirectly increase the C_T.