Artificial intelligence reveals past climate extremes by reconstructing historical records

Domain: earth_science
Published: 2024-10-24
Editorial outcome: minor_revision

Abstract
The understanding of recent climate extremes and the characterization of climate risk require examining these extremes within a historical context. However, the existing datasets of observed extremes generally exhibit spatial gaps and inaccuracies due to inadequate spatial extrapolation. This problem arises from traditional statistical methods used to account for the lack of measurements, particularly prevalent before the mid-20th century. In this work, we use artificial intelligence to reconstruct observations of European climate extremes (warm and cold days and nights) by leveraging Earth system model data from CMIP6 through transfer learning. Our method surpasses conventional statistical techniques and diffusion models, showcasing its ability to reconstruct past extreme events and reveal spatial trends across an extensive timespan (1901-2018) that is not covered by most reanalysis datasets. Providing our dataset to the climate community will improve the characterization of climate extremes, resulting in better risk management and policies.

Introduction
In 2023, global temperatures reached unprecedented levels in climatic history, marking the warmest year on record by a substantial margin1. Correspondingly, Europe experienced its second warmest year on record, surpassed only by the unprecedented temperatures of 2020. The preceding year, 2022, registered as the third warmest for Europe, was notably characterized by the warmest summer observed in the region2. Due to the persistent warmth, the summer of 2022 exhibited numerous extreme events, notably featuring intense and prolonged heat waves with considerable regional differences. These heatwaves, in conjunction with persistently low levels of rainfall and other contributing factors, led to widespread drought conditions, impacting ecosystems, agriculture, and water resources3. A significant and alarming consequence of these climatic conditions is the escalation of wildfire activity. Again, the European summer of 2022 witnessed the highest total wildfire emissions in the last 15 years and the second-highest burnt area ever4. These wildfires not only underscore the immediate impacts of extreme weather events but also reflect a broader environmental crisis, driven by climate change.

To understand these significant changes in the European climate, it is crucial to gain a temporal and spatial perspective by questioning: How do these recent extremes compare to past events throughout the 20th Century, and what are the regional variations of the long-term trends?

The Working Group I contribution to the Sixth Assessment (AR6) Report from the IPCC reached the conclusion that Europe is one of the regions of the world with the highest increase in the intensity and frequency of hot extremes since 19605. These changes are assessed by examining trends of extreme indices such as the ones recommended by the former WMO joint Expert Team on Climate Change Detection and Indices (ETCCDI) from the Commission for Climatology/World Climate Research Program/Commission for Oceanography and Marine Meteorology6. These form part of a now larger set of WMO-recommended indices, some of which have been developed with specific sectors in mind7. In the present study, we are focusing on four ETCCDI indices, which measure the frequency of hot and cold temperature extremes, as defined in Table 1.

These indices are part of the HadEX3 dataset which contains 29 land-based ETCCDI indices derived from daily station observations of temperature and precipitation data8 as part of the full Climpact set of indices. On the one hand, the use of extreme indices derived from station data offers an extended historical perspective into the evolution of the climate, which is free from the biases inherent to climate models. On the other hand, a large amount of missing values, which vary in both time and space, are present in the observational data and are especially prominent and problematic in the early half of the 20th Century.

To circumvent the lack of information, it is common to infill the gaps using statistical interpolation methods such as Inverse Distance Weighting (IDW), Angular Distance Weighting (ADW, 9), the thin plate spline interpolation10 or Kriging11. However, these methods suffer from well-known limitations that reduce their applicability, especially when the data is very scarce. For instance, ADW only considers the spatial and angular proximity between points. By neglecting more complex relationships within the data, the ADW tends to produce overly smooth interpolations for many climate variables. Kriging is better at capturing spatial variability but requires the definition of an appropriate variogram model by setting various parameters. Moreover, it is also computationally demanding and sensitive to outliers. In recent years, deep learning-based inpainting techniques have emerged as a groundbreaking approach to reconstructing missing data in climate datasets12–18. As epitomized by the work of ref. 15 on the reconstruction of the HadCRUT4 dataset19, artificial neural networks can significantly outperform traditional infilling methods on relevant metrics, such as the root mean square error or the spatial correlation. Furthermore, they show a remarkable proficiency in capturing intricate spatial climatic patterns even in the context of a large proportion of missing data. Depending on the infilling task, either CNN (Convolutional Neural Network)-based or GAN (Generative Adversarial Network)-based approaches are generally adopted. As shown in ref. 14, the choice of one approach or the other hinges on whether pixel-level accuracy (CNN) or physical realism (GAN) is more relevant to the task.

In this study, we are opting for a CNN-based approach20. It is applied to the reconstruction of extreme indices from an intermediate product of the HadEX3 dataset with monthly frequency. This intermediate product, called hereafter HadEX-CAM, is an non-infilled version of HadEX3 obtained before the aggregation of the station data using the ADW gridding method (for more details, see “Methods”, Section HadEX3 and HadEX-CAM). This study focuses on a smaller region from the global HadEX-CAM dataset, namely the European continent. This choice is motivated by the relative abundance of data in the recent period, as well as the temperate but diverse climate of the region. Because these conditions are favorable to traditional statistical methods, this choice aims to provide a challenging benchmark for our AI method.

Results
Comparative evaluation of reconstruction methods
The deep-learning technique used to reconstruct the extreme indices, hereafter referred to as Climate Reconstruction AI (CRAI), is based on a U-Net made of partial convolutions (see “Methods”, Section Deep learning based reconstruction method). The AI models are trained using historical simulations with Earth System Models (ESMs) from the CMIP6 archive (see “Methods”, Section Dataset up for the training and the evaluation). To evaluate these models, we use data from three types of datasets that were not included in the model training: a simulation dataset (from CMIP6 models), a reanalysis dataset (ERA5), and an observational dataset (HadEX-CAM). An overview of the datasets considered in this study and their utilization is given in Supplementary Table S1. The evaluation is carried out using the root mean square error (RMSE), the Spearman rank-order correlation coefficient (SROCC), the Wasserstein distance (WD, 21), and the coefficient of determination (R2) calculated on the reconstructed values only. The results of this evaluation are systematically compared with other interpolation methods, namely IDW, Kriging, and diffusion models (see “Methods”, Section Comparative methods).

The test dataset used for the evaluation is formed by selecting a spatial sample from each month within the original pool of simulation data taken from the CMIP6 archive (see “Methods”, Section Dataset up for the training and the evaluation). The remaining data is subsequently used for the creation of the training and validation datasets. To ensure the heterogeneity of the test dataset, each monthly sample is selected randomly among all models and realizations. A mask of missing values extracted from the corresponding month of the HadEX-CAM dataset is applied to each sample in order to create artificial missing values. The resulting masked dataset is then reconstructed using the trained models and compared with the original complete data. Table 2 shows the RMSE and the SROCC for the reconstruction of the test dataset. It reveals that CRAI outperforms IDW and Kriging with a relatively constant improvement for all metrics and extreme indices. This trend is also evident in other metrics (see Supplementary Table S4), such as the coefficient of determination (R2) and the Wasserstein distance (WD), which is a measure of the distance between the probability distributions of the reference and the reconstructed data. The performance of the diffusion models closely approaches that of CRAI; however, CRAI remains superior in the majority of cases. This notable performance, achieved with fewer computational resources, indicates that the partial convolutional layers used in CRAI, but not in diffusion models, are particularly well-suited for reconstructing such climate data. All methods exhibit lower RMSE for the TX90p and TN90p indices than the TX10p and TN10p. As shown in the Supplementary Table S3, it is likely attributable to the lower mean and standard deviation of the TX90p and TN90p indices. For the SROCC, it is the other way around, i.e., a higher correlation is obtained for TX10p and TN10p compared to TX90p and TN90p.

Despite the successful results of the CRAI models applied to the test dataset, it is unclear how the models trained on simulation data perform on observational data, such as HadEX-CAM. To estimate the ability of the models for generalization, it is convenient to perform an evaluation on reanalysis datasets. Reanalysis datasets are particularly suitable for cross-validation as they are created by assimilating historical observational data and do not contain missing values. By comparing the standard deviation and the SROCC of different reanalysis datasets with HadEX3 and HadEX-CAM (see Supplementary Fig. S2), we have determined that ERA5 presents a strong similarity with the target data, although it spans a shorter period of time (1940–2018). For the sake of simplicity, we have then restricted our evaluation to the reconstruction of ERA5 only.

The evaluation of ERA5 has been performed following the same procedure as described for the test dataset. The results of this evaluation are shown in Table 3 and Supplementary Table S5. The RMSE, SROCC, WD, and R2 scores are globally better for all methods compared to the evaluation on the test dataset (Table 2 and Supplementary Table S4). This is likely due to the fact that the reconstruction of ERA5 is evaluated on a more recent period of time (1940–2018), hence with more valid values, i.e., more statistical information available for the reconstruction. This also means that the results of each evaluation metric are more uniform across indices. The CRAI models remarkably outperform all the other methods for all metrics and indices, except for three results that equal those of diffusion models. The strong performance of CRAI on different types of data highlights its ability for generalization, particularly when compared to diffusion models.

The temporal evolution and the spatial distribution of the RMSE values of the reconstruction of ERA5 using CRAI are shown in Supplementary Fig. S3. As expected, the values are larger for time periods (see Supplementary Fig. S3a) and regions (see Supplementary Fig. S3b) with a higher prevalence of missing values in the dataset (as shown in Supplementary Fig. S4). Nevertheless, the comparison shown in Fig. 1 between the RMSEs of CRAI and the best statistical method (Kriging) reveals that our method achieves better reconstructions for the whole time period and for most regions of Europe. More particularly, the time series in Fig. 1a, obtained by computing the mean of the RMSE difference across the spatial domain, exhibits stronger variations and slightly lower values for 1940–1960 and 2000–2018. It demonstrates that CRAI is particularly efficient for larger quantities of missing data in comparison with Kriging. This finding aligns with the spatial distribution of valid values in Supplementary Fig. S4b, revealing that regions characterized by a higher prevalence of missing values (e.g., northern Africa) exhibit an enhanced performance of the CRAI reconstruction with respect to Kriging (Fig. 1b). For completeness, a similar analysis is shown in Supplementary Fig. S5 where the agreement of the CRAI reconstruction with ERA5 is even larger compared to HadEX3.

