Temperature and radiative responses to anthropogenic aerosols over the Mediterranean Basin based on CMIP6 Earth system models

Here, we assess the amplification of near-surface warming in the Mediterranean (MED) resulting from global anthropogenic aerosol (AA) reductions, based on simulations from CMIP6 Earth system models (ESMs). The effective radiative forcing (ERF) and near-surface temperature (TAS) exhibit decreasing trends until around 1980 followed by increasing trends, driven by air pollution control policies. The annual mean ERF at the top-of-atmosphere over the MED changes by $\mathrm{2.37}\pm \mathrm{1.06}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ between the peak AA period (1970–1979) and the near-present period (2005–2014). During this interval, the annual mean TAS increases by 0.67±0.37 °C. Overall, the multi-model ensemble shows a robust amplification of warming over the MED on annual scale resulting from global AA reductions from 1970–1979 to 2005–2014, in good agreement with observational datasets over land. The model simulations indicate that AAs are responsible for 49 % (39 %) of the annual (summer) warming between the two periods. In the winter, ESMs produce an overestimated warming of 1.19 °C, with AAs contributing 60 % to this warming. Finally, we show that circulation changes caused by AA reductions can play an additional role in the redistribution of regional temperature changes apart from the radiative effects per se. Our results reveal a strong link between the recent acceleration of MED warming and global AA decreases, which unmask additional greenhouse gas-driven warming. This study highlights the sensitivity of the MED to recent global AA emission changes and the need for climate policies that couple air quality improvements with rapid greenhouse gas mitigation.

1 Introduction

The Mediterranean Basin (MED) region is considered a hotspot of climate change (Cos et al., 2022; Lionello and Scarascia, 2018) that has experienced enhanced warming since the late 20th century, particularly during the summer (Doblas-Reyes et al., 2021; Gutiérrez et al., 2021). This trend cannot be fully explained by internal climate variability alone and can be attributed to external anthropogenic forcing (Allen et al., 2006; Hodnebrog et al., 2024; Rashid, 2021). The decrease of anthropogenic aerosols (AAs) over Europe, driven by European air pollution control policies (Turnock et al., 2016), has been linked to temperature increases in the Euro-Mediterranean region (Boé et al., 2020; Dong et al., 2017; Glantz et al., 2022; Nabat et al., 2014, 2025; Urdiales-Flores et al., 2023; Zanis et al., 2012). However, other parts of the world such as North America (primarily the U.S. with the Clean Air Act) have also reduced their AA emissions, which could also influence the European climate (Elkins et al., 2025; Yang et al., 2020). Global concentrations of AAs in the atmosphere have decreased since around the 1980s, as suggested by changes in aerosol optical depth (AOD) over Europe and Eastern North America (Cherian and Quaas, 2020; Gulev et al., 2021; Klimont et al., 2013; Wei et al., 2019). Sulfate atmospheric loads have largely reduced during the last decades (Aas et al., 2019; Gulev et al., 2021), with subsequent radiative changes (Kalisoras et al., 2024; Myhre et al., 2017a; Szopa et al., 2021).

Global AA decrease has caused changes in the Earth's radiative budget, which are expressed in terms of the effective radiative forcing (ERF) (Forster et al., 2021; Hansen et al., 2005; Sherwood et al., 2015). According to the Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report (AR6) WGI (Forster et al., 2021; Gulev et al., 2021; Szopa et al., 2021), the ERF caused by AAs at the top-of-atmosphere (TOA) was negative throughout the industrial period over Europe, with a shift by the late 1970s–early 1980s towards less negative values and a weakening of the negative ERF in the first decades of the 21st century. Parallel to ERF weakening, an increase in shortwave radiation reaching the surface, a phenomenon referred to as the “brightening effect” (Wild, 2009, 2012), has been observed over Europe and the MED (Alexandri et al., 2017; Galanaki et al., 2021; MedECC, 2020; Nabat et al., 2014; Sanchez-Lorenzo et al., 2015; Wild, 2012). Arguably this can be attributed to reductions in sulfate aerosols (Allen et al., 2013; Nabat et al., 2014), which scatter incoming solar shortwave radiation and interact with clouds (Bellouin et al., 2020; Kasoar et al., 2016; Myhre et al., 2013).

A number of recent studies used observational datasets, reanalysis products, and global and regional climate models (GCMs and RCMs, respectively) to examine the enhanced surface warming and increase in temperature extremes over the Euro-Mediterranean region due to reductions in global AAs, which earlier had masked the continuous warming from greenhouse gases until around the 1980s (Boé et al., 2020; Cos et al., 2022; Dong et al., 2017; Drugé et al., 2021; Feng and Qian, 2024; Glantz et al., 2022; Hodnebrog et al., 2024; Nabat et al., 2014, 2025; Navarro et al., 2017; Samset et al., 2018; Tian et al., 2020; Urdiales-Flores et al., 2023). Decreases in AAs caused radiative imbalances over Europe and a shift from dimming to brightening, which was linked to increased warming at the surface (Folini and Wild, 2011; Glantz et al., 2022; Hodnebrog et al., 2024; Nabat et al., 2014, 2025; Sanchez-Lorenzo et al., 2015). Dong et al. (2017) showed that near-surface temperature and temperature extremes have increased since the 1990s in Western Europe during the summer, with AAs playing a key factor. Glantz et al. (2022) found an increase of about 1 °C attributed to AA decreases between the months April and September during the period 1979–2020 over Central and Eastern Europe. According to Nabat et al. (2025), each 0.1 reduction in AOD over Central Europe led to a 0.3 °C rise in near-surface temperature. In an earlier study, Nabat et al. (2014) showed that AA changes explain 81±16 % of the brightening and 23±5 % of the surface warming simulated for the period 1980–2012 over the European region based on reanalysis-driven coupled regional climate model. Turnock et al. (2016) estimated a mean surface temperature increase of around 0.45±0.11 °C over Europe between 1970 and 2010 from European AA reduction measures. Urdiales-Flores et al. (2023) carried out an attribution analysis to understand the drivers of accelerated warming since the 1980s over the MED (based on observations, the median trend for 1981–2020 exceeded 0.4 °C per decade). They calculated a surface warming of 0.10 °C per decade over the MED for the period 1981–2020 attributed to AOD reductions, while they attributed 0.18 °C per decade to radiative forcing from well-mixed greenhouse gases and 0.03 °C per decade to soil moisture decrease. Boé et al. (2020) investigated summer temperature changes over Europe using global and regional climate models, and assessed the influence of AAs in surface warming with sensitivity experiments (also discussed in Boé, 2016). Yang et al. (2020) examined the trends of AAs over Europe during 1980–2018 and found a warming of 2.0 W m⁻² in Europe due to local decreases in sulfate loading, with 12 % attributed to changes in non-European emissions. However, circulation changes attributed to reduced AA emissions also play an important role in near-surface temperature changes over Europe (Delworth et al., 2022; Dong and Sutton, 2021; Nabat et al., 2020; Zanis, 2009).

Despite the plethora of studies investigating aerosol-driven changes over the Euro-Mediterranean region, there still remain several knowledge gaps: (a) how do global climate models simulate the Mediterranean response to global anthropogenic aerosol emission changes on basin-wide scale, (b) how do changes in aerosol optical depth, the radiative budget, and atmospheric circulation translate to horizontal and vertical temperature changes over the Mediterranean, (c) how well do climate models perform when simulating the observed land temperature amplification on annual and seasonal scales?

