Summer and winter maxima of daily maximum temperature are obtained from a 14-member CMIP5 ensemble of GCMs as well as three reanalysis surrogate observation datasets. Likewise, summer and winter minima of daily minimum temperature are obtained from the same datasets. Summer is defined as June-July-August (JJA) in the Northern Hemisphere and December-January-February (DJF) in the Southern Hemisphere and winter is defined as DJF and JJA in the Northern and Southern Hemispheres, respectively. All these data are obtained for historical (1970–1999) and, for GCMs only, future (2070–2099) periods. Future GCM data is obtained for the Relative Concentration Pathway 4.5 (RCP4.5) moderate greenhouse gas emissions trajectory scenario14. The climate models and reanalysis datasets are summarized in SI: Table S7.

GEV and GPD models45 are fit separately to each of the four classes (summer maxima, summer minima, winter maxima and winter minima) of extrema at each grid point for all GCMs and reanalyses. All minima are negated to enable the model-fitting process and transformed back afterward45. To assess ensemble skill in emulating historical statistical attributes of extremes, historical parameters from GEV and GPD fits are averaged over latitudinal bands, for land (Figure 5 and Figure S2) and ocean (Figure 6 and Figure S3). Shape parameters are almost always significantly less than 0 or not statistically differentiable from 0 and almost without exception never greater than 0, for both models (Figures S1 and S4).

We perform a bootstrap resampling analysis similar to recent work17,46 with seasonal extrema data to obtain robust estimates of changes in the distributions of extrema. Specifically, at each grid cell, for each GCM and separately for each of summer and winter maxima and minima (generally stated, extrema), we treat future extrema data (29 data points) as a pool of approximately independent data to iteratively resample from with replacement. We draw 29 samples with replacement from this pool 100 times. In each of those 100 random samples of size 29, we fit GEV parameters via maximum likelihood45 and subsequently simulate a random realization from the fitted parameters. As such, we expect to obtain robust estimates of extrema percentiles where sufficient realizations exist. Rarely, unstable (i.e., large) estimates of large shape parameters led to unrealistic bootstrap realizations. In these cases, where resultant temperature was simulated to be above 375 Kelvins for maxima or below 180 Kelvins for minima, these realizations were discarded and replacement realizations were simulated. We note that the asymmetry results described next do not change materially if the bootstrap analysis is not employed; the bootstrap enhances confidence in the robustness of the particular extrema percentiles used subsequently.

Separately for each of summer and winter maxima and minima, we estimate spatial fields of historical climatology Z by averaging historical seasonal extrema at each grid cell over time. Then, for each grid cell, we subtract Z from each T = 1,…t,…29 historical seasonal extrema spatial fields H t , which yields 29 spatial fields ΔH t = H t − Z. We perform the same operation for the S = 1,…s,…100 bootstrapped future spatial fields of extrema, F s , which yields 100 spatial spatial-temporal fields ΔF s = F s − Z. ΔF and ΔH represent the three-dimensional concatenation of 100 future and 29 historical spatial fields, respectively and are considered deviations from the historical climatology Z, in Kelvins.

We define four factors, with their levels in parentheses: X 1 - Tail type (minima or maxima), X 2 - Season (summer or winter), X 3 - Terrain type (land or ocean), X 4 - Region (SP, SHEX, TX, NHEX, or NP). Over each possible combination of the four and for each of the 14 GCMs, denoted X 5 , we obtain percentiles (0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95) from each of ΔF and ΔH, defined as P F and P H , respectively, where subscripts F and H denote “Future” and “Historical”, respectively. As percentiles increase, events always increase in temperature (i.e., 95th percentile of winter minima is hotter than the 5th). These percentiles are subtracted from each other to yield ΔP = P F − P H that shows projected changes in the distribution of extrema deviations (Figures 2–3). We note that percentiles for each combination of factors may be obtained from effective samples sizes that differ substantially and as such especially the more extreme percentiles (0.05 and 0.95) may be relatively more robust with larger sample sizes. For example, closer to the poles, size of grid cells becomes smaller and as a result the effective sample size of those regions shrinks. Because of this and since the bootstrap also implies that estimates of more extreme percentiles (e.g., 0.01 and 0.99 and beyond) are less robust, we do not venture that far into the tails in the following asymmetry analysis.

Next, we describe relationships between the X 1 through X 5 and the asymmetry inferred visually from Figures 2 and 3. “Asymmetry”, or for short, Y, in extrema projections is constructed by subtracting the change in the 5th percentile of extrema (ΔP05) from the change in the 95th percentile (ΔP95), i.e., Y = ΔP95 − ΔP05. We fit a linear mixed effect model47 to obtain maximum likelihood estimates relationships between four factors, X 1 through X 4 , as well as their two way interactions, with Y and to estimate how much variability in Y can be explained by systematic GCM differences, X 5 . Those estimates, as well as their variability in terms of inter-GCM differences, are portrayed for each possible combination of X 1 through X 4 in Figure 4. More details on the linear model implementation and results can be found in the SI.

To further test the robustness of the asymmetry insight, we select three GCMs and repeat the above entire GEV bootstrap-driven asymmetry analysis but first remove a historical daily seasonal cycle from each grid cell of each GCM's historical and future output. This analysis ensures that the findings in this study are not an artifact of seasonal timing (e.g., seasonal maxima or minima that will tend to be repeatedly extracted from very specific times of the year). Indeed we find this analysis strongly implies the robustness of the asymmetric projections to seasonal timing of extremes (SI: Figure S5 and Table S5).

Similar analyses are performed using the GPD. When using the GPD, we choose grid-wise 99th percentiles as the location parameters and fit scale and shape parameters to data exceeding those thresholds. With the GPD, the location parameters for both 20th century and projected data are estimated directly as percentiles and therefore they are not associated with statistical model-fitting uncertainty. Of primary interest here is whether the GPD analysis reveals similar patterns compared to the GEV model, especially in terms of evaluating historical GCM runs against reanalysis statistics. Outputs related to the GPD can be found in the SI: Figures S2–S4.