To evaluate our model in a context that is closer to our intended use case, we have adopted an original approach consisting of evaluating the model using the HadEX-CAM dataset itself. For this purpose, we have created additional missing values in the original HadEX-CAM dataset by retaining for each timestep only the data points corresponding to the valid values of January 1901 (timestep with the highest prevalence of missing values). Supplementary Table S6 summarizes the evaluation of the masked HadEX-CAM dataset reconstruction. As expected, the overall performance is poorer for all methods and for all metrics compared to the previous evaluations (test dataset and ERA5). This can be attributed to the reduced amount of valid values across the full-time period, as shown by the deterioration of the evaluation metrics in Supplementary Table S7 when applying the same mask to ERA5. The results shown in Supplementary Tables S6 and S7 exhibit significantly reduced discrepancies between methods in comparison with the previous evaluations (test and ERA5 datasets). Nevertheless, the overall improved performance of CRAI compared to the other methods remains evident. It is important to note that consistency across all the evaluations presented in this study has been ensured by employing the same CRAI and diffusion models, which were trained using multiple masks of missing values. Applying these models to datasets masked with a single mask of missing values (as in Supplementary Tables S6 and S7), therefore, constitutes a challenging assessment of their performances.

Regional spatial mean
Having established that our AI trained model (CRAI) is outperforming various interpolation methods (IDW, Kriging, diffusion) in the reconstruction of three types of dataset, we now apply it to the reconstruction of the full field of the HadEX-CAM dataset over the European domain.

The differences between HadEX3, HadEX-CAM, ERA5 and the CRAI reconstruction of HadEX-CAM can be observed in the temporal domain by looking at the spatial mean of the extreme indices (Fig. 2). At the European scale, the four datasets present a good overall agreement for both indices. All results align with the conclusions of the AR6 IPCC report5 by exhibiting a clear increase (decrease) in the frequency of warm days (cool nights). Discrepancies across datasets are larger for TN10p, especially before 1960, where HadEX-CAM, CRAI, and ERA5 predict lower values compared to HadEX3. It corresponds to the period of time where the amount of missing values is the largest (see Supplementary Fig. S4a). For the same reason, the magnitude of discrepancies between HadEX3 and the other datasets is higher for the Mediterranean region (MED) that is the region containing the largest amount of missing values (see Supplementary Fig. S4b). A similar assessment stands for the TN90p and TX10p indices (see Supplementary Fig. S6). As expected, the ensemble spread of the CRAI reconstructions obtained from the twenty trained models (see Methods, Section Deep learning based reconstruction method) also reflects the evolution and the distribution of missing values in the HadEX-CAM dataset. It is barely distinguishable from the ensemble average, except for the first half of the 20th century and for the Mediterranean region, where the ensemble spread is the largest. However, the ensemble spread generally remains lower than the difference between the ensemble average and the other datasets.

Interestingly, ERA5 and CRAI greatly reduce the magnitude of the TX90p hot extremes in the 1945–1950 period. Using November 1947 as an example (Fig. 3), we can see that both ERA5 and the CRAI reconstruction show a warm spell in the Iberian peninsula, as well as western Morocco and southwest Ireland. In contrast, HadEX3 shows the warm days extending across all of Morocco and northern Algeria. For this date, only 31 valid grid point values are available in HadEX-CAM in the Mediterranean region, mostly concentrated in the Iberian peninsula and the Balkans. Based on this patchy information, the ADW gridding method in HadEX3 extrapolates the high temperatures measured in Spain into northern Africa and hence overestimates the TX90p values in the Mediterranean region. In addition, HadEX3 also presents smoother fields of warm and cold days across the continent, while ERA5 and the CRAI reconstruction show more intricate structures. It is quantified by the global Moran’s I values22,23 shown in Fig. 3 that indicate a higher spatial autocorrelation for HadEX3 compared to ERA5 and the CRAI reconstruction. This over-smoothing is a common feature of the ADW method (see refs. 24, 25) and leads to a clear underestimation of the high TX90p values in the Iberian Peninsula. Similar effects are observed for other dates and regions in the dataset. They illustrate one of the main limitations of the traditional interpolation methods that aim at reconstructing missing values based solely on the statistical information available in a rather limited spatiotemporal context. For instance, the monthly TX90p values in HadEX3 are regridded by incorporating station data from the same month within a search radius. This radius is determined by the decorrelation length scale calculated for that calendar month from all years and for stations within a wide latitude band8. In contrast, the CRAI models are trained using a large amount of physically consistent data (grounded in the specific CMIP6 model physics) that spans all months from 1901–2014. Hence, it has the ability to learn climatic patterns specific to the European climate, which are accounted for in the predictions of the gridded fields. This methodological distinction explains the larger differences found between the CRAI reconstruction and HadEX3 (Fig. 2 and Supplementary Fig. S6) for a period of time (e.g., 1901–1950) and regions (e.g., Mediterranean) where the scarcity of valid data points compromises the applicability of purely statistical methods. The same issues are also detected in the Kriging results (Fig. 1b), where the largest RMSEs are found in the Mediterranean region.

Trend analysis
Most of the published articles reporting spatial trend analysis of ETCCDI indices that include European regions26–31 are focusing on recent trends. This is partly related to the temporal limitation of the available observational datasets (e.g., E-OBS32 spanning 1950–2023). Our AI methodology offers the possibility to extend the scope of the spatial trend analysis to a longer time period with a greater precision compared to HadEX3.

Despite similar spatial mean values at the scale of the continent (Fig. 2), strong differences emerge between the HadEX3 and CRAI spatial trend analysis. Figures 4 and 5 show the long-term linear trend calculated for the entire time period (1901–2018) using the median of pairwise slopes estimator33, 34. We acknowledge that changes over this period are not linear (Fig. 2) but use this approach to summarize changes over time without losing the spatial information. The reconstruction of the HadEX-CAM dataset using CRAI exhibits a higher level of spatial heterogeneity compared to HadEX3, with complex patterns and pronounced spatial variability. The CRAI reconstruction indicates a rather limited TX90p trend for north Africa and along the coasts of southern Turkey and Syria while showing an increase in the frequency of warm days in central Europe and the Baltic Sea (Fig. 4). For TN10p however, both HadEX3 and CRAI predict a high decrease in the frequency of cool nights in north Africa and western Europe (Fig. 5). Across the continent, the CRAI reconstruction shows a more contrasted picture with regional variations, including notable negative trend values (e.g., in Ukraine and Romania) that are not accounted for in HadEX3. Both datasets also disagree in the Middle East, where the CRAI reconstruction indicates a smaller decrease of TN10p compared to HadEX3.

As depicted in Fig. 5 and Supplementary Fig. S7, the HadEX3 TN10p and TX10p spatial trends exhibit a strong resemblance. In contrast, for CRAI, TX10p diverges from TN10p in many regions, displaying a modest reduction in the occurrence of cool days across North Africa and the Middle East, while accentuating the trend notably along the Baltic Sea and southern UK. HadEX3 and the CRAI reconstruction present similar spatial patterns for TN90p (see Supplementary Fig. S8). In particular, both datasets show a notable increase in the frequency of warm nights for a large area that extends from the west of Norway to Georgia and which is not found in the TX90p maps.

Results for more recent trends (1980–2018) are also shown (see Supplementary Figs. S9, S10), including ERA5 results. Here again, the reconstruction using CRAI presents rather detailed spatial structures compared to HadEX3, similar to those observed in ERA5 and quantified by the global Moran’s I values shown in the figures. The resemblance between the CRAI and the ERA5 trends is particularly notable in certain regions with local extreme values such as the Black-Sea coasts.

Early cases of heatwaves and cold waves
The CRAI models allow for a detailed regional investigation of extreme events, especially in the early time period when existing observational datasets are based on fewer observations and have more missing values. The summer of 1911 was notable for an intense and deadly heatwave that affected various parts of Europe, especially Metropolitan France. Many French regions experienced temperatures well above their normal averages from July to September. As shown in Fig. 6 and Supplementary Fig. S11, hot extremes in southern France were still detected in September 1911 by CRAI, extending the hot summer into fall. Despite the lack of data (visible in the central panel), values over 40% of the days exceeding the 90th percentile of the maximum temperatures are being predicted by the CRAI models. The standard deviations derived from the twenty models trained for TX90p and TN90p (see Supplementary Fig. S12) provide a quantification of the models’ uncertainty. In the region affected by the heatwave, the standard deviation is consistently low (below 3%), suggesting a high level of confidence in the prediction of these hot extremes. A similar spatial pattern across western-central Europe is confirmed by the indices derived from the 20th Century Reanalysis Version 335 in Supplementary Fig. S13. HadEX3 has some indication of increased hot extremes in this region, however it does not show the dramatic heat.