This study aims to highlight the sensitivity of the Mediterranean region to global emission changes of anthropogenic aerosols by building on existing knowledge. We address the aforementioned knowledge gaps using high-complexity, state-of-the-art Earth system models (ESMs) endorsed by the sixth phase of the Coupled Model Intercomparison Project (CMIP6) (Eyring et al., 2016) to investigate temperature changes over the Mediterranean Basin in terms of radiative and circulation changes throughout the historical period, focusing on two periods of interest: the decade that AAs peak over Europe (1970–1979) and the near-present decade (2005–2014). To our knowledge, this is the first study that calculates changes in TOA and surface radiative fluxes to quantify fast and transient (i.e., a mixture of fast and incomplete slow responses from coupled-ocean simulations) radiative responses resulting from global changes in AA emissions and tries to connect them with changes in temperature both at the surface and throughout the troposphere on annual and seasonal scales. Furthermore, AA-induced changes to atmospheric circulation are examined in an effort to explain temperature changes not related to AA radiative effects. Finally, observational datasets are used to evaluate the temperature changes over the MED simulated by ESMs on both annual and seasonal scales. The rest of the paper is structured as follows. Details about the CMIP6 ESMs and simulations used, along with a description of our methodology, are given in Sect. 2. Our results are discussed in Sect. 3, while in Sect. 4 the main conclusions of our work are summarized.

2 Data and methodology

2.1 Data description

In this work, historical simulations from eight ESMs were used (Table 1), which were performed within the framework of the Aerosol Chemistry Model Intercomparison Project (AerChemMIP; Collins et al., 2017) and supported by CMIP6 (Eyring et al., 2016). Anthropogenic emissions of aerosols as well as precursors of aerosols and ozone (excluding methane) used by climate models are taken from van Marle et al. (2017) and Hoesly et al. (2018), while each model uses its own natural aerosol (e.g., dust and sea salt) emissions. All models include both aerosol indirect effects (i.e., the cloud albedo effect and the cloud lifetime effect) (Albrecht, 1989; Twomey, 1974, 1977), except for CNRM-ESM2-1 which only represents the cloud albedo effect (and hence is lacking a cloud lifetime effect). Additionally, MRI-ESM2-0 and NorESM2-LM parameterize aerosol effect on ice clouds, thus impacting the LW ERF (e.g., Smith et al., 2020; Zelinka et al., 2023). A short description of each ESM is given in the Appendix.

Table 1 Information on model resolution (horizontal and vertical), variant label, and references for each ESM used in this study. Each experiment (see Table 2) has a
variant label r_ai_bp_cf_d, where a is the realization index, b is the initialization index, c is the physics index, and d is the forcing index.

^a The hist-piAer simulation is identical to the hist-piNTCF simulation as CNRM-ESM2-1 has no tropospheric chemistry, and therefore no ozone precursors, which means that the two configurations (hist-piAer and hist-piNTCF) are identical.^b The 0.7°×0.7° is approximate for the TL255 grid of IFS. Aerosols and atmospheric chemistry are simulated with the Tracer Model version 5 (TM5) with horizontal resolution 3°×2° (lon × lat), with 34 levels in the vertical and a model top at about 0.1 hPa.^c GFDL-ESM4 uses a cubed-sphere grid with ∼100 km resolution (1o nominal horizontal resolution). Results are regridded to a 1.25°×1° (lon × lat) grid for analysis.^d The 1.125°×1.125° is approximate for the MRI-AGCM3.5. Aerosols are simulated with the Model of Aerosol Species in the Global Atmosphere mark-2 revision 4-climate (MASINGAR mk-2r4c) with horizontal resolution 1.875°×1.875° (lon × lat), 80 vertical levels and a model top at 0.01 hPa.^e The distinction between p1 and p2 exists only for prescribed-SST simulations. These versions differ slightly in the treatment of the fluxes of energy, moisture, and momentum at the surface between the ocean and the atmosphere. The treatment of these fluxes in the p2 prescribed-SST simulations is closest to the fully-coupled (including active ocean) simulations. The p1 prescribed-SST simulations have a slightly different formulation of these fluxes (the flux calculation here corresponds more with the standard flux calculation in CESM2). Hence, the p2 version is closest related to the fully-coupled simulation (Dirk Olivié, personal communication, 2023).

Within the framework of CMIP6 and AerChemMIP ESMs performed two sets of transient historical experiments for the period between 1850 and 2014 (Table 2). The first set includes the CMIP6 historical (Eyring et al., 2016) and hist-piAer coupled-ocean experiments, and the second set includes the histSST and histSST-piAer experiments with prescribed sea surface temperatures (SSTs) and sea ice (Collins et al., 2017). The historical experiment branches from the piControl experiment and is forced by time-evolving, externally imposed anthropogenic and natural forcings that are based on observations, such as solar variability, volcanic aerosols, and changes in land use and greenhouse gas and aerosol emissions or concentrations. The histSST experiment shares the same forcings as the historical and SSTs and sea ice are prescribed to the monthly mean time-evolving values from the corresponding historical simulation of each model. Including evolving SSTs, rather than using a fixed repeating climatology throughout the simulation, means that the underlying climate state in histSST is consistent with the historical experiment. The hist-piAer and histSST-piAer experiments are identical to historical and histSST, respectively, but they use fixed pre-industrial anthropogenic aerosol and aerosol precursor emissions corresponding to the year 1850, while all other forcings remain the same.

Table 2 List of transient historical simulations covering the historical period (1850–2014) used in this work. The histSST and histSST-piAer experiments use prescribed SSTs and sea ice, whereas in the historical and hist-piAer experiments SSTs are allowed to evolve. The year indicates that the emissions or concentrations are fixed to that year, while “Hist” means that the concentrations or emissions evolve as for the historical experiment (more information in Collins et al., 2017; Eyring et al., 2016).

