The Earth’s climate is an open system, and the energy entering or leaving the Earth determines the Earth Energy Imbalance (EEI). The EEI is defined as the net radiative flux at the top of the atmosphere, i.e. as the difference between the incoming solar shortwave radiation and the outgoing infrared longwave radiation, and is considered the most fundamental metric for climate change. Anthropogenic greenhouse gas emissions caused a radiative imbalance at the top of the atmosphere that led to excessive heat within the Earth climate system (Trenberth and Stepaniak, 2004; von Schuckmann et al., 2016; Von Schuckmann et al., 2020; Trenberth et al., 2009; Hansen et al., 2011; Trenberth et al., 2014).

The World Ocean accounts for the uptake of 89% (90%) of the EEI over the period from 1971 to 2018 (2010-2018, respectively) in the upper 700 m (Levitus et al., 2012; Rhein et al., 2013; von Schuckmann et al., 2016; Von Schuckmann et al., 2020). The heat gain over land accounts for the uptake of 6% (5%) over these periods, while 4% (3%) go into the melting of grounded and floating ice and 1% (2%) go into atmospheric warming (Von Schuckmann et al., 2020). This quantification of the EEI – clearly shows that the ocean plays a major role in the heat uptake and storage. Therefore, the ocean heat content trend is a key component for the quantification of the EEI (Knutti and Rugenstein, 2015; Von Schuckmann et al., 2020; Cheng et al., 2023) and gives an estimate of the net climate forcing (Hansen et al., 2005; Levitus et al., 2005). It dominates the energy uptake with most of the warming absorbed by the upper layer (Rhein et al., 2013) and therefore acts as a strong buffer to climate change.

Due to the high heat capacity of seawater (Millero et al., 1973), small changes in seawater temperature lead to a large heat storage and therefore a strong increase in Ocean Heat Content (OHC). The global ocean warming for the period from 1955 to 2010 in the upper 700 m was found to be 0.18°C, corresponding to an increase in OHC of 0.27 Wm^-2 (per unit area of the World Ocean; Levitus et al., 2012). The OHC is also directly linked to sea level rise through thermal expansion of the warmer waters, causing a thermosteric sea level rise of 0.41 mmyr^-1 for the upper 700 m layer of the World Ocean from 1955 to 2010 (Levitus et al., 2012; Church and White, 2011).

The region considered for our study is the Mediterranean Sea (MS), which is defined as one of the climate hotspots and plays an important role for climate change prediction (Giorgi, 2006; Bernstein et al., 2008; Lionello and Scarascia, 2018; Menna et al., 2022). Due to its reduced dimension, semi-enclosed morphology, and anti-estuarine circulation (evaporation exceeds over precipitation and river runoff), it responds faster to climatic changes than the global ocean (Giorgi, 2006; Schroeder et al., 2017). The MS occupies about 0.32% of the total volume of the global ocean, yet its unique geological history has resulted in high biodiversity with 7-10% of all known marine species (Bianchi and Morri, 2000; Coll et al., 2010). During the twentieth century it has undergone rapid changes with accelerating trends in temperature and salinity (Béthoux et al., 1998; Ozer et al., 2017; Schroeder et al., 2017) and climate model projections show a further acceleration in the future (Giorgi, 2006; Giorgi and Lionello, 2008; Kirtman et al., 2013). The MS is showing an increase not only in temperature and salinity, but also in sea level and freshwater fluxes during the period 1993-2020 (Menna et al., 2022). Many species are struggling to adapt to the rapidly warming waters, and more frequent marine heat waves are causing mass mortality events (Garrabou et al., 2022; Frölicher et al., 2018; Oliver et al., 2021; Smith et al., 2022).