The extreme heat had a significant impact on various aspects of European societies, including agriculture, health, and daily life. In fact, the effects of the heatwave on the population are even apparent in demographic data, especially in the death rate. For instance, a demographic study36 estimates this heatwave to have caused more than 40,000 deaths in France, mostly infants and seniors. While not entirely correlated, it is possible to use these data as a proxy to uncover indirect evidence of a large number of warm days and nights. The French map of increased senior deaths for 1911 obtained from 37, 38 and shown in Supplementary Figs. S14, S15 exhibits strong similarities with the high TX90p and TN90p values predicted by the CRAI models in this region. The SROCC calculated using the normalized maps indicate a significant spatial correlation between the CRAI reconstruction and the mortality increase for both indices: SROCC = 0.39 for TX90p and SROCC = 0.51 for TN90p. That is not the case for the other datasets, which present inconsistent spatial correlations depending on the index (TX90p or TN90p).

The coldwave of 1929 was a severe winter weather event that impacted various parts of the world, including North America and Europe. The cold spell occurred during the winter of 1928–1929 and was particularly intense during January and February of 1929. For instance, in Fig. 7, the CRAI reconstruction shows that many central European countries experienced daily minimum temperatures below the 10th percentile for most of the days in February (over 80% of the days). In contrast, normal temperature conditions (with respect to the base period) are observed for northern Africa, Spain, Ireland and northern Norway. Weather reports of the Meteorological Magazine in UK39 describe significant variations of the cold wave in the British Isles: “In many parts of England, February 1929, will be remembered as the coldest February experienced since 1895: in Scotland the conditions were less extreme, while in western Ireland the mean temperature was slightly above normal”. This is more accurately reflected in the spatial distribution of TN10p values in the CRAI reconstruction and the 20th Century Reanalysis (see Supplementary Fig. S18) than in HadEX3. These regional differences are also consistent with the distribution of mean sea-level temperatures published in the Monthly Weather Report of the Met Office UK for the year 1929 (see Supplementary Fig. S19). It is also possible to find evidence of the cold winter spell in non-scientific journals such as the Thurgauer Jahrbuch40, which describes the impact of the cold temperatures in Europe and, more particularly, on the population from the canton of Thurgau (Switzerland). Interestingly, the journal mentions an unusual warmth in Scandinavia, where “farmers were able to cultivate their fields”. This statement corroborates the low TN10p and TX10p values reconstructed by the CRAI models in Scandinavia (Fig. 7 and Supplementary Fig. S16) that are not clearly visible in HadEX3.

Methods
HadEX3 and HadEX-CAM
The HadEX37,8 dataset presents indices calculated from daily precipitation accumulations as well as maximum and minimum temperatures derived from over 30,000 weather stations. Index values are calculated for each station, with the data subsequently undergoing quality control checks before being interpolated onto a regular latitude-longitude grid using the ADW method9. For full details of the construction of HadEX3, see ref. 8 and references therein.

The ADW gridding applied in HadEX3 is useful as the method can extrapolate the sparse station data to provide relatively contiguous fields of these index values. However, for the purposes of investigating the use of AI to perform infilling, the available HadEX3 data files are not appropriate, having already been infilled using the ADW approach. Therefore an alternative dataset, which uses the same input data but a different gridding method that does not interpolate, was created in order to be reconstructed by the AI algorithm. The parameters of the 1.875°×1.25° longitude-latitude grid were kept the same as in HadEX3.

The four percentile-based temperature indices (Table 1) are calculated from the in situ data using the Climpact2 software following the HadEX3 methodology42, 43. However, to convert the in situ indices into a sparse-gridded form without doing any interpolation, we use the Climate Anomaly Method (CAM)44 after calculating the indices at each station. This gridding method was used by ref. 45 as one of their selection of gridding methods when assessing structural uncertainties in a dataset of gridded extremes indices, and has a long history of being used in the creation of gridded datasets of global surface temperature anomalies (e.g., refs. 19, 44).

The CAM requires station time series in anomaly form (i.e., the difference from a climatology calculated over a reference period) and then at each timestep, performs a simple average of all stations within the grid box with valid data. Valid data from only a single station within a grid box is necessary to result in a value being set for the grid box. The indices in this work already use a reference period when calculating the percentile thresholds, so we use the same reference period (1981–2010) to calculate climate anomalies.

After applying the CAM, the set of gridded anomalies needs to be converted back to actual values by adding on the climatology, both for the AI algorithm and to be able to compare the final grids with e.g., HadEX3 or ERA546, 47. During the early stages of this work, adding the grid box climatology as calculated from the station data series themselves resulted in values fractionally <0% or >100% times for some grid boxes, which are clearly unrealistic. These arose because of data completeness requirements that are imposed on the calculation of these indices for the station data by Climpact2, which result in missing months in the station time series of indices if more than 3 daily values are missing. When calculating these 4 indices, the exceedance above the percentile threshold is derived from the daily data. However, when missing data criteria result in missing months during the reference period, the climatological value can be slightly different from the theoretically expected value of 10% (36.5 days). In cases where there is a single station in a grid box, this causes no problems. But for grid boxes where multiple stations are combined and for months where one or more of these stations have missing values, if the remaining have particularly negative anomalies in that month, adding the box average climatology can result in values <0%. Similarly for cases where particularly high anomalies are present, resulting in values >100%. Therefore, we use the theoretical value of 10% as the climatology in all cases, and furthermore, mask values that fall outside of the range 0≤v≤100% in case of any residual floating point issues (though none were found during testing).

Further quality controls were conducted subsequent to the computation of the indices. These checks facilitated the identification of anomalous data attributed to values measured in Lugano’s station (number 100750 in ECAD, 48, 49, the source of most of the data over the European region). In order to preserve the integrity of HadEX-CAM, the data associated with this station have been excluded across the entire temporal range.

Deep learning-based reconstruction method
The methodology employed to reconstruct the extreme indices is based on the work by ref. 15 and is referred to as CRAI. CRAI is an implementation of a deep-learning algorithm that makes use of a U-Net architecture50 and partial convolution layers20. As shown by ref. 20, partial convolution layers are more suitable for the reconstruction of large and irregular regions of missing values (as in the HadEX-CAM dataset) than standard convolution layers.

The neural network takes two types of inputs: the maps of extremes indices (e.g., from HadEX-CAM), and the masks of missing values which are extracted from the HadEX-CAM dataset. For the training, we used the stochastic gradient descent approach by loading batches of 16 samples taken from simulation data (see Section Dataset up for the training and the evaluation) with randomly selected masks of missing values. The masks from HadEX-CAM are used to create artificial missing values in the training dataset, which is originally complete. The inputs are propagated through the neural network following the mask-update procedure described in ref. 20, and the mean-absolute-error loss function is applied to the regions of missing values (land only) by comparing the output with the original complete data. In order to constrain the output to be in the interval of permitted values (from 0% to 100%), we apply a rescaled sigmoid function after the last layer. The weights and biases of the layers are updated after each iteration by calculating the gradients through backpropagation.

For all the calculations, the maximum number of iterations has been set to 1 million and the learning rate to 5×10−5. The optimal bias-variance trade-off of the selected models is confirmed by analyzing the training and validation loss values over iterations. The number of encoding/decoding layers and the number of output channels for each layer have been determined by performing a hyperparameter search on a regular grid of 18 hyperparameter values for the TX90p index. To account for the inherent stochasticity of our method, we have trained ten models for each configuration of hyperparameters and selected the optimal configuration based on the best validation RMSE average across these 10 runs. For the sake of simplicity and given that such a hyperparameter search is very demanding in terms of computing resources, we have chosen the same optimal configuration of hyperparameters for all extreme indices (see Supplementary Fig. S1). To improve predictive accuracy, twenty separate models with identical configurations were trained for each extreme index. The resulting predictions were then averaged to produce a single outcome. As shown in Fig. 2 and Supplementary Figs. S6, S12 and S17, this ensemble of predictions can be used as well to estimate the uncertainty of the model.

Dataset up for the training and the evaluation
The reconstruction of the four monthly extreme indices under consideration can be done by following two strategies:
• the reconstruction of the daily maximum (minimum) temperature fields prior to the calculation of the TX90p and TX10p (TN90p and TN10p) indices.
• the reconstruction of the monthly TX90p and TX10p (TN90p and TN10p) indices calculated from the incomplete daily maximum (minimum) temperature fields.

In this study, we have opted for the second strategy based on two main considerations. Firstly, the daily frequency of the data can result in highly sparse spatial distribution, potentially complicating the AI model’s reconstruction task. Secondly, the HadEX3 dataset is derived from a mixture of daily station data and pre-calculated indices. Therefore, reconstructing the extreme indices ensures alignment with the data employed to create the HadEX3 dataset, maintaining consistency in our approach.

Given the limited number of samples and the incomplete nature of the monthly HadEX-CAM dataset, it is required to employ a transfer learning methodology that incorporates data from Earth system models. For the training, the validation, and the evaluation of the CRAI and diffusion models, we are using the data from 45 historical forcing simulations from 8 CMIP6 models51–58 at a monthly temporal frequency spanning the 1901–2014 period. From all the models of the CMIP6 archive, only those whose spatial resolution is higher than the one of the HadEX-CAM dataset (1.875°×1.25° longitude-latitude) have been retained. The details of the models and members used in this study are given in Supplementary Table S2. The extreme indices are calculated for each member of each model by using the ETCCDI implementation in CDO59. Following the recommendation of ref. 60, the calculation of the extreme indices is performed first, and the remapping to the HadEX-CAM grid is performed in a second step. The resulting remapped data is randomly split into a training, a validation, and a test set. The random split approach is justified by the relatively low temporal autocorrelation between consecutive months of extreme indices while allowing the model to learn from the entire timespan, encompassing periods with and without significant climate change effects. The total number of samples (number of samples per month) is respectively 50616 (37), 9576 (7), and 1368 (1).