To evaluate the near-surface air temperature evolution in CMIP6 coupled-ocean experiments, three gridded observational datasets (Schneider et al., 2013) were employed: (a) the Berkeley Earth (https://berkeleyearth.org/data/, last access: 13 June 2025) dataset (Rohde and Hausfather, 2020), which provides high-resolution (1°×1°) monthly global land-ocean temperature data starting from the year 1850, (b) the Climatic Research Unit of the University of East Anglia (https://crudata.uea.ac.uk/cru/data/hrg/, last access: 13 June 2025) gridded Time Series (CRU TS v4.08) (Harris et al., 2020), which offers a high-resolution (0.5°×0.5°) monthly grid of land-based observations going back to 1901, and (c) the University of Delaware land temperature monthly gridded data (v4.01) at high resolution (0.5°×0.5°) beginning in 1900 and provided by the National Oceanic and Atmospheric Administration (NOAA) Physical Sciences Laboratory (PSL), Boulder, Colorado, USA, from their website (https://psl.noaa.gov, last access: 13 June 2025). The CRU TS and University of Delaware datasets contained absolute temperature values (in °C). The Berkley Earth dataset contained temperature anomalies relative to the 1951–1980 baseline climatology, which were converted to absolute temperature values (in °C).

2.2 Calculation methodology

In this work, monthly outputs from the ESM simulations were used (summarized in Table 3). All data were first regridded to a common spatial grid (2.8125°×2.8125°) by applying bilinear interpolation. The effects of changes in global anthropogenic emissions of aerosols on AOD, near-surface air temperature (TAS), and circulation (sea level pressure, geopotential height, and winds) were calculated by subtracting hist-piAer from the historical simulation (i.e., the coupled experiments).

Table 3 Description of the CMIP6 variables used in this study.

The transient aerosol ERF was calculated over the historical period using experiments with prescribed SSTs and sea ice. By fixing SSTs and sea ice, all other parts of the system are allowed to respond until reaching steady state (Hansen et al., 2005; Sherwood et al., 2015). While not technically an ERF (since SSTs and sea ice are evolving), the impact of transient SSTs and sea ice on ERF is regarded to be small (Collins et al., 2017; Forster et al., 2016). Here, ERF was estimated at the top-of-atmosphere (ERF_TOA) and at the surface (ERF_SURF) following Forster et al. (2016). ERF_TOA was calculated as the difference between the incoming and outgoing radiative fluxes at the top-of-atmosphere (Eq. 1), while ERF_SURF as the difference between the surface downwelling and upwelling radiative fluxes (Eq. 2):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{\text{ERF}}_{\text{TOA}}=\mathrm{\Delta }{F}_{\text{TOA}}=\mathrm{\Delta }[\text{rsdt}-(\text{rsut}+\text{rlut}\left)\right],\text{(2)}& {\displaystyle }& {\displaystyle }{\text{ERF}}_{\text{SURF}}=\mathrm{\Delta }{F}_{\text{SURF}}=\mathrm{\Delta }\left[\right(\text{rsds}+\text{rlds})-(\text{rsus}+\text{rlus}\left)\right],\end{array}$$

where F_TOA and F_SURF are the net (downward minus upward) radiative fluxes at the top-of-atmosphere and at the surface, respectively, and Δ denotes the difference between the perturbation and the control experiment (histSST minus histSST-piAer).

Supplementarily with ERF, we calculated the AA-induced change in net radiative fluxes at the top-of-atmosphere and the surface (denoted as $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ and $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$, respectively) using the same equations as above and by taking the difference between historical and hist-piAer this time. While ERF depicts the radiative perturbation caused by AAs after accounting for rapid adjustments and fast climate responses (since SSTs and sea ice are prescribed), the net radiative flux differences calculated from the coupled-ocean experiments ($\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ and $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$) quantify a quasi–radiative forcing caused by AAs, accounting for a mixture of fast and incomplete slow responses (due to climate feedbacks and ocean heat uptake) under the evolving transient historical conditions. Hence, the temperature differences from the coupled-ocean experiments represent transient climate responses being far from an equilibrium state (in which the initial radiative aerosol perturbation will be fully absorbed by the climate system response), since the ocean does not have enough time to fully equilibrate to the evolving historical forcings.

The statistical significance of each ESM's results was tested at the 95 % confidence level using a two-tailed t-test. The robustness of the multi-model ensemble results was estimated using a comprehensive method based on statistical significance and the agreement between ESMs on the sign of change. In cases where at least 80 % (i.e., at least 6 out of 8) of the ESMs show statistically significant change and at least 80 % of the ESMs agree with the multi-model mean on the sign of change, then these results are considered robust. If less than 80 % (i.e., 5 or less out of 8) of the ESMs show not statistically significant change or less than 80 % of the models do not agree with the multi-model mean on the sign of change, then they are characterized as non-robust or conflicting, respectively.

Here, we give emphasis on inter-model variability, and the robustness of the multi-model ensemble results was estimated using a comprehensive method based on statistical significance and the agreement between ESMs on the sign of change. This is a common approach of quantifying uncertainty, and particularly inter-model uncertainty (e.g., Liu et al., 2018; Myhre et al., 2017b; Samset et al., 2018; Szopa et al., 2021). Calculating the multi-model ensemble mean reduces internal model variability, as the averaging suppresses the random/internal component of each individual model (since internal variability is largely uncorrelated across the different models) and enhances the common forced signal. It should be stressed that, although the multi-model mean dampens internal variability, there is still need of further analysis using long-term ensemble of piControl experiments (or ensemble spread) for each model to evaluate whether the remaining signal is statistically distinguishable from noise (Collins et al., 2017). The accurate assessment of individual model uncertainty related to model internal variability would require initial condition large ensemble members from single models (Jones et al., 2024; Parsons et al., 2020).

To study the effects caused by global AA reduction, two periods were selected: (a) 1970–1979, which corresponds to the peak aerosol decade over Europe, and (b) 2005–2014 as the last decade available in CMIP6 historical simulations which represents the near-present. When examining a single period (1970–1979 or 2005–2014), changes reflect the effect of AAs relative to the pre-industrial period, whereas the difference between the two periods (2005–2014 minus 1970–1979) represents the effects of AA decreases since around 1980. In the next section, the multi-model ensemble results are presented for the study region (10° W–40° E, 30–45° N). Area-weighted mean values refer to the multi-model means unless stated otherwise.

3 Results

3.1 Temporal evolution of AOD and radiative changes

AA concentrations peaked in the late 1970s to early 1980s in Europe and the MED. Consequently, AOD change reached a maximum during that period (Fig. 1, left column). The AOD maximum is larger in magnitude during the boreal summer (JJA) than the winter (DJF).

Figure 1 Timeseries of the AA-induced evolution of aerosol optical depth (ΔAOD; left column), effective radiative forcing (ERF; middle column), and radiative flux change (ΔF_coupled; right column) over the Mediterranean region (10° W–40° E, 30–45° N) for the historical period (1850–2014) on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). Changes in AOD and ΔF_coupled were calculated from the difference historical minus hist-piAer, whereas ERF from the difference histSST minus histSST-piAer. The lines represent the 10 year moving average, while shading expresses the inter-model variability (one standard deviation).

Likewise, the transient annual mean ERF both at TOA and the surface attained its most negative values during the late 1970s to early 1980s, dominated by the evolution of ERF in JJA (Fig. 1, middle column). In DJF both ERFs (ERF_TOA and ERF_SURF) are weaker in magnitude than in JJA resulting from less incoming radiation. In the winter ERF_SURF becomes weaker than ERF_TOA during 1975–1985, which is not the case for the summer or annual ERFs. Inter-model variability is large in DJF, with values reaching up to around $\pm \mathrm{3}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ during the negative peak.

In a similar manner, the annual mean ΔF_coupled both at TOA and the surface show a trend towards more negative values until around the mid-1980s followed by a reverse positive trend afterwards (Fig. 1, right column). The annual mean $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ is weaker in magnitude than ERF_TOA as the climate system has partly adjusted to the AA-induced radiative perturbation through ongoing slow responses and changes in land and ocean temperatures. A weakening in magnitude can be seen towards the end of the historical period, driven by the evolution of their shortwave (SW) rather than their longwave (LW) component (not shown). $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ and $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$ show smaller changes throughout the historical period during the winter than during the summer. In JJA, $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ attains its most negative values in the late 1970s – early 1980s, while $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$ exhibits the most negative peak in the mid-1980s, slightly later than ERF_SURF. Overall, the SW component practically drives the negative ΔF_coupled and its temporal evolution. This can be explained by the evolution of the AOD of sulfate aerosols over the MED (Fig. S1 in the Supplement), which is the dominant AA type and closely follows the trend of total AOD (Fig. 1). Changes in the AOD of other types of AAs are nearly an order of magnitude smaller than those of sulfates, while dust aerosols show little to no change (see also Haugvaldstad et al., 2025). This reveals that sulfates are the main driver of AOD change, consequently affecting ERF and ΔF_coupled. However, interactions between aerosols and clouds could also play a role in the simulated radiative changes at TOA.

The temporal evolution of the above quantities is in line with previous studies suggesting that reductions in (mainly sulfate) AOD led to increased SW radiation reaching the surface (enhanced brightening effect) over the Euro-Mediterranean region since the 1980s (e.g., Glantz et al., 2022; Nabat et al., 2014, 2025; Sanchez-Lorenzo et al., 2015). Global reductions in sulfate emissions have a significant effect over Europe and cause a positive radiative forcing, which is linked to surface warming, in agreement with the findings of Yang et al. (2020).

3.2 Spatial pattern of AOD changes

