As of 2023, global mean temperature has risen by about 1.45 +/- 0.12 degrees C with respect to the 1850-1900 pre-industrial (PI) baseline according to the World Meteorological Organization. This rise constitutes the first period of substantial global warming since the Last Deglaciation, when global temperatures rose over several millennia by about 4.0-7.0 degrees C according to proxy reconstructions. Similar levels of warming could be reached in the coming centuries considering current and possible future emissions. Such warming causes widespread changes in the climate system, of which the mean state provides only an incomplete picture. Instead, fluctuations around the mean and in higher-order statistics need to be considered. Indeed, climate's variability and the distributions of climate variables change with warming, impacting, for example, ecosystems and the frequency and intensity of extremes. However, previous investigations of climate variability focus mostly on measures such as variance, or standard deviation, and on quasi-equilibrium states such as the Holocene or Last Glacial Maximum (LGM). Changes in the tails of distributions of climate variables and transition periods such as the Last Deglaciation remain largely unexplored.Therefore, we investigate changes of climate variability on annual to millennial timescales in 15 transient climate model simulations of the Last Deglaciation. This ensemble consists of models of varying complexity, from an energy balance model to Earth system models (ESMs), and includes sensitivity experiments, which differ only in terms of their underlying ice sheet reconstruction, meltwater protocol, or consideration of volcanic forcing. The ensemble simulates an increase in global mean temperature of 3.0-6.6 degrees C between the LGM and Holocene. Against this backdrop, we examine whether common patterns of variability emerge in the ensemble. To this end, we compare the variability in surface climate during the LGM, Deglaciation, and Holocene by estimating and analyzing the distributions and power spectra of surface temperature and precipitation. For analyzing the distribution shapes, we turn to the higher-order moments of variance, skewness, and kurtosis. These show that the distributions cannot be assumed to be normal, a precondition for commonly used statistical methods. During the LGM and Holocene, they further reveal significant differences, as most simulations feature larger temperature variance during the LGM than the Holocene, in line with results from reconstructions.As a transition period, the Deglaciation stands out as a time of high variance in surface temperature and precipitation, especially on decadal and longer timescales. In general, this dependency on the mean state increases with model complexity, although there is a large spread between models of similar complexity. Some of that spread can be explained by differences in ice sheet, meltwater, and volcanic forcings, revealing the impact of simulation protocols on simulated variability. The forcings affect variability not only on their characteristic timescales. Rather, we find that they impact variability on all timescales from annual to millennial. The different forcing protocols further have a stronger imprint on the distributions of temperature than precipitation. A reanalysis of the LGM exhibits similar global mean variability to most of the ensemble, but spatial patterns vary. However, paleoclimate data assimilation combines model and proxy data information using a Kalman-filter-based algorithm. More research is needed to disentangle their relative impact on reconstructed levels of variability. As such, uncertainty around the models' abilities to capture climate variability likewise remains, affecting simulations of all time periods: past, present, and future. Decreasing this uncertainty warrants a systematic model-data comparison of simulated variability during periods of warming.