To assess the quality of the trained models, we are also making use of two reanalysis datasets: ERA546, 47 and 20CRv335. The extreme indices are derived from these two datasets using the same procedure as for the CMIP6 models. For the quantitative evaluation, we primarily rely on the extreme indices calculated from the ERA5 dataset. This choice is motivated by the Taylor diagrams61 presented in Supplementary Fig. S2 that compare the standard deviation and SROCC of four reanalysis datasets (ERA5, 20CRv3, JRA-5562, and MERRA-263) with the HadEX3 and HadEX-CAM datasets. These diagrams were prepared as outlined in ref. 64, which also presents a detailed comparison of reanalysis datasets with HadEX3. Supplementary Fig. S2 reveal that ERA5 aligns most closely with the target data, making it the most suitable reanalysis dataset for evaluating the reconstruction of the HadEX-CAM dataset. While ERA5 is one of the most advanced atmospheric reanalyzes currently available, it has limited temporal coverage (1940–2018). To assess the CRAI models for earlier time periods, we then resort to 20CRv3, which is the only reanalysis dataset analyzed here, covering the early 20th Century (see Supplementary Table S1).

Comparative methods
To assess the quality of the CRAI reconstructions, it is advisable to establish a comparison with other approaches. As a baseline, we have opted for two statistical methods commonly employed for the interpolation of geospatial data: Inverse Distance Weighting (IDW) and ordinary Kriging. The IDW technique stands out as the simplest method, as it estimates the missing values by weighting the contribution of the surrounding valid values with weights that are inversely proportional to their distances. The implementation of this method is straightforward and requires defining a single parameter, the power parameter, which has been set to 4.

The ordinary Kriging technique65 is more sophisticated as it makes use of a variogram analysis to incorporate more spatial correlation in the reconstruction. For the present study, we have created our own implementation of the statistical method, based on the PyKrige Python library66. A common set of hyperparameters has been used for all the reconstructions presented in this article. The variogram model (exponential) and the number of averaging bins (200) have been determined by performing a hyperparameter search on the test dataset.