High emissions of AAs and their precursors in the period 1970–1979 over Europe have led to an increase in AOD mainly in Central Europe and its surrounding regions relative to the pre-industrial period, thus increasing AOD in the northern parts of the MED (Fig. 2). The annual mean AOD change due to AAs during 1970–1979 over the MED is 0.11±0.04 (inter-model variability; 1 standard deviation), with changes being slightly larger (smaller) in magnitude in JJA (DJF). Changes in AA emissions and concentrations have led to reductions in AOD in 2005–2014 (with respect to 1970–1979) in Central European regions and consequently in the MED of $-\mathrm{0.05}\pm \mathrm{0.02}$ annually, with the largest AOD reductions shown in JJA ($-\mathrm{0.07}\pm \mathrm{0.02}$). Spatially, Southeastern Europe (Balkan Peninsula) shows the largest annual and seasonal AOD changes between the periods 1970–1979 and 2005–2014 within the MED region. These results are in line with the findings of Aas et al. (2019) and Cherian and Quaas (2020) over the Euro-Mediterranean region. Our multi-model results are robust for most of the MED (except the southern and western parts) on both annual and seasonal scales. The MPI-ESM-1-2-HAM and MRI-ESM2-0 models show the largest AOD change between the two periods (−0.086 and −0.087, respectively; Table 4), with the former exhibiting the largest change in DJF (−0.115; Table 5) and the latter in JJA (−0.127; Table 6).

Figure 2 Changes in AOD due to anthropogenic aerosols over the Mediterranean (calculated from the difference historical minus hist-piAer). Spatial patterns are presented for the periods 1970–1979 (left column), 2005–2014 (middle column) and their difference (2005–2014 minus 1970–1979; right column) on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). The multi-model ensemble means over the Mediterranean region (shown as a box) are shown along with the inter-model variability (one standard deviation) at the top of each panel. Colored areas without markings indicate robust changes, while hatched ($/$) and cross-hatched (×) areas indicate non-robust changes and conflicting signals, respectively.

Table 4 Annual mean values of AOD change (ΔAOD; unitless), ERF at TOA and at the surface (in W m⁻²), ΔF_coupled at TOA and at the surface (in W m⁻²), and near-surface air temperature change (ΔTAS; in °C) over the Mediterranean for the periods 1970–1979, 2005–2014 and their difference (2005–2014 minus 1970–1979). Values are presented for each model, along with the multi-model ensemble mean and the inter-model variability (one standard deviation; SD).

Table 5 As in Table 4, but for the winter (DJF).

Table 6 As in Table 4, but for the summer (JJA).

3.3 Spatial pattern of ERF changes

ERF_TOA spatial patterns are presented in Fig. 3, while ERF_SURF spatial patterns are illustrated in Fig. S2 in the Supplement. AAs exert a negative ERF_TOA over Central Europe in 1970–1979 relative to the pre-industrial period on annual scale (Fig. 3). ERF_TOA greatly increases in magnitude in JJA and spreads mainly over the eastern and partly over the western and northern parts of Europe. However, it weakens in DJF and relocates over the MED. During the summer the largely negative ERF_TOA over regions north of the MED is offset by the positive ERF_TOA over the Sahara Desert, leading to a multi-model mean ERF_TOA of $-\mathrm{4.16}\pm \mathrm{1.74}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ over the MED. This is slightly weaker than the annual mean of $-\mathrm{4.43}\pm \mathrm{1.55}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. ERF_SURF exhibits a spatial pattern that is very similar to ERF_TOA, but it is stronger (i.e., more negative) due to the more negative ERF_SURF over Africa (Fig. S2). Unlike ERF_TOA, ERF_SURF in JJA is larger in magnitude than the annual mean ERF_SURF.

Figure 3 ERF_TOA due to anthropogenic aerosols over the Mediterranean (calculated from the difference histSST minus histSST-piAer). Spatial patterns are presented for the periods 1970–1979 (left column), 2005–2014 (middle column) and their difference (2005–2014 minus 1970–1979; right column) on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). The multi-model ensemble means over the Mediterranean region (shown as a box) are shown along with the inter-model variability (one standard deviation) at the top of each panel. Colored areas without markings indicate robust changes, while hatched ($/$) and cross-hatched (×) areas indicate non-robust changes and conflicting signals, respectively.

The weakening of the forcing due to AA reduction leads to a smaller annual mean ERF_TOA of $-\mathrm{2.06}\pm \mathrm{1.00}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ over MED in 2005–2014. ERF_TOA becomes even smaller in JJA ($-\mathrm{0.99}\pm \mathrm{1.70}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) because of the weaker radiative forcing gradient between Central Europe and Northern Africa. However, the multi-model results are robust only over the Eastern MED on an annual and winter scales. The large inter-model variability in the summer during the period 2005–2014 stems from MRI-ESM2-0 (Table 4), which simulates a strongly negative ERF_TOA ($-\mathrm{4.77}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$). The MRI-ESM2-0 model produces the most negative ERF_TOA and ERF_SURF over the MED in all cases except in the winter of 2005–2014, where MPI-ESM-1-2-HAM takes the lead in both ERFs (Table 5).

In general, the ensemble differences of ERF_TOA and ERF_SURF between the periods 2005–2014 and 1970–1979 indicate that all ESMs point towards an ERF weakening (less negative ERF) that is more prominent in JJA. The summer ERF change dominates the spatial pattern of the annual mean change in ERF (mind that the differences in DJF are weaker and non-robust). The annual and summer ERF weakening is confined to the European part of the MED, where the multi-model results are robust. The largest changes in ERF_TOA and ERF_SURF between the two examined periods are seen for EC-Earth3-AerChem, MPI-ESM-1-2-HAM, MRI-ESM2-0, and UKESM1-0-LL, which differ in magnitude depending on the season. Our results regarding ERF_TOA over the MED are in agreement with the findings of studies examining regional ERFs (e.g., Kalisoras et al., 2024; Szopa et al., 2021).

3.4 Spatial pattern of ΔF_coupled changes

Changes in the spatial pattern of $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ are shown in Fig. 4, while the respective changes in the spatial distribution of $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$ are presented in Fig. S3 in the Supplement. Overall, the spatial pattern of $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ is nearly identical to that of ERF_TOA but weaker in magnitude. The annual mean $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ in 1970–1979 is $-\mathrm{2.75}\pm \mathrm{1.78}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ over the MED and is largest over Central and Eastern Europe (Fig. 4). As in the case of ERF_TOA, the annual spatial pattern of $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ in the period 1970–1979 is dominated by its respective summer pattern. The negative $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ in JJA is spread over the European part of the MED region and has a multi-model mean value of $-\mathrm{1.97}\pm \mathrm{2.42}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ over the MED. On the other hand, over the southern and southeastern parts of MED, $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ becomes positive (albeit not statistically significant), with a magnitude greater than the respective ERF_TOA magnitude over the same areas. In the winter of 1970–1979, $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ is generally weaker ($-\mathrm{1.91}\pm \mathrm{1.56}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) and robust only over the Balkan Peninsula. After global AA emission reductions, the annual mean $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ increases by $\mathrm{1.95}\pm \mathrm{1.09}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ but remains negative. Similarly to ERF_TOA, the largest increases are seen over Central Europe. In the winter, there is a non-robust increase in $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ over the MED, mainly due to increases over its northern parts. Conversely, the summer $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ undergoes a significant increase and becomes positive in 2005–2014 ($\mathrm{0.27}\pm \mathrm{1.96}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$). Although there are conflicting signals on the summer mean $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ in 2005–2014, all ESMs simulate a trend towards more positive values in the transition from 1970–1979 to 2005–2014.

