2.1. Model and observational data

We use all available CMIP6 models which provide daily data for 2 m surface air temperature in the historical period. In total, this amounts to 43 models, which can be broadly grouped into 22 families by developing institutions (see Supplementary Table S3 for a full list). To represent the observations, we use four datasets selected to cover some of the diversity in observational products: (a) the fifth generation of the European Centre for Medium-Range Weather Forecasts (ECMWF) Retrospective Analysis (ERA5; Hersbach et al., Reference Hersbach, Bell, Berrisford, Hirahara, Horányi, Muñoz-Sabater, Nicolas, Peubey, Radu, Schepers, Simmons, Soci, Abdalla, Abellan, Balsamo, Bechtold, Biavati, Bidlot, Bonavita, Chiara, Dahlgren, Dee, Diamantakis, Dragani, Flemming, Forbes, Fuentes, Geer, Haimberger, Healy, Hogan, Hólm, Janisková, Keeley, Laloyaux, Lopez, Lupu, Radnoti, Rosnay, Rozum, Vamborg, Villaume and Thépaut2020); (b) the Modern-Era Retrospective analysis for Research and Applications version 2 (MERRA2; GMAO, 2015; Gelaro et al., Reference Gelaro, McCarty, Suárez, Todling, Molod, Takacs, Randles, Darmenov, Bosilovich, Reichle, Wargan, Coy, Cullather, Draper, Akella, Buchard, Conaty, da Silva, Gu, Kim, Koster, Lucchesi, Merkova, Nielsen, Partyka, Pawson, Putman, Rienecker, Schubert, Sienkiewicz and Zhao2017); (c) the Twentieth Century Reanalysis version 3 (20CR; Slivinski et al., Reference Slivinski, Compo, Sardeshmukh, Whitaker, McColl, Allan, Brohan, Yin, Smith, Spencer, Vose, Rohrer, Conroy, Schuster, Kennedy, Ashcroft, Brönnimann, Brunet, Camuffo, Cornes, Cram, Domínguez-Castro, Freeman, Gergis, Hawkins, Jones, Kubota, Lee, Lorrey, Luterbacher, Mock, Przybylak, Pudmenzky, Slonosky, Tinz, Trewin, Wang, Wilkinson, Wood and Wyszynski2021), which assimilates only surface pressure observations and, hence, relies considerably more on modeling for its output; (d) the Daily Optimum Interpolation Sea Surface Temperature dataset version 2 (DOISST; Huang et al., Reference Huang, Liu, Banzon, Freeman, Graham, Hankins, Smith and Zhang2021), which could be considered as more direct observations compared to the other datasets. The first three observational datasets provide 2 m surface air temperature also over the oceans equivalent to the models, while DOISST is based on interpolated in situ and remote sensing measurements of sea surface temperature, which is a slightly different but closely related measure of temperature.

The time period from 1982 to 2014 is used to match the first availability of all observational datasets and the end of the historical forcing period in CMIP6, respectively, and all daily temperature fields (including the high-resolution simulations) are regridded to 2.5° × 2.5° resulting in a total of 10,368 grid cells (72 latitudes × 144 longitudes). Since DOISST does not provide data on land, we apply a common land-sea mask to all datasets, using only the 6,888 ocean grid cells. Sensitivity analysis (omitting DOISST) shows that results are similar if the analysis is based on all grid cells.

