Transient global hindcast simulations were performed with the three-dimensional intermediate complexity Earth system model LOVECLIM (9, 50) to investigate the millennial-scale climate variations of the last glacial [50 to 30 ka ago (ka B.P.)]. Initial conditions for the transient run were based on an equilibrium spin-up simulation with an atmospheric CO2 content of 207.5 ppmv (parts per million by volume) and an orbital forcing and estimated ice sheet orography for 50 ka B.P. (9). In the transient run, these glacial boundary conditions were continuously updated and a time-varying anomalous freshwater flux to the North Atlantic (55°-10°W, 50°-65°N) was applied as forcing (9). The time scale of the model data was originally placed on the GICC05 (B.P.) chronology (9). For the present study, we used simulated stadial-interstadial anomalies of annual mean sea surface temperature, 10% spring sea ice extent, and sea ice thickness for the North Atlantic. To calculate stadial-interstadial anomalies displayed in Fig. 1, we averaged stadial conditions over 39.9 to 38.46, 36.34 to 35.5, 34.5 to 33.7, and 33.1 to 32.45 ka B.P., and interstadial conditions over 40.2 to 40.0, 38.34 to 37.9, 35.35 to 34.9, 33.6 to 33.4, and 32.38 to 32.16 ka B.P. For the model-proxy data comparison, we used simulated time series of spring (April-May-June) sea ice cover averaged over 60°N-62.5°N, 1.25°W-1.25°E and 60°N-80°N, 15°W-20°E. The time scale of simulated time series of sea ice cover shown in Fig. 2F was transferred to the GICC05 (b2k) chronology (51, 52).