Figure 4 ΔF_coupled at TOA due to anthropogenic aerosols over the Mediterranean (calculated from the difference historical minus hist-piAer). Spatial patterns are presented for the periods 1970–1979 (left column), 2005–2014 (middle column) and their difference (2005–2014 minus 1970–1979; right column) on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). The multi-model ensemble means over the Mediterranean region (shown as a box) are shown along with the inter-model variability (one standard deviation) at the top of each panel. Colored areas without markings indicate robust changes, while hatched ($/$) and cross-hatched (×) areas indicate non-robust changes and conflicting signals, respectively.

The spatial distribution of $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$ (Fig. S3) resembles that of ERF_SURF (Fig. S2) but, unlike the above, their annual and seasonal means over the MED are of similar magnitude. A notable paradox is the seasonal change of $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$ after the decrease in global AA emissions (Fig. S3; right column). The DJF $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$ shows a small but robust increase ($\mathrm{0.95}\pm \mathrm{1.17}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$), especially over continental Europe. This winter increase is not robust over such a large area for the rest of the metrics of radiative flux change (ERF_TOA, ERF_SURF, $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$). On the other hand, the change in the JJA $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$ is larger in magnitude ($\mathrm{2.71}\pm \mathrm{1.26}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) than in DJF, but it is not a robust feature among ESMs.

After examination of the SW and LW components of ΔF_coupled at TOA and the surface (not shown), it becomes clear that the SW component drives the pattern and evolution of ΔF_coupled (indicating SW radiative loss at TOA and the surface), while the LW component partly masks the ΔF_coupled change. The dominance of SW aerosol forcing is consistent with earlier studies emphasizing the role of aerosols in shaping the radiative balance over the Euro-Mediterranean region (e.g., Boé et al., 2020; Glantz et al., 2022; Nabat et al., 2014). The LW $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ and $\mathrm{\Delta }{F}_{\text{SURF}}^{\text{coupled}}$ are positive during both periods on annual and seasonal scales (Tables S1–S3 in the Supplement). This indicates LW radiative energy gain at TOA and the surface, but both weaken after the AA reduction. It should be noted that the annual spatial pattern is driven by the JJA pattern in both cases (SW and LW). After the AA decrease, changes in SW and LW ΔF_coupled are generally more prominent over Central and Eastern Europe on annual and summer scales.

3.5 Temporal evolution of temperature changes – CMIP6 vs. Observations

The evolution of TAS in the MED region throughout the historical period imposed by AAs (Fig. 5) follows the evolution of transient ERF and ΔF_coupled. The annual mean TAS exhibits a negative peak in the late 1970s to early 1980s, which is stronger in JJA compared to DJF and shows an increasing trend towards more positive values by the end of the historical period.

Figure 5 Timeseries of the AA-induced near-surface air temperature (TAS) evolution over the Mediterranean region (10° W–40° E, 30–45° N) for the historical period (1850–2014) on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). TAS changes were calculated from the difference historical minus hist-piAer. The lines represent the 10 year moving average, while shading indicates the inter-model variability expressed as one standard deviation.

To gain confidence in the AA-induced TAS changes presented in Fig. 5, the ability of CMIP6 ESMs to capture the historical trends of TAS in the MED is examined. The TAS evolution over land from the historical simulation is compared to those from three gridded observational datasets (Fig. 6). On annual and summer scales, observations and the historical simulation agree on the evolution of land TAS in MED, with a cooling trend after the 1940s that peaks in the mid-1970s, followed by a warming trend afterwards. The absolute temperature levels in historical simulation are in good agreement with the University of Delaware and CRU TS datasets, especially in JJA, whereas there is a distinct warm bias in JJA and cold bias in DJF with respect to the Berkley Earth dataset.

Figure 6 Timeseries of land TAS evolution over the Mediterranean region (10° W–40° E, 30–45° N) for the historical period (1850–2014) as captured by CMIP6 ESMs and observations on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). The lines represent 10 year moving averages. In the cases of the CMIP6 experiments, shading indicates the inter-model variability expressed as one standard deviation.

On an annual scale over land, the three observational datasets show a warming ranging from 1.14–1.20 °C between the periods 2005–2014 and 1970–1979. The ensemble from the historical simulations indicates a relatively similar warming of 1.34 °C with a small overestimation in the range of 12 %–18 % (Table 7). The ensemble from the hist-piAer simulations shows a smaller warming of 0.68 °C, thus indicating that AA reductions contribute roughly 0.66 °C (49 %) to the warming in transient historical simulations between 2005–2014 and 1970–1979.

Table 7 Annual and seasonal means for near-surface air temperature (TAS; in °C) over the Mediterranean region for the periods 1970–1979, 2005–2014 and their difference (2005–2014 minus 1970–1979). Values are presented for the three observational datasets and the multi-model ensembles for the CMIP6 historical and hist-piAer experiments. Only gridpoints over land were considered here.

In the summer, the three observational datasets show a stronger land warming in the range of 1.68–1.81 °C between 2005–2014 and 1970–1979, and the ensemble from the historical simulations indicates a similar warming of 1.67 °C with a slight underestimation in the range of 1 %–8 %. The ensemble from the hist-piAer simulations shows a smaller warming of 1.02 °C, which indicates that reductions in AAs contribute roughly 0.65 °C (39 %) to the warming in CMIP6 historical simulations between the periods 2005–2014 and 1970–1979.

In the winter, air temperature over land evolves differently in the historical simulation compared to observations. While ESMs exhibit a similar, yet weaker cooling–warming pattern which is also seen in JJA, observational datasets show a general temperature increase since the 1970s. For example, in DJF over land, the three observational datasets show a slight warming in the range of 0.21–0.34 °C between 2005–2014 and 1970–1979. However, the ensemble from the historical simulations indicates an overestimated warming of 1.19 °C, with hist-piAer simulations revealing that AA decrease contributes roughly 0.71 °C (60 %) to this warming.

CMIP6 models produce a distinct cooling during the period 1960–1990 in DJF followed by an excessive warming trend from around 1995 onwards, which can arguably be attributed to aerosol concentration biases in the ESMs (Zhang et al., 2021b). In addition, the relatively weaker radiative effects from AAs in DJF than in JJA as well as atmospheric circulation, which is generally more disturbed in DJF than in JJA in the midlatitudes (Akritidis et al., 2021), influence the differences in DJF land temperature evolution between ESMs and observations. It has to be highlighted here that all three observational datasets used here show significant discrepancies between each other in representing the seasonal land TAS evolution, especially in DJF. In some cases, these discrepancies are larger between the observational datasets than between the observational data and the CMIP6 historical ensemble.

3.6 Spatial pattern of horizontal and vertical temperature changes based on CMIP6

Changes in the spatial pattern of TAS due to AAs as simulated by CMIP6 models are presented in Fig. 7. AAs cause a surface cooling relative to the pre-industrial era, which is most intense during 1970–1979 with an annual mean of $-\mathrm{1.49}\pm \mathrm{0.52}\mathrm{°}\mathrm{C}$ over the MED. The central and eastern parts of Europe experience the largest TAS decrease on both annual and seasonal scales. The surface cooling in the MED for the period 1970–1979 is more enhanced in the summer ($-\mathrm{1.89}\pm \mathrm{0.52}\mathrm{°}\mathrm{C}$) than in the winter ($-\mathrm{1.14}\pm \mathrm{0.70}\mathrm{°}\mathrm{C}$). The same applies to both Europe and Northern Africa, a feature that is robust among ESMs. After AA emission reduction, towards the period 2005–2014, mainly Central Europe and the Eastern MED undergo a surface warming both annually and seasonally (as can be noted in Fig. 7 from the difference between the periods 2005–2014 and 1970–1979). Specifically, in the MED the mean increase in TAS is of the same magnitude (around 0.7 °C) on an annual scale and during the winter and summer months. However, the annual change is robust in most of the MED, in agreement with the findings of Urdiales-Flores et al. (2023), whereas in DJF it is not statistically significant (despite the agreement on the sign of change between ESMs). In JJA, the robust results are confined in the Eastern MED, which is highlighted as a hotspot of summer warming due to AA reduction in all ESMs used in this work. The largest surface warming from 1970–1979 to 2005–2014 is produced by the EC-Earth3-AerChem and the UKESM1-0-LL models on both annual and seasonal scales, while the smallest by the BCC-ESM1 model.