In addition, we use temperature projections for the end of the century (2091–2100) driven by the high emission scenario SSP5-8.5 (Meinshausen et al., Reference Meinshausen, Nicholls, Lewis, Gidden, Vogel, Freund, Beyerle, Gessner, Nauels, Bauer, Canadell, Daniel, John, Krummel, Luderer, Meinshausen, Montzka, Rayner, Reimann, Smith, van den Berg, Velders, Vollmer and Wang2020), where available, to test the robustness of our approach under severe warming (see Supplementary Table S3). Finally, we also draw on 1 year (February 2020 to January 2021) of prototype data from cycle one of the NextGEMS project and use global temperature fields from an ICON-Sapphire run with an atmospheric resolution of 5 km (experiment ID: dpp066; Hohenegger et al., Reference Hohenegger, Korn, Linardakis, Redler, Schnur, Adamidis, Bao, Bastin, Behravesh, Bergemann, Biercamp, Bockelmann, Brokopf, Brüggemann, Casaroli, Chegini, Datseris, Esch, George, Giorgetta, Gutjahr, Haak, Hanke, Ilyina, Jahns, Jungclaus, Kern, Klocke, Kluft, Kölling, Kornblueh, Kosukhin, Kroll, Lee, Mauritsen, Mehlmann, Mieslinger, Naumann, Paccini, Peinado, Praturi, Putrasahan, Rast, Riddick, Roeber, Schmidt, Schulzweida, Schütte, Segura, Shevchenko, Singh, Specht, Stephan, Vogel, Wengel, Winkler, Ziemen, Marotzke and Stevens2022).

2.2. Statistical and machine learning classifiers

We use two different statistical and machine learning methods to separate models from observations and from each other. First, logistic regression, which allows insights into the learned coefficients but has the limitation of being a linear method. Second, a convolutional neural network (CNN) which represents rather the other end of the complexity spectrum, being able to learn nonlinear spatial relations between features but lacking the easy interpretability of logistic regression.

Logistic regression is linear in its parameters and takes an $ M\times N $ matrix as input, where $ N $ is the number of samples (days) and $ M $ is the number of features (ocean grid cells). For training an additional vector of length $ N $ is provided containing the true labels $ {y}_n $ (i.e., 0 for “model” and 1 for “observation”) for each sample. A continuous predicted probability $ {\hat{p}}_n $ is then assigned to each sample $ {X}_n $ based on its features $ {x}_{n,m} $ and the regression parameters $ {w}_m $ via the logistic function:

$$ {\hat{p}}_n=\frac{1}{1+\exp \left(-\left({w}_0+{\sum}_{m=1}^M{w}_m{x}_{n,m}\right)\right)} $$

The outcome is in $ \left[0,1\right] $ with 1 indicating the highest chance that a given day belongs to the “observation” category and at the same time the lowest chance that belongs to the “model” category since the categories are complementary:

$$ \hat{p}\left({y}_n=0|{X}_n\right)=1-\hat{p}\left({y}_n=1|{X}_n\right) $$

The predicted category $ {\hat{y}}_n $ is based on this probability $ {\hat{p}}_n $ with a decision threshold of 0.5. The logistic regression classifier is constrained by a L₂ regularization (i.e., a penalty on the squared sum of regression coefficients $ {w}_m $ ) to avoid overfitting and to ensure smooth coefficients in space. Overall coefficients $ {w}_m $ are searched to minimize the expression:

$$ {\min}_w\left\{-C\sum \limits_{n=1}^N\left[{y}_n\log \left({\hat{p}}_n\right)+\left(1-{y}_n\right)\log \left(1-{\hat{p}}_n\right)\right]+\sum \limits_{m=0}^M{w}_m^2\right\} $$

Note that in the implementation based on scikit-learnFootnote ² used in this paper, the regularization parameter $ C $ is multiplied with the residuals rather than with the coefficients, meaning that smaller values of $ C $ lead to stronger regularization. $ C $ is optimized using 5-fold cross-validation. More general details on regularized logistic regression can be found, for instance, in Hastie et al. (Reference Hastie, Tibshirani and Friedman2009).

To complement the logistic regression classifier, we use a CNN as a second method. Deep neural networks, such as CNNs, can have considerably more trainable parameters (often organized in multiple layers) and can be less interpretable than traditional methods. This means that they can be, on the one hand, more prone to overfitting but, one the other hand, can also learn more complex, nonlinear relationships in the data. Therefore, their use in different scientific disciplines, not least in climate sciences, has been rapidly increasing in recent years (Kashinath et al., Reference Kashinath, Mustafa, Albert, Wu, Jiang, Esmaeilzadeh, Azizzadenesheli, Wang, Chattopadhyay, Singh, Manepalli, Chirila, Yu, Walters, White, Xiao, Tchelepi, Marcus, Anandkumar, Hassanzadeh and Prabhat2021; Hsieh, Reference Hsieh2022). Due to their complexity, many different design choices in the exact layout of the network are possible; here, we use an out-of-the-box setup for image classification without hyperparameter tuning but adjusted to the resolution of the daily temperature maps used in the input layer. Overall, the CNN consists of an input layer, eight hidden layers, and an output layer. See Supplementary Table S2 for details about the layout.