In addition to these statistical methods, we have also reconstructed the fields of extreme temperature indices using diffusion models67, 68, which are state-of-the-art generative AI models. For this purpose, diffusion models have been trained for the reconstruction of each extreme index using a modified version of the guided diffusion code released at https://github.com/openai/guided-diffusion. Inspired by techniques from non-equilibrium thermodynamics69, diffusion models involve a forward diffusion process where noise is progressively added to the input data until it becomes pure noise. A U-Net is then trained to denoise the noisy data, integrating time embeddings to condition the model on specific timesteps that represent the level of noise added in the forward diffusion process. During the training, the model minimizes the difference between predicted denoised maps and the original maps using a L1 loss function. The results presented here are obtained with models trained using the same configuration of hyperparameters for all extreme indices. In particular, the models are using 2000 diffusion steps, a U-Net made of 3 encoding and 3 decoding layers with 128, 256, 512, 256, 128, 1 output channels respectively, a batch size of 16 samples, a learning rate of 10−5 and a total of 500,000 iterations. This configuration was found to give the best RMSE and SROCC results among 12 configurations evaluated during a hyperparameter search on the test datasets. The trained models were applied twenty times on each evaluated dataset in order to obtain a comprehensive sample of predictions. These predictions were then averaged to produce a single outcome for each index and dataset.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "I have attached my review as a PDF.", "reviewer_2": "This paper introduces a novel AI-based methodology for reconstructing historical observations of climatic extremes (warm and cold days and nights) over Europe from 1901 to 2018. The AI model leverages transfer learning using Earth system model data from CMIP6 to enhance the coverage and accuracy of the observational dataset HadEX-CAM.\nThe key points, based on my understanding, are:\n1) Compared to existing interpolation method, such as Angular Distance Weighting and Kriging, the AI model has higher accuracy in terms of regional variability and event detection, and shows promise in identifying complex spatial trends over the extensive 1901-2018 period.\n2) Compared to reanalysis, such as ERA5, the AI model is cheaper, and can extent to early 20th century when reanalysis is missing. It helps to reveal the 1911 heatwave event and the 1929 cold spell event with more clarity.\nThe drawback of this manuscript, based on my understanding, are:\n1) The partial convolution neural network methodology is not new by the Year 2024. Many works, including the first one from a contributing author here, have explored this methodology for reconstructing historical climate dataset. The authors declared that this methodology is superior to GAN. First, GANs can use partial CNN as backbone. These two are not comparable in this sense. Second, GANs are probabilistic method, which can yield uncertainty of the inpainting results. In principle, it is superior to the methodology proposed here. If the authors would like to stress the superiority of their methodology, Angular Distance Weighting and Kriging are too weak benchmarks. GANs could be good. Probabilistic diffusion models are now state-of-the-art for inpainting tasks.\n2) The problem setup here is too convoluted. If the authors would like to make use of the underlying low-dimensional statistical structure of climate data for the inpainting task, it is better to apply the raw climte variables, instead of derived statistics, as the former enjoys more spatial consistency.\n3) The presenting style here is more like a technical report, which might not benefit a broader readers of Nature Communications.\nIn summary, this AI-based climate reconstruction method is not very new, and does not offer significant contribution to the community. Some minor comments: for the title, we can not reconstruct gap, but reconstruct records. A line-by-line comment is not given, as there is no line number or truly exciting findings in the manuscript.", "reviewer_3": "The paper describes an AI method to gap fill the HadEX3 data set for extreme indices over Europe. The authors use an inpainting approach, which they apply in a transfer learning configuration. As is to be expected (and as is found in many other papers), the AI approach is statistically better than the alternative statistical approaches (Kriging) studied. However, I would appreciate a bit more discussion on what the actual goal and use of the data set is rather than just statistical arguments. Only then can the AI gap filling be appreciated. Otherwise the paper is well written. Below are some comments.\nMajor comments\nThe paper shows different analyses that could be done with the gap-filled data set. For instance, trends can be calculated form complete records rather than from incomplete ones. Specific heatwaves or cold waves can be studied (1911, 1929), or climate impacts can be attributed (Fig. S11). I agree that calculating trends from infilled data sets is more convenient (and, since it is infilled anyway, it might be better done using a good method). Here one has to have the confidence that AI infilling does not change the statistical properties, which of course holds also for other methods. So in this case I can see the benefit of AI infilling. However, for the heatwaves or cold waves shown, I would, as a scientist, prefer reanalysis products over a gap-filled HadEX3 data set. And the same holds for the attribution of mortality – it is dangerous to assume that the (unobserved) pattern of the heat wave should be close to that of mortality. Also, I am not sure we can really get a better physical understanding of climate extremes from the AI infilled data set, as implied in the abstract (no physical constraints are built in, as far as I understand).\nAs I understand, inpainting preserves the available original grid cell values, while HadEX3 does not. This should be discussed. Why does HadEX3 not preserve the grid cell value? Probably there are good reasons for that, but they should be mentioned (the paper is only submitted). This also means that there are different strategies in terms of what the data set should be (not just the method). This of course also raises questions of comparability (this may not be a big issue, the authors have quite a good concept by comparing their infilled data also with a Kriging approach and with reanalyses). It would just be nice to read why the data set should or should not preserve the observations.\nOn the approach: I am not an expert in AI, so cannot comment much. However, what would be relevant is to get a feeling of the number of gaps. This is in Fig. S2, but this should also be mentioned in the text. How does the model performance vary as a function of gaps?\nCan confidence intervals be calculated for the infilled data?\nFig. S11 is dangerous – what do we expect? Higher mortality in areas with more heat, but also in more urbanized areas? And is it really the days above the 90th percentile or would an absolute threshold be more appropriate? Why do the author show TX90p and not TN90p? This is a clear case for which I would prefer a reanalysis.\nVisually it becomes obvious that the AI infilling produces fields that are less smooth then HadEX3, and sometimes they even appear patchy. The authors speculate that this is more realistic and refer to ERA5 in the supplement. However, looking at the supplement, at least visually, it seems like ERA5 is clearly smoother than the infilled product. In any case, rather than just speculating whether or not the spatial noise is realistic: Could you quantify that?\nIn order to get a realistic spatial variability, could you include the covariance matrix into the loss function such as to train the correct smoothness?\nMinor\nIs there any way AI could deal with observation errors?\nOn p. 2 the authors write “the choice of one approach or the other hinges on whether pixel-level accuracy (CNN) or physical realism (GAN) is more relevant for the task. In this study, we are opting for a CNN-based approach” Does this mean that the authors aim for high pixel-level accuracy? Why?\nSeveral figures in the Supplement show 20CRv3, please mention this also in the text (the reader keeps wondering why 20CRv3 is not used in these cases – there is no mention in the text).\nIn the conclusions you mention “higher spatial variance”? Is that good or bad?\nIn the conclusions, please also mention reanalyses.\nTable 1: These are statistics based on monthly values, and the unit of RMSE is %, right? Can you give some additional info (e.g., the standard deviation) that would allow the reader to judge whether this is good or bad.\nCan you add ERA5 to Fig. 3?\np. 3: “Finally, unlike HadEX3 where calendar months and annual values are kept separate, here the monthly indices are interleaved into chronological order to form HadEX-CAM”. I do not understand this sentence.\np. 3: I am not sure what you then really do about values <0% and >100%. The authors argue that this has to do with missing data and say that “for these indices, the exceedance above the percentile threshold is calculated from the daily data.” But this is true for all indices.\nFig. S1a,b,c,d: say which is which"}, "author_response": "To all reviewers\nWe thank the reviewers for critically assessing our manuscript and providing\nconstructive comments for improvements. We are pleased that all reviewers see\nmerit in our study, while also identifying some issues that affected the clarity or\ncompleteness of our manuscript. We carefully considered all the comments\nprovided and diligently revised our manuscript in response to these comments. The\nrevisions made since the previous version of the manuscript have been highlighted\nin blue in the revised main manuscript.\nTo increase the clarity of our point-by-point response to each comment from the\nreviewers (see below), we adopted the following color-coding system:\n● The original comments and remarks from the three reviewers are shown in\ndark grey\n● Our answers are shown in black\n● The changes introduced in the revised version of the manuscript are shown in\nblue\nOur more substantial revisions can be summarised as follows:\nIn response to a comment by Reviewer #2, we have implemented a new AI method\nemploying diffusion models (see section entitled “Comparative methods” in the\nrevised manuscript). The corresponding results have been added to Table 2, 3, S4,\nS5, S6, and S7. To avoid confusion with our original AI method frequently\nreferenced as 'AI model' or 'AI reconstructions' throughout our original manuscript,\nwe have chosen to refer to our original method as ‘CRAI’ (standing for Climate\nReconstruction AI), aligning it with the name of the implementing code. We have\nupdated the naming consistently throughout the revised manuscript.\nWe have chosen to utilise the reviewer's feedback as an opportunity to expand our\nstudy by increasing the number of trained CRAI models per index from 12 to 20.\nThe results remain nearly identical across all tables and figures, underscoring the\nrobustness of our method. Nonetheless, we opted to update all tables and figures in\nthe revised manuscript. We would like to highlight that this change has no impact\non the discussion, which remains unchanged in our manuscript.\nReviewer #1\nMajor Comments:\n1) The data section should be edited to improve clarity. In particular:\na) -After a couple reads through Section 2.1, I am still not entirely sure how the\nHadEX-CAM data set was constructed and some of the choices made when\nconstructing the data set seem poorly motivated. Please consider re-writing this to\nimprove clarity assuming the reader doesn’t have familiarity with HadEX3 etc.\nResponse: We have revised this section, and hope that our changes have made our\nsteps and motivation clearer.\nb) This section is really bogged down with jargon. Perhaps make a clear list of the\nvarious data sets and software you are referencing here (e.g. HadEX2, HadEX3,\nHadCRUT3, HadEX-CAM, Climpact2, ERA5).\nResponse: The inclusion of HadCRUT3 & HadCRUT4 where to indicate the heritage\nof our gridding approach, and are not used anywhere in this study. And similarly\nHadEX2 was mentioned only to show where the CAM gridding method was used in\nan investigation into structural uncertainties. We have rephrased this part to remove\nexplicit mention of these datasets in the revised manuscript but retain the\nreferences.\nc) Section 2.1, paragraph 3, “to do this, monthly indices are calculated…”: are these\nmonthly indices of temperature or percentile indices?\nResponse: We have clarified this statement with “The four percentile-based\ntemperature indices…”\nd) Section 2.1, paragraph 3, “in-situ indices into a gridded form”: perhaps use\n“sparse gridded form”?\nResponse: We thank the reviewer for suggesting this change. It has been integrated\ninto the revised manuscript.\ne) Section 2.1 paragraph 4: I am confused by the last sentence, please clarify\nwhat you mean.\nResponse: Upon review, we found the corresponding sentence was not essential\nfor the understanding of the work and datasets. Therefore, we have removed it to\nenhance the manuscript's conciseness and readability.\nTo provide clarification for the reviewer: the HadEX netCDF data files have 13 fields\nfor indices available at monthly and annual timescales. One contains annual values,\nand the remaining 12 are for the calendar months. So that the January fields are\nJanuary 1901, January 1902….January 2018. As the AI algorithm required\nchronological data fields, these were constructed by interleaving the 12 fields\ntogether.\nf) Section 2.1 paragraph 5, 1st sentence. What is an “actual” value?\nResponse: In short, not an anomaly (actual = anomaly + climatology). We now\nstate this more clearly in the sentence: “These need to be converted back to\n\"actual\" values by adding on the climatology as calculated over the reference\nperiod…”\ng) Section 2.1 paragraph 5, \"during the early stages of this work…\": I am confused\nabout why you would do this. Shouldn't anomalies be computed on the temperature\ndata and percentile indices calculated on that? Then the indices would be bounded\nto 0-100%?\nResponse: These four indices have been defined by the WMO to be calculated\nfrom the daily temperatures, rather than the anomalies. Hence, the Climpact2 code,\ndeveloped by the WMO ET-SCI takes in daily values rather than anomalies.\nIn the case of the four indices in this work they have a theoretically expected value\nof on average 36.5days/year during the reference period (as the 90th/10th percentile\nwould be exceeded/subceeded on 10% of the days on average). In practice this is\nnot exactly achieved at a grid box level where grid boxes contain more than one\nstation because of missing data criteria.\nThe 10th percentile thresholds are set at the daily level for each station, and\ncalculated using a 5 day window centred on the day in question and a 30-year\nreference period. However, if there are more than 3 days missing in a specific\nmonth (including during the reference period), then the monthly value for the index\nis not calculated on a specific station. This can result in different numbers of\nstations contributing to the grid box values at different time steps.\nHence, for a station with complete data over the reference period, the climatology\nover this period will be close to 36.5 days (10%*365). However, in the case of a few\nmonths missing because of missing daily data (but still sufficient for the thresholds\nto be calculate), and the missing months were cooler than average, then the\nclimatology as calculated from the remaining months for this station will be biased\nwarmer and hence no longer equal to 36.5 days (the exact direction will depend on\nwhether it is an index from the 10th or 90th percentile). When combined with other\nstations in the grid box, the resulting grid box average will also be affected if the\nclimatology is calculated from the values of the monthly indices themselves.\nAs an example below:\nStation 1: No missing days in the record, climatology of 36.5d, Value for example\nmonth, 0d\nStation 2: No missing days in the record, climatology of 36.5d, Value for example\nmonth, 0d\nStation 3: Missing days in the record: Missing months in the climatology, hence\nvalue of 35d, Value for example month: Missing\nAverage climatology for grid box = (36.5+36.5+35)/3 = 36d\nStation 1 anomaly = -36.5\nStation 2 anomaly = -36.5\nStation 3 anomaly = N/A\nAverage anomaly for grid box = -36.5\nHence grid box actual value = -0.5d\nh) Please note the total size of the data set (how many monthly realizations are there\nin total and what is the date range?) I think this is in the abstract but it should\nprobably be in the data section too.\nResponse: The number of monthly realizations for the training, validation and test\nsets is given in the section entitled “Data setup for the training and the evaluation”:\n“The total number of samples (number of samples per month) is respectively 50616\n(37), 9576 (7), 1368 (1).” Following the suggestion of the reviewer, the total number\nof samples and the time span for each dataset considered in this study has also\nbeen added to Table S1 of the Supplementary information.\ni) The ordering of introducing the application data set, then model description, then\ntraining data set, and then introducing ERA5 as an additional testing data set but\nnot in the data section, is a bit clunky.\nResponse: The structure of the sections has been changed in the revised\nmanuscript and the reanalysis datasets are now introduced in the Section “Data\nsetup for the training and the evaluation”.\nj) It could be helpful to add a quick road-map to the start of section 2 outlining each\nof the main data sets and their roles in the study before jumping in to detailed\ndescriptions (HadEX-CAM: data set the CNN will be deployed on, CMIP6: training\nvalidation and testing data, ERA5: additional testing data set).\nResponse: We thank the reviewer for suggesting this modification. We have\nrestructured the revised manuscript to clarify the description of the employed\ndatasets and added Table S1 in the supplementary information to give an overview\nof their use in the study.\n2) Regarding the CMIP6 data sets. Did you make an effort to account for temporal\nauto-correlation when selecting your train/validation/test split? Selecting samples\nrandomly can lead to validation samples with very similar counterparts in the\ntraining set for a data set like this. Ultimately, the model does very well on unseen\nERA5 testing data so it does not seem that an over-fit has occurred, but perhaps\nthis is something to address in the text.\nResponse: We thank the reviewer for bringing this important point to our attention.\nWe acknowledge the potential issue of temporal auto-correlation when creating the\ntrain, validation, and test datasets by randomly selecting the simulation data from\nthe CMIP6 archive. As described in the methodological section, for each iteration of\nthe training process, we are using 16 randomly selected training samples and\npairing each one with a randomly selected mask of missing values. By doing so, we\nensure that the input data used for training are most likely different from those used\nto compute the validation loss, further reducing the risk of overfitting.\n3) How did you choose the hyper-parameters of your U-Net model? It appears to\nhave far fewer channels then a typical U-Net.\nResponse: We thank the reviewer for giving us the opportunity to clarify this part.\nThe hyper-parameter search has been performed using a fixed learning rate (5x10-5)\nand a maximum number of iterations of 500000 with the aim to find the best\nconfiguration of hyperparameters in terms of RMSE. During this search, a total of 18\nconfigurations of hyperparameters have been explored (a typographical error stating\n180 configurations in the original manuscript has been corrected in the revised\nmanuscript). These configurations correspond to all possible combinations of the\nfollowing hyperparameter values: encoding layers (3, 4, 5), additional pooling layers\n(0, 3), number of output channels in the bottleneck layer (128, 256, 512). Due to the\nnon-deterministic nature of our AI method, the calculation has been repeated 10\ntimes for each configuration of hyperparameters. The best configuration has been\nselected by looking at the best RMSE averaged across the 10 runs. Such a\nhyperparameter search is very demanding in terms of computing resources and has\nbeen carried out for the TX90p index only (it has been clarified and added in the\nrevised manuscript). Therefore, we assumed that this configuration of\nhyperparameters is also suitable for the other indices under consideration.\nMinor comments:\n1) Section 2.4: Why not show the results using the inverse distance weighting\ntechnique as well. Isn’t this just an additional row in Tables 1 and 2?\nResponse: We have added the results of the inverse distance weighting method in\nall the tables (Table 2, 3, S4, S5, S6, S7)\n2) Tables 1 and 2: What are the units for RMSE?\nResponse: We have added the units (%) of the RMSE in the caption of the tables\n3) Figure 1: This needs a much better caption. Please try to provide more detail\nabout the network architecture here. Also the shading in the diagram is not\nexplained.\nResponse: We have extended the caption of Fig. S1 (previously Fig. 1) to give more\ndetails about the network architecture\n4) Section 2.2: “high bias and variance… … as a function of iterations.” is a very\nawkward sentence and “discarded” may not be the intended meaning. Please\nreword.\nResponse: This sentence has been rephrased in the revised manuscript: “The\noptimal bias-variance trade-off of the selected models is confirmed by analysing the\ntraining and validation loss values over iterations.”\n5) Page 6, last paragraph, second sentence. Please be more specific that you are\nreferring to a time-series of difference in RMSE between the two methods.\nResponse: This sentence has been completed to avoid any possible confusions:\n“...the time series in Fig. 1a, obtained by computing the mean of the RMSE\ndifference across the spatial domain...”\n6) Page 7, last paragraph, first sentence is a bit awkward. Perhaps change to: “To\nevaluate our model in a context approximating our intended use case,…”. Mostly,\n“target data” is probably not the best term to use here.\nResponse: We thank the reviewer for the suggestion. We have modified the\nsentence using the suggested change.\n7) Regarding Fig 3.: You might note somewhere that the spatial average can cause a\nsignificant portion of the pixel level error to cancel, which leads to fairly high\nagreement between these lines. Said differently: the mean squared error is larger\nthan the squared difference in means.\nResponse: We agree with the point outlined by the reviewer and we expect the\nreader to know that spatial averaging can mitigate some of the pixel-level\ndiscrepancies. This aspect is taken into consideration in the subsequent sections of\nour manuscript as we extend our analysis to spatial analysis.\n8) Figures 5 and 6: The color bar is labeled with “days/10 years” while the caption is\n“days/decade.” Best to choose one.\nResponse: We thank the reviewer for noticing this mismatch. The captions of Fig. 4\nand 5 (previously Fig. 5 and 6) have been changed accordingly\n9) Page 12, paragraph 1, “casualties” has different meanings depending on field.\nResponse: “casualties” have been replaced by “deaths” in the revised manuscript\n10) Throughout the manuscript opening quotes are facing the wrong way.\nResponse: We thank the reviewer for noticing these typos. We corrected all of them\nin the revised manuscript\n11) Section 2.2, paragraph 2, “see section 2.2.3”: there is no section 2.2.3\nResponse: We thank the reviewer for noticing this mistake. It has been corrected in\nthe revised manuscript.\nReviewer #2\nThe drawback of this manuscript, based on my understanding, are:\n1) The partial convolution neural network methodology is not new by the Year 2024.\nMany works, including the first one from a contributing author here, have explored\nthis methodology for reconstructing historical climate dataset. The authors declared\nthat this methodology is superior to GAN. First, GANs can use partial CNN as\nbackbone. These two are not comparable in this sense. Second, GANs are\nprobabilistic method, which can yield uncertainty of the inpainting results. In\nprinciple, it is superior to the methodology proposed here. If the authors would like\nto stress the superiority of their methodology, Angular Distance Weighting and\nKriging are too weak benchmarks. GANs could be good. Probabilistic diffusion\nmodels are now state-of-the-art for inpainting tasks.\nResponse: We thank the reviewer for thoughtful comments regarding the\nmethodology used in our study and the comparison with alternative approaches.\nWe acknowledge that partial CNNs have been used since 2018 and that it is not the\nmost recent deep-learning methodology available. However, we do not think that\nthe novelty of the approach is necessarily correlated with its efficiency for the\nintended use case in this study. In fact, many of these methods (such as the GANs\nor the probabilistic diffusion models mentioned by the reviewer) have been\ndeveloped in the field of computer vision for tasks which are quite different from the\nreconstruction of missing values in climate data. In particular, a great effort has\nbeen made in the last years to develop generative methods able to produce diverse\nand realistic images for applications such as text-to-image generation. While\ngenerating visually appealing outputs, these methods have not been developed to\nbetter represent the underlying physical processes associated with climate data and\nwe are not expecting them to produce better RMSEs than partial CNNs. That is the\nmeaning of the sentence “the choice of one approach or the other hinges on\nwhether pixel-level accuracy (CNN) or physical realism (GAN) is more relevant for\nthe task” from the original manuscript and that is based on the conclusions of Geiss\net al., 2021 that applied these methods to climate data. However, we would like to\nclear any doubts in this regard, and for that reason we have now also adapted a\nguided-diffusion code released by OpenAI in order to train diffusion models to\nreconstruct the climate datasets presented in our study. The methodology and the\nresults have been added to the revised manuscript. The diffusion models employed\nin this study make use of a common configuration of hyperparameters that was\nfound to give the best RMSE and SCORR results among 12 configurations\nevaluated on the test datasets for TX90p. The evaluation presented in Table 2, 3, S4,\nS5, S6, S7 shows that the diffusion models are performing worse compared to the\npartial CNNs for most of the datasets, extreme indices and evaluation metrics,\nincluding the new ones (R2 score and Wasserstein distance).\n2) The problem setup here is too convoluted. If the authors would like to make use\nof the underlying low-dimensional statistical structure of climate data for the\ninpainting task, it is better to apply the raw climte variables, instead of derived\nstatistics, as the former enjoys more spatial consistency.\nResponse: Reconstructing the raw climate variables as suggested by the reviewer\nis a possible strategy that we have considered at the beginning of the work and that\nis mentioned as a possible extension of the work in the conclusions of the original\nmanuscript. However, we have opted for another strategy that consists in\nreconstructing the extreme indices derived from the incomplete temperature fields.\nThis decision is based on two main considerations detailed in the revised\nmanuscript:\n“Firstly, the daily frequency of the data can result in highly sparse spatial\ndistribution, potentially complicating the AI model's reconstruction task. Secondly,\nthe HadEX3 dataset is derived from a mixture of daily station data and\npre-calculated indices. Therefore, reconstructing the extreme indices ensures\nalignment with the data employed to create the HadEX3 dataset, maintaining\nconsistency in our approach.”\n3) The presenting style here is more like a technical report, which might not benefit a\nbroader readers of Nature Communications.\nResponse: We thank the reviewer for giving us the opportunity to improve the\nformat of the article. We have reordered the sections and adapted the text in the\nrevised manuscript to make it more suitable for the broader readership of Nature\nCommunications.\nSome minor comments: for the title, we can not reconstruct gap, but reconstruct\nrecords.\nResponse: We thank the reviewer for bringing this wording issue to our attention.\nWe have corrected the title following the suggestion of the reviewer: “Artificial\nIntelligence Reveals Past Climate Extremes by Reconstructing Historical Records”\nReviewer #3\nMajor comments\nThe paper shows different analyses that could be done with the gap-filled data set.\nFor instance, trends can be calculated form complete records rather than from\nincomplete ones. Specific heatwaves or cold waves can be studied (1911, 1929), or\nclimate impacts can be attributed (Fig. S11). I agree that calculating trends from\ninfilled data sets is more convenient (and, since it is infilled anyway, it might be\nbetter done using a good method). Here one has to have the confidence that AI\ninfilling does not change the statistical properties, which of course holds also for\nother methods. So in this case I can see the benefit of AI infilling. However, for the\nheatwaves or cold waves shown, I would, as a scientist, prefer reanalysis products\nover a gap-filled HadEX3 data set. And the same holds for the attribution of\nmortality – it is dangerous to assume that the (unobserved) pattern of the heat wave\nshould be close to that of mortality.\nResponse: We thank the reviewer for the detailed comments. In particular, we\nappreciate the understanding of both the benefits and the limitations of infilling\nclimate data with AI methods. We agree that while AI infilling can offer significant\nadvantages compared to traditional statistical methods, it is crucial to remain\ncautious about potential changes in statistical properties. To further demonstrate\nthe robustness of our AI method in this regard, we have included new tables in the\nrevised manuscript (Table S4, S5, S6, and S7). These tables present results from\nadditional evaluation metrics that measure the variance (R2 score) and the\ndistribution (Wasserstein distance) of the data. They show that our AI method\nperformed better overall compared to the other methods for these new metrics as\nwell.\nWe understand that reanalysis datasets may be preferred to gap-filled HadEX3 or\nHadEX-CAM dataset for specific analysis. However, reanalyses also have their\nlimitations (even in the modern period, depending on which ECV is being studied).\nThey are also reliant on the amount of available input data, which is limited to in situ\nand sonde measurements at times earlier than the satellite period. Furthermore, only\nfew reanalysis products cover the early 20th Century (20CRv3 included in this work\nis one of a handful), and these assimilate typically only surface pressure\nobservations and surface winds (in addition to using reconstructions of monthly\nmean sea-surface temperatures) - so that the temperature and precipitation\nextremes are primarily an output of the (circulation-constrained) atmospheric model\nused in the reanalysis. Some earlier work has exposed substantial uncertainties also\nin these long-term reanalysis datasets for extremes, in particular in the first half of\nthe 20th century (see e.g. Donat et al 2016: https://doi.org/10.1002/2016JD025480).\nPlease note that our intention is not to claim that our infilled dataset surpasses\nexisting reanalysis products for the analysis of heatwave or coldwave events.\nRather, we consider our infilled dataset as a complementary resource that can be\nused by the community in conjunction with other existing products. We believe that\nthe novelty of our approach, with its benefits (e.g., preserving the original\nobservations, transferring the learning from ESMs) and limitations, can provide new\ninsights in specific scenarios (e.g., in time spans that are not typically encompassed\nby reanalysis datasets). Furthermore, the present work aims at demonstrating the\ncapabilities of the method that is computationally more efficient than data\nassimilation processes and can be easily scaled to reconstruct climate data at a\nhigh temporal and spatial resolution. For instance, once the AI models are trained,\nonly 1-2 seconds and a small amount of RAM on a standard laptop are needed to\nreconstruct an entire dataset such as HadEX-CAM.\nAlso, I am not sure we can really get a better physical understanding of climate\nextremes from the AI infilled data set, as implied in the abstract (no physical\nconstraints are built in, as far as I understand).\nResponse: As pointed out by the reviewer, no physical constraints such as the ones\nused in Physics-informed neural networks are implemented in our model.\nWhile incorporating specific physical constraints could be a valuable approach for\nfuture research, we have determined that it is not essential for the scope of the\npresent study. We think that properly configured and trained neural networks should\nbe capable of independently inferring some of the underlying physical principles\nfrom the training data (which are physically consistent according to the climate\nmodel formulations). An example of this ability of our method is shown in Fig. 3.\nDespite a large number of missing values in the Mediterranean region, our trained\nmodel successfully reconstructed a complex spatial pattern in northern Africa,\nclosely matching the pattern predicted by ERA5, which is grounded in physical\nlaws. It is very different from the predictions of statistical methods such as HadEX3\nor Kriging that tends to extrapolate the high TX90p values from the Iberic peninsula\nbased solely on statistical considerations.\nAs I understand, inpainting preserves the available original grid cell values, while\nHadEX3 does not. This should be discussed. Why does HadEX3 not preserve the\ngrid cell value? Probably there are good reasons for that, but they should be\nmentioned (the paper is only submitted). This also means that there are different\nstrategies in terms of what the data set should be (not just the method). This of\ncourse also raises questions of comparability (this may not be a big issue, the\nauthors have quite a good concept by comparing their infilled data also with a\nKriging approach and with reanalyses). It would just be nice to read why the data\nset should or should not preserve the observations.\nResponse: We thank the reviewer for giving us the opportunity to further clarify this\npoint. As indicated in the section entitled “HadEX3 and HadEX-CAM”, both datasets\nimply the aggregation of station data but using different weighting methods.\nThe HadEX3 ADW routine uses inverse-distance weighting along with an angular\ncomponent when calculating the values for the grid boxes. If there are multiple\nstations contributing to a grid box, the value will necessarily be a blend of these\nstations. In HadEX3, using ADW, there need to be sufficient (>=3) stations to be\nwithin the search radius to calculate a valid grid box. The decorrelation length scale\nused with this radius, is derived from the correlation structure of the station time\nseries. For these 4 indices, the decorrelation length scale used in HadEX3 is of\norder 1000km (note that the TX10/90p / TN10/90p indices have the largest\ndecorrelation length scales across all indices within HadEX3). There are no grid box\nvalues to preserve when building HadEX3, as the creation of the dataset from the\nstation time series sets the grid box values as a weighted mean of all stations within\nthe decorrelation length scale. However, stations at some distance can influence\nthe values at a particular location, especially and most importantly, when there are\nno values in the grid box at that time, where stations up to ~1000km away can\ncontribute to the grid box. This has been the construction method for HadEX,\nHadEX2 and HadEX3 and so is an intrinsic part of these versions of this dataset.\nHadEX3 is already an interpolated product and cannot be reconstructed using our\nAI models. For this reason, we have created the HadEX-CAM that uses the CAM\nmethod for gridding, which does not interpolate (see Section “HadEX3 and\nHadEX-CAM). To maximise data availability, a single station is sufficient to set a\nvalid grid box value, but if there are more stations within a grid box (likely towards\nthe end of the HadEX3 record), these are blended with a simple average.\nWhen the AI models are applied to the HadEX-CAM dataset, they give predictions\nfor all the gridboxes, including those corresponding to the original values. These\npredictions are bias adjusted values that can be used to assess observation errors.\nHowever, in the current version of the manuscript, we have chosen to retain the\noriginal valid values in our dataset for two main reasons. Firstly, these values are\nexpected to reflect the real-world conditions more closely, ensuring more relevance\nfor climate research. Secondly, these values have been obtained using a novel\napproach and stand out as valuable information on their own that could bring new\ninsights to the community for future research. If the reviewer finds it useful, we\ncould provide the reconstructed HadEX-CAM datasets containing the bias-adjusted\nvalues as well.\nThe reviewer raises a legitimate question by mentioning the comparability of the\ndatasets. To circumvent this problem, all the metrics presented in Table 2, 3, S4,\nS5, S6, and S7 are computed for the reconstructed values only.\nOn the approach: I am not an expert in AI, so cannot comment much. However,\nwhat would be relevant is to get a feeling of the number of gaps. This is in Fig. S2,\nbut this should also be mentioned in the text. How does the model performance\nvary as a function of gaps?\nResponse: To help the reader to have a better understanding of the number of\ngap/valid values in the dataset, we have added a sentence to the caption of Fig. S4\n(previously Fig. S2) specifying the total number of valid values in the reconstructed\nHadEX-CAM dataset. We have also added a new figure to the supplementary\ninformation (Fig. S3) showing specifically the performance of the model (RMSE) as a\nfunction of time (Fig. S3a) and space (Fig. S3b). Two sentences have been added to\nthe revised manuscript to summarise the most relevant findings of the figure: “As\nexpected, the values are larger for time periods (see Supplementary Fig. S3a) and\nregions (see Supplementary Fig. S3b) with a higher prevalence of missing values in\nthe dataset (as shown in Supplementary Fig. S4).”\nCan confidence intervals be calculated for the infilled data?\nResponse: The two AI methods we have employed in this study (CRAI and diffusion\nmodels) allow for the calculation of uncertainties. In the case of CRAI, they are\nobtained by training separately multiple models using the same configuration of\nhyperparameters. As the initial conditions (random weights initialization) and\nbatches of training samples (maps of extreme indices and mask of missing values)\nare different for each calculation, we obtain slightly different results for each trained\nmodel. These differences characterise the epistemic uncertainty and give an\nindication of the stability and reliability of our models. In the original manuscript, we\ndecided to exclude model uncertainties from the figures and tables to avoid\noverloading them with information. However, in response to the reviewer's insightful\ncomment, we have reconsidered this approach and incorporated the uncertainty\nmeasures into four key figures in our revised manuscript (Fig. 2, Fig. S6, Fig. S12,\nFig. S17). In particular, Fig. S12 and Fig. S17 are two new figures that illustrate the\ndistribution of the epistemic uncertainty across the spatial domain for two indices of\nthe heatwave (Fig. 6) and coldwave (Fig. 7) events respectively.\nFig. S11 is dangerous – what do we expect? Higher mortality in areas with more\nheat, but also in more urbanized areas?\nThe purpose of the comparison presented in Fig. S14 (previously Fig. S11) is to give\nan indication on the capability of the different datasets to describe the analysed\nheatwave event or its possible impacts. Given the absence of better alternative\nmethodologies, we view this comparison as a pragmatic approach that offers\nindirect (proxy) evidence rather than conclusive proof. To underscore these\nlimitations in the interpretability of the results, we have included an additional\nstatement in the revised manuscript: “While not entirely correlated, it is possible to\nuse these data as a proxy to uncover indirect evidence of a large number of warm\ndays and nights.“\nAnd is it really the days above the 90th percentile or would an absolute threshold be\nmore appropriate?\nResponse: We agree with the reviewer that using absolute temperature thresholds\ncould be more appropriate for the comparison with the mortality increase shown in\nFig. S14 (previously Fig. S11). The TX90p and TN90p indices indicate the\noccurrence of unusually warm temperatures (relative to the climate conditions at a\nspecific location or time). An advantage of using relative threshold is that these\nindices are not affected by possible model biases or inconsistencies in orography\nacross datasets. Other indices (e.g. those measuring the intensity of the hottest\nconditions, like TXx) could be considered in the future, but these may pose different\nchallenges to the AI-based reconstruction method. We also would like to\nunderscore that the purpose of this comparison is not to establish causality\nbetween the specific temperature index values and the mortality increase, but rather\nto use the mortality increase as a proxy that seems to confirm that during this time\nthere were factors leading to increased mortality, which coincided with the hot\nconditions.\nWhy do the author show TX90p and not TN90p? This is a clear case for which I\nwould prefer a reanalysis.\nResponse: We acknowledge that both warm days (TX90p) and warm nights\n(TN90p) can lead to severe health impacts. To address this, we have added a similar\nanalysis for the TN90p index (Fig. S15) to the revised supplementary information. As\ndetailed in the revised manuscript, the findings for TN90p corroborate those for\nTX90p, revealing a strong spatial correlation between the AI-reconstructed\nHadEX-CAM dataset and the demographic data. Interestingly, the 20cr reanalysis\ndataset exhibits a low spatial correlation between TN90p and the mortality increase,\nin contrast to the high spatial correlation observed for TX90p. Regarding the use of\nreanalysis data, please see our earlier response highlighting the large uncertainties\nalso in reanalysis of these early periods of the 20th century.\nVisually it becomes obvious that the AI infilling produces fields that are less smooth\nthen HadEX3, and sometimes they even appear patchy. The authors speculate that\nthis is more realistic and refer to ERA5 in the supplement. However, looking at the\nsupplement, at least visually, it seems like ERA5 is clearly smoother than the infilled\nproduct. In any case, rather than just speculating whether or not the spatial noise is\nrealistic: Could you quantify that?\nResponse: We thank the reviewer for giving us the opportunity to clarify this point in\nthe revised manuscript. To quantify the smoothness of the fields, we have employed\nthe global Moran’s I analysis that measures the spatial autocorrelation and\ndetermines the degree of clusterization and dispersion in the data. The results of\nthis analysis have been added to Fig.3, Fig. S9 and Fig. S10 in the revised\nmanuscript and commented in the text. The values exhibit systematically closer\nalignment between ERA5 and CRAI compared to HadEX3, confirming the structural\nsimilarities between the two datasets, as outlined in the manuscript. Moreover,\nHadEX3 consistently shows the highest values, indicating smoother spatial patterns\ncompared to ERA5 and CRAI.\nIn order to get a realistic spatial variability, could you include the covariance matrix\ninto the loss function such as to train the correct smoothness?\nResponse: We acknowledge the suggestion to integrate the covariance matrix into\nthe loss function to potentially refine the model's guidance toward accurate spatial\nvariability. Based on our assessment, we conclude that such integration is feasible\nbut is however not imperative for the present study. Indeed, our current CRAI\nmodels already demonstrate robust capability in capturing spatial variability, as\nshown by the high R2 scores presented in tables S4, S5, S6, and S7 added to the\nrevised Supplementary information.\nMinor\nIs there any way AI could deal with observation errors?\nResponse: Our method also generates predictions for the original valid values of\nthe input data, which can be potentially employed to assess the observation errors.\nWhile the reviewer raises an interesting point, we believe that addressing\nobservation errors specifically would necessitate an in-depth investigation beyond\nthe scope of the current study.\nOn p. 2 the authors write “the choice of one approach or the other hinges on\nwhether pixel-level accuracy (CNN) or physical realism (GAN) is more relevant for\nthe task. In this study, we are opting for a CNN-based approach” Does this mean\nthat the authors aim for high pixel-level accuracy? Why?\nResponse: As underscored by Geiss et al (2021), the choice of one approach or the\nother depends on the intended task. While producing visually appealing\nreconstructions that can be used for illustrative purposes, GANs tend to hallucinate\nand predict spatial patterns that are climatically incoherent. For rigorous scientific\nanalysis, it is therefore preferable to use CNNs that are less prompt to\nhallucinations. As a drawback, it is often mentioned that CNNs can produce\nunrealistically smooth results. However, the global Moran’s I analysis depicted in\nFig. 3, Fig. S9 and Fig. S10 suggests a comparable spatial autocorrelation between\nERA5 and CRAI, indicating that our reconstructions do not suffer from pronounced\nover-smoothing effects.\nSeveral figures in the Supplement show 20CRv3, please mention this also in the text\n(the reader keeps wondering why 20CRv3 is not used in these cases – there is no\nmention in the text).\nResponse: We thank the reviewer for bringing this point to our attention. We have\nadded a new paragraph in Section “Data setup” which presents all the reanalysis\nproducts used in this work.\nIn the conclusions you mention “higher spatial variance”? Is that good or bad? In\nthe conclusions, please also mention reanalyses.\nResponse: To clarify our former statement, we have modified the corresponding\nsentence in the conclusion of the revised manuscript: “In addition, it maintains the\nmean field accuracy of HadEX3 while giving a more detailed representation of local\nclimate conditions, similar to modern reanalysis products such as ERA5 and\n20CRv3.”\nTable 1: These are statistics based on monthly values, and the unit of RMSE is %,\nright? Can you give some additional info (e.g., the standard deviation) that would\nallow the reader to judge whether this is good or bad.\nResponse: We thank the reviewer for bringing to our attention that additional\nstatistical information would be beneficial to assess the RMSE (and other metrics)\npresented in the manuscript. To provide additional context, we have included a new\ntable (S3) in the revised manuscript with the mean values and standard deviation of\nthe reference datasets. We have also added the unit of RMSE (%) to the caption of\nTable 2, 3 and S6.\nCan you add ERA5 to Fig. 3?\nResponse: Following the suggestion from the reviewer, we have added to the ERA5\nresults to Fig. 2 (previously Fig.3) in the revised manuscript.\np. 3: “Finally, unlike HadEX3 where calendar months and annual values are kept\nseparate, here the monthly indices are interleaved into chronological order to form\nHadEX-CAM”. I do not understand this sentence.\nResponse: We thank the reviewer for pointing out the lack of clarity in that\nsentence. It referred to how data is structured in the datasets. In HadEX3, the global\nfields for each calendar month are presented separately (fields for all Januaries, all\nFebruaries...etc), as this allows the decorrelation length scale to be calculated for\neach calendar month and to adapt to the seasonal cycle (Dunn et al. 2020).\nHowever, for the AI algorithm, a series of chronological monthly fields is required,\nand is created from the separate monthly grids (the annual field is not used in this\nwork). Upon review, we realised that the sentence in question was not essential to\nthe main focus of our study. Consequently, we have decided to remove it from the\nmanuscript to improve clarity and conciseness.\np. 3: I am not sure what you then really do about values <0% and >100%. The\nauthors argue that this has to do with missing data and say that “for these indices,\nthe exceedance above the percentile threshold is calculated from the daily data.”\nBut this is true for all indices.\nResponse: In light of the insightful comments from Reviewer #1 and Reviewer #3,\nwe have rephrased this section in the revised manuscript.\nWhat we intended to highlight with the second clause mentioned by the reviewer is\nthat in these 4 indices the threshold and exceedances are calculated from the daily\ndata, but when we process these to form climate anomalies, this is done on the\nmonthly indices, which are masked by the missing data criteria. We hope this is\nclearer in our revised text, and our reasons for masking the final grid box values\nwhich are <0% and >100%.\nFig. S1a,b,c,d: say which is which\nResponse: We thank the reviewer for pointing out this omission. The labels of each\nextreme index have been added to Fig. S2 in the revised manuscript."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thank you for the revised manuscript and detailed responses. Overall, I think the study has been greatly improved through the addition of the IDW and the diffusion model for comparison with the proposed method, and recommend publication.\nI have a couple minor suggestions you can implement if you wish:\nThe ordering of the revised manuscript is still unusual. Why not follow the typical structure of a scientific paper? (abstract, introduction, data, methods, results, discussion, conclusions).\nThe first section isn’t labeled.\nRegarding the train/validation/test split: were the CMIP6 dataset the only one used for model evaluation I would still be concerned about the ML dataset construction. This dataset is likely time-autocorrelated so randomly selecting individual samples to use for testing will lead to information leakage between the training and testing data. Holding out a specific set of years for testing or holding out one of the 45 historical forcing simulations in its entirety, for example, is probably a more appropriate way to construct a testing set for this type of problem. The authors’ response to this in the rebuttal document is that unique random dropouts of input data were used during training, which undoubtedly does help the model learn to generalize, but my main concern here is with the target data. All this said, the fact that the model has comparable performance on the ERA data that it did not see during training I think is sufficient evidence that an overfit has not occurred. The authors might consider explicitly treating this topic in the text.", "reviewer_3": "This is the revised version of a paper I have reviewed previously. The authors have done a great job and have responded to my comment adequately. For instance, they have quantified the “smoothness” of the fields, which was one of their (unsupported) arguments, and they now explain several steps of the procedure in more detail so that I now understand the procedure. The revised version of the paper is an interesting contribution to the field, which will see more AI-supported products in the future. I therefore recommend publication.\nPerhaps one comment: The IPCC gold standard is multiple lines of independent evidence. The reanalysis data set 20CRv3 and their HadEX infilled data set are, as far as I can tell, 100% independent. There are no common input data and no common methods or assumptions (reanalyses are only used in evaluation, not training, as I understand). This is the strength I see in this paper (compared also to other approaches that use anything for training), perhaps this complementary aspect is not yet highlighted enough in this paper. We need multiple data sets."}, "author_response": "To all reviewers\nWe thank the reviewers for their constructive comments and positive assessment of\nour work. We carefully considered the final suggestions provided and diligently\nrevised our manuscript in response to these suggestions.\nTo increase the clarity of our point-by-point response to each comment from the\nreviewers (see below), we adopted the following color-coding system:\n● The original comments and remarks from the three reviewers are shown in\ndark grey\n● Our answers are shown in black\n● The changes introduced in the revised version of the manuscript are shown in\nblue\nReviewer #1\nThe ordering of the revised manuscript is still unusual. Why not follow the typical\nstructure of a scientific paper? (abstract, introduction, data, methods, results,\ndiscussion, conclusions). The first section isn’t labeled.\nResponse: The ordering of the manuscript is following the Nature Communications\nguidelines and cannot be changed.\nHolding out a specific set of years for testing or holding out one of the 45 historical\nforcing simulations in its entirety, for example, is probably a more appropriate way\nto construct a testing set for this type of problem. The authors’ response to this in\nthe rebuttal document is that unique random dropouts of input data were used\nduring training, which undoubtedly does help the model learn to generalize, but my\nmain concern here is with the target data. All this said, the fact that the model has\ncomparable performance on the ERA data that it did not see during training I think is\nsufficient evidence that an overfit has not occurred. The authors might consider\nexplicitly treating this topic in the text.\nResponse: We thank the reviewer for highlighting once again this issue. We have\nadded a new sentence to the revised manuscript in order to address this point: “The\nrandom split approach is justified by the relatively low temporal autocorrelation\nbetween consecutive months of extreme indices while allowing the model to learn\nfrom the entire time span, encompassing periods with and without significant\nclimate change effects.”\nReviewer #3\nThe reanalysis data set 20CRv3 and their HadEX infilled data set are, as far as I can\ntell, 100% independent. There are no common input data and no common methods\nor assumptions (reanalyses are only used in evaluation, not training, as I\nunderstand). This is the strength I see in this paper (compared also to other\napproaches that use anything for training), perhaps this complementary aspect is\nnot yet highlighted enough in this paper. We need multiple data sets.\nResponse: We thank the reviewer for suggesting this modification. We have added\na new sentence to the revised manuscript in order to strengthen the complementary\naspect of our AI dataset: “In fact, given that our dataset is based on distinct\nassumptions and methodology, it can serve as an independent and complementary\nresource to reanalysis products.”"}]

Decision letter