Figure 7 Changes in TAS due to anthropogenic aerosols over the Mediterranean (calculated from the difference historical minus hist-piAer). Spatial patterns are presented for the periods 1970–1979 (left column), 2005–2014 (middle column) and their difference (2005–2014 minus 1970–1979; right column) on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). The multi-model ensemble means over the Mediterranean region (shown as a box) are shown along with the inter-model variability (one standard deviation) at the top of each panel. Colored areas without markings indicate robust changes, while hatched ($/$ and cross-hatched (×) areas indicate non-robust changes and conflicting signals, respectively.

The magnitude and seasonal asymmetry of the warming simulated by CMIP6 ESMs are consistent with previous observational and modeling studies reporting accelerated warming over the Euro-Mediterranean region following the decline of anthropogenic aerosol emissions, particularly during summer months (e.g., Dong et al., 2017; Nabat et al., 2014; Urdiales-Flores et al., 2023).

The warming caused by the decrease in AA emissions and concentrations is also apparent at higher levels in the troposphere. In Fig. 8 (9) the zonal (meridional) temperature means from the surface to 50 mb are presented considering only the MED longitudes (latitudes). The warming extends up to 200 hPa over the MED region (defined by the vertical dashed lines) as part of the Northern Hemisphere (NH) extratropical warming from 1970–1979 to 2005–2014 (Fig. 8). In the summer, the temperature increase affects a larger part of the NH troposphere, reaching up to 90° N, whereas it is more geographically confined in the winter. The epicenter of the annual mean warming is located at the MED longitudes, with an exception in JJA, during which the warming extends to around 60° E (Fig. 9). The magnitude of the warming decreases with height, except in the case of the NH midlatitudes in DJF (Fig. 8f).

Figure 8 Changes in zonal mean temperature due to AAs (calculated from the difference historical minus hist-piAer) for the MED longitudes (10° W–40° E). Cross sections (pressure level over latitude) are presented for the periods 1970–1979 (left column), 2005–2014 (middle column) and their difference (2005–2014 minus 1970–1979; right column) on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). Dashed vertical lines indicate the latitudinal boundaries of the MED region.

Figure 9 Changes in meridional mean temperature due to AAs (calculated from the difference historical minus hist-piAer) for the MED latitudes (30–45° N). Cross sections (pressure level over longitude) are presented for the periods 1970–1979 (left column), 2005–2014 (middle column) and their difference (2005–2014 minus 1970–1979; right column) on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). Dashed vertical lines indicate the longitudinal boundaries of the MED region.

3.7 Spatial pattern of circulation changes

Based on the results presented above, AA-induced TAS spatial changes cannot be fully explained by radiative changes alone. ERF_TOA and $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ changes between the periods 1970–1979 and 2005–2014 are not always spatially collocated with the respective TAS changes in the Euro-Mediterranean region (20–65° N, 25° W–55° E). This can be shown by calculating their spatial correlation coefficients on an annual scale (0.65 and 0.55, respectively), in DJF (0.52 and 0.49, respectively), and in JJA (0.69 and 0.60, respectively). This reveals that apart from the radiative effects and aerosol–cloud interactions (Kalisoras et al., 2024), atmospheric circulation changes induced by AAs may modify horizontal and vertical heat transport, thus playing an additional role in the redistribution of TAS changes.

To examine the thermal transfer through atmospheric circulation due to AAs, circulation changes at 850 hPa and sea-level pressure (SLP) relative to the pre-industrial period and due to decreases in AAs after around 1980 are shown in Fig. 10 and Fig. S4 in the Supplement, respectively.

Figure 10 Changes in geopotential height and winds at 850 hPa due to anthropogenic aerosols over the Mediterranean (calculated from the difference historical minus hist-piAer). Spatial patterns are presented for the periods 1970–1979 (left column), 2005–2014 (middle column) and their difference (2005–2014 minus 1970–1979; right column) on an annual scale (top row) and in boreal winter (DJF; middle row) and summer (JJA; bottom row). The multi-model ensemble means over the Mediterranean region (shown as a box) are shown along with the inter-model variability (one standard deviation) at the top of each panel. Colored areas without markings indicate robust changes, while hatched ($/$) and cross-hatched (×) areas indicate non-robust changes and conflicting signals, respectively.

On an annual scale, AA reduction (Figs. 10c and S4c) causes a cyclonic anomaly over the North Atlantic and Western Europe, which extends towards the central and eastern parts of Europe. This atmospheric circulation anomaly at 850 hPa and SLP is responsible for an enhancement of the westerly/southwesterly flow over the MED latitudes, along with an inflow from the northwestern parts of Africa towards Southern Europe. This indicates a warm advection anomaly which could contribute to the annual local TAS maximum over Southeastern Europe and Eastern MED.

During the winter, an anticyclonic anomaly (relative to the pre-industrial era) over Scandinavia in 1970–1979 (Fig. 10d), which could induce a cold advection anomaly towards Central and Southeastern Europe, strongly weakens or even disappears in the period 2005–2014 (Fig. 10e; also seen in Fig. S4e). Consequently, after the decrease in AAs from 1970–1979 to 2005–2014, a cyclonic anomaly extending from the North Atlantic to Eastern/Northeastern Europe (Figs. 10f and S4f) intensifies the westerly/southwesterly flow over the European midlatitudes (as well as the European part of the MED) and could possibly contribute to the warming pattern over Central and Eastern Europe, along with the warming due to radiative effects.

On the other hand, in JJA there are cyclonic anomalies (relative to the pre-industrial period) caused by AAs that extend from the North Atlantic to Northwestern Africa and Southeastern Europe during both 1970–1979 and 2005–2014. These anomalies could impose a southerly warm advection over the Eastern MED that may partly mask the radiative cooling effect of AAs in the region for both periods (Fig. 10g and h, respectively). The circulation anomalies, resulting from the difference between 2005–2014 and 1970–1979 following the reduction of AAs, are generally weak (Fig. 10i). It could be mentioned that the emergence of a non-robust cyclonic anomaly over the Eastern MED and the Middle East (Figs. 10i and S4i) may contribute to the summer warming amplification over the Eastern MED possibly through a weakening of the Etesian wind system (Dafka et al., 2016) and a warm advection anomaly from Northeastern Africa. However, it should be noted that ESMs show large differences when simulating circulation changes attributed to AA reduction, as the large inter-model variability and the lack of statistical significance suggest.

4 Conclusions