The 2-dimensional temperature fields from each day can directly be interpreted by the CNN equivalent to a image classification task; therefore, the input layer takes a $ K\times L\times N\times 1 $ matrix where $ N $ is, again, the number of days, $ K $ is the number of latitudes, $ L $ is the number of longitudes, and $ 1 $ is a single color channel representing the temperature values. Note that for this case land grid cells are included but set to a constant fill value. The convolutional layers use the rectified linear unit as activation function, while the output layer uses the softmax operator to assign probabilities to each of classification categories $ {z}_n $ :

$$ {\hat{p}}_n=\frac{e^{z_n}}{\sum_n^N{e}^{z_n}} $$

This means that, in contrast to the logistic regression, the CNN is set up to separate also between multiple output categories (multi-class case) and a given input is identified as belonging to the class with the highest probability. For the binary case with only the “model” and “observation” classes, discussed in the first part of this paper, this means that the probabilities are complementary and the threshold for assignment is 0.5 as for the logistic regression. To avoid overfitting for the CNN, part of the training data is used for validation during the training process and an early stopping criterion is applied on the validation loss. The evolution of loss and accuracy during the training epochs are shown in Supplementary Figures S2 and S3. All classifiers are mostly well-calibrated with the CNN showing a tendency for overconfidence, in particular, for the multi-class cases (Supplementary Figures S4–S6).

The performance of the classifiers is, on the one hand, assessed using the overall accuracy which is defined as the number of correct predictions divided by the number of total predictions. On the other hand, also the confidence in the predictions is assessed based on the probabilities assigned to each test sample from each dataset individually.

2.3. Out-of-sample frameworks and preprocessing

We use daily, land-masked temperature fields as samples, resulting in 6,888 grid cells being used as features for the classification. The 2-dimensional daily fields are either used directly in the case of the CNN classifier or flattened to a 1-dimensional feature array in the case of logistic regression. Note that the changing grid cell area with latitude due to longitude convergence is not explicitly accounted for in this study but may be implicitly learned by the classifiers.

Training and validation days are drawn from the 20-year period 1982–2001, as test data all days from the two temporally out-of-sample 10-year periods 2005–2014 and 2091–2100 are used. Therefore, all results are based on test samples that were never seen in training.

For the separation of “model” versus “observation” in Section 3, we use an additional dataset out-of-sample framework to explore whether the learned parameters can be generalized to unseen datasets: a separate classifier is trained for each dataset, where all samples from the dataset in question are withheld from training. To be even stricter, for climate models not only the model in question is withheld in training, but also all closely related models (see Supplementary Table S3 for a list of the model groups). We chose this approach because it has been shown that closely related models (i.e., models developed at the same institution) can have very similar output (Brunner et al., Reference Brunner, Pendergrass, Lehner, Merrifield, Lorenz and Knutti2020). The only exceptions to this are the regression coefficients shown in Figure 1 and the bootstrap tests for which all datasets were used in training.

Figure 1. (a) Logistic regression coefficients learned from 17,200 randomly drawn daily samples in the period 1982–2001 to separate models and observations. (b) Climatological mean, multi-model mean temperature difference to the mean over the four observational datasets in the period 2005–2014. See Supplementary Figure S7 for corresponding maps of the individual models. Coefficients and climatologies are calculated from daily data with the global mean removed.

In both machine learning and climate modeling, there are different preprocessing and bias-correction options and corresponding terminologies. These include more general approaches such as feature normalization and more domain-specific approaches such as bias correction. The latter are based on physical system understanding, such as the insight that the mean temperature bias is not a relevant predictor of a model’s ability to simulate changes in the climate system (Giorgi and Coppola, Reference Giorgi and Coppola2010). Here, we opt for two domain-specific bias-correction approaches from climate sciences: (a) from each daily temperature field, the global mean over that field is removed and (b) from each daily temperature field, the seasonal mean field is removed in addition to (a).

