Attribution of the unprecedented heat event of August 2023 in Barcelona (Spain) to observed and projected global warming

Urban Climate 56 (2024) 102019 Available online 24 June 2024 2212-0955/© 2024 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/). Attribution of the unprecedented heat event of August 2023 in Barcelona (Spain) to observed and projected global warming Marc Lemus-Canovas a , b , * , Eduard Montesinos-Ciur´ o c , Tania Cearreta-Innocenti d , Roberto Serrano-Notivoli e , Dominic Roy´ e f , g a CRETUS, Non-linear Physics Group, Universidade de Santiago de Compostela, Galicia, Spain b Center for Climate Change and Transformation, Eurac Research, Viale Druso 1, 39100 Bolzano, Italy c Vice-Rectorate for Research, University of Barcelona, Gran Via de les Corts Catalanes 585, 08007 Barcelona, Spain d Department of Geography, University of Barcelona, Montalegre 6, 08001, Barcelona, Spain e Department of Geography and Regional Planning, Environmental Sciences Institute, University of Zaragoza, Pedro Cerbuna 12, 50009, Zaragoza, Spain f Climate Research Foundation (FIC), 28013 Madrid, Spain g Valencian International University, School of Engineering, Science and Technology - Robotics and Industrials, 46002, Valencia, Spain ARTICLE INFO Keywords: Extreme event attribution Heatwave Urban climate Global warming Barcelona ABSTRACT The study analyses observed and numerical simulations of daily maximum and minimum temperature from 1920 onwards and to investigate the unprecedented heat event that occurred in 21–23 August 2023 in Barcelona. The historical changes in the intensity of such events, their expected future changes under scenarios of +1.5 ◦C, +2 ◦C, and +3 ◦C, and the future exposure of populations to such kind of events are examined using the flow analogues approach. The findings indicate a significant increase in observed temperatures for similar heatwaves to those occurred in August 2023. The study also emphasises the impact of global warming on the intensification of heat events over the impact of urbanization. Additionally, after examining the role of natural variability in temperature changes, we concluded that global warming is the primary factor driving the increase in heatwave intensity. In terms of the frequency of such events, we found that extreme heat events, such as the August 2023 heatwave, will become 2 and 5 times more likely with a global summer warming of 2 ◦C and 3 ◦C, respectively. This will expose a large portion of the population to dangerous heat levels highlighting the importance of limiting global warming to 1.5 ◦C to mitigate the impacts on urban populations. 1. Introduction The intensity and impact of heatwaves are escalating around the globe (Perkins-Kirkpatrick and Lewis, 2020), presenting significant challenges within the context of already observed and projected global warming. In particular, urban areas are even more vulnerable to heat extremes due to the high population densities and the rapid urban growth they are experiencing (Liu et al., 2022; Lenton et al., 2023; Masselot et al., 2023). Moreover, the urban heat island (UHI) effect further intensifies these conditions, making cities warmer than their surrounding rural areas (Oke, 1973). On the other hand, the Mediterranean region, recognized as a climate change hotspot (Cramer et al., 2018; Lionello and Scarascia, 2018), is warming at a faster rate than the global average, with the most * Corresponding author at: CRETUS, Non-linear Physics Group, Universidade de Santiago de Compostela, Galicia, Spain E-mail address: [email protected] (M. Lemus-Canovas). Contents lists available at ScienceDirect Urban Climate journal homepage: www.elsevier.com/locate/uclim https://doi.org/10.1016/j.uclim.2024.102019 Received 14 March 2024; Received in revised form 27 May 2024; Accepted 14 June 2024 Urban Climate 56 (2024) 102019 2 pronounced increases during the summer months (Cos et al., 2022). This accelerated warming is marked by a rising frequency of heatwave patterns impacting southern Europe (Fischer and Sch¨ ar, 2010; Serrano-Notivoli et al., 2022) alongside prolonged droughts (Faranda et al., 2023; Garrido-Perez et al., 2024) and increased forest fires (Turco et al., 2018), among other climate related risks. As documented in recent studies of heatwave events in other regions of the world (Barriopedro et al., 2020; Xu et al., 2021; White et al., 2023), extreme temperatures are linked to a significant deviation from the expected large-scale circulation pattern. In southwestern Europe, where Barcelona is located, high temperature episodes or heat waves are usually driven by subtropical high-pressure systems (García-Herrera et al., 2005), commonly called subtropical ridges. These manifest as relatively narrow bands of positive geopotential height anomalies extending from subtropical latitudes towards southern Europe (Sousa et al., 2018). Elevated geopotential values sustained over time favor subsidence, hence strong adiabatic warming, and also an increase in surface radiative fluxes (Zschenderlein et al., 2019; Domeisen et al., 2023). In addition to these atmospheric dynamics, the observed global warming is the primary driver behind the consistent upward trend in urban temperatures (Liu et al., 2022). Notably, Mediterranean cities have observed sharp increases in both minimum (nighttime) and maximum (daytime) temperatures, with the rises in nighttime temperatures often surpassing those during the day (Founda et al., 2019). Urbanization, including urban growth, further contributes to local warming, enhancing the effects of global warming, as highlighted in various recent studies (Sun et al., 2016; Nogueira and Soares, 2019; Estrada and Perron, 2021; Liu et al., 2022). This interplay between urbanization and global warming is pivotal in determining the future climate of Mediterranean cities (Zhou et al., 2022). However, compared to the influence of global warming, the contribution of urbanization is generally minor (Liu et al., 2022; Zhou et al., 2022). Despite these findings, there remains a significant gap in understanding how global warming influences the intensity and frequency of specific heatwave events, or how such events can be attributed to anthropogenic climate change. Studies on severe heatwaves, such as those in Europe in 2003 (Stott et al., 2004), Russia in 2010 (Otto et al., 2012), and the US/Canada in 2021 (Philip et al., 2022; Zhang et al., 2023), have primarily focused on large-scale impacts, often extending beyond national borders, without offering insights tailored to local contexts. Therefore, conducting attribution studies at the urban scale could be particularly valuable for designing adaptation and mitigation strategies at the municipal level, especially given the significant impacts of urban heat on human health (Mitchell et al., 2016; Vicedo-Cabrera et al., 2021; Ballester et al., 2023; Iungman et al., 2023; Huang et al., 2023). Furthermore, accurately understanding the observed and projected frequency of specific events is crucial for quantifying how much of the population (Wang et al., 2023) is likely to be exposed to dangerous heatwaves. However, conducting research in urban areas often means relying on in-situ observations, which can present challenges due to various factors. Van Oldenborgh et al. (2022) highlight that local conditions, such as the design of the thermometer screen and its immediate environment, can affect the measurement of maximum temperatures. This underscores the complexity of gathering accurate climate data in urban settings, where numerous microclimatic variables can influence temperature (Peng et al., 2018; Yue et al., 2019) and, by extension, the assessment of heatwave risks and their potential impact on urban populations. Natural variability, particularly teleconnection patterns, have been studied for several decades since they play a significant role in influencing the temperature of Southern Europe. Studying such patterns is relevant when attributing extreme events to isolate them from the role of anthropogenic global warming. For instance, the El Ni˜ no-Southern Oscillation (ENSO) has been seen to play a partial role in influencing the climate variability of Europe, especially precipitation (Br¨ onnimann, 2007; Cagnazzo and Manzini, 2009). The role of ENSO in driving temperature variability in Western Europe is less clear since its impact is nonlinear and nonstationary (Br¨ onnimann et al., 2007; Martija-Díez et al., 2021). In Western Europe, Martija-Díez et al. (2021) noted that warm summer events are related to La Ni˜ na. Precisely, Faranda et al., 2023 used the ENSO to account for the role of natural variability when attributing the intensity of an extreme Mediterranean summer heatwave to anthropogenic climate change. On the other hand, the Atlantic Multidecadal Oscillation (AMO) is a natural mode of variability of the Northern Atlantic Sea Surface Temperature (SST), characterised by a main period of about 65–70 years (Schlesinger and Ramankutty, 1994; Delworth and Mann, 2000) also influencing temperature variability in Southern Europe. For instance, Zampieri et al. (2017) observed that during warm phases of the AMO, Southern Europe experienced increased temperatures, contributing to heatwaves and dry conditions. This study represents a significant step forward in comprehending the effects of global warming on Mediterranean urban zones, focusing particularly on the densely populated city of Barcelona (Spain), by conducting a thorough analysis of the record-breaking heat event of August 2023. This research not only highlights the exacerbated heat extremes due to global warming but also underscores the future risk of population exposure to similar heatwaves. This is particularly crucial for advocating adaptation strategies at the municipal level (Cuthbert et al., 2022; Goonesekera and Olazabal, 2022; Pietrapertosa et al., 2023). To this end, we address three major questions through our study: •How has the intensity of an extreme heat event, like the one experienced in August 2023, changed historically and how is it expected to change in the future, especially in scenarios of a +1.5 ◦C, +2 ◦C, and +3 ◦C warmer world? •How much will the frequency and likelihood of events, similar to the August 2023 heatwave, increase? •What portion of the population will be exposed to such events in the future? By answering these questions, our findings aim to draw societal attention to this pressing issue and serve as a valuable resource for policymakers in devising effective climate adaptation and mitigation strategies. The remainder of this paper is structured as follows: Section 2 discusses the materials and methods used and presents the study area. In section 3, we first present results on the local and synoptic characterisation of the August 2023 heat event. Afterwards we detail the contribution of global warming and natural variability to both observed and future intensification events like the 2023. Finally on this section, we examine the future frequency and probability of events as much or more extreme than the 2023 one, and the projected population exposed to such events. Finally, section M. Lemus-Canovas et al. Urban Climate 56 (2024) 102019 3 4 presents discussion and conclusions. 2. Materials and methods 2.1. Climate and urban setting The municipality of Barcelona is located in the northwest of the Mediterranean basin, between the arid and warm northern African climate and the humid and mild European climate (Cramer et al., 2018). The contrast between them is partly explained by the influence of the surrounding oceans, their interaction with the land surface, and the general atmospheric circulation characteristics in the mid-latitudes (Bo´ e and Terray, 2014). Barcelona presents a Mediterranean climate, with warm summers (average temperature in July and August of 24.5 ◦C) and mild winters (average temperature in January of 10 ◦C). The seasonal precipitation regimes present equinoctial maxima, with the absolute in autumn, while the summer is characterised by the typical period of Mediterranean aridity (De Luis et al., 2010). The city of Barcelona occupies a strategic geographical position on the Iberian Peninsula, boasting significant connections with the rest of Europe. With a population of 1,660,122 inhabitants (INE 2023), it is the second most populated municipality in Spain, also owning a high-density rate of 16,339 inhabitants/km 2 . It covers an area of approximately 100 km 2 , forming the centre of one of the densest metropolitan areas in Europe and the Western world (Fig. 1a). This fact is reflected in its urban morphology, characterised by the density of buildings and the lack of large open or green spaces (Fig. 1b). Barcelona is mostly urbanised in its flatter areas. Despite having hosted thousands of industrial buildings since the 19th century, following a process of industrial displacement in the mid-20th century, the most common land use today is residential, and those dedicated to service economy facilities. The high population density and intricate land cover pose significant challenges in terms of urban planning, mobility, and resource management. This situation can lead to congestion in urban infrastructure, making access to basic services difficult and increasing pressure on natural and environmental resources. 2.2. Data collection The Meteorological Service of Catalonia (SMC) directly supplied the observed time series of maximum (TX) and minimum (TN) temperatures for the Barcelona observatory, starting in 1914, with a daily temporal resolution. The SMC itself has conducted various quality controls on these data (Prohom et al., 2023). Moreover, to historically contextualize the summer of 2023, we used a series of reconstructed average monthly temperatures for Barcelona from 1780 to the present (Prohom et al., 2016), making it the oldest series available on the Iberian Peninsula to date. For the synoptic variables in this study, we analysed the geopotential height at 500 hPa (Z500) and the temperature at 850 hPa (T850) within the domain [10◦W–12◦E–34◦N–50◦N] from the ERA–5 reanalyses (Hersbach et al., 2020) and 20thCRv3. We averaged the ERA-5 hourly data on a daily scale from 1940 to 2023. We directly downloaded the 20thCRv3 data at a daily resolution from 1914 to 2015. We resampled both datasets to the same spatial resolution (1◦x1◦). Section 2.3 details how we combined these datasets into a single series from 1914 to 2023. Extending this time series was crucial to set three distinct periods for analyzing the warming escalation (see Section 2.5), which allowed for a more comprehensive examination of Fig. 1. Urban Structure and Local Climate Zones of Barcelona, Spain. (a) Orthoimage of Barcelona, with an inset map highlighting the city's location in southwestern Europe. (b) A colour-coded representation of Barcelona's local climate zones (LCZs), each colour corresponding to a different LCZ that characterizes the urban and natural environments within the city. M. Lemus-Canovas et al. Urban Climate 56 (2024) 102019 4 historical changes. Regarding numerical simulation data, this study employed Z500 and urban TX/TN for the historical forcing of CMIP5 from 1920 to 2005 and the radiative forcing of RCP 8.5 from 2006 to 2100, using 40 ensemble members of the Community Earth System Model version 1 (CESM1) Large Ensemble Project (Kay et al., 2015). Since historical simulations extend only to 2005, we bridged the historical simulation to 2023 by appending data from the RCP8.5 experiment for the period 2006–2023. We obtained these data with a daily temporal resolution and a spatial resolution of 1◦, matching the geographical domain of the reanalysis data. The nearest grid point to the coordinates of Barcelona was extracted to perform the analysis. To assess the percentage contributions of Northern Hemisphere warming and urban growth to the annual 3-day maximum temperatures (detailed in Section 2.4), we employed two different datasets. Firstly, we used the HadCRUT5 dataset (Morice et al., 2021), which contains monthly mean temperature anomalies for the Northern Hemisphere covering the period from 1850 to present, with anomalies computed using the 1961–1990 period as a reference. Secondly, we utilized the HISDAC dataset (Uhl et al., 2023), offering a series that tracks the number of buildings constructed in Barcelona from 1900 to 2020. To jointly analyse both datasets with the observed temperature time series we established a common period set in 1925–2010 (see section 2.4). We employed Eurostat's annual population projections at the NUTS3 level to ascertain the population's degree of exposure to extreme temperature events. These projections were constructed using deterministic “what-if” scenarios that model hypothetical developments in population size and structure, based on assumptions about fertility, mortality, and migration rates (Eurostat, 2024). The projection begins with the population as of January 1, 2019, applies mortality rates to calculate the expected number of deaths, and adjusts for immigration and emigration to estimate population changes. For the purposes of this study, we focused on the NUTS 3 region of Barcelona, which encompasses the entire population of the province of Barcelona, totaling 5,575,204 inhabitants as of 2019. To specifically estimate the population of the city of Barcelona, which stood at 1,650,358 inhabitants in 2019, we employed a stationary approximation. This approach assumes that the ratio between the city's population and that of the broader NUTS3 region remains constant over time. This assumption is based on the relationship observed at a specific point in time (2019), facilitating the extrapolation of the city's population into future scenarios under study. POPBCN,year =POPBCN NUTS3,year–(POPBCN NUTS3,2019–POPBCN,2019)(1) where POP BCN, year is the projected population of the city of Barcelona for a specific year; POP BCN NUTS3, year is the projected population for the NUTS3 region of Barcelona for that year; POP BCN NUTS3, 2019 is the observed population in 2019 for the NUTS3 region of Barcelona; POP BCN, 2019 is the observed population in 2019 for the city of Barcelona. Finally, and just for illustrative purposes in the study area section, we used an orthoimage provided by the Cartographic and Geological Institute of Catalonia (ICGC) captured by the Sentinel-2 satellites (winter 2019–2020) and without clouds at 10 m resolution. Additionally, we employed the global map of local climate zones (Demuzere et al., 2022) masked for the Barcelona administrative border, to introduce the main land uses with climate implications at the urban scale. The abovementioned administrative limits were retrieved from the Open Data BCN portal maintained by the Barcelona City Council. 2.3. Bias correction of ERA-5, 20thCR and CESM1 datasets Here, we separate between 1) the combination of ERA-5 and 20thCRv3 reanalysis for the creation of the Z500 and T850 series and 2) the bias correction of the CESM1 numerical simulations of Z500 and TX/TN. We combined the two reanalysis datasets, ERA-5 and 20thCRv3, using the methodology described in Faranda et al. (2023), aiming to eliminate systematic differences between them and thus prevent biased results in the analyses discussed in sections 2.4 and 2.5. To correct biases between 20thCRv3 and ERA5, we employed a linear rescaling technique over their common period (1940–2015). We designated ERA5 as the reference dataset and 20thCRv3 as the dataset to be rescaled. The recalibration formula we developed is as follows: Yʹ= μ ERA5 + (Y− μ 20thCRv3) • c(2) Y ′ refers to the adjusted version of original 20thCRv3 data (Y) to align with ERA-5 standards for the period 1914–1939. The terms μ ERA5 and μ 20thCRv3 represent the mean values of the ERA5 and 20thCRv3 datasets, respectively, calculated over the shared period from 1940 to 2015. The recalibration factor c is determined through a linear variance-based method, calculating it as the ratio of the standard deviations ( σ ) of ERA5 to 20thCRv3 for the same common period. This ratio is expressed as: c= σ ERA5/ σ 20thCRv3 (3) This approach applies a bias correction to the 20thCRv3 data for the period 1914–1939, which is not covered by ERA-5, thereby aligning it more closely with the ERA5 dataset. Nevertheless, when we provided results for both the reanalysis combined dataset and the CESM1 simulations, we employed 1920 as starting year in both series since it is the starting year available in CESM1 simulations we employed. For the bias correction of the CESM1 numerical simulations, we tailored our approaches to suit the variables Z500 and TX/ TN individually. To maintain consistency with the adjusted ERA-5/20thCRv3 dataset, we applied the same bias correction technique to Z500. This adjustment is grounded on the overlapping period between ERA-5 and the CESM1 historical simulations, spanning from 1940 to 2005. This ensures that the Z500 data from CESM1 align with the reanalysis data, enhancing the coherence of our climate analysis across different datasets. For the variables TX and TN, we adopted a different strategy, employing the non-parametric quantile mapping method. This approach involves aligning the empirical quantiles of the CESM1 distributions for TX and TN with those derived M. Lemus-Canovas et al. Urban Climate 56 (2024) 102019 5 from in-situ observations in Barcelona. Such a methodology has been widely applied to correct for bias between in-situ observations and models in the past (Thrasher et al., 2012; Tong et al., 2021; Enayati et al., 2021). Details of this methodology can be found in Gudmundsson et al. (2012), and its application was performed through the R package qmap (Gudmundsson, 2016). 2.4. Heat event detection and its characterisation To identify the event, we first calculated the 3-day running mean of the maximum (TXm3) and minimum daily (TNm3) temperature anomalies for Barcelona Fabra Observatory (hereafter, Barcelona) series. We then computed the annual absolute maximums for TXm3 (TNm3), resulting in TXx3 (TNx3). Notably, the all-time maximums for TXx3 (TNx3) occurred in the most recent year, 2023. Consequently, we defined the three dates corresponding to TXx3 (TNx3) in 2023 as the dates of the heat event under investigation. We computed these anomalies relative to their corresponding calendar dates for the period 1952–1981. This period was chosen after calculating a 30-year running mean of the monthly mean summer temperature anomalies (June–July–August) from 1780 to 2023, compared to the preindustrial period of 1850–1900, as shown in Fig. S1. This analysis established a thirty-year period (1952–1981) with an anomaly close to 0 ◦C (pre-industrial level) within the data available for the Barcelona series (1914–2023). For the synoptic characterisation of the event, we calculated temperature anomalies at 850 hPa and geopotential heights at 500 hPa for the period 1952–1981, as well as we computed the percentiles of these variables for the entire period available. We also estimated the degree of contribution from the hemispheric mean temperature and the changes in the number of buildings to the series of 10-year running mean TXx3 (TNx3) anomalies. We performed this calculation using a multiple linear regression (eq. 4) from the R package relaimpo (Groemping, 2007), which allows for the decomposition of the model's explained variance into nonnegative contributions. To avoid overprediction due to their increasing trends, we detrended both predictors before fitting the regression model. The purpose of this analysis was to determine the extent to which urban growth and hemispheric warming have influenced TXx3 (TNx3) values. TXx3 =β₀+β₁(NHtemp)+β₂(NBuildings) + ε (4) where TXx3 is the annual summer (June–July–August) maxima of the 3-day running mean of the daily maximum temperature anomaly; NHtemp is the hemispheric mean temperature anomaly series for summer, while NBuildings is the series of number of buildings constructed in the city. Note that for this calculation, the common period was 1920–2015 and eq.4 was also computed for TNx3 as a predictand. 2.5. Attribution analysis and population exposure to future warming To attribute the three-day heat event, our study used the flow analogues approach, as outlined by J´ ez´ equel et al. (2018). This method aims to infer the behavior of a surface variable, in a past (counterfactual) climate, as well as in a present or future (factual) climate. The key criterion for this inference is that the atmospheric circulation characteristics during these periods must be identical or similar to those observed during the target heat event. Specifically, this approach focuses on the atmospheric circulation patterns at 500 hPa associated with the event under study. By constraining the analysis to circulation patterns, particularly at the 500 hPa level, we can better understand how similar atmospheric conditions might influence surface temperatures in different climatic periods. To identify days or analogues similar to each of the three days of the heat event, we calculated the root mean square differences (RMSD) against the observed Z500 anomaly field for each event day across the specified synoptic domain (10◦W–12◦E–34◦N–50◦N). This process was undertaken for each day of the event, with the search confined to the boreal summer months (June, July, August) and excluding the event's occurrence year (2023). We conducted the analogue search over three historical periods: 1920–1949, 1952–1981, and 1993–2022. This approach allowed us to estimate the degree of thermal amplification of such events in the current era compared to the two earlier periods. For projecting the future amplification of maximum (TX) and minimum (TN) temperatures for an event akin to that of 2023, the search was not constrained by specific periods but was based on three warming thresholds: 1.5 ◦C, 2 ◦C, and 3 ◦C above pre-industrial levels. We calculated these warming levels following the methodology proposed by Hauser et al. (2022), which involves calculating a 20-year centered moving average of global temperature anomalies relative to the baseline period of 1850–1900. The onset of a 20-year period averaging above one of the specified warming thresholds marks the commencement of a global warming level. These future warming levels were computed for each of the 40 CESM1 simulations employed in our study, with variability in the starting and ending years of each warming level illustrated in Fig. S2. For each analogue search, we used the respective warming level periods defined for each CESM1 simulation to estimate the temperature increase of TX and TN in events like the one under study. Once the analogue days have been identified, we reconstruct the probability distribution of daily TX (maximum temperature) anomalies, conditioned on the circulation pattern of the 2023 event, using the method proposed by Zhang et al. (2023). This process unfolds as follows: 1) For each day, one of the top20 analogues –with lower RMSD– is randomly selected, thus obtaining a random flow-conditioned simulation of the 2023 event; 2) The average of the reconstructed TX anomalies of the event is calculated, thus obtaining the TXm3; 3) By repeating the above steps 1000 times, the probability distribution of the circulation-conditioned TXm3 anomalies of the 2023 heat event is obtained. For comparison, the distribution of TXm3 anomalies is randomly reconstructed from a series of days of the same length as the event, but in this case not conditional on circulation, thus providing a control simulation. In all cases, the anomalies were defined with respect to each date in the calendar for the period 1952–1981. To incorporate the effect of natural interdecadal variability into our analysis, we adopt the method proposed by Faranda et al. M. Lemus-Canovas et al. Urban Climate 56 (2024) 102019 6 (2023), using monthly indices derived from NOAA/ERSSTv5 data, as accessed through KNMI's Climate Explorer. We apply the El Ni˜ no South Oscillation index in its region 3.4 (hereafter referred to as ENSO), as defined by Huang et al. (2017), and calculate the AMO index following the method described by Trenberth and Shea (2006). A positive ENSO indicates El Ni˜ no conditions, whereas a negative value denotes La Ni˜ na. To evaluate the significance of differences between factual and counterfactual distributions, we conducted a two-sided Cram´ er-von Mises test at the 0.05 significance level in all instances. A p-value lower than 0.05 leads to the rejection of the null hypothesis (H =0), asserting that the two samples originate from identical distributions (Anderson, 1962). Thus, if we reject the null hypothesis that there is no difference between the two distributions across periods, we cannot discount that the thermodynamic or dynamical distinctions among the analogues might partially arise from these natural variability modes rather than anthropogenic forcing. Conversely, inability to reject the null hypothesis of equivalent distributions suggests that the observed changes in the analogues are attributable to human activities. Using the 2023 TXx3 (TNx3) as a threshold, we calculate the frequency as the number of times this threshold is exceeded over an average of 3 consecutive days each summer. By considering both climate and population projections, we calculate the population exposure as the number of summer events (i.e., the frequency) multiplied by the projected number of people in the city of Barcelona. Finally, the impacts of extreme heat events that are avoided at 1.5◦C compared with 2◦C warmer climates are investigated using the eq. (5) proposed by Li et al. (2018): AI =C2.0−C1.5 C2.0⋅100 (5) Fig. 2. (a)Evolution of maximum temperature (TX, salmon line) and minimum temperature (TN, blue line) for the summer 2023 in Barcelona. The top-center inset plot in (a) shows the 3-day running mean TX (TN) anomalies (relative to 1952–1981) for Barcelona in August. The segment colored in magenta indicates the event under study (21–23 of August) (b) Evolution of the anomalies of TXx3 (salmon line) and TNx3 (blue line) in Barcelona during 1914–2023. TXx3 (TNx3) is defined as the annual summer (June–August) maxima of the 3-day centered running mean of the daily maximum (minimum) temperature anomaly in Barcelona. The red (blue) dotted series show the 10-years running mean for TXx3 (TNx3). The red (blue) horizontal dashed lines indicates the TXx3 (TNx3) anomalies (21–23 August) for 2023. (c) 3-days running mean anomalies for geopotential height (shading) at 500 hPa and temperature at 850 hPa for the 21–23 August 2023. (d) The probability density functions (PDF) of the summer mean temperature anomalies (relative to 1850–1900) reconstructed time series (see Data and Methods) for Barcelona during 1780–2023. The black dashed vertical line indicates the summer mean temperature (TM) for 2023. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) M. Lemus-Canovas et al. Urban Climate 56 (2024) 102019 7 where AI is avoided impacts and C1.5 and C2.0 are the changes in frequency for the 1.5◦C and 2◦C warming climate levels. The avoided impacts of 1.5◦C warmer climates compared with 3◦C warmer climates are also calculated by Eq. (5), but C2.0 is replaced by C3.0, representing the changes in the 3◦C warmer world threshold. 3. Results 3.1. Observed characteristics of the 2023 heat event in Barcelona We first analyse the observed characteristics of the event, both in the in-situ records and in the synoptic fields at altitude. On 23 August a TX of 38.8 ◦C was recorded (Fig. 2a), this value was the second warmest in August since the beginning of the temperature series in 1914. Although the absolute record for this TX series is 39.8 ◦C in 1982, the daily mean (TM) and minimum (TN) temperature recorded on the 23rd of August of 2023 were the highest ever recorded for the last 110 years. Precisely, if we take the TXm3 and TNm3 anomalies for the 21–23 of August 2023, both are the warmest absolute values of the whole series (Fig. 2b), with mean anomalies for days 21–23 of +11.3 ◦C and +10.8 ◦C, respectively (Fig. 2a,b). Surface temperatures had their main driver at altitude, where anomalies at 500 hPa exceeded +20 dm and those at 850 hPa reached 12 ◦C in the north of the Iberian Peninsula, while they were close to +9 ◦C over the Barcelona vertical (Fig. 2c). All these anomalies exceeded the 99.5th percentile over large parts of western Europe (Fig. S3). Despite such temperature anomalies for the days from 21 to 23 and, overall, for the whole month of August 2023, the mean temperature of that summer was not the warmest on record since 1780 (Fig. 2d), but it is among the 10 warmest summers of the last 243 years, i.e. in the 96th percentile. 3.2. Influence of observed and projected warming on a heatwave-like 2023 3.2.1. Observed changes One of the primary goals of this research was to ascertain the extent to which global warming has intensified a heat event like that of August 2023. To this end, we reconstructed the 21–23 August event using analogous conditions from two past periods with minimal warming signals and a more recent period that exhibits a distinct global warming trend. As shown in Fig. 3, the simulations for the average temperature of three analogue days, conditioned by a circulation similar to that of the 2023 event, always show higher temperature anomalies than those of random analogue days. Normally, in the absence of a warming signal, the median of these non- Fig. 3. The distribution of the 3-day mean maximum (a-b, TXm3) and minimum (c-d, TNm3) temperature anomalies for flow analog events (analogs, colored) generated using daily 500 hPa geopotential height from 21 to 23 August and randomly picked days (control, grey) using an in-situ time series (a-c) and CESM1 simulations (b-d). The red dashed horizontal line indicates the maximum (minimum) temperature anomalies of the average of 21–23 August 2023 from the observed time series. Boxes show the 25th–75th percentile ranges and whiskers span the 1st–99th percentiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) M. Lemus-Canovas et al. Urban Climate 56 (2024) 102019 8 circulation-conditioned simulations aligns closely with a 0 ◦C anomaly, a pattern observed in the periods 1920–1949 and 1952–1981, for both observed data and CESM1 numerical simulations. However, in the recent period of 1993–2022, which is already influenced by global warming (IPCC, 2023), both circulation-conditioned and non-conditioned simulations exhibit a marked increase in the median temperature. This rise is evident across both datasets and variables. Specifically, for observed TXm3 (Fig. 2a), the general increase in thermal anomalies for a set of three random days is approximately +2.1 ◦C (Table 1), while CESM1 simulations (Fig. 3b) reveal a positive anomaly of about +1.4 ◦C (Table 1). Notably, for simulations conditioned on the circulation of the 2023 event in this latest period, the increase in temperature anomaly is significantly higher in the numerical simulations - with a median of the simulations around +7.5 ◦C - compared to the observed data simulations - with a median around +6.6 ◦C. This discrepancy suggests a potential bias in the CESM1 data results, despite prior bias corrections. For TNm3 results, both observed (Fig. 3c) and CESM1 (Fig. 3d), the outcomes are consistent with those mentioned for TXm3. However, the thermal anomaly in the 1993–2022 period simulations is about 0.9 ◦C higher for TX than TN in the observed estimation whereas it is around 1.4 ◦C in CESM1 outcomes, highlighting the differential impact of warming on daytime versus nighttime temperatures. To evaluate if natural variability, such as ENSO and AMO changes, contributed to the observed thermal increase in the most recent period, we analysed their impacts (Fig. 4). The distributions of both indices during the summer months show no significant shifts across the three study periods, allowing us to largely discount a direct impact of ENSO and AMO variability on the observed temperature increases. Specifically, the AMO index shows the most variation (Fig. 4a), transitioning from weakly positive values in the first period (1920–1949) to negative in the second (1952–1981), and returning to positive in the latest period (1993–2022). These fluctuations might be expected to influence the thermal intensity of Barcelona's heat events across these periods. However, as Fig. 3 illustrates, an increase in event intensity is only evident in the latest period, despite experiencing different AMO phases in the earlier periods, where the temperature impact of such heatwaves remained nearly constant. For ENSO, the shift is from neutral values in the initial two periods to slightly negative values in the most recent period (Fig. 4b), hinting at a minor influence on the exacerbation of recent heatwaves. Yet, these shifts in the ENSO index are not statistically significant when compared to earlier periods. Thus, our analysis concludes a non-significant role of teleconnections in the thermal intensification of recent heat events similar to the one in August 2023. Finally, to test whether urban growth in Barcelona influenced this observed warming, the contribution of hemispheric mean temperature was separated from the growth in the number of buildings in the same city (Fig. S4). Such an analysis showed that the weight of hemispheric thermal variability in driving temperature extremes in Barcelona is clearly higher (73%) than the weight of urban modification (15%). This reinforces the idea that the role of global warming is behind the recent warming of observed heat waves such as the one in 2023. 3.2.2. Projected changes In our research, we aimed not only to understand the observed role of warming but also to estimate the future amplification of such events in a warmer climate. To achieve this, we conducted a second experiment, searching for analogues in the RCP8.5 scenario from the CESM1 simulations across different warming levels: +1.5 ◦C, +2 ◦C, and +3 ◦C (Fig. 5a). Notably, in a future warmer climate by these degrees, with atmospheric circulation similar to the 2023 event, the exceedance of the observed temperature thresholds for both TXm3 and TNm3 will become increasingly common. A particularly striking finding is that the rate of increase for TNm3 values in such events will outpace that of TXm3. In a world warming by +3 ◦C, the increase for both variables averages +9.3 ◦C (Table 1), which represents an increase of +1.8 ◦C for TXm3 and +3.2 ◦C for TNm3 compared to the 1993–2022 period. This results in a heightened likelihood of experiencing events with TXm3 ≥11.3 ◦C or TNm3 ≥10.8 ◦C in a warmer climate (Fig. 5b). These probabilities and temperature increases are further elucidated in Table 1. From the analysis, it's clear that the probability of surpassing the TXx and TNx thresholds identified by the 2023 event was non-existent up until 1981. However, during the period 1993–2022, there was a slight increase in probability, albeit only for TXm3. This probability becomes significant for both variables upon reaching the +1.5 ◦C warming level, suggesting a likelihood of 5% or higher that these temperature thresholds will be exceeded. As the warming level Table 1 Summary of the observed and projected warming flow-conditioned to the event of 21–23 of August, showcasing temperature increases in terms of maximum (TX) and minimum (TN) values. The table presents observed data (OBS) and projections from the Community Earth System Model (CESM1) across three distinct periods: 1920–1949, 1952–1981, and 1993–2022, alongside warming levels of 1.5 ◦C, 2 ◦C, and 3 ◦C. Mean anomalies and standard deviations (Std) are provided, with anomalies relative to the 1850–1900 pre-industrial baseline in parentheses. Period/Warming level TX (◦C) TN (◦C) OBS CESM1 OBS CESM1 Mean Std Mean Std Mean Std Mean Std 1920–1949 +4.0 (+0.4) 0.1 (3.1) +5.7 (+0.4) 1.8 (3.3) +3.6 (+0.43) 0.1 (2.5) +4.0 (+0.4) 1.6 (2.5) 1952–1981 +4.5 (+0.1) 0.2 (3.0) +5.5 (+0.3) 1.6 (3.2) +4.1 (−0.0) 0.1 (2.3) +4.0 (+0.3) 1.5 (2.5) 1993–2022 +6.6 (+2.1) 0.5 (3.1) +7.5 (+1.4) 1.9 (3.4) +5.7 (+1.8) 0.3 (2.5) +6.1 (+1.3) 1.7 (2.7) 1.5 ◦C +8.6 (+2.9) 1.9 (3.4) +7.7 (+2.8) 1.7 (2.9) 2 ◦C +8.7 (+3.9) 2.0 (3.6) +8.3 (+4.0) 1.8 (3.2) 3 ◦C +9.3 (+5.4) 2.1 (4.1) +9.3 (+5.7) 1.8 (3.65) M. Lemus-Canovas et al. Urban Climate 56 (2024) 102019 9 reaches +3 ◦C, the chances of breaking these records again ascend to 24% for TN and 19% for TX. Consequently, such events, which were either non-existent (1920–1981) or highly unlikely (1993–2022), will transition to being unusual yet probable occurrences in a climate warmed by +1.5 ◦C or more. 3.2.3. Future frequency and exposition of heat extremes like 2023 The recorded observations of the event in August 2023 indicated anomalies around +11.3 ◦C for TXm3 and +10.8 ◦C for TNm3. Fig. 4. Analysis of the interannual and interdecadal variability. Violin plots for past periods (blue; 1920–1949 and 1952–1981) and most recent period (red; 1993–2022) periods for AMO (a) and ENSO (b) values corresponding to the analogue's months (June–July-August). Black lines in violin plots indicate mean values. Labels in violin plots report the results H of the two-sided Cram´ er-von Mises test at the 0.05 significance level with the corresponding p-values (see Data and methods). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 5. a) The distribution of the mean maximum (red box) and minimum (yellow box) temperature anomalies (unit: ◦C) for flow analog events in given levels of projected global warming (1.5 ◦C, 2 ◦C and 3 ◦C) by means of CESM1. The dashed red (yellow) line corresponds to the TXm3 (TNm3) value reached in the event of 21–23 August 2023. Boxes show the 25th–75th percentile ranges and whiskers span the 1st–99th percentiles. b) probability of exceeding TXx3 (red) 2023 and TNx3 (yellow) in the historical scenario periods (1920–1949; 1952–1981 and 1993–2022) and in the projected warming levels (1.5,2, and 3 ◦C) in CESM1 ensemble. The inner point indicates the mean values of the simulations, whereas whiskers span the 1st–99th percentiles. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) M. Lemus-Canovas et al.