This work quantifies the contribution of global AA emission changes to recent Mediterranean warming based on simulations from eight CMIP6 ESMs including historical and hist-piAer coupled-ocean experiments, as well as histSST and histSST-piAer experiments with prescribed SSTs and sea ice. Reductions in AA emissions due to air pollution control policies cause a decrease in AOD over the MED which is mainly driven by the sulfate emissions decrease. The negative peaks in ERF and ΔF_coupled (either at TOA or surface) coincide with the positive AOD peak around the 1980s over the Euro-Mediterranean region, and afterwards a positive trend follows until the second decade of the 21st century, especially on annual scale and JJA. The AA-induced changes in ERF and ΔF_coupled for the peak aerosol period (1970–1979) and the near-present period (2005–2014), are robust over the northern half of the MED and are largest at the surface than at TOA on annual scale and JJA. The JJA patterns of ERF and ΔF_coupled dominate their respective annual spatial patterns in all cases. Changes in the spatial patterns of ERF and ΔF_coupled are dominated by their SW component on both annual and seasonal scales. The midlatitudes of Europe undergo the largest and most robust change in $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$, with JJA changes playing a pivotal role in the changes seen on an annual scale. The AA-induced changes in $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ are smaller in magnitude than ERF_TOA. This can arguably be attributed to the partial adjustment of the climate system through fast and incomplete slow responses after changes in land and ocean surface temperatures.

The evolution of AA-induced TAS changes is similar to that of ERF and ΔF_coupled, with a decreasing trend until approximately 1980 followed by an increasing trend toward the near present, which is strongest during JJA and robust over the central and eastern parts of the MED. The AA-induced temperature changes represent transient climate responses (fast and incomplete slow responses) under evolving historical conditions, which are far from being in an equilibrium state because the ocean does not have enough time to fully equilibrate to the evolving historical forcings.

On annual scale as well as during JJA, observations and the multi-model ensemble agree on the evolution of land TAS in MED, with a cooling trend after the 1940s that peaked in the mid-1970s and a warming trend afterwards. Thus, CMIP6 models are capable of capturing the enhanced MED surface warming due to AA reductions and generally agree on the temperature range with observations. From the peak aerosol period 1970–1979 towards the recent decade 2005–2014, the ensemble of CMIP6 simulations indicate a warming of 1.34 °C over land in the MED region on an annual scale and 1.67 °C during JJA, in relatively good agreement with three observational datasets. The model simulations also indicate that AA reductions between the two decades contribute to this warming roughly by 0.66 °C (49 %) on an annual scale and 0.65 °C (39 %) in JJA. In DJF, there are differences between ESMs and observations with regard to the temperature evolution. While ESMs exhibit a weak cooling–warming pattern before and after the early 1980s, observational datasets show a general monotonic temperature increase. Furthermore, the ensemble of the historical model simulations indicates strongly overestimated warming for DJF (from 1970–1979 towards the recent decade 2005–2014) in comparison to the observational datasets. This could arguably be attributed to overestimated aerosol concentrations in ESMs responsible for the excessive cooling during the period 1960–1990 as well as the different climate model responses to the relatively weaker AA radiative effects and more disturbed midlatitude atmospheric circulation conditions in DJF than in JJA. A warming due to AA decrease is also apparent at higher levels in the troposphere over the MED region, up to around 200 hPa.

ERF_TOA and $\mathrm{\Delta }{F}_{\text{TOA}}^{\text{coupled}}$ changes (from 1970–1979 towards the recent decade 2005–2014) are not always spatially collocated with the respective TAS changes in the Euro-Mediterranean region, indicating that apart from the radiative effects and aerosol–cloud interactions, atmospheric circulation changes induced by AAs may modify horizontal and vertical heat transport, thus playing an additional role in the redistribution of regional TAS changes and being an important factor in the MED warming.

Cyclonic anomalies are induced by the AA decrease both at sea-level and at 850 hPa on an annual scale and during DJF. These anomalies extend from the North Atlantic to Central and Eastern Europe and could be responsible for warm advection from Northern Africa. During JJA, there is a southerly warm advection towards Southeastern Europe from the northern parts of Africa, which may slightly mitigate the radiative cooling from AAs relative to the pre-industrial era. When AAs are decreased, JJA cyclonic anomalies emerge over the Eastern MED that could contribute to the amplification of the warming over the Eastern MED and Southeastern Europe through weakening of the Etesian wind system and warm advection anomaly from Northeastern Africa. Nevertheless, ESMs show large differences when simulating circulation changes attributed to AA reduction, as revealed by the large inter-model variability and the lack of statistical significance.

Overall, our results demonstrate that recent past acceleration of Mediterranean warming is strongly linked to global aerosol reductions, which have unmasked additional greenhouse gas–driven warming. However, due to the complex nature of aerosol–radiation–cloud interactions resolved in ESMs, which affect the amount of radiation reaching the surface and being trapped within the atmosphere, the robust quantification of AAs contribution in the MED surface warming demands further work, especially when considering the lack of model agreement regarding AA-induced changes in atmospheric circulation that impact heat transfer and temperature anomalies. The design of the AerChemMIP historical experiments used here does not allow attribution of regional AA emission changes to the MED climate. Furthermore, climate model parameterizations limit their ability to simulate certain climate feedbacks, sometimes resulting in large discrepancies between different ESMs.

Appendix A

In this section, the eight CMIP6 ESMs used in this work are described, focusing on their atmospheric general circulation and atmospheric chemistry models.

The Beijing Climate Center Earth System Model version 1 (BCC-ESM1) (Wu et al., 2020; Zhang et al., 2021a) uses the BCC Atmospheric General Model version 3 (BCC-AGCM3-Chem) (Wu et al., 2019) with interactive atmospheric chemistry. The 63 out of 66 chemical species considered in BCC-AGCM3-Chem are the same as those in the “standard version” of the Model for Ozone and Related chemical Tracers version 2 (MOZART2) (Horowitz et al., 2003). Sulfates are considered as prognostic and their treatment follows the Community Atmosphere Model (CAM) version 4 interactive atmospheric chemistry scheme (CAM4-chem) (Lamarque et al., 2012). Two types of organic and black carbon are treated in BCC-AGCM3-Chem, one water-insoluble and one water-soluble tracer. Sea salt and dust aerosols are classified into four size bins each. The mass mixing ratios of bulk aerosols are treated as prognostic variables following the NCAR Community Atmosphere Model (CAM3) (Collins et al., 2004), directly affecting the radiative transfer in the atmosphere. BCC-AGCM3-Chem uses prognostic aerosol masses to estimate the liquid cloud droplet number concentration and considers indirect effects of aerosols (Wu et al., 2019).

The CNRM-ESM2-1 model (Michou et al., 2020; Séférian et al., 2019) developed by CNRM-CERFACS incorporates the Reactive Processes Ruling the Ozone Budget in the Stratosphere Version 2 (REPRO-BUS-C_v2) atmospheric chemistry scheme for 63 chemistry species, in which chemical evolution is calculated only above the 560 hPa level (Michou et al., 2011; Morgenstern et al., 2017). Below that level, the concentrations of certain species are relaxed either toward the yearly evolving global mean abundances (Meinshausen et al., 2017) and for other species toward the 560 hPa value. The Tropospheric Aerosols for ClimaTe In CNRM (TACTIC_v2) interactive tropospheric aerosol scheme has been adapted for use in CNRM-ESM2-1 (Michou et al., 2015), which implements a sectional representation of sea-salt and desert dust aerosols (3 bins each), black carbon and organic matter (2 bins each to distinguish hydrophilic and hydrophobic particles), and sulfates (1 bin). Sulfate precursors evolve in sulfate aerosols with dependence on latitude. These 12 species of tropospheric aerosols are considered as prognostic variables and are affected by transport, dry deposition, and sub-cloud and in-cloud scavenging (Séférian et al., 2019). The cloud droplet number concentration is determined on the concentrations of sea salt, sulfate aerosols, and organic matter as in Menon et al. (2002), thus representing only the cloud albedo effect.

