The transient simulation was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR) coupled to the ice-sheet model mPISM and the solid-earth model VILMA. MPI-ESM includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. Embedded into MPIOM1.6 is an Eulerian iceberg model (Erokhina and Mikolajewicz 2024). mPISM is based on PISM0.7.3. and was employed using a polar stereographic grid with a resolution of 10 km in the northern hemisphere and 15 km for Antarctica. VILMA was used in its 1D configuration and discretized on a Gaussian F128 grid with a nominal resolution of 0.7°. The applied setup was introduced in detail in Mikolajewicz et al. (2024). Following an asynchronous spin-up from 46 ka to 26ka during which MPI-ESM was run with an acceleration factor of 10, the model was integrated without acceleration from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). This simulation used climatological volcanic forcing representing PI conditions. Forcing fields between MPI-ESM, mPISM, and VILMA as well as changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask were updated every 10 years. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations, bug fixes and volcanic boundary forcing. These differences are described according to the CMIP6 nomenclature, where r denotes the realization, i the initialization method, p differences in the physics and f in the forcing. The ensemble is described in detail in Mikolajewicz et al. (2024). This particular simulation corresponds to the simulation D1.1 in Mikolajewicz et al. (2024). Note that the time runs with an offset of +25001 years, meaning that the year range 1-25000 of the data set equals the years 25000-1 BP (Before Present). All data are global datasets with the exception of output from the ice sheet model where regional northern hemispheric and southern hemispheric domains are used.

The transient simulation was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR) coupled to the ice-sheet model mPISM and the solid-earth model VILMA. MPI-ESM includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. Embedded into MPIOM1.6 is an Eulerian iceberg model (Erokhina and Mikolajewicz 2024). mPISM is based on PISM0.7.3. and was employed using a polar stereographic grid with a resolution of 10 km in the northern hemisphere and 15 km for Antarctica. VILMA was used in its 1D configuration and discretized on a Gaussian F128 grid with a nominal resolution of 0.7°. The applied setup was introduced in detail in Mikolajewicz et al. (2024). Following an asynchronous spin-up from 46 ka to 26ka during which MPI-ESM was run with an acceleration factor of 10, the model was integrated without acceleration from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). This simulation used annualy varying volcanic forcing (Schindlbeck-Belo et al., 2024). Forcing fields between MPI-ESM, mPISM, and VILMA as well as changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask were updated every 10 years. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations, bug fixes and volcanic boundary forcing. These differences are described according to the CMIP6 nomenclature, where r denotes the realization, i the initialization method, p differences in the physics and f in the forcing. The ensemble is described in detail in Mikolajewicz et al. (2024). This particular simulation corresponds to the simulation D2.1 in Mikolajewicz et al. (2024). Note that the time runs with an offset of +25001 years, meaning that the year range 1-25000 of the data set equals the years 25000-1 BP (Before Present). All data are global datasets with the exception of output from the ice sheet model where regional northern hemispheric and southern hemispheric domains are used.

The transient simulation was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR) coupled to the ice-sheet model mPISM and the solid-earth model VILMA. MPI-ESM includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. Embedded into MPIOM1.6 is an Eulerian iceberg model (Erokhina and Mikolajewicz 2024). mPISM is based on PISM0.7.3. and was employed using a polar stereographic grid with a resolution of 10 km in the northern hemisphere and 15 km for Antarctica. VILMA was used in its 1D configuration and discretized on a Gaussian F128 grid with a nominal resolution of 0.7°. The applied setup was introduced in detail in Mikolajewicz et al. (2024). Following an asynchronous spin-up from 46 ka to 26ka during which MPI-ESM was run with an acceleration factor of 10, the model was integrated without acceleration from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). This simulation used climatological volcanic forcing representing PI conditions. Forcing fields between MPI-ESM, mPISM, and VILMA as well as changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask were updated every 10 years. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations, bug fixes and volcanic boundary forcing. These differences are described according to the CMIP6 nomenclature, where r denotes the realization, i the initialization method, p differences in the physics and f in the forcing. The ensemble is described in detail in Mikolajewicz et al. (2024). This particular simulation corresponds to the simulation D1.2 in Mikolajewicz et al. (2024). Note that the time runs with an offset of +25001 years, meaning that the year range 1-25000 of the data set equals the years 25000-1 BP (Before Present). All data are global datasets with the exception of output from the ice sheet model where regional northern hemispheric and southern hemispheric domains are used.