Since the late 1980s, many studies have been conducted to quantify the trends of temperature and OHC in the MS (Vargas-Yáñez et al., 2008; Vargas-Yáñez et al., 2017; Iona et al., 2018). The time series of the averaged heat content anomaly from 1993 to 2019 (0- 700 m depth layer) shows a yearly trend of 1.4 ± 0.3 Wm^-2 (OHC trend from Copernicus Marine Service; https://marine.copernicus.eu/access-data/ocean-monitoring-indicators). The yearly trend of temperature and sea level which was assessed from 1993 to 2021 is of 0.035 ± 0.002°C and 2.70 ± 0.83 mm, respectively (Copernicus Temperature and Sea Level Trend; https://marine.copernicus.eu/access-data/ocean-monitoring-indicators). The sea level rise is due to melting glaciers (Galassi and Spada, 2014; Marcos and Tsimplis, 2007; Lambeck and Bard, 2000) and thermal expansion (Criado-Aldeanueva et al., 2008; Tsimplis et al., 2008). In the surface layer, Pisano et al. (2020) found a yearly temperature trend of 0.041 ± 0.006°C for the whole MS from 1982 to 2018, using a satellite-based dataset. The trend shows uneven spatial patterns, with yearly trends of 0.036 ± 0.006°C and 0.048 ± 0.006°C for the Western Mediterranean and the Levantine-Aegean basin, respectively. In the period 1982-2020, Juza and Tintoré (2021) estimated trend values of 0.032 ± 0.002°Cyr^-1 and 0.044 ± 0.002°Cyr^-1 as averages in the Western and Eastern Mediterranean, respectively. Ibrahim et al. (2021) defined a trend of 0.033 ± 0.004°Cyr^-1 in the Ionian Sea and Levantine basin. The rapid sea surface warming in the last two decades was associated with a strong increase in marine heatwave days (Bensoussan et al., 2019; Ibrahim et al., 2021; Juza and Tintoré 2021; Hamdeno and Alvera-Azcarate, 2023; Dayan et al., 2023).

The warming of the surface layer is transported toward deeper layers in the dense water formation areas of the MS: the Gulf of Lion (GoL; Grignon et al., 2010; Durrieu de Madron et al., 2013; Tamburini et al., 2013; Estournel et al., 2016; Conan et al., 2018), the Southern Adriatic Pit (SAP; Vilibić and Supić, 2005; Artegiani et al., 1989; Klein et al., 2000; Gačić et al., 2001; Manca et al., 2002; Cardin et al., 2011; Rubino et al., 2012; Janeković et al., 2014), the Aegean Sea and the Rhodes Gyre (Theocharis and Georgopoulos, 1993; Lascaratos et al., 1999; Malanotte-Rizzoli et al., 2003; Nittis et al., 2003; Theocharis et al., 2014; Roether et al., 2007; Kubin et al., 2019). Moreover, boundary currents also play an important role in transporting the changing water mass characteristics to deeper layers (Waldman et al., 2018; Kubin et al., 2019; Pinardi et al., 2019).

The trends of the MS can also influence the thermohaline cell of the North Atlantic (Potter and Béthoux, 2004; Béthoux et al., 1998; Rahmstorf, 1998). Through the Strait of Gibraltar, Mediterranean Intermediate Waters are exported to the Atlantic, therefore monitoring the temperature and OHC variations in the MS can contribute to the understanding of changes at global level.

This study aims to better understand the temperature changes and the heat storage in the MS. We used sub-surface temperature measurements by Argo floats, which is considered the most practical method to quantify the changes in temperature and OHC (Hansen et al., 2011; Von Schuckmann and Le Traon, 2011; von Schuckmann et al., 2016). We present the first climatology of OHC from in situ data in the MS (2001 to 2020). Trends of temperature and OHC for surface, intermediate, and deeper layers were calculated from 2005 to 2020 for the entire Mediterranean, as well as for specified sub-basins (Western and Eastern Mediterranean, the WWestMed, the WWestMed North, the WWestMedSouth and the deep dense water formation areas such as the Gulf of Lion and the South Adriatic Gyre; see Figure 1 ).

Argo floats are autonomous quasi-Lagrangian instruments that move with the ocean currents. They are considered expendable platforms but the small dimension of the MS offers a good chance of potential recoveries (Gonzalez-Santana et al., 2023). With the help of an external bladder, they modify their volume and are able to ascend and descend in the water column. They are programmed to be parked at a specific drifting depth (mostly at 350 or 1000 m) where they stay for a specified period (1-10 days, mostly 5 days in the MS; Poulain et al., 2007). Then, they descend to greater depths - up to 2000 m and 4000 m in the case of Core or Deep Argo floats. During their ascent back to the surface, they measure temperature and salinity (and eventually other biogeochemical parameters) throughout the water column (Wong et al., 2020). The typical vertical resolution of the profiles are 2 m (for 0-100 m), 10 m (for 100-700 m) and 25 m (for 700-2000 m). At the surface, Argo floats transmit the data to satellites and descend again to repeat their diving cycle.

Quality controlled Argo float data were downloaded from the Coriolis Global Data Assembly Center - GDAC (ftp://ftp.ifremer.fr/ifremer/argo/dac/coriolis; Cabanes et al., 2021; Argo float data and metadata from Global Data Assembly Centre (Argo GDAC), SEANOE, https://doi.org/10.17882/42182, 2020.). For this study, only data with the best quality control flags (qc=1, “good data”) were considered for each float profile and, to additionally ensure the correctness of the data, a visual inspection of all the temperature profiles was performed. The quality controlled data were subsampled over time intervals of 5 days in order to obtain a homogenous data set, i.e. one Argo float profile every 5 day interval.

The total number of Argo float profiles subsampled every 5 days exceeds 38.000 ( Figure 2A ). The temporal and geographical distribution of these profiles show a strong increase starting in 2012. The geographical distribution (Western and Eastern Mediterranean) is similar for the different seasons ( Figure 2B ).

There are different methods to derive the OHC: it can be inferred by the surface heat flux from space, the thermal expansion of the sea (thermosteric sea level rise) from space or directly by using in situ temperature profiles. These subsurface temperature measurements are the most practical way to quantify the excess heat (Kirtman et al., 2013). In this study, we therefore use the available Argo float data to assess the climatology and the warming trends of temperature and OHC in the MS. For the OHC computation, different formulas were tested (Equation 1-3; ϱ z : depth dependent density; C p : specific heat capacity; θ : potential temperature; S: salinity; T: temperature) and the results were found to be consistent with one another.

For this study, we have chosen Equation 3, which has been used for recent estimates of OHC for the global ocean as well as for the MS (both for a depth range of 0-700 m and a time range from 2005 to 2019 and 1993 to 2019, respectively; https://marine.copernicus.eu/access-data/ocean-monitoring-indicators?category=105&region=all&search). The mean OHC per Argo float profile 〈 O H C 〉 is proportional to the mean temperature 〈 T 〉 within a certain depth interval Δ z = z 2 − z 1 , with a reference density of ϱ_{0 =} 1030 kgm^-3 and a specific heat capacity of C_p =3980 Jkg^-1°C^-1 (Von Schuckmann et al., 2009). Units of OHC are Jm^-2.

Float profiles available in the period 2001-2020 were organized in bins of 1°x1°. The climatologies of the 〈 O H C 〉 were estimated for each bin within different vertical layers (different   Δ z ) : the whole layer measured by the Argo floats (5-2000 m), the surface and intermediate layers considered together (5-700 m), the surface (5-150 m), intermediate (150-700 m) and deep (700-2000 m) layers separately. The OHC distribution in the deepest part of the MS is not yet reported in the literature, so the deepest layer measured by floats (1500-2000 m) was also considered for the climatologies.

The trends of temperature and OHC were estimated in the period from 2005 to 2020, due to a limited number of float profiles available before 2005. For the trend calculation, we considered surface and intermediate layers together (5-700 m), the deeper layer (700-2000 m) and surface (5-150 m) and intermediate (150-700 m) layers separately.

Trends were computed using a linear least squares method combined with an autocorrelation correction, following Mieruch et al. (2008) and Weatherhead et al. (1998). The reason for using such a correction is that environmental data are typically autocorrelated. Thus the typical assumption of uncorrelated (white noise) fit residuals in estimating linear trends cannot be justified, and would yield to an underestimation of the trend uncertainties, the errors. Our trend analysis and the interpretation is based on annual data, therefore we can use the simple trend model:

where Y_t are the annual parameters (e.g. temperature), μ is the ordinate, ω represents the trend, X_t is the time (years) and N_t is the autocorrelated noise and modelled as an autoregressive process of order 1 (AR[1]):

where φ is the autocorrelation at lag one (-1<φ<1) and ϵ_t is considered as the independent non-correlated noise component. For the details of the estimation of the trend and its error we refer to Mieruch et al. (2008). Finally, according to probability theory, an estimated trend ω is considered statistically significant with a probability of 95% if it is twice as large as its error:

The trend analysis of temperature and OHC was done for the whole Mediterranean Sea as well as for specified sub-basins and was expressed in °C per year [°C·yr^-1] and Watt per square metre [W·m^-2] per year, respectively.

This work adds new insights on the OHC distribution and variability in different depth layers of the MS, focusing on the differences among the sub-basins and on the behavior of deep dense water convection sites. The 1°x1° climatologies from 2001 to 2020 show an increase in OHC from west to east ( Figures 4 , 5 ). The dense water formation sites such as the Gulf of Lion, the SAP, the Aegean Sea and the Rhodes Gyre show the expected low amount of OHC ( Figure 4A ). Surface layers (5-150 m; Figure 4B ) show a high amount of OHC along the boundary currents within the Levantine Sea. Intermediate layers (150-700 m; Figure 4C ) also exhibit a high amount of OHC within the Levantine, the Aegean and within parts of the Ionian Sea. The deepest layer measured by the Argo floats (1500-2000 m; Figure 5C ) reveals four main different sub-basins regarding the OHC distribution, with an increase in OHC from west to east: the West Western Mediterranean, the Tyrrhenian Sea, the Ionian Sea and the Levantine Sea.

Time series of temperature and OHC from 2005 to 2020 revealed significant positive trends in the upper 700 m of the water column, of 0.041± 0.012°C·yr^-1 and of 3.59 ± 1.02 W·m^-2·yr^-1, respectively ( Figure 6 ; Table 1 ).

Compared to the annual OHC trend of 1.5 ± 1.02 Wm^-2, estimated in the period from 1993 to 2018 by Von Schuckmann et al. (2018; see Equation 3 and for dataset description: https://marine.copernicus.eu/access-data/ocean-monitoring-indicators?category=105&region=all&search), the annual trend derived from Argo float data in the period 2005 - 2020 appears more than doubled. Although this result comes out from two different datasets, it clearly indicates a significant increase of the heat accumulation rate over the last 15 years.

The upper 700 m of the Western Mediterranean are warming fastest with an increase in temperature of 0.070± 0.015°C·yr^-1, corresponding to an increase in OHC of 5.72 ± 1.28 W·m^-2·yr^-1 ( Table 1 ). These results pave the way for a new “view” of Mediterranean warming trends than those derived from SST-based estimations (Pisano et al., 2020; Juza and Tintore, 2021), which defined an uneven spatial pattern with increasing trends going from west to east. However, according to the float dataset, the Western basin holds less heat (lower OHC in Figure 4 ) than the Eastern basin (higher OHC in Figure 4 ), but with a steep warming trend ( Tables 1 , 3 ).

The southern part of the Western Mediterranean Sea (WWestMed South in Figure 1 ) exhibits the strongest sub-sub-basin warming trend of the MS for the upper 700 m, which is related to the warming of the resident Mediterranean Intermediate Water (MIW). The MIW is a water mass composed of a mixture of waters coming from the Eastern basin: the Levantine Intermediate (LIW) and the Cretan Intermediate Water (CIW, https://ciesm.org/MWM_Acronyms/MedWaterMassAcronyms.pdf). The heat gain in the MIW (0.019 ± 0.007°Cyr^-1, see Table 3 ) and in particular in the Western Intermediate Water (WIW, https://ciesm.org/MWM_Acronyms/MedWaterMassAcronyms.pdf) of 0.019 ± 0.007°Cyr^-1 also provides information about the warming of the Mediterranean Outflow Waters (MOW), the saline and warm water mass located in the intermediate depths of the North Atlantic and produced by the outflow of Mediterranean water through the Strait of Gibraltar (Millot et al., 2006; Naranjo et al., 2015; de Pascual-Collar et al., 2019). This water mass is one of the most important intermediate water masses in the North Atlantic; it is involved in the North Atlantic Deep Water formation processes and its properties and variability have a significant impact on global climate (Aldama-Campino and Döös, 2020).

Significant warming trends are evident in the intermediate and deep layers ( Figure 7 ) of the two deep convection sites in the MS (Gulf of Lion, South Adriatic; Figure 1 ), with an exceptionally strong warming trend in the South Adriatic in 2013-2020 ( Figure 7B , right panel; Table 1 ). This result is consistent with recent studies showing a steady temperature increase in the deep layer of this area (Civitarese et al., 2005; Cardin et al., 2011; Cardin et al., 2020). Thermohaline properties in the SAP are closely related to the periodic reversal (from anticyclonic to cyclonic and vice-versa) of the Northern Ionian Gyre (NIG; Gačić et al., 2010; Menna et al., 2019; Rubino et al., 2020; Gačić et al., 2021) on a quasi-decadal temporal scale (Civitarese et al., 2010; Mihanović et al., 2021). NIG reversals affect the water mass distribution among the Eastern Mediterranean sub-basins (Gačić et al., 2011; Bessières et al., 2013; Reale et al., 2017; von Schuckmann et al., 2019), in turn influencing the thermohaline properties of the whole Mediterranean on decadal and multi-decadal scales (Gačić et al., 2013; Schroeder et al., 2017; Placenti et al., 2022). The effect of the NIG reversal on the SAP intermediate layer results in a clear oscillation in the temperature time series ( Figure 7A , right panel). Larger values are observed during the cyclonic phase (2014-2015 and 2020-2021), which favors the inflow in the SAP of water of Levantine origin (warmer and saltier), and lower values during the anticyclonic phase, which favors the inflow of water of Atlantic origin (colder and less salty). Although less pronounced, a temperature oscillation is also observed in the deep layer that is in phase with that of the intermediate layer, emphasizing the role of convective mixing and winter convection in the vertical distribution of properties and associated long-term variations in the water column.

The Gulf of Lion also shows quasi-decadal variations in the intermediate layer ( Figure 7A , left panel), which are probably also related to the LIW temperature variations associated with the NIG reversals. Gačić et al. (2013) defined a travel time of the LIW between the Sicily Channel and the Gulf of Lion of ~ 10-12 yr, while Placenti et al. (2022) found a lag of ~ 9-10 yr between the NIG reversals and their effect in the Sicily Channel. Combining these two results, it can be argued that the temperature increase (decrease) in the intermediate layer of the Gulf of Lion ( Figure 7A , left panel) in 2014-2017 (2018-2019) can be related to the temperature increase (decrease) observed in the Sicily Channel in 1998-2006 (2007-2010), which in turn is related to the anticyclonic (cyclonic) NIG of 1988-1996 (1998-2005) (see Figure 6 in Placenti et al., 2022). Quasi-decadal temperature fluctuations in the deep layer in the Gulf of Lion ( Figure 7B , right panel) appear to be out of phase with the intermediate layer ( Figure 7A , right panel), revealing a weak communication between the two layers. This result can be a consequence of the absence of bottom-reaching convection in the Western Mediterranean since 2014 (Josey and Schroeder, 2013; Margirier et al., 2020). Under these conditions, the intermediate waters in the Gulf of Lion became warmer and saltier throughout the basin (Margirier et al., 2020), supporting the steep trend observed in Figure 7A (left panel) compared to the previous period ( Table 3 ).

The advection and diffusion of the temperature and the OHC leads to a general warming of the deeper sea. The sinking of surface water masses does not only take place within the deep dense water formation sites, but also along boundary currents. This has been shown by theoretical (Waldman et al., 2018; Pinardi et al., 2019) and experimental (Kubin et al., 2019) studies within the MS: the net sinking along the boundary currents is due to the conservation of potential vorticity and serves to balance the friction along the boundary currents. As a consequence, the warming of deeper layers might be more intense than previously thought and can potentially shift deeper ecosystems within the deep dense water formation sites (Tamburini et al., 2013; Coma et al., 2009).

In fact, the high amount of OHC in the surface layer along the boundary currents in the Levantine Sea ( Figure 4B ) is transported to intermediate (and eventually also to deeper) layers. The warm Levantine Intermediate Water (LIW) enters the Aegean Sea through the Cretan Straits as Modified Levantine Intermediate Water (MLIW; Velaoras and Lascaratos, 2010), leading to warm Aegean intermediate waters and as a consequence also to warm intermediate waters exhibiting a high amount of OHC within the Ionian Sea (outflow of the Aegean Sea; Figure 4C ).

The warming of these surface and intermediate water masses also influences the preconditioning for the deep convection site in the Gulf of Lion (Grignon et al., 2010; Estournel et al., 2016).

The subduction and spreading of surface water anomalies and the consequent warming of intermediate and deep water masses are already showing their feedback on ocean dynamics and the atmosphere (air-sea fluxes) and will continue to do so in the coming years, decades, or even centuries as these warming water masses spread or re-emerge (Artale et al., 2018; Lo Bue et al., 2021). Marine heatwaves have also become and will become more frequent (Frölicher et al., 2018; Darmaraki et al., 2019a; Darmaraki et al., 2019b; Oliver et al., 2021; Juza et al., 2022; Hamdeno and Alvera-Azcarate, 2023). These unprecedented changes will stress ecosystems and accelerate the loss of biodiversity and the extinction of several marine species who cannot adapt to such rapid temperature changes (Garrabou et al., 2022; Smith et al., 2022; Fragkopoulou et al., 2023). Moreover, impacts of climate change also threaten indigenous cultures and their knowledge about the ocean and ecosystems (Bindoff et al., 2019).

Therefore, this study should act as another wake up call for policy makers and society. The stabilization of the climate - agreed by the UNFCCC in 1992 and the Paris agreement in 2015 - would require that the EEI is reduced to approximately zero in order to achieve the Earth’s system quasi-equilibrium, which corresponds to a decrease of atmospheric CO₂ from 410 ppm to 353 ppm (Von Schuckmann et al., 2020). This can only happen if social, economic and ecological development are no longer seen as separate parts, but as interconnected (Berkes, 2017; Von Schuckmann et al., 2020).

The 2021-2030 UN Decade of Ocean Science for Sustainable Development aims to create a more holistic and integrated approach, with an emphasis also on indigenous people and the traditional knowledge of local people to achieve a truly sustainable approach and not just green- or bluewashing projects. Blue economy and blue growth is often doing more harm than good, because it is still based on exploitation and on the concept of economic growth on a finite planet (Ehlers, 2016), while indigenous people and their traditional conservation and management are based on the respect of nature and taking care of the land (Minerbi, 1999; Kealiikanakaoleohaililani and Giardina, 2016; Berkes, 2017; Witte and Xué, 2018; UNESCO Man and Biosphere Programme, Reed, 2019). Such holistic ways of understanding the environment offer alternatives to the prevailing consumption-oriented values of Western societies (Berkes and Turner, 2006; Kimmerer, 2012; Brondízio et al., 2021).

Therefore, possible solutions cannot only be of technological nature, but require an urgent and strong shift of our way of thinking and of our entire worldview to make sure that future generations can experience healthy, living oceans and ecosystems.

Publicly available datasets were analyzed in this study. This data can be found here: Argo float data and metadata from Global Data Assembly Centre (Argo GDAC), SEANOE, https://doi.org/10.17882/42182, 2020.

EK: Conceptualization, Data curation, Investigation, Methodology, Writing – original draft, Writing – review & editing, Software, Visualization. MM: Funding acquisition, Investigation, Visualization, Writing – original draft, Writing – review & editing. EM: Funding acquisition, Writing – review & editing. GN: Data curation, Writing – review & editing. SM: Investigation, Software, Writing – original draft, Writing – review & editing. PP: Conceptualization, Funding acquisition, Methodology, Writing – review & editing.