The EC-Earth3-AerChem model (van Noije et al., 2021) is an extended version of EC-Earth3 (Döscher et al., 2022) with an added component to simulate tropospheric aerosols, methane, ozone, and atmospheric chemistry. It uses the McRad radiation package to treat the radiative transfer in clouds (Morcrette et al., 2008). Cloud droplet formation is described by the Abdul-Razzak and Ghan (2000) aerosol activation scheme. Atmospheric chemistry and aerosols are simulated with the Tracer Model version 5 release 3.0 (TM5-mp 3.0), which includes sulfate, black carbon, organic aerosols, mineral dust, and sea salt described by the modal aerosol microphysical scheme M7 (Vignati et al., 2004). Four water-soluble modes (nucleation, Aitken, accumulation, and coarse) and three insoluble modes (Aitken, accumulation, and coarse) are considered in M7, with a lognormal size distribution that has a fixed geometric standard deviation describing each mode. M7 describes the evolution of the total mass and total particle number of each species for each mode, and explains water uptake, new particle formation, and aging through coalescence and condensation (Vignati et al., 2004).

The GFDL-ESM4 model (Dunne et al., 2020; Horowitz et al., 2020) consists of the Geophysical Fluid Dynamics Laboratory (GFDL) Atmosphere Model version 4.1 (AM4.1), which includes an interactive tropospheric and stratospheric gas-phase and aerosol chemistry scheme. It is similar to its previous version (AM4.0) (Zhao et al., 2018a, b), but with some updates: it treats nitrate and ammonium aerosols explicitly, the aging of black and organic carbon aerosols from hydrophobic to hydrophilic forms depends on the calculated concentrations of hydroxyl radical (Liu et al., 2011), and the oxidation of sulfur dioxide and dimethyl sulfide to produce sulfate aerosols is driven by the gas-phase oxidant concentrations and cloud pH simulated by the online chemistry scheme (Paulot et al., 2016). The AM4.1 has an online representation of gas-phase tropospheric and stratospheric chemistry (Horowitz et al., 2020). Aerosols are represented as bulk concentrations of nitrate, ammonium, sulfate, and hydrophilic and hydrophobic black and organic carbon particles, plus five bins each for sea salt and mineral dust aerosols. Sulfate and hydrophilic black carbon aerosols are regarded as being internally mixed by the radiation code (Larry W. Horowitz, personal communication, 2024).

The MPI-ESM-1-2-HAM model consists of the Max Planck Institute for Meteorology Earth System Model (MPI-ESM1.2), which uses the atmospheric general circulation model ECHAM6.3 (Mauritsen et al., 2019), coupled with the Hamburg Aerosol Model version 2.3 (HAM2.3) (Tegen et al., 2019). The ECHAM6.3-HAM2.3 (Lohmann and Neubauer, 2018; Mauritsen et al., 2019; Neubauer et al., 2019e; Tegen et al., 2019) uses a two-moment cloud microphysics scheme to examine aerosol–cloud interactions and to improve cloud simulation (Lohmann and Neubauer, 2018). Aerosol–cloud interactions are simulated in liquid, mixed-phase, and ice clouds (Neubauer et al., 2019e). In ECHAM6.3 the radiative effects of anthropogenic aerosols in the radiation calculation are treated using the MACv2-SP parameterization (Fiedler et al., 2017; Stevens et al., 2017). The aerosol microphysics module HAM (Stier et al., 2005; Zhang et al., 2012) calculates the evolution of aerosol particles considering the black and organic carbon, sulfate, sea salt, and mineral dust. In its default version, HAM simulates the aerosol size spectrum using the M7 (Vignati et al., 2004) scheme (Tegen et al., 2019). The updated version (HAM2.3) uses the scheme by Abdul-Razzak and Ghan (2000) for aerosol activation, along with other modifications for aerosol–cloud interactions, and updated emissions of aerosols and aerosol precursors from anthropogenic and natural sources (Neubauer et al., 2019e; Tegen et al., 2019).

The MRI-ESM2.0 model (Kawai et al., 2019; Oshima et al., 2020; Yukimoto et al., 2019c) consists of the atmospheric general circulation model with land processes MRI-AGCM3.5, coupled with the MRI Chemistry Climate Model version 2.1 (MRI-CCM2.1) atmospheric chemistry model, and the Model of Aerosol Species in the Global Atmosphere mark-2 revision 4-climate (MASINGAR mk-2r4c) aerosol model (Kawai et al., 2019; Yukimoto et al., 2019c). The MRI-CCM2.1 model simulates the evolution and distribution of ozone and other trace gases in the troposphere and middle atmosphere (Yukimoto et al., 2019c). The MASINGAR mk-2r4c model treats physical and chemical processes of atmospheric aerosols and includes black and organic carbon, non-sea-salt sulfate, mineral dust, sea salt, and aerosol precursor gases. Sea salt and mineral dust size distributions are divided into 10 discrete bins, whereas lognormal size distributions represent the sizes of the other aerosols (Oshima et al., 2020). External mixing is assumed for all aerosol species in this model, except for the radiation process in MRI-AGCM3.5, in which hydrophilic BC is assumed to be internally mixed with sulfate (Oshima et al., 2020; Yukimoto et al., 2019c). The conversion rate of hydrophobic to hydrophilic black carbon is based on the Oshima and Koike (2013) parameterization and depends on the rate at which condensable materials cover hydrophobic black carbon (Oshima et al., 2020). The aerosol activation into cloud droplets is based on the Abdul-Razzak et al. (1998), Abdul-Razzak and Ghan (2000), and Takemura et al. (2005) parameterizations.

The NorESM2-LM model (Seland et al., 2020) is the “low-resolution” version of the Norwegian Earth System Model version 2 (NorESM2), which is based on the Community Earth System Model (CESM2.1) (Danabasoglu et al., 2020) but with several modifications. NorESM2 employs the CAM6-Nor atmosphere model, which uses modified (compared to CESM2.1) parameterization schemes for aerosols and aerosol–radiation–cloud interactions (Kirkevåg et al., 2013, 2018), and the OsloAero6 atmospheric aerosol module (Kirkevåg et al., 2018), which describes the formation and evolution of sulfates, black carbon, organic matter, mineral dust, sea salt, and secondary organic aerosols. The major difference between OsloAero6 and other aerosol modules is the division of tracers into “background” and “process” tracers. The “background” tracers are mainly primary emitted particles that form lognormal modes and contribute to the aerosol number concentration. The “process” tracers change the chemical composition and shape of the initially lognormal background modes (Kirkevåg et al., 2018).

The UKESM1 model (Mulcahy et al., 2020; Sellar et al., 2019, 2020), with HadGEM3-GC3.1 (Kuhlbrodt et al., 2018; Williams et al., 2018) as its physical core, uses the UK Chemistry and Aerosol (UKCA) interactive stratosphere–troposphere chemistry scheme (UKCA StratTrop) (Archibald et al., 2020) and the GLOMAP microphysical aerosol scheme (Mann et al., 2010, 2012). The UKCA StratTrop model employs 84 chemical species and represents the chemistry of 81 of these (oxygen, nitrogen, and carbon dioxide are not treated as chemically active species) (Archibald et al., 2020). GLOMAP is a two-moment modal aerosol microphysics scheme which treats the evolution, the sources and sinks of sulfates, black carbon, organic matter, and sea salt across five lognormal size modes. Five aerosol size modes are represented in GLOMAP (nucleation, soluble Aitken, insoluble Aitken, accumulation, and coarse). Mineral dust, on the other hand, is simulated independently using the CLASSIC dust scheme (Bellouin et al., 2011) and can be regarded to be externally mixed with the GLOMAP aerosols (Mulcahy et al., 2020).