The transient simulation was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR) coupled to the ice-sheet model mPISM and the solid-earth model VILMA. MPI-ESM includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. Embedded into MPIOM1.6 is an Eulerian iceberg model (Erokhina and Mikolajewicz 2024). mPISM is based on PISM0.7.3. and was employed using a polar stereographic grid with a resolution of 10 km in the northern hemisphere and 15 km for Antarctica. VILMA was used in its 1D configuration and discretized on a Gaussian F128 grid with a nominal resolution of 0.7°. The applied setup was introduced in detail in Mikolajewicz et al. (2024). Following an asynchronous spin-up from 46 ka to 26ka during which MPI-ESM was run with an acceleration factor of 10, the model was integrated without acceleration from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). This simulation used climatological volcanic forcing representing PI conditions. Forcing fields between MPI-ESM, mPISM, and VILMA as well as changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask were updated every 10 years. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations, bug fixes and volcanic boundary forcing. These differences are described according to the CMIP6 nomenclature, where r denotes the realization, i the initialization method, p differences in the physics and f in the forcing. The ensemble is described in detail in Mikolajewicz et al. (2024). This particular simulation corresponds to the simulation D1.3 in Mikolajewicz et al. (2024). Note that the time runs with an offset of +25001 years, meaning that the year range 1-25000 of the data set equals the years 25000-1 BP (Before Present). All data are global datasets with the exception of output from the ice sheet model where regional northern hemispheric and southern hemispheric domains are used. Please note that 10 years of model output from MPI-ESM, spanning the time from 381-370 BP are missing. The corresponding files only contain missing values.

The transient simulation was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR) coupled to the ice-sheet model mPISM and the solid-earth model VILMA. MPI-ESM includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. Embedded into MPIOM1.6 is an Eulerian iceberg model (Erokhina and Mikolajewicz 2024). mPISM is based on PISM0.7.3. and was employed using a polar stereographic grid with a resolution of 10 km in the northern hemisphere and 15 km for Antarctica. VILMA was used in its 1D configuration and discretized on a Gaussian F128 grid with a nominal resolution of 0.7°. The applied setup was introduced in detail in Mikolajewicz et al. (2024). Following an asynchronous spin-up from 46 ka to 26ka during which MPI-ESM was run with an acceleration factor of 10, the model was integrated without acceleration from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). This simulation used annualy varying volcanic forcing (Schindlbeck-Belo et al., 2024). Forcing fields between MPI-ESM, mPISM, and VILMA as well as changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask were updated every 10 years. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations, bug fixes and volcanic boundary forcing. These differences are described according to the CMIP6 nomenclature, where r denotes the realization, i the initialization method, p differences in the physics and f in the forcing. The ensemble is described in detail in Mikolajewicz et al. (2024). This particular simulation corresponds to the simulation D2.2 in Mikolajewicz et al. (2024). Note that the time runs with an offset of +25001 years, meaning that the year range 1-25000 of the data set equals the years 25000-1 BP (Before Present). All data are global datasets with the exception of output from the ice sheet model where regional northern hemispheric and southern hemispheric domains are used.

The transient simulation was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR) coupled to the ice-sheet model mPISM and the solid-earth model VILMA. MPI-ESM includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. Embedded into MPIOM1.6 is an Eulerian iceberg model (Erokhina and Mikolajewicz 2024). mPISM is based on PISM0.7.3. and was employed using a polar stereographic grid with a resolution of 10 km in the northern hemisphere and 15 km for Antarctica. VILMA was used in its 1D configuration and discretized on a Gaussian F128 grid with a nominal resolution of 0.7°. The applied setup was introduced in detail in Mikolajewicz et al. (2024). Following an asynchronous spin-up from 46 ka to 26ka during which MPI-ESM was run with an acceleration factor of 10, the model was integrated without acceleration from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). This simulation used annualy varying volcanic forcing (Schindlbeck-Belo et al., 2024). Forcing fields between MPI-ESM, mPISM, and VILMA as well as changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask were updated every 10 years. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations, bug fixes and volcanic boundary forcing. These differences are described according to the CMIP6 nomenclature, where r denotes the realization, i the initialization method, p differences in the physics and f in the forcing. The ensemble is described in detail in Mikolajewicz et al. (2024). This particular simulation corresponds to the simulation D2.3 in Mikolajewicz et al. (2024). Note that the time runs with an offset of +25001 years, meaning that the year range 1-25000 of the data set equals the years 25000-1 BP (Before Present). All data are global datasets with the exception of output from the ice sheet model where regional northern hemispheric and southern hemispheric domains are used.

The transient simulation was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR) coupled to the ice-sheet model mPISM and the solid-earth model VILMA. MPI-ESM includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. Embedded into MPIOM1.6 is an Eulerian iceberg model (Erokhina and Mikolajewicz 2024). mPISM is based on PISM0.7.3. and was employed using a polar stereographic grid with a resolution of 10 km in the northern hemisphere and 15 km for Antarctica. VILMA was used in its 1D configuration and discretized on a Gaussian F128 grid with a nominal resolution of 0.7°. The applied setup was introduced in detail in Mikolajewicz et al. (2024). Following an asynchronous spin-up from 46 ka to 26ka during which MPI-ESM was run with an acceleration factor of 10, the model was integrated without acceleration from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). This simulation used annualy varying volcanic forcing (Schindlbeck-Belo et al., 2024). Forcing fields between MPI-ESM, mPISM, and VILMA as well as changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask were updated every 10 years. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations, bug fixes and volcanic boundary forcing. These differences are described according to the CMIP6 nomenclature, where r denotes the realization, i the initialization method, p differences in the physics and f in the forcing. The ensemble is described in detail in Mikolajewicz et al. (2024). This particular simulation corresponds to the simulation D2.4 in Mikolajewicz et al. (2024). Note that the time runs with an offset of +25001 years, meaning that the year range 1-25000 of the data set equals the years 25000-1 BP (Before Present). All data are global datasets with the exception of output from the ice sheet model where regional northern hemispheric and southern hemispheric domains are used.