For (b), the seasonal mean temperature field is calculated individually for each dataset and only based on the training data. For each day-of-the-year and grid cell, it is calculated as the average over $ \pm 15 $ days around that day-of-the-year and over all training years (i.e., a mean over 31 days × 20 years = 620 values). For the ICON-Sapphire model, the seasonal cycle is estimated based on only 1 year (i.e., 31 values centered around each day of the year).

Finally, each sample is associated with one of two label types to be predicted by the machine learning classifiers: either “model”/“observation” (binary case) or the name of the dataset (multi-class case; see the dataset ID column in Supplementary Table S3). A summary of all classifiers, bias corrections, and out-of-sample strategies employed in this work can be found in Supplementary Table S1. The living code that implements the general method is available on GitHub (https://github.com/lukasbrunner/model_learning) and from Brunner and Sippel (Reference Brunner and Sippel2023). Supplementary Figures S7 and S8 show an example test day from each of the datasets and for the two preprocessing cases, respectively. While, in the first case, latitudinal temperature gradients due to differing solar insolation clearly dominate the temperature pattern, in the second case patches of regionally cooler or warmer temperatures emerge, related to the synoptic atmospheric conditions on that day. Note that on a given day (March 21st 2010 in Supplementary Figures S7 and S8) models are not expected to simulate identical weather patterns compared to the observations. The models are free-running and, therefore, simulate different synoptic situations and related temperature patterns while the four observation-based datasets assimilate measurements from the same day leading to very similar surface temperature patterns.

3.1. Regression coefficient maps to separate models and observations

First, we train a single logistic regression classifier on all available datasets to establish the regularization parameter and examine the learned coefficients. For this case, we use data with only the global mean removed. We use 200 different, randomly drawn training days from each of the 43 models resulting in 8,600 training samples labeled “model” which are matched by 8,600 random days labeled “observation” (2,150 from each of the four observational datasets). The classifier manages to correctly identify the vast majority of test samples for this setup so that the number of training samples was limited to this amount to save computational resources. It can be assumed that increasing the number of training samples would slowly improve the classification skill even further. The results are also robust to using different random training samples (see Section S3 in the Supplementary material for results from a 100-member bootstrap test).

Since temperatures between neighboring grid cells, used as features, are not independent, the 5-fold cross-validation yields a strong L₂ regularization parameter of about $ C=0.002 $ , which is used for all logistic regression classifiers in this section. The spatial patterns in the features that are important for separating models and observations are reflected in the regression coefficients learned by the classifier and are shown in Figure 1a. The distribution of the coefficients identifies areas important for the separation of climate models and observations on a daily scale. In addition, the sign of the coefficients can be interpreted physically, with positive values indicating regions where models tend to be warmer than the observations and vice-versa.

The most prominent region of negative coefficients is found in the North Atlantic, near the so-called North Atlantic warming hole (Chemke et al., Reference Chemke, Zanna and Polvani2020; Keil et al., Reference Keil, Mauritsen, Jungclaus, Hedemann, Olonscheck and Ghosh2020). Here models appear to systematically underestimate temperatures on a daily basis compared with observations and relative to the global mean. In contrast, there are regions of high coefficients at the eastern edges of the Pacific and Atlantic ocean basins. These regions correspond to persistent model biases in the representation of clouds and their radiative effects which are known to occur on all timescales from daily to decadal (Williams et al., Reference Williams, Bodas-Salcedo, Déqué, Fermepin, Medeiros, Watanabe, Jakob, Klein, Senior and Williamson2013; Hsi Yen Ma et al., Reference Hsi Yen Ma, Klein, Williams, Boyle, Bony, Douville, Fermepin, Medeiros, Tyteca, Watanabe and Williamson2014; Brient et al., Reference Brient, Roehrig and Voldoire2019; Bock et al., Reference Bock, Lauer, Schlund, Barreiro, Bellouin, Jones, Meehl, Predoi, Roberts and Eyring2020; Chen et al., Reference Chen, Wang, Bao and Li2022). In the equatorial Pacific, known as a region with notorious climate model biases that typically show too cold and too narrow equatorial cold tongues, negative coefficients are also found by the logistic regression classifier. This is accompanied by warm biases to the north and south (shown as positive regression coefficients) associated with the models’ representation of the intertropical convergence zone (Hirota et al., Reference Hirota, Takayabu, Watanabe and Kimoto2011; Li and Xie, Reference Li and Xie2014; Tian and Dong, Reference Tian and Dong2020). Overall, there are notable consistencies in several of the large-scale patterns between the logistic regression coefficients and the climatological mean, multi-model mean biases (Figure 1a,b).

However, there are also several regions where the two do not match, such as high northern latitudes, parts of the southern ocean, and the Antarctic coast. This can be an indication that the corresponding features are not shared across all (or at least most) potential days drawn from different models and from across the seasonal cycle. The high northern latitudes are briefly discussed here as one example for such a case with the clear climatological cold bias in the multi-model mean (Figure 1b) not being reflected in corresponding patterns in the regression coefficients (Figure 1a). This is probably caused by a combination of reasons as this region is known for its large (climatological) model spread (see, e.g., Notz and Community, Reference Notz and Community2020 and Supplementary Figure S9), its seasonally varying biases with cold biases mainly coming from the winter season (Davy and Outten, Reference Davy and Outten2020), and large internal variability which reduces the size of the regression coefficients there (Barnes et al., Reference Barnes, Hurrell, Ebert-Uphoff, Anderson and Anderson2019; Sippel et al., Reference Sippel, Meinshausen, Fischer, Székely and Knutti2020). In general, it is, therefore, not a priori clear if patterns important for the classification of individual days and long-term climatological biases (averaged over multiple models) would match at all, but our findings show that for several regions this is the case.

3.2. Logistic regression classification of out-of-sample datasets

In the rest of this section, we use the dataset out-of-sample approach described in Section 2.3 to show whether the classifiers can be generalized to datasets unseen in training. The probabilities assigned to the test samples by each corresponding classifier are aggregated in Figure 2a. For the vast majority of cases, the logistic regression classifiers assign the correct category with close to 100% probability leading to an overall accuracy of 99.4% (excluding ICON-Sapphire discussed below). Some of this skill may be due to remaining dependencies across families, which were quite loosely defined based on institutions in this study. However, it could also be an indication for the existence of persistent distinguishing features that can be transferred between models even in presence of the large internal variability on daily timescales.

Figure 2. Distribution of predicted probabilities for the dataset out-of-sample test days: for each dataset, the probabilities are estimated by a classifier which has not been trained on this dataset. The vertical dotted line at 0.5 marks the decision threshold between the two categories. ICON-Sapphire is never used in training and has only 1 year of data available. (a) Results for logistic regression classifiers using data with the daily global mean removed. (b) Same as (a) but for the convolutional neural network. (c) Same as (b) but using data with the seasonal cycle removed in addition.

For several models, a fraction of test samples is predicted with less certainty (boxes and whiskers emerging from the zero-line in Figure 2a and reliability diagram in Supplementary Figure S4). The model families most prone to get confused with observations are CMCC, CNRM, EC-Earth3, GFDL, and INM. An intuitive interpretation of this behavior might be that these models provide the “best” representation of “true” temperatures on a daily basis as they are most similar to the observations (at least viewed through the lens of logistic regression). However, this is not conclusive and would require additional evidence. First, there are only very few samples with a nonzero probability of belonging into the observation category for any given model, so the confusion could be due to chance (i.e., a model sample could be classified as observations, but for the wrong reasons so that no conclusion should be drawn about the overall performance of the model). A second important consideration is possible codebase overlap also between models and the observational datasets, which could be picked up by the classifier as discussed in more detail in Section 4. Nevertheless, the approach of intentionally exploiting misclassifications of models as observations is still a promising avenue for future research as recently also highlighted by work from Labe and Barnes (Reference Labe and Barnes2022), who find that for the Arctic test samples observations are most prone to get confused for certain models.

Focusing on model families that provide high and low resolution variants of the same model (AWI, CMCC, CNRM, HadGEM, MPI, and NorESM2), one can speculate about some resolution dependence in the probability to be misclassified. Several finer-resolution model variants appear to have a higher chance to be misclassified compared to their coarser-resolution siblings (model names including HR, MR/MM, or LR/LL indicating relatively higher to lower resolution in Figure 2a). Such a behavior would be consistent with studies based on climatological timescales such as the findings of Bock et al. (Reference Bock, Lauer, Schlund, Barreiro, Bellouin, Jones, Meehl, Predoi, Roberts and Eyring2020), who note that long-standing regional model biases are smaller in higher resolution versions of the same model.

3.3. Classifying samples from the kilometer-scale ICON-Sapphire model

To investigate this further and to highlight a potential application of our approach, we include preliminary results from the NextGEMS project and predict samples from 1 year of data from a global, storm-resolving (atmospheric resolution 5 km) simulation using ICON-Sapphire (Hohenegger et al., Reference Hohenegger, Korn, Linardakis, Redler, Schnur, Adamidis, Bao, Bastin, Behravesh, Bergemann, Biercamp, Bockelmann, Brokopf, Brüggemann, Casaroli, Chegini, Datseris, Esch, George, Giorgetta, Gutjahr, Haak, Hanke, Ilyina, Jahns, Jungclaus, Kern, Klocke, Kluft, Kölling, Kornblueh, Kosukhin, Kroll, Lee, Mauritsen, Mehlmann, Mieslinger, Naumann, Paccini, Peinado, Praturi, Putrasahan, Rast, Riddick, Roeber, Schmidt, Schulzweida, Schütte, Segura, Shevchenko, Singh, Specht, Stephan, Vogel, Wengel, Winkler, Ziemen, Marotzke and Stevens2022). Due to the high resolution, processes that need to be parameterized in CMIP6-type models can be explicitly resolved in ICON-Sapphire, which can be expected to have considerable impacts also on the daily temperature fields. However, they do not obviously show in the climatological mean difference to the observations which is comparable to coarser resolved models and shows similar patterns, although we note that this could be a coincidence given that only a single year is used (see Supplementary Figure S10). This means that based on the comparison of climatologies alone, one might assume that the logistic regression classifier should be able to unambiguously identify ICON-Sapphire as a model.

To investigate this, we classify test samples from ICON-Sapphire using the logistic regression classifier trained on all other datasets. However, despite the similar climatologies, the classifier is unable to clearly identify ICON-Sapphire as either model or observation with about half of the samples predicted to be in either category. This indicates that the classification is (at least partly) based on more complex relationships in the high-dimensional feature space than can be easily assessed by the comparison shown in Figure 1. The explicit resolution of processes which are parameterized in coarser models seems to lead to differences in the daily temperature fields in ICON-Sapphire that prevent a correct classification, even though they do not clearly emerge in the climatological mean. This points to a potentially highly encouraging emergent behavior of the ICON-Sapphire simulations, but in principle some effect of compensating errors cannot be excluded, and hence we warrant a careful interpretation of this result, which will be verified as soon as a longer simulation with this model is available. Should such a relationship hold in future research, it will enable innovative, new ways of model evaluation based on much shorter timescales than the 20+ years typically used.

To test whether the results for the ICON-Sapphire model are merely an artifact of the logistic regression, we also use a more complex, but less interpretable CNN. The CNN classifiers achieve a similar overall accuracy of 99.7% for the same dataset out-of-sample framework and samples with the global mean removed. In general, the pattern of models with samples that get misclassified is quite consistent between the logistic regression and the CNN (Figure 2a,b). While ICON-Sapphire still stands out as the model most likely to be confused for an observation, confirming the results from the logistic regression classifier, the CNN is able to correctly identify the majority (75.1%) of test samples, suggesting that it learns some more fundamental model properties that persist also at high resolution. Compared to the rather simple “model fingerprint” shown for the logistic regression in Figure 1a, the classifications from the CNN are, thus, likely to be based on more complex relationships. We plan to investigate these relationships and the importance of different grid cells for the skill of the CNN in future work, drawing on techniques that, for example, aim to reveal regions of higher and lower importance in the temperature maps used as input (Bach et al., Reference Bach, Binder, Montavon, Klauschen, Müller and Samek2015).

3.4. CNN classification of out-of-sample datasets without climatological bias

Based on these results, we, next, test if models and observations can still be separated in the absence of any climatological biases. To do this, we now also remove the mean seasonal cycle from each sample (see Section 2.3 for methodological details, Supplementary Figure S11 for the resulting multi-model mean bias equivalent to Figure 1b, and Supplementary Figure S12 for a breakdown by individual models). This means that any dataset-specific persistent regional biases that might have served as a basis for separation so far are now removed, along with any biases in the equator-pole gradient or between the hemispheres. Therefore, the classifiers can now only train on the spatial relationships of the remaining internal variability (which could be interpreted as daily weather), making the classification task considerably harder.

For this case, logistic regression no longer has any skill as the only remaining sources of information are nonlinear relations between the spatial structures of daily global weather and the test samples are all centered around the decision threshold (not shown). In contrast, the CNN achieves an overall accuracy of about 94.2% demonstrating the power of this nonlinear method. Figure 2c shows the corresponding breakdown of predicted probabilities revealing that now almost all models get confused for observations a number of times but still all are classified correctly for the vast majority of test samples, with the sole exception being the ICON-Sapphire model.

For ICON-Sapphire, the seasonal cycle has been estimated using only the 31-day running window as only 1 year is available. This differs from the other datasets where the seasonal cycle has been calculated over the full 20 years of the training period (Section 2.3). Therefore, the results for ICON-Sapphire need to be interpreted with care and should be revisited once more data are available. Nevertheless, we show these preliminary findings here, to highlight that, based on these results the structure of the remaining internal variability in ICON-Sapphire is recognized as more closely resembling observations than CMIP6 generation models. In fact, ICON-Sapphire is more frequently classified as observation than three of the four observational datasets (ERA5, MERRA2, and 20CR) by the respective classifiers trained with the tested dataset withheld. While this result is preliminary and not conclusive due to the limited amount of data available, it, again, points toward very encouraging properties of kilometer-scale models that warrant closer investigation as more data become available.

In turn, the DOISST dataset is identified correctly with perfect accuracy in all three cases shown in Figure 2 with only individual misclassifications appearing even when bootstrapping the training data (Supplementary Figure S1). Optimistically interpreted, this could mean that the classifiers are picking up on the fact that DOISST is the dataset with the least amount of model included (see Section 2.1) and that they have indeed learned some fundamental distinguishing features. This interpretation is supported by the fact that 20CR is most prone to get confused for a model, while also being the observational dataset that relies most heavily on a model for its output.

The overall skill of this binary classification of out-of-sample datasets is notable, in particular, when considering that bias correcting each model by subtracting the mean seasonal cycle and the daily global mean effectively removes the entire time-persistent regional bias as well as any global mean offset between the datasets. This means that the only remaining sources of information to learn from are amplitude and spatial dependencies of the remaining daily temperature variability. One remaining source of model-observation differences for this case could arise from the coupling of atmosphere and ocean and the resulting surface energy balance in the models.

In a planned follow-up study, we will build on the results presented here and analyze the origin of this skill in more detail, using explainable machine learning techniques (e.g., Bach et al., Reference Bach, Binder, Montavon, Klauschen, Müller and Samek2015) as well as more specific approaches drawing on domain knowledge from climate sciences. The latter will include, for example, targeted masking of certain regions where climate models are known to perform particularly well/poorly to systematically investigate which areas of the globe are essential for skill and how this might depend on the models used. Such a combination of general and domain-specific approaches to classification skill is important to account for special properties of the temperature maps used compared to more general image classification. These include fundamental properties of the climate system such as the imprint of topography, temperature gradients, and circulation patterns, which are to some extent common to all datasets.