Detection and attribution of heat waves with the Multivariate Autoencoder Flow-Analogue Method (MvAE-AM)

Contents lists available at ScienceDirect Atmospheric Research journal homepage: www.elsevier.com/locate/atmosres Detection and attribution of heat waves with the Multivariate Autoencoder Flow-Analogue Method (MvAE-AM) Cosmin M. Marina a,∗, Jorge Pérez-Aracil a, Ronan McAdam d, Eugenio Lorente-Ramos a, Niklas Luther b, Eduardo Zorita c, Enrico Scoccimarro d, Jürg Luterbacher b, Elena Xoplaki b,d, Sancho Salcedo-Sanz a aDepartment of Signal Processing and Communications, Universidad de Alcalá, 28805 Alcalá de Henares, Madrid, Spain bDepartment of Geography, Climatology, Climate Dynamixs and Climate Change, Justus Liebig University, Giessen, Germany cHelmholtz-Zentrum Hereon, Germany dCMCC Foundation - Euro-Mediterranean Center on Climate Change, Italy A R T I C L E I N F O Dataset link: Climate Data Store Keywords: Heat waves Attribution Analogue method Autoencoders Deep learning Explainable AI A B S T R A C T Heat waves (HWs) are complex, multivariate, extreme weather events that cause significant harm to human health, ecosystems, and economies. Correct detection and attribution of HWs to anthropogenic climate change is important to better understand the underlying mechanisms and to improve predictions. In this work, we address this issue and propose a multivariate version of a hybrid approach to reconstruct heat waves, consisting of the AM and deep Autoencoders (MvEA-AM algorithm), improving existing less effective methods used until now, such as the multivariate Analogue Method (MvAM). The proposed hybrid approach produces a more reliable representation of the event than the classical MvAM for reconstructing and attributing HWs in Europe. The explainable and interpretable analysis of the obtained results is based on leveraging the SHapley Additive exPlanations (SHAP) method to explain deep learning algorithms, a capability that is not achievable with the MvAM. This explainability analysis shows that our model learns useful features during the training of the algorithm, which are aligned with the Physics of the problem, and employs the correct features during reconstruction and attribution analysis of the HWs considered. 1. Introduction Extreme heat has become more frequent and more intense globally and is expected to worsen under anthropogenic climate change (Perkins-Kirkpatrick and Lewis, 2020). Heat waves (HWs) are among the deadliest meteorological hazards (Barriopedro et al., 2023; Amengual et al., 2014; Campbell et al., 2018), significantly affecting infrastructure, livelihoods, society, and natural systems. Climate change has introduced unprecedented regimes of extreme heat, traditionally considered part of the seasonal cycle in mid-latitude to subtropical regions. In high latitudes, cumulative heat disrupts ecosystems by melting permafrost, while new extremes in mid-latitudes affect energy and health systems. Dangerous extreme heat events can now occur anywhere, leading to unexpected and highly disruptive impacts. Thus, the development of new models for the analysis, prediction, or attribution of heat waves has become a significant topic in atmospheric and extreme event research (Faranda et al., 2022; Salcedo-Sanz et al., 2023; Xu et al., 2021). ∗Corresponding author. E-mail address: [email protected] (C.M. Marina). There are several methods to address the detection, prediction, and attribution of heat waves (Jézéquel et al., 2018; Ma et al., 2024). The most widely used approaches are based on probabilistic reconstruction, such as the flow-analogues-based methods, in particular, the standard Analogue Method (AM) (Zorita and Von Storch, 1999). The AM is a classical statistical downscaling and reconstruction technique whose primary objective is to predict local meteorological variables based on large and often synoptic scale predictors. That is, large scale variables such as pressure or geopotential height. It can also be applied to reconstruct large-scale fields from local observations (Zorita and Von Storch, 1999), and has been used in extreme event prediction and attribution problems (Cattiaux et al., 2013; Ren et al., 2020; Om et al., 2024). By reconstruction of events, we refer to the process of estimating or simulating meteorological events, such as heat waves, through the identification of atmospheric states (analogues) that closely resemble the target event. Specifically, the approach involves leveraging multivariate predictors to estimate the spatial and temporal characteristics of extreme temperatures. https://doi.org/10.1016/j.atmosres.2025.108409 Received 9 May 2025; Received in revised form 28 July 2025; Accepted 5 August 2025 Atmospheric Research 328 (2026) 108409 Available online 13 August 2025 0169-8095/© 2025 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC license ( http://creativecommons.org/licenses/bync/4.0/ ). C.M. Marina et al. In its more basic version, the AM is based on the K-nearest neighbour (KNN) algorithm (Syriopoulos et al., 2023), based on the hypothesis that two similar large-scale states of the atmosphere lead to similar local effects (Lorenz, 1969). More specifically, two atmospheric states are considered analogues when they exhibit resemblance in terms of a similarity criterion and objective variables. Hence, AM involves searching through a meteorological archive for specific situations in the past that exhibit properties similar to the target situation based on a chosen similarity criterion. Previous studies have aimed to enhance the classical Analogue Method by integrating it with deep learning approaches, such as Autoencoder (AE) models—neural networks designed for efficient data compression and reconstruction. Yang and Grooms (2021), Grooms (2021) have proposed a novel use of AM with Variational AutoEncoders (VAEs) to create an optimal interpolation ensemble in a data assimilation problem. This methodology involves training an AE (Kingma et al., 2019) on the available data to generate a latent space (encoder part of the AE) from which an ensemble is constructed. The decoder part of the AE is then used to create an analogue-based ensemble of reconstructions with improved performance compared to previous approaches, such as the classical AM. Another approach involving AM with AEs was presented in Miloshevich et al. (2023), where a method based on VAEs, convolutional networks, and AM was proposed to obtain the probabilities of prolonged heat waves in France and Scandinavia. Pérez-Aracil et al. (2024) put forward a novel hybrid method that combines AM and AEs. The approach combines AEs and AM that delivers a probabilistic reconstruction of any meteorological event. The key idea was that the latent space of an AE enables a compact representation of the data, so that the AM can find more suitable analogues by focusing on this optimized latent space of the predictor rather than on the original spatially resolved fields. In PérezAracil et al. (2024), the AE-AM method was applied to temperature reconstruction during heat waves from flow-analogues, as a univariate method (using only the surface pressure field or Z500 as input). In this paper, we present a multivariate version of the AE-AM by Pérez-Aracil et al. (2024) – the MvEA-AM algorithm – which improves the results of the AE-enhanced AM by considering multiple input fields to characterize a given meteorological event. We show that the MvEA-AM can produce better results than the classical multivariate AM, and AE-AM (the MVEA-AM univariate counterpart) in problems of detection (reconstruction) and attribution of heat waves in Europe. An explainable and interpretable analysis of the results has been done using the SHapley Additive exPlanations (SHAP) method. This method allows for explaining the output of the model, enabling the interpretation of the patterns captured by the model during the training process. In the context of this work, the attribution focuses on quantifying anthropogenic climate change to the altered occurrence probability, intensity, and duration of heat waves. Specifically, for top European HWs, as shown in Section 3. This involves comparing reconstructed events under recent dynamics (with anthropogenic influence) and past dynamics (natural variability-only or with minimal anthropogenic influence) climate conditions that isolate the role of human influence in observed trends. This study focused on two key research questions: (1) Can a hybrid approach that combines the Analogue Method (AM) with deep AutoEncoders improve the reconstruction and attribution of HWs compared to classical methods? (2) How can explainability techniques enhance our understanding of the meaningful information extraction process of the AutoEncoder in the context of HW analysis? The rest of the paper has been structured in the following way: the next Section presents the methods considered in this work, including a description of the classical AM, multivariate AM, and MvAE-AM approaches. We also present here the SHAP method for the explainability of deep learning algorithms, which will be used in the experimental part of the paper to verify that the MvAE-AM aligns to the physics of the problem. Section 3 describes the data sources and case studies considered, focused on reconstructing and attributing different European HWs. Section 4 presents the results obtained using the MvAE-AM for the reconstruction and attribution of the selected HWs, and compares its performance with the classical multivariate AM and the univariate AE-AM approaches. Finally, Section 5 summarizes the main findings and offers final reflection on the research undertaken. 2. Methods This section details the classical multivariate Analogue Method (MvAM), the proposed Multivariate Autoencoder Flow-Analogue Method (MvAE-AM), and the SHAP explainability framework. First, we outline the theoretical foundations of MvAM and its limitations. Next, we introduce the MvAE-AM architecture, which combines autoencoders with the analogue method to improve heat wave reconstruction. It will be applied and evaluated in the experimental section of the paper. Finally, we describe the SHAP methodology used to interpret the model’s outputs, ensuring alignment with physical processes. 2.1. Multivariate Analogue Method (MvAM) The original AM is based on the principle that, in a stationary climate, two similar synoptic states of the atmosphere should lead to similar local effects (Lorenz, 1969). Given a dataset of an independent variable 𝑿, consisting of spatial maps with its own time stamp, and an objective map for a specific time moment (in this case describing a HW), 𝑥ℎ𝑤, we can define the Univariate Analogue Search (original AM method) as obtaining the synoptic event 𝑿, in terms of space and time, with the minimum distance to the map simultaneous to the target (HW), i.e. 𝑇 min 𝑡=0 √ √ √ √ 𝑚 ∑ 𝑖=1 (𝑥ℎ𝑤,𝑖 −𝑥𝑖,𝑡)2,(1) where 𝑥𝑖,𝑡 is the value of the 𝑖th grid point of the spatial map on the 𝑡th time, 𝑥ℎ𝑤,𝑖 is the value of the 𝑖th grid point of the target map, and 𝑚 is the total number of grid points on the map. As shown, we define the distance metric using the Euclidean one (with 𝑝= 2). Note that this can be generalized to any Minkowski distance, where 𝑝 stands for the order of the distance. Nevertheless, there is no unique way to perform AM in a multivariate scenario, where several independent variables are used to identify states similar to the HW event. In fact, the literature presents more than one approach to perform a Multivariate AM, with two main variants (Zorita and Von Storch, 1999; Caillouet et al., 2019) being most commonly used, each with own drawbacks. The first approach is a direct extension of the univariate AM for multiple features. Given 𝑛 input variables, the Analogue Search (1) is carried out in all 𝑛 dimensions (Zorita and Von Storch, 1999). In this case, the function to be minimized is given by Eq. (2). 𝑇 min 𝑡=0 √ √ √ √ 𝒏×𝑚 ∑ 𝑗=1 (𝑥ℎ𝑤,𝑗 −𝑥𝑗,𝑡)2(2) It can be seen as a concatenation of maps in a higher dimension or a flattening of all dimensions in a vector of 𝒏×𝑚. This approach is simpler to implement directly, however, it incurs high computational costs and the high dimensionality undermines the accuracy of the analogues closeness metric. The second approach is focused on building an ensemble of hierarchical searches (Caillouet et al., 2019) with 𝑛 searches being performed on maps of size 𝑚. For the specific case of the variables used here (see Section 3 for details), an arbitrary order of the 4 variables is needed. The first search found the best analogues on the first variable (e.g. 5000 analogues). These analogues are then used as the domain for the second analogues search on the second variable. Subsequently, this process is Atmospheric Research 328 (2026) 108409 2 C.M. Marina et al. Fig. 1. MvAE-AM model architecture. repeated until all searches are fulfilled. The analogues obtained after the last search are the output of this multivariate method. This approach has lower computational cost, however, a higher time demand. Besides, it is necessary to set a hierarchy between variables, which could introduce bias in the search. Note that, for example, identifying the closest states in MSL and subsequently searching for the nearest states in PEva, yields significantly different results compared to conducting this search in the reverse order. When dealing with a limited number of predictor variables, one can perform a correlation analysis or leverage existing knowledge regarding the significance of certain variables in relation to the problem at hand. However, establishing a hierarchy becomes challenging when scaled up, and even when a hierarchy is determined through correlation analysis or prior insights, its status as the optimal hierarchy is not ensured. Therefore, in order to avoid self-introduced bias by the method, we will use the first approach by Zorita and Von Storch (1999) as a benchmark to compare in this work. 2.2. Multivariate AutoEncoder Flow-Analogues (MvAE-AM) An AutoEncoder (AE) is a Deep Learning (DL) technique to reduce data dimensionality (Goodfellow et al., 2016; Liu et al., 2023), which has been widely used in atmospheric science (Singh and Goyal, 2023). The model comprises three main parts: the encoder, the latent space, and the decoder (see Fig. 1 for a graphical representation of the architecture). The primary function of the encoder is to compress the essential input data into a lower-dimension representation known as the latent space. The decoder, on the other side, serves as an auxiliary component, typically less complex than the encoder, and aims to reconstruct the original input from the codified latent space. This design compels the encoder to represent the input in a way that facilitates reconstruction by a simpler network. An AE-based model has been successfully used as an alternative DL-based AM for HW reconstruction in a univariate case by PérezAracil et al. (2024), known as AE-AM. Further, multichannel AEs have been used to deal with multivariate or multimodal problems in the literature (Nikisins et al., 2019). In this work, we propose extending the AE-AM model to multivariate problems through a multi-channel approach. As shown in Fig. 1, each variable, 𝑥1, 𝑥2,…, 𝑥𝑛, enters a different channel. Then, the AE extracts the meaningful information from each channel and creates the univariate latent space. Thereafter, the decoder tries to reconstruct the inputs to 𝑥′ 1, 𝑥′ 2,…, 𝑥′ 𝑛. We call this new model Multivariate AutoEncoder Flow-Analogues Method (MvAE-AM). The main advantage of this procedure is that it can mix information from different variables at different scales and regions in a univariate virtual space. Then, the Analogue Search is conducted in the latent space: 𝑇 min 𝑡=0 √ √ √ √ √ 𝑘 ∑ 𝑗=1 (ℎ𝑤,𝑗 −𝑗,𝑡)2,(3) where 𝑘 is the reduced new dimension of search that corresponds to the size of the latent space, . This redefinition of the Analogue Search allows for dealing with the problem of the loss of accuracy of the similarity metric. That is, we are in a smaller univariate space, and the noise from the different input variables is removed. One of the advantages of the MvAE-AM model is the possibility to apply explainability techniques (XAI techniques). In this study, we undertake a detailed investigation of the model’s explainability. We thus seek to ascertain the regions of the variables on which the model focuses to encode latent space. This interpretability cannot be obtained with the classical MvAM method. These model capabilities open up a whole horizon of possibilities and analysis, mainly for the attribution of extreme events. 2.3. Shapley Additive exPlanations (SHAP) Interpreting model predictions in deep learning (DL) is crucial in various domains, with particular importance in weather and climate sciences (Bommer et al., 2024). In the last few years, different works introduced XAI techniques to climate-related applications, with the objective of leveraging explainable AI to enhance the interpretability of outcomes derived from deep learning (DL) methodologies (Dikshit and Pradhan, 2021). In particular, SHAP is a theoretical approach used to interpret the outcomes of ML and DL models (Lundberg and Lee, 2017). It constructs an additive feature attribution model, assigning a specific predicted importance value to each variable, considering its dimensionality. The SHAP values measure the effect of the input on the output of the model. SHAP unifies different alternative explanation methods, including Local Interpretable Model-Agnostic Explanations (LIME), Deep Learning Important FeaTures (DeepLIFT) and Layer-wise Relevance Propagation (LRP). SHAP ensures that the ML/DL model learns useful features during training, according to the physics of the problem, and employs the correct features during the reconstruction of the considered event, the HW (Lundberg and Lee, 2017). In this work, we use a DeepLIFT specifically designed for DL models and multidimensional tensor data (Shrikumar et al., 2017). In particular, it uses background samples of the dataset to approximate feature importance values integrating throughout the samples. 3. Data and case studies For this work, ERA5 Reanalysis data (Hersbach et al., 2020) from 1940 to 2022, with a 2◦ grid resolution in latitude and longitude have been used. A coarser resolution than the original ERA5 (0.25◦ x0.25◦ ◦ ) is used here due to GPU memory constraints. To assess the impact of the coarser resolution, we performed a sensitivity analysis on the resolution value for the France 2003 HW. Fig. A.23 shows that the results remain largely consistent, while the execution at 2◦ resolution required 10 times less computational time compared to 0.25◦. The maximum daily Atmospheric Research 328 (2026) 108409 3 C.M. Marina et al. temperature at 2 metres (𝑇𝑚𝑎𝑥) is the target (𝐗∗) variable. As for the predictors variables (𝐗), MSL (mean sea level pressure), SM (soil moisture, as volumetric soil water), PEva (potential evaporation), and Z500 (geopotential height at 500hPa) are used. About PEva, negative values means evaporation and positive values corresponds to condensation. As indicated previously, the set of all four predictor variables has been denoted as ZPMS (Z500 - PEva - MSL - SM) in this work. For the identification of HWs, we define a HW as a period of at least 3 consecutive days during which the 𝑇𝑚𝑎𝑥 exceeds a specific threshold. The threshold is set to the 90th percentile (P90) of the 𝑇𝑚𝑎𝑥 calculated over a reference period from 1980 to 2010 using a 31-day moving window for smoothing (see Russo et al. (2015) for details). This is the most common choice in literature (Sánchez-Benítez et al., 2020) for HW definition, but it is not the only option (Barriopedro et al., 2023; Fenner et al., 2018; N et al., 2025). Note that any of the most common HW definitions in the literature could be used in the proposed methodology, without significant changes in the results obtained. For the comparison between the MvAE-AM and MvAM methods, as well as for the attribution analysis, we consider eight of the most intense, long-lasting and well documented European HWs (see Table 1). The duration of the HWs varies from one week to one month, and the spatial domain regions for ZPMS and 𝑇𝑚𝑎𝑥 are listed in Table 1. The Skill-Score function used to quantify the improvement of the MvAE-AM method against the MvAM method is defined as follows: 𝑆𝑆 =(1 − 𝑀𝐴𝐸𝑀𝑣𝐴𝐸−𝐴𝑀 𝑀𝐴𝐸𝑀𝑣𝐴𝑀 )× 100.(4) For the explainability discussion, ten alternative different case studies are used (see Table 2). We consider two types of scenarios: extreme events (HWs), and not extreme events (no-HW). Five events of each type are considered approximately centred over France, to enable both individual comparisons and an analysis of the average behaviour of these two types of situation. Moreover, since each of the HWs has a different persistence and duration we have selected 10 days from each event. That is, to avoid giving more weight to one event over another and to make the comparison of scenarios as fair as possible. Finally, for the attribution analysis, recent and past dynamics periods need to be defined. The recent dynamics period, between 1981 and 2022, is called from now on the ‘‘Post’’ period, whereas the past dynamics period, between 1940 and 1980, will be, hereafter, denoted the ‘‘Pre’’ period. The attribution that we perform in this work is a conditional probabilistic attribution, subject to a specific extreme event (HW in this case). This means that all analysis carried out is focused on specific case studies in a conditional way. Note that the objective of the conditional probabilistic attribution is not to prove causality between events or relationships, as dynamical attribution does. Further discussions on this topic can be found in Keellings and Hernández Ayala (2019), Risser and Wehner (2017). 4. Results In this section, a comparison between the classical MvAM method and the proposed MvAE-AM method is conducted, in terms of the reconstruction of the eight large European heat wave events (see Table 1) with four different predictor variables (see also Section 3). Taking advantage of the SHAP XAI-technique, an explainability analysis of the reconstruction obtained is then performed to identify the most important regions of each predictor field for HWs and no-HW (not extreme) case studies (see Table 2). Finally, on the reconstruction results of the HW events and the important regions, an attribution analysis is conducted. 4.1. Comparison of multivariate heat waves reconstruction between MvAEAM and MvAM methods First, the performance of the MvAE-AM and the MvAM methods to reconstruct the HWs is compared. Analogues of the original ZPMS fields and the encoded latent space are retrieved from the dataset, and the corresponding 𝑇𝑚𝑎𝑥 field of each analogue are retained to reconstruct the HW events. Fig. 2 shows the distributions of the 𝑇𝑚𝑎𝑥 reconstructions obtained by the MvAM and MvAE-AM methods with a latent space dimension of 1200, where the red lines correspond to the average 𝑇𝑚𝑎𝑥 of the target event from ERA5. Tables A.3–A.6 contain a sensitivity analysis of the latent space values, from which 1200 is selected as the most stable. Note that the proposed MvAE-AM method outperforms the classical MvAM in four of the considered HWs, namely France 2003, Spain 1995, Greece 1987 and Germany 2006 with a Skill-Score of 10.25%, 7.17%, 7.72% and 26.89%, respectively (see Tables A.3 and A.4). In the case of Russia 1954 HW, the distributions are almost identical, with a slight improvement of the MvAE-AM. As for the other three case studies, the classical MvAM provides slightly better distributions and Skill-Score (see Tables A.5 and A.6). It is interesting to note that only in two of the eight HWs considered (Germany 2006 and Russia 1954), the reconstructions clearly reach the target 𝑇𝑚𝑎𝑥. For the other six cases of HWs, both methods clearly underestimate the extreme temperatures reached during the heat waves. Additionally, for two heat waves (Greece 1987 and Balkans 2007), both methods reach the target 𝑇𝑚𝑎𝑥, but only on a few analogues. In the other HWs studied, the reconstruction does not come close to the target 𝑇𝑚𝑎𝑥. It should be noted that an underestimation of the observed HWs magnitude is common when applying AM-based methods to reconstruct extreme events. This can be partially explained by: (1) a biased selection of the heat wave events towards the regions and intervals of maximum severity; (2) the fact that some of the analysed HWs were record-breaking events over their respective regions of occurrence (i.e., with few historical analogues of the given severity); and (3) the lack of consideration of other amplification factors that could have contributed to exacerbating the magnitude of the HW events (e.g. land–atmosphere coupling; see e.g. Barriopedro et al. (2023) and references therein). The previous analysis is further confirmed by Tables A.3–A.6 (see the Appendix), where details of the results obtained by the MvAM and the MvAE-AM are provided. Focusing on the Skill-Score results (see Eq. (4)), the proposed MvAE-AM method outperforms the classical MvAM in France 2003, Spain 1995, Greece 1987, and Germany 2006 HWs. The improvement range varies from 7.17% to 26.89%, while the difference in the 𝑇𝑚𝑎𝑥 field from 0.12 ◦C to 0.24 ◦C. For the France 2003 HW, the proposed method surpasses the classical MvAM by the greatest temperature difference (0.24 ◦C), and a Skill-Score of 10.24%. On the other hand, the Germany 2006 HW presents the largest SkillScore (26.89%), with a 𝑇𝑚𝑎𝑥 difference of 0.17 ◦C. For the scenarios of the Poland 1994, Balkans 2007, and Russia 2010 HWs, the classical MvAM stands out over the MvAE-AM, especially for the extreme case of the mega-HW of Russia 2010. For the Poland 1994 and Balkans 2007 HWs the Skill-Score is −5.72% and −5.02% and the difference in terms of 𝑇𝑚𝑎𝑥 is 0.12 ◦C and 0.09 ◦C, respectively. For Russia 2010 HW, the Skill-Score is −12.17% and the difference in 𝑇𝑚𝑎𝑥 is 0.26 ◦C. Finally, for Russia 1954 HW, no method clearly outperforms the other. For latent space dimensions of 1200 (half of the original dimension), the MvAE-AM surpasses the performance of the MvAM, while for lower dimensionality, the classical method obtains better reconstructions. In this case, a trade-off between performance and dimension is required. Figs. 3and 4 show the France 2003 HW in ERA5 and the two analogues reconstructions, illustrating the predictors and target field 𝑇𝑚𝑎𝑥. Differences between the reconstructions and ERA5 are expressed in terms of the standard deviation of those differences. Hatched areas in each figure indicate regions where the difference between the event and the analogue reconstructions is less than one standard deviation (diff < Atmospheric Research 328 (2026) 108409 4 C.M. Marina et al. Table 1 Summary information of the heat waves considered for reconstruction and attribution. Heat wave Duration ZPMS grid (X) 𝑇𝑚𝑎𝑥 grid (X∗) France 2003 01 Aug.–19 Aug. 32◦N–70◦N 42◦N–50◦N 28◦W–30◦E 6◦W–8◦E Spain 1995 16 Jul.–24 Jul. 32◦N–70◦N 34◦N–42◦N 28◦W–30◦E 10◦W–4◦E Greece 1987 18 Jul.–27 Jul. 28◦N–66◦N 34◦N–44◦N 8◦W–50◦E 18◦E–32◦E Germany 2006 09 Jul.–31 Jul. 24◦N–72◦N 44◦N–54◦N 28◦W–30◦E 4◦W–16◦E Poland 1994 21 Jul.–11 Aug. 32◦N–70◦N 48◦N–56◦N 18◦W–40◦E 14◦E–26◦E Balkans 2007 15 Aug.–28 Aug. 32◦N–70◦N 40◦N–52◦N 8◦W–50◦E 18◦E–42◦E Russia 2010 16 Jul.–19 Aug. 32◦N–70◦N 38◦N–60◦N 22◦E–80◦E 40◦E–60◦E Russia 1954 01 Jul.–12 Jul. 32◦N–70◦N 44◦N–60◦N 8◦W–50◦E 28◦E–48◦E Fig. 2. Comparison of the reconstructed 𝑇𝑚𝑎𝑥 distributions obtained by MvAE-AM and AM, for latent space dimension of 1200 in the HWs of France 2003, Spain 1995, Greece 1987, Germany 2006, Poland 1994, Balkans 2007, Russia 2010 and Russia 1954, where the target is the average 𝑇𝑚𝑎𝑥 of the events. Atmospheric Research 328 (2026) 108409 5 C.M. Marina et al. Fig. 3. The France 2003 HW MSL, Z500 and 𝑇𝑚𝑎𝑥 in ERA5 and the MvAE-AM and MvAM reconstructions,. Hatched areas denote regions where differences between the 2003 event and the analogue reconstructions are lower than one standard deviation of the differences (diff < 1 std). Table 2 Summary information of the case studies considered for the XAI analysis. ZPMS and 𝑇𝑚𝑎𝑥 grid are the same as France 2003 from Table 1. Type event Duration Heat wave 01 Aug.–10 Aug. 2003 Heat wave 14 Jul.–23 Jul. 2013 Heat wave 09 Jun.–18 Jun. 2017 Heat wave 27 Jun.–06 Jul. 2019 Heat wave 20 Aug.–29 Aug. 2022 Not extreme 01 Jul.–10 Jul 2004 Not extreme 01 Jul.–10 Jul. 2012 Not extreme 01 Aug.–10 Aug. 2014 Not extreme 06 Jun.–15 Jun. 2019 Not extreme 21 Jun.–30 Jun. 2021 1 std), highlighting the areas of best agreement. The black boxes on the 𝑇𝑚𝑎𝑥 panels present the areas of interest, where the HWs occurred and thus the reconstruction focus. Furthermore, the differences between the analogues reconstructions for all variables are depicted in Fig. 5. Regarding the target field 𝑇𝑚𝑎𝑥, both methods capture the main patterns of the France 2003 HW (see Fig. 3). The high 𝑇𝑚𝑎𝑥 characterizes the central and southern part of the Iberian Peninsula, with lower temperatures over northern Europe and the Alpine region, an anomalously high temperature pole in France, and the extended area around Paris. As Fig. 3 shows, both methods provide analogues colder than 𝑇𝑚𝑎𝑥 of the event, except for the hatched regions, where the difference is less than one standard deviation. Then MvAE-AM results in a more precise reconstruction, while except for Norway and the British Isles, MvAE-AM provides hotter analogues and is closer to the target event (Fig. 5). For MSL, the HW event is characterized by a high-pressure region in the centre of Europe, surrounded by two low-pressure areas in the north-west and south-est of the considered region (see Fig. 3). This can also be noted in the Z500 field, where an omega-blocking or ridges systems is observed (Sousa et al., 2021). These meteorological events are usually associated with HWs and extreme temperatures (Neal et al., 2022; Fazel-Rastgar and Sivakumar, 2024). Regarding the performance of both approaches, the proposed MvAE-AM is more accurate in reconstructing the high-pressure and north-west low-pressure areas, as can be seen on the hatched areas of Fig. 3. Both analogues can correctly reconstruct the low-pressure area of the South-East. Fig. 5 shows that the MvAE-AM reconstructs a lower pressure than MvAM on the low-pressure region of the North-West and a higher pressure on the high-pressure area on the Centre-West. While the MvAE-AM outperforms in the reconstruction of the omega-block (or ridge) itself, the MvAM seems to have better performance in the top region of the Z500 map, above the British Isles (see hetched areas in Fig. 3). SM and PEva reconstruction differences shown in Fig. 5 are highly connected in our interest region: the central Europe, from the Iberian peninsula to the Balkans. The event seems to be dominated by a low SM (less than 0.1𝑚3𝑠3) and a high evaporation (where negative PEva means evaporation and positive PEva means condensation), due to the dry characteristics of a HW scenario. The MvAE-AM outperforms the traditional method throughout the map, especially in the region of interest. In this case, negative values in Fig. 5 means that MvAE-AM found drier analogues than MvAM. Atmospheric Research 328 (2026) 108409 6 C.M. Marina et al. Fig. 4. The France 2003 HW SM and PEva in ERA5 and the MvAE-AM and MvAM reconstructions. Hatched areas denote regions where differences between the event and the analogue reconstructions are lower than one standard deviation of the differences (diff < 1 std). Fig. 5. France 2003 HW analogues reconstructions differences, MvAE-AM and MvAM, for the MSL, Z500, 𝑇𝑚𝑎𝑥 SM and PEva. Hatched areas denote regions where differences between the analogues reconstructions are lower than one standard deviation of the differences (diff < 1 std). Atmospheric Research 328 (2026) 108409 7 C.M. Marina et al. Fig. 6. Comparison of SHAP importance region for the MSL input channel of the AE (France 2003 HW); (a) Stands for the mean importance map; (b) (d) (f) (h) (j) represent the importance map of the five HW case studies in France (see Section 3); (c) (e) (g) (i) (k) show the difference between the case studies and the mean importance map, where lined part represents regions with the difference less than the standard deviation (diff < std). 4.2. Explainability analysis and discussion with SHAP It is possible to see that the SHAP technique offers explainability to the MvAE-AM model, which is a differencing aspect of the proposed methodology versus the traditional MvAM. As stated in Section 2.3, SHAP provides a shapely metric of the given input for specific case studies, that can be interpreted as the regions of fields where the model focuses to codify the most relevant information. The SHAP technique adds explainability to the MvAE-AM model, distinguishing the proposed methodology from the traditional MvAM. As described in Section 2.3, SHAP quantifies each input’s contribution to specific case studies, highlighting the regions within fields that the model prioritizes to encode the most significant information. This analysis is performed for each input field (Z500, PEva, SM, MSL), and the discussion is along which field regions are used by the AE to construct the latent space. Note that for this explainability analysis we consider 5 HWs occurred in France, and other 5 periods of not HW (not extreme event) in the same region, so we can appreciate differences in the output of MvAE-AM model for both cases. Table 2 in Section 3 shows the periods considered in this analysis. Fig. 6 shows a comparison of the shapely important regions for the case studies of the 5 HWs considered, for MSL input variable. Note that, in this case, the regions of Britain, northern of France, the North of Germany and Denmark, the Baltic and North Sea are the most relevant zones of information for the algorithm (in terms of SHAP values). As can be seen in Fig. 6, each HW has its own particularities, but all of them agree on the mentioned pattern of importance. The red region, i.e., from North of France to North of Germany and the North Sea, has a direct positive impact on the latent space codification. In contrast, the blue region of Britain and the Baltic Sea has an inverse (or negative) contribution to the latent space of the encoder. Fig. 7 shows the SHAP analysis for the no-HW cases (also focused on MSL input variable). In this case there is an agreement between cases, where the Atlantic and North regions are the most relevant ones, and the maps visibly different patterns compared to the HW cases. Fig. 8 shows the differences between the mean importance map for the HWs and no-HWs cases (MSL variable). This plot emphasizes the previous analysis, where extreme (HW) and non-extreme (no-HW) events have differentiated characteristic patterns. Note that all of this importance analysis is carried out in terms of the regions in which the AE is focusing to codifying the field in the latent space. Figs. A.12–A.14, in the Appendix section, show the results for the SHAP analysis in SM, PEva and Z500. Figs. A.15–A.17 show the analysis for the non-extreme events (no-HW). Figs. A.18–A.20 show the difference between extreme and non-extreme important regions. 4.3. Attribution analysis As stated in Section 3, we defined two periods to perform the attribution analysis: Pre and Post periods. The first corresponds to the past dynamics period, from 1940 to 1980, while the second one represents the recent dynamics period, from 1981 to 2022. We trained two different AEs separately, each one for each period. AE trained in Pre period has not seen extreme temperature events, while the AE trained on Post knows the dynamics of HW events. Then, we can compare how each AE performs when they have to reconstruct a HW using analogues. This comparison, along with the expected differences between a method limited to a past dynamics period and one restricted to a recent dynamics period, allows for the attribution of the anthropogenic effect. The preliminary assumption is that the Pre-AE does not know how to encode an HW. Although these events may have occurred in the Pre period, they have been so few times, that the AE might treat them as noise. Therefore, when reconstructing the HW, this AE will find very distant events. In contrast, the Post-AE is in a more advantageous position. It has seen a variety of HWs and has been able to learn how to encode them correctly. If we add to this the fact that in the Post period there are more cases of HWs, we expect to see large differences (of the order of 1 ◦C) between the reconstructions of each AE. Fig. 9 presents the comparison between the distribution of reconstructed 𝑇𝑚𝑎𝑥 by each AE for the eight large European HWs considered, where the orange distributions correspond to the Post-AE and the blue distributions to the Pre-AE. We can see that for Germany 2006 and Russia 1954, the Post-AE is able to reach the actual values of 𝑇𝑚𝑎𝑥 Atmospheric Research 328 (2026) 108409 8 C.M. Marina et al. Fig. 7. Comparison of SHAP importance region for the MSL input channel of the AE; (a) Stands for the mean importance map, (b) (d) (f) (h) (j) represent the importance map of the five not extreme (no-HW) case studies in France (see Section 3); (c) (e) (g) (i) (k) show the difference between the case studies and the mean importance map, where lined part represents regions with the difference less than the standard deviation (diff < std). Fig. 8. Difference between the mean importance map of the five France HW cases and the mean importance map of the five France no-HW cases for variable MSL, where lined part represents regions with the difference less than the standard deviation (diff < std). of the HWs. Moreover, the model demonstrates some capability in reconstructing the temperature patterns of Greece 1987 and Balkans 2007 HWs. However, this reconstruction is not as definitive as in the earlier ones, relying on select analogues from the distribution tails. Nevertheless, for France 2003, Spain 2005, Poland 1994 and Russia 2010, while Post-AE is closer to the target than Pre-AE, none of the AEs find similar enough 𝑇𝑚𝑎𝑥 to the HWs. A similar pattern emerges when the variability between the obtained analogues is examined. Where the target values are reached, the distributions of analogues found by Post-AE have a higher standard deviation (std) than those found by Pre-AE. If we subtract the Post-AE std from the Pre-AE std, we obtain: +0.09 for Greece 1987, +0.08 for Germany 2006, +0.05 for Balkans 2007, and +0.11 for Russia 1954. On the other hand, if we extend this to HWs where the target values of 𝑇𝑚𝑎𝑥 are not reached we get: −0.02 for France 2003, −0.04 for Spain 1995, −0.09 for Poland 1994, and −0.02 for Russia 2010. We next evaluate the temperature intensity increases between Pre and Post-based reconstruction, Only in the case of Greece 1987 the difference is less than 1 ◦C, namely 0.6 ◦C. For Russia 2010, the increase is about 1 ◦C. Then, for France 2003, Spain 1995, and Germany 2006 the increase is between 1.1 ◦C and 1.4 ◦C. Subsequently, the increase of 𝑇𝑚𝑎𝑥 is higher than 1.5 ◦C for Poland 1994, Balkans 2007 and Russia 1954. To contextualize these differences, extend the SHAP-based explainability framework (Section 2.3) to compare variable importance between Preand Post-period models. Since we have a trained model for each period, we can apply SHAP to understand which regions of which variables are influential in encoding an HW for each period. This approach enables us to detect any notable shifts in the dynamics of the variables when comparing the actual dynamics period to the past dynamics period. For instance, we have conducted this comparison on the variable MSL. Atmospheric Research 328 (2026) 108409 9 C.M. Marina et al. Fig. A.18. Difference between the mean importance map of the five France HW cases and the mean importance map of the five France not extremes (no-HW) cases for variable SM, where lined part represents regions with the difference less than the standard deviation (diff < std). Fig. A.19. Difference between the mean importance map of the five France HW cases and the mean importance map of the five France not extremes (no-HW) cases for variable PEva, where lined part represents regions with the difference less than the standard deviation (diff < std). Atmospheric Research 328 (2026) 108409 16 C.M. Marina et al. Fig. A.20. Difference between the mean importance map of the five France HW cases and the mean importance map of the five France not extremes cases (no-HW) for variable Z500, where lined part represents regions with the difference less than the standard deviation (diff < std). Fig. A.21. Comparison of SHAP importance region for the Z500 input channel of the AE trained only on the Pre period. (a) is the mean importance map, (b) (d) (f) (h) (j) are the importance map of the five France heat waves case studies (see Section 3), (c) (e) (g) (i) (k) are the difference between the case studies and the mean importance map, where lined part represents regions with the difference less than the standard deviation (diff < std). Atmospheric Research 328 (2026) 108409 17 C.M. Marina et al. Table A.4 Average MvAE-AM results for different dimensions of the latent space (heat waves of Greece 1987 and Germany 2006). Greece 1987 Latent dim. Avg. ZPMS Diff. Avg. 𝑇𝑚𝑎𝑥 Diff. Avg. 𝑇𝑚𝑎𝑥 Std. 𝑇𝑚𝑎𝑥 SS Target – – – 30.1892 – – MvAM – 0.1277 1.6913 28.4979 0.3034 – MvAE-AM 400 0.0190 1.8024 28.3868 0.3295 −6.57 % 600 0.0189 1.5608 28.6284 0.5680 7.72% 800 0.0214 1.6384 28.5508 0.5064 3.12% 1200 0.0253 1.5786 28.6106 0.5380 6.66% Germany 2006 Latent dim. Avg. ZPMS Diff. Avg. 𝑇𝑚𝑎𝑥 Diff. Avg. 𝑇𝑚𝑎𝑥 Std. 𝑇𝑚𝑎𝑥 SS Target – – – 26.1944 – – MvAM – 0.1121 0.6589 25.5355 0.3499 – MvAE-AM 400 0.0158 0.5556 25.6388 0.3612 15.68% 600 0.0174 0.4817 25.7127 0.3621 26.89% 800 0.0215 0.4870 25.7074 0.3585 26.08% 1200 0.0205 0.5160 25.6784 0.3584 21.69% Table A.5 Average MvAE-AM results for different dimensions of the latent space (heat waves of Poland 1994 and Balkans 2007). Poland 1994 Latent dim. Avg. ZPMS Diff. Avg. 𝑇𝑚𝑎𝑥 Diff. Avg. 𝑇𝑚𝑎𝑥 Std. 𝑇𝑚𝑎𝑥 SS Target – – – 28.7656 – – MvAM – 0.1362 2.1121 26.6535 0.3628 – MvAE-AM 400 0.0182 2.3125 26.4531 0.3963 −9.49 % 600 0.0222 2.2428 26.5228 0.3834 −6.19 % 800 0.0217 2.2330 26.5326 0.3870 −𝟓.𝟕𝟐% 1200 0.0264 2.2353 26.5303 0.3782 −5.83 % Balkans 2007 Latent dim. Avg. ZPMS Diff. Avg. 𝑇𝑚𝑎𝑥 Diff. Avg. 𝑇𝑚𝑎𝑥 Std. 𝑇𝑚𝑎𝑥 SS Target – – – 29.2392 – – MvAM – 0.1346 1.3042 27.9350 0.4398 – MvAE-AM 400 0.0171 1.7779 27.4613 0.2642 −36.32 % 600 0.0205 1.3697 27.8695 0.4541 −𝟓.𝟎𝟐% 800 0.0197 1.3940 27.8452 0.4306 −6.89 % 1200 0.0229 1.3983 27.8409 0.4264 −7.22 % Table A.6 Average MvAE-AM results for different dimensions of the latent space (heat waves of Russia 2010 and Russia 1954). Russia 2010 Latent dim. Avg. ZPMS Diff. Avg. 𝑇𝑚𝑎𝑥 Diff. Avg. 𝑇𝑚𝑎𝑥 Std. 𝑇𝑚𝑎𝑥 SS Target – – – 32.6340 – – MvAM – 0.1693 2.0693 30.5647 0.2375 – MvAE-AM 400 0.0128 2.3517 30.2823 0.2413 −13.65 % 600 0.0156 2.3811 30.2529 0.2853 −15.07 % 800 0.0163 2.3809 30.2531 0.2791 −15.06 % 1200 0.0168 2.3212 30.3128 0.2746 −𝟏𝟐.𝟏𝟕% Russia 1954 Latent dim. Avg. ZPMS Diff. Avg. 𝑇𝑚𝑎𝑥 Diff. Avg. 𝑇𝑚𝑎𝑥 Std. 𝑇𝑚𝑎𝑥 SS Target – – – 28.8298 – – MvAM – 0.1346 0.9010 27.9288 0.5425 – MvAE-AM 400 0.0168 1.1177 27.7121 0.5474 −24.05 % 600 0.0189 0.9332 27.8966 0.5628 −3.57 % 800 0.0186 1.0940 27.7358 0.5567 −21.42 % 1200 0.0212 0.8565 27.9733 0.5580 4.94% Atmospheric Research 328 (2026) 108409 18 C.M. Marina et al. Fig. A.22. Comparison of SHAP importance region for the Z500 input channel of the AE trained only on the Post period. (a) is the mean importance map, (b) (d) (f) (h) (j) are the importance map of the five France heat waves case studies (see Section 3), (c) (e) (g) (i) (k) are the difference between the case studies and the mean importance map, where lined part represents regions with the difference less than the standard deviation (diff < std). Fig. A.23. Reconstruction of 𝑇𝑚𝑎𝑥 distribution obtained by MvAE-AM and MvAM for France 2003, with 0.25◦ grid resolution. Data availability Data is completely open through the Copernicus ClimateDataStore. References Amengual, A., Homar, V., Romero, R., Brooks, H., Ramis, C., Gordaliza, M., Alonso, S., 2014. Projections of heat waves with high impact on human health in Europe. Glob. Planet. Chang. 119, 71–84. Barriopedro, D., García-Herrera, R., Ordóñez, C., Miralles, D., Salcedo-Sanz, S., 2023. Heat waves: Physical understanding and scientific challenges. Rev. Geophys. e2022RG000780. Bommer, P., Kretschmer, M., Hedström, A., Bareeva, D., Höhne, M., 2024. Finding the right XAI method—a guide for the evaluation and ranking of explainable AI methods in climate science. Artif. Intell. Earth Syst. 3, e230074. Caillouet, L., Vidal, J., Sauquet, E., Graff, B., Soubeyroux, J., 2019. SCOPE climate: a 142-year daily high-resolution ensemble meteorological reconstruction dataset over France. Earth Syst. Sci. Data. Campbell, S., Remenyi, T., White, C., Johnston, F., 2018. Heatwave and health impact research: A global review. Heal. Place 53, 210–218. Cattiaux, J., et al., 2013. US heat waves of spring and summer 2012 from the flow-analogue perspective. Bull. Am. Meteorol. Soc. 94, S10–S13. Dikshit, A., Pradhan, B., 2021. Interpretable and explainable AI (XAI) model for spatial drought prediction. Sci. Total. Environ. 801, 149797. Faranda, D., Bourdin, S., Ginesta, M., Krouma, M., Noyelle, R., Pons, F., Yiou, P., Messori, G., 2022. A climate-change attribution retrospective of some impactful weather extremes of 2021. Weather. Clim. Dyn. 3, 1311–1340. Fazel-Rastgar, F., Sivakumar, V., 2024. A case study on more recent heat wave occurred in South Africa, based on background weather synoptic and dynamic characteristics analysis. Bull. Atmos. Sci. Technol.. Fenner, D., Holtmann, A., Krug, A., Scherer, D., 2018. Heat waves in Berlin and Potsdam, Germany – long-term trends and comparison of heat wave definitions from 1893 to 2017. Int. J. Climatol. 39, 2422–2437. Goodfellow, I., Bengio, Y., Courville, A., 2016. Deep Learning. MIT Press. Grooms, I., 2021. Analog ensemble data assimilation and a method for constructing analogs with variational autoencoders. Q. J. R. Meteorol. Soc. 147, 139–149. Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., et al., 2020. The ERA5 global reanalysis. Q. J. R. Meteorol. Soc. 146, 1999–2049. Jézéquel, A., Yiou, P., Radanovics, S., 2018. Role of circulation in European heatwaves using flow analogues. Clim. Dyn. 50, 1145–1159. Keellings, D., Hernández Ayala, J., 2019. Extreme rainfall associated with hurricane maria over puerto rico and its connections to climate variability and change. Geophys. Res. Lett. 46, 2964–2973. Kingma, D., Welling, M., et al., 2019. An introduction to variational autoencoders. Found. Trends Mach. Learn. 12, 307–392. Liu, T., Wang, J., Liu, Q., Alibhai, S., Lu, T., He, X., 2023. High-ratio lossy compression: Exploring the autoencoder to compress scientific data. IEEE Trans. Big Data 9, 22–36. Lorenz, E., 1969. Atmospheric predictability as revealed by naturally occurring analogues. J. Atmos. Sci. 26, 636–646. Lundberg, S., Lee, S., 2017. A unified approach to interpreting model predictions. arXiv Preprint arXiv:1705.07874. Atmospheric Research 328 (2026) 108409 19 C.M. Marina et al. Ma, K., Gong, H., Wang, L., 2024. Attribution of the concurrent extreme heatwaves in Northern Europe and Northeast Asia in 2018. Atmos. Res. 107506. Miloshevich, G., Lucente, D., Yiou, P., Bouchet, F., 2023. Extreme heatwave sampling and prediction with analog markov chain and comparisons with deep learning. arXiv Preprint arXiv:2307.09060. N, S., Gunther, S., Kjellstrom, T., Lee, J., 2025. Advancing heat wave definitions: a policy review towards prioritizing health impacts of extreme heat. Environ. Res. Lett. 20. Neal, E., Huang, C., Nakamura, N., 2022. The 2021 Pacific northwest heat wave and associated blocking: Meteorology and the role of an upstream cyclone as a diabatic source of wave activity. Geophys. Res. Lett. 49. Nikisins, O., George, A., Marcel, S., 2019. Domain adaptation in multi-channel autoencoder based features for robust face anti-spoofing. In: 2019 International Conference On Biometrics. ICB, pp. 1–8. Om, K., GuoYu, R., Jong, S., Hyon-Ok, O., Kim, S., Kang-Chol, O., 2024. Extreme rainfall events variation during 1860–1909 in the Korean Peninsula: Investigation of the possible circulation mechanism by a method of analogue. Atmos. Res. 298, 107120. Pérez-Aracil, J., Marina, C., Zorita, E., Barriopedro, D., Zaninelli, P., Giuliani, M., Castelletti, A., Gutiérrez, P., Salcedo-Sanz, S., 2024. Autoencoder-based flowanalogue probabilistic reconstruction of heat waves from pressure fields. Ann. N. Y. Acad. Sci. 1541, 230–242. Perkins-Kirkpatrick, S., Lewis, S., 2020. Increasing trends in regional heatwaves. Nat. Commun. 11. Ren, L., Zhou, T., Zhang, W., 2020. Attribution of the record-breaking heat event over Northeast Asia in summer 2018: the role of circulation. Environ. Res. Lett. 15, 054018. Risser, M., Wehner, M., 2017. Attributable human-induced changes in the likelihood and magnitude of the observed extreme precipitation during hurricane harvey. Geophys. Res. Lett. 44 (12), 457-412. 464. Russo, S., Sillmann, J., Fischer, E., 2015. Top ten European heatwaves since 1950 and their occurrence in the coming decades. Environ. Res. Lett. 10, 124003. Salcedo-Sanz, S., Pérez-Aracil, J., Ascenso, G., Ser, J.Del., Casillas-Pérez, D., Kadow, C., Fister, D., Barriopedro, D., García-Herrera, R., Giuliani, M., et al., 2023. Analysis, characterization, prediction, and attribution of extreme atmospheric events with machine learning and deep learning techniques: A review. Theor. Appl. Climatol. 1, 44. Sánchez-Benítez, A., Barriopedro, D., García-Herrera, R., 2020. Tracking iberian heatwaves from a new perspective. Weather. Clim. Extrem. 28, 100238. Shrikumar, A., Greenside, P., Kundaje, A., 2017. Learning important features through propagating activation differences. CoRR. abs/1704.02685. Singh, S., Goyal, M., 2023. An innovative approach to predict atmospheric rivers: Exploring convolutional autoencoder. Atmos. Res. 289, 106754. Sousa, P., Barriopedro, D., García-Herrera, R., Woollings, T., Trigo, R., 2021. A new combined detection algorithm for blocking and subtropical ridges. J. Clim. 34, 7735–7758. Syriopoulos, P., Kalampalikis, N., Kotsiantis, S., Vrahatis, M., 2023. KNN classification: a review. Ann. Math. Artif. Intell. 1, 33. Xu, P., Wang, L., Huang, P., Chen, W., 2021. Disentangling dynamical and thermodynamical contributions to the record-breaking heatwave over Central Europe in 2019. Atmos. Res. 252, 105446. Yang, L., Grooms, I., 2021. Machine learning techniques to construct patched analog ensembles for data assimilation. J. Comput. Phys. 443, 110532. Zorita, E., Von Storch, H., 1999. The analog method as a simple statistical downscaling technique: Comparison with more complicated methods. J. Clim. 12, 2474–2489. Atmospheric Research 328 (2026) 108409 20