The transient simulation was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR) coupled to the ice-sheet model mPISM and the solid-earth model VILMA. MPI-ESM includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. Embedded into MPIOM1.6 is an Eulerian iceberg model (Erokhina and Mikolajewicz 2024). mPISM is based on PISM0.7.3. and was employed using a polar stereographic grid with a resolution of 10 km in the northern hemisphere and 15 km for Antarctica. VILMA was used in its 1D configuration and discretized on a Gaussian F128 grid with a nominal resolution of 0.7°. The applied setup was introduced in detail in Mikolajewicz et al. (2024). Following an asynchronous spin-up from 46 ka to 26ka during which MPI-ESM was run with an acceleration factor of 10, the model was integrated without acceleration from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). This simulation used annualy varying volcanic forcing (Schindlbeck-Belo et al., 2024). Forcing fields between MPI-ESM, mPISM, and VILMA as well as changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask were updated every 10 years. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations, bug fixes and volcanic boundary forcing. These differences are described according to the CMIP6 nomenclature, where r denotes the realization, i the initialization method, p differences in the physics and f in the forcing. The ensemble is described in detail in Mikolajewicz et al. (2024). This particular simulation corresponds to the simulation D1.4 in Mikolajewicz et al. (2024). Note that the time runs with an offset of +24051 years, meaning that the year range 50-25000 of the data set equals the years 25000-50 BP (Before Present). All data are global datasets with the exception of output from the ice sheet model where regional northern hemispheric and southern hemispheric domains are used.

For transient-deglaciation-prescribed-glac1d-methane (r1i1p1f1): The transient experiment was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR). The model includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. The model has been extended with a full methane cycle, including a simplified atmospheric chemistry, as described in Kleinen et al. (2020), Kleinen et al. (2021) and Kleinen et al. (2022). For the experiment, the model was integrated from a glacial state at 23 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). Ice sheets and surface topographies were prescribed from GLAC-1D (Tarasov et al., 2012) reconstructions. All forcing fields are updated every 10 years of the simulations and initiate changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask. For transient-deglaciation-prescribed-glac1d-methane (r1i1p1f2): This experiment is derived from the r1i1p1f1 experiment. It is branched off at 18 ka BP, and starting at 15.2 ka BP, meltwater from the Laurentide ice sheet was removed from the system and stored. The accumulated meltwater was released over a period of 1200 years starting at 12.8 ka BP. This induced a collapse of the AMOC, leading to climatic conditions very similar to the Younger Dryas cold period seen in climate reconstructions. After the end of the meltwater release, the circulation recovered quickly and climatic conditions converged with the r1i1p1f1 experiment. Important: Please be aware that due to CMORization constraints data set time values concerning to "Before Present(BP)" time got a year offset (+25001). Means for example year range 2001 to 25000 model time is equal to 23000 to 1 BP.

The transient simulation was performed with the Max Planck Institute for Meteorology Earth System Model version 1.2 in coarse resolution (MPI-ESM-CR). The model includes the spectral atmospheric model ECHAM6.3 at T31 horizontal resolution (approx. 3.75°) and 31 vertical levels, the land surface vegetation model JSBACH3.2, and the primitive equation ocean model MPIOM1.6 with a nominal resolution of 3°. The applied setup was introduced in detail in Kapsch et al. (2021 and 2022). For the experiment, the model was integrated from a glacial state at 26 ka to the year 1950 with prescribed atmospheric greenhouse gas concentrations (Köhler et al., 2017) and insolation (Berger & Loutre, 1991). Ice sheets and surface topographies were prescribed from ICE-6G (Peltier et al., 2015) reconstructions. All forcing fields are updated every 10 years of the simulations and initiate changes in the topography, glacier mask, river pathways, ocean bathymetry, and land-sea mask. Meltwater from ice sheets is calculated as the temporal derivative of ice thickness at grid points covered by grounded ice sheets. The derived meltwater is then distributed by the hydrological discharge model and finally released into the ocean as freshwater. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations and bug fixes. The ensemble is described in detail in the Supplementary Material of Kapsch et al. (2022). The current simulation refers to ICE6G_P2 in Kapsch et al. (2022) and is significantly colder during the glacial than model version P1, indicating a stronger sensitivity to greenhouse gas changes. The glacial cooling over the North Atlantic is more realistic in P2 and P3 than in P1. P2 and P3 are also more sensitive to small changes in the meltwater forcing than P1. The experiment was performed as part of a model ensemble that contains simulations differing in terms of their tuning parameters, parameterizations, bug fixes and forcing. These differences are described according to the CMIP6 nomenclature, where r denotes the realization, i the initialization method, p differences in the physics and f in the forcing. The ensemble is described in detail in the Supplementary Material of Kapsch et al. (2022).