The PMIP4 contribution to CMIP6 Part 4: Scientific objectives and experimental design of the PMIP4-CMIP6 Last Glacial Maximum experiments and PMIP4 sensitivity experiments. Geoscientific

. The Last Glacial Maximum (LGM, 21 000 years ago) is one of the suite of paleoclimate simulations included in the current phase of the Coupled Model Intercomparison Project (CMIP6). It is an interval when insolation was similar to the present, but global ice volume was at a maximum, eustatic sea level was at or close to a minimum, greenhouse gas concentrations were lower, atmospheric aerosol loadings were higher than today, and vegetation and land-surface characteristics were different from today. The LGM has been a focus for the Paleoclimate Modelling Intercomparison Project (PMIP) since its inception, and thus many of the problems that might be associated with simulating such a radically different climate are well documented. The LGM state provides an ideal case study for evaluating climate model performance because the changes in forcing and temperature between the LGM and pre-industrial are of the same order of magnitude as those projected for the end of the 21st century. Thus, the CMIP6 LGM experiment could provide additional information that can be used to constrain es-timates of climate sensitivity. The design of the Tier 1 LGM experiment ( lgm ) includes an assessment of uncertainties in boundary conditions, in particular through the use of different reconstructions of the ice sheets and of the change in dust forcing. Additional (Tier 2) sensitivity experiments have been designed to quantify feedbacks associated with land-surface


Introduction
The Last Glacial Maximum (LGM), dated ∼ 21 000 years BP, is the last period during which the global ice volume was at its maximum, and eustatic sea level at or near to its minimum, ∼ 115 to 130 m below the present sea level (Lambeck et al., 2014;Peltier and Fairbanks, 2006). It has been defined as a relatively stable climatic period, in between two major intervals of iceberg discharge into the North Atlantic, Heinrich events 1 and 2 (Mix et al., 2001). In addition to ex-panded Greenland and Antarctic ice sheets, there were large ice sheets over northern North America and northern Europe. They caused large perturbations to the atmospheric radiative balance due to their albedo, and to atmospheric circulation because they were several kilometres high and therefore acted as large topographic barriers to the atmospheric flow. They also caused changes in coastlines and bathymetry due to the change in sea level and the mass load of the ice sheets. The atmospheric radiative budget was different at the LGM from the pre-industrial state due to much lower atmospheric greenhouse gas (GHG) concentrations (e.g. Bereiter et al., 2015, for CO 2 ;Loulergue et al., 2008, for CH 4 ). Both the change in ice sheets and in GHG concentrations are negative radiative forcings and contribute, with impacts of similar orders of magnitude, to a climate much colder than today (e.g. Yoshimori et al., 2009;Brady et al., 2013). They are the main drivers of differences in the LGM atmosphere compared to present or pre-industrial conditions. The ocean, continental surface, and carbon cycle respond and feed back to the atmosphere: the ocean circulation is affected by changes in the atmosphere as well as in coastlines and bathymetry; atmospheric and vegetation changes alter the atmospheric chemistry and aerosol loads; climate changes as well as CO 2 lowering modify the distribution and productivity of vegetation.
The LGM is extensively documented by continental, ice, and marine indicators. Sea surface temperature reconstructions from different indicators (MARGO Project Members, 2009) indicate a cooling from a few • C in the tropics to more than 10 • C at higher latitudes. Tracers of ocean circulation (e.g. δ 13 C, Lynch-Stieglitz et al., 2007;Pa/Th, ε Nd , Böhm et al., 2015) indicate convection in the North Atlantic, producing intermediate waters (the so-called Glacial North Atlantic Intermediate Waters, or GNAIW, Lynch-Stieglitz et al., 2007) rather than deep waters (North Atlantic Deep Water, NADW) characteristic of the modern ocean. Pollen and plant macrofossil records indicate that LGM vegetation patterns were very different from today, with expansion of steppe and tundra in Eurasia, and reduced cover of moist forests in the tropics (Prentice et al., 2000. Pollen-based climate reconstructions (e.g.  generally show a cooling compared to the present, which can reach more than 10 • C for mean annual temperature at some locations. Dry conditions and the reduction in vegetation cover led to major changes in dust emission, recorded in ice cores, marine sediments, and loess/paleosol deposits. Based on global compilations of these records, it has been estimated that the LGM was 2-4 times dustier than the Holocene on global average (Kohfeld and Harrison, 2001;Maher et al., 2010). However, the spatial variability of changes in dust deposition rates is very large, with a 20-fold increase shown in polar ice cores (Lambert et al., 2008;Steffensen et al., 2008). These changes in dust loading reflect changes in surface characteristics, winds, and precipitation. They also represent an important feedback from the climate system onto atmospheric radiative properties, which include direct and indirect effects on the atmospheric radiation budget through scattering and absorption of radiation and dust-cloud interactions (Boucher et al., 2013), which can alter regional climates (Claquin et al., 2003;Mahowald et al., 2006;Takemura et al., 2009;Hopcroft et al., 2015). Dust deposition changes can also impact the global carbon cycle, in particular because of the potential fertilization effect that dust-borne iron may exert on the Southern Ocean marine ecosystems and carbon sequestration in the deep ocean (Martin et al., 1990;Bopp et al., 2003;Kohfeld et al., 2005).
Modelling the LGM climate has been a focus for the Paleoclimate Modelling Intercomparison Project (PMIP) since its beginning (Joussaume and Taylor, 1995), progressing from simulations with atmospheric general circulation models (AGCMs), using prescribed ocean conditions or coupled to slab ocean models, to simulations using fully coupled atmosphere-ocean general circulation models (AOGCMs), some of which included vegetation dynamics, in the second phase of PMIP (PMIP2: Braconnot et al., 2007), and Earth system models (ESMs) with interactive carbon cycles in PMIP's third phase (PMIP3: Braconnot et al., 2012). The progression from AGCMs to AOGCMs has allowed oceanic reconstructions to be used for model evaluation and analysis of the physical consistency (as represented by models) of continental and oceanic reconstructions (e.g. . At each phase in this evolution, PMIP has taken into account new knowledge about boundary conditions, in particular for the form of the ice sheets, as well as the new capabilities of climate models. This paper describes the experimental set-up for the LGM experiments for PMIP4-CMIP6. Compared to the previous phases of PMIP, the new aspects of the PMIP4-CMIP6 simulations are the inclusion of dust forcing, either by using models in which the dust cycle is interactive or by prescribing atmospheric dust concentrations, so as to consider the interactions between dust and radiation. This is expected to cause significant differences in simulated regional climates and to have impacts on ocean biogeochemistry through a more realistic representation of dust input at the ocean surface; and explicit consideration of the uncertainties in ice-sheet reconstructions and the impact of different reconstructions of ice-sheet elevation on simulated climate. Consideration of uncertainties in boundary conditions is particularly important when comparing the model results to paleoclimatic reconstructions and drawing conclusions about the capabilities of the state-of-the-art models that are used for future climate projections.
This paper provides guidelines on the implementation of the PMIP4 LGM experiment in the CMIP6 climate models. This LGM experiment is a Tier 1 CMIP6 experiment (as defined in Eyring et al., 2016, i.e. Tier 1 defines experiments with the highest priority) and one of the two possible entry cards for PMIP4 . It is also a reference experiment for additional sensitivity experiments, which are considered Tier 2 and described here. These additional experiments are designed to improve our understanding of the simulated LGM climate. Section 2 presents how the LGM experiment will address CMIP6 questions. Section 3 describes the LGM PMIP4-CMIP6 experiments and the PMIP4 sensitivity experiments that were designed to address these questions. Section 4 details the implementation of LGM simulations. Section 5 finally outlines the analysis plan of the LGM experiments. There are two companion papers which document the other PMIP4-CMIP6 experiments: the last interglacial and mid-Holocene  and the last millennium (Jungclaus et al., 2017). In addition, Kageyama et al. (2016) provide an overview of the PMIP4-CMIP6 project.
2 The relevance of the LGM experiment for CMIP6 The LGM experiments are directly relevant to CMIP6 questions 1 and 2 (Eyring et al, 2016): "How does the Earth System respond to forcing?", and "What are the origins and consequences of systematic model biases?".
1. What are the responses of the Earth system to the LGM forcings?
In the following, we use the word "forcing" from the point of view of the CMIP6-type climate models. We include GHG and ice sheets in this term as these are prescribed in the CMIP6-PMIP4 LGM experiments, even though these are interactive components of the full climate system. Our current understanding of the LGM climate is then based on the response of the Earth system to the following forcings: decreased atmospheric GHG concentrations, and impacts of the ice sheets and associated changes in topography, bathymetry, coastlines, and Earth surface types on the atmosphere and the ocean. The change in GHG is well constrained, but there are non-negligible differences in ice-sheet reconstructions and a major goal in PMIP4-CMIP6 is to explore the impact of these differences on climate. Differences between the ice sheets are expected to cause differences in climate above and around the ice sheets (e.g. Löfverström et al., 2014Löfverström et al., , 2016, but also at a larger scale if the changes in large-scale circulation are sufficiently large to have an impact on the North Atlantic Ocean circulation (e.g. Roberts et al., 2014;Ullman et al., 2014;Zhang et al., 2014;Beghin et al., 2016). Several studies have shown that changes in vegetation cover and increases in dust loading affect LGM climates (e.g. Maher et al., 2010;Albani et al., 2014). The design of the PMIP4-CMIP6 simulations allows the impact of vegetation and dust forcing to be explored systematically. The Tier 2 PMIP4 sensitivity experiments will separate the influence of individual forcings (GHG and ice sheets) on the LGM climate. Thus, the PMIP4-CMIP6 Tier 1 LGM experiment, and the associated Tier 2 sensitivity experiments, will help to understand the response to multiple forcings, the sensitivity to individual forcings, and how the responses to individual features and forcings combine to produce the full LGM response.
2. Can models represent the reconstructed climatic and environmental changes for the LGM?
Model evaluation based on LGM climate or environmental reconstructions has been an ongoing activity since the beginning of PMIP (Braconnot et al., 2012;Harrison et al., 2014Harrison et al., , 2015Annan and Hargreaves, 2015). Model-data comparisons have been performed at data sites and this has helped identify discrepancies in the LGM experimental set-up (e.g. for the eastward extension of the Fennoscandian ice sheet which had a strong impact on summer temperatures, ). Data-model comparison has helped to establish the realism of large-scale climatic features, such as polar amplification, land-sea contrast and precipitation scaling with temperature (Masson-Delmotte et al., 2006;Izumi et al., 2013;Li et al., 2013;Lambert et al., 2013;Schmidt et al., 2014), and the ocean state and deep circulation (e.g. Otto-Bliesner et al., 2007;Muglia and Schmittner, 2015;Zhang et al., 2013). Benchmarking in comparison to paleoclimatic surface reconstructions (land and ocean) has shown there has been little improvement from PMIP2 to PMIP3, especially at the regional scale Annan and Hargreaves, 2015). However, in PMIP4, given improvements in the climate models themselves, the inclusion of additional boundary conditions (dust, vegetation) and updates to pre-existing boundary conditions (e.g. ice sheets, river routing, GHGs) in line with latest knowledge, the simulations of regional climate should be more realistic. In addition, models now explicitly represent processes or climate system components such as marine biogeochemistry, oxygen and carbon isotopes, dust emission and transport, and vegetation dynamics, making it possible to make direct comparisons with environmental records and reducing the uncertainties re-sulting from the interpretation of these records in terms of climate signals in model-data comparisons. An important aspect of the data-model comparisons will be to determine whether there is sufficient data to characterize and quantify differences in regional climates resulting from the uncertainties in the imposed boundary conditions (i.e. different ice sheets, different representations of vegetation and/or dust forcing).
3. What are the roles of each component of the climate system, or of specific processes within the climate system, in producing the LGM climate?
The LGM climate is the result of a combined set of forcings and feedbacks. In particular, decreased GHG and increased dust act on the atmospheric radiative forcings and feedbacks; changes in sea ice provide a feedback to atmospheric radiation, atmosphere-ocean exchanges, and ocean circulation (deep water formation); the ice sheets and vegetation changes act on the albedo and surface energy fluxes; ice-sheet topography, decreased sea level, and modified bathymetry act on the atmospheric and oceanic circulations; the decreased atmospheric CO 2 concentration acts on vegetation and the way it exchanges water and CO 2 with the atmosphere via changes in water-use efficiency. Thus, much can be learnt about the respective role and magnitude of key feedbacks affecting Earth's energetics by analysing the PMIP4-CMIP6 LGM experiments, as well as the PMIP4 sensitivity experiments. These were developed following a few studies led with single models (e.g. Klockmann et al., 2016;Brady et al., 2013;Pausata et al., 2011). We also expect analyses of the impacts of these LGM forcings to strongly benefit from diagnostics developed by the Modelling Intercomparison Projects (MIPs) dedicated to these components and processes, such as OMIP for the ocean (Griffies et al., 2016;Orr et al., 2017), SIMIP for sea-ice processes (Notz et al., 2016), LS3MIP for the land surface (van den Hurk et al., 2016), AerChemMIP for dust (Collins et al., 2017), CFMIP for clouds (Webb et al., 2017), and RFMIP for radiative forcing diagnostics (Pincus et al., 2016).

Can the LGM climate constrain climate sensitivity?
The amplitude of the temperature change from the LGM to the pre-industrial state is of the same order of magnitude as climate warming projected for the end of the 21st century. The potential of the LGM reconstructions for constraining climate sensitivity has been shown, with climate models of intermediate complexity (Schneider von Deimling et al., 2006;Schmittner et al., 2011) as well as with CMIP-type models Hargreaves et al., 2012;Annan and Hargreaves, 2015). These studies, as well as Schmidt et al. (2014), point to the LGM tropical SSTs in particular. This would strongly benefit from progress on re-constructions of these tropical SSTs, which have been strongly debated since CLIMAP (1981) and still show discrepancies (as summarized in e.g. Annan and Hargreaves, 2015). On the other hand, the studies using CMIP-type models have shown that more individual simulations than presently available are required to establish statistically significant relationships. Analysis of the processes involved in the temperature response to the forcings (i.e. GHG for current-to-future warming, and ice sheets and GHG for the LGM-to-pre-industrial warming) are essential for this investigation, because while some feedbacks appear to work in a similar manner for LGM-to-pre-industrial and for future warming, feedbacks such as the cloud radiative feedback do not (Yoshimori et al., 2009). The relative magnitudes of the different feedbacks also vary between those two climates, so that the relationship with climate sensitivity is not straightforward . Changes in vegetation and dust, which produce changes in regional climate, also need to be taken into account when regional reconstructions (such as over the tropical oceans) are used to constrain climate sensitivity (Hopcroft and Valdes, 2015a). By increasing the number of simulations available, including important regional forcings, and focusing on uncertainties in these forcings, the LGM PMIP4-CMIP6 experiments will provide a much better data set to re-examine climate sensitivity.

PMIP4-CMIP6 experiments and PMIP4 sensitivity experiments
This section describes the PMIP4-CMIP6 Tier 1 LGM climate experiment, termed "lgm", as well as complementary PMIP4 Tier 2 sensitivity experiments. These are also summarized in Table 1. Section 4 describes how to implement the associated boundary conditions.

The Tier 1 PMIP4-CMIP6 lgm experiment
The lgm simulation is a CMIP6 Tier 1 experiment, as well as one of the two possible PMIP4 entry cards (i.e. one of the two simulations that must be performed by modelling groups wishing to officially take part in PMIP4). The lgm simulation will be compared to the CMIP DECK (Diagnostic, Evaluation and Characterization of Klima) pre-industrial control (piControl) for 1850 CE and the CMIP6 historical experiment (Eyring et al., 2016) and must therefore be run using the same version (including level of complexity and the interactive feedbacks) and resolution of the model and following the same protocols for implementing external forcings as in these two reference simulations. The minimum set of changes that must be made for the lgm simulation, compared to the set-up of the piControl, are the insolation, GHG, and ice-sheet forcings (see Sect. 4 for the implementation of these changes). There are several plausible alternatives for the ice-sheet forcing and modelling groups can choose between one of three options ( Fig. 1): the ice-sheet reconstruction produced for PMIP3 (Abe-Ouchi et al., 2015), the ICE-6G_C reconstruction Argus et al., 2014), or the GLAC-1D (Tarasov et al., 2012;Briggs et al., 2014;Ivanovic et al., 2016). However, if the PMIP4 transient last deglaciation experiment is run , the modelling groups should ensure consistency between the LGM simulation and subsequent transient phase of the experiment, when possible.
The dust and vegetation forcing in the Tier 1 lgm experiment must be imposed in a manner that is consistent with the DECK simulations. Models that include interactive dust, for example, should allow interactive emissions at the LGM. For this purpose, two alternative reconstructions of LGM dust emission regions are provided for models without dynamic vegetation (Hopcroft et al., 2015;Albani et al., 2016: see the PMIP4 website) and modelling groups are free to choose either one of these. If dust-enabled models do not include dynamical vegetation, then vegetation should be changed in the LGM dust emission regions so that dust emission can occur (e.g. by imposing bare soil or a fractional grass cover). Both dust data sets provide atmospheric mass concentrations, which could alternatively be used to compute a corresponding radiative forcing in a consistent manner as for the reference simulations. Modelling groups can also use a climatology of atmospheric dust mass concentrations produced offline by their own dust model, using dust emission regions and vegetation as above. Otherwise the lgm simulation should be run using the same forcing as for the DECK and historical runs (i.e. with no increase in dust). Unless a model includes dynamic vegetation or interactive dust, the vegetation should be prescribed to be the same as in the DECK and historical runs.
The relative flexibility of the set-up summarized above reflects the range of model configurations foreseen for CMIP6 and PMIP4, in particular in terms of the representation of vegetation and dust. In the case of the ice sheets, it also reflects the uncertainties in our knowledge of the boundary conditions. Taking this uncertainty into account is new to the PMIP LGM experimental design but is essential for evaluating the CMIP6-type models. The differences between the reconstructed ice-sheet altitude can be as large as several hundred metres (e.g. over the North American ice sheet), which can be enough to induce differences in the Atlantic jet stream (Beghin et al., 2016).

Sensitivity to vegetation and dust
Experiments designed to test the sensitivity of the LGM climate to vegetation and dust, run with model versions or set- LGM simulation with interactive vegetation, with no corresponding DECK simulation, then this is also considered a PMIP4 sensitivity run, which should be named lgm_v2 (2 for PMIP4 only).
Simulations with or without changes in dust are already included in the PMIP4-CMIP6 protocol, so the sensitivity to dust can be analysed through these simulations. However, if a modelling group runs an LGM experiment with interactive dust but with no corresponding DECK simulation, this simulation would be a PMIP4 sensitivity experiment, named lgm_d2n, with n varying according to the data used to set the emission regions (see Sect. 4.11). Sensitivity experiments with vegetation and dust different from the PMIP4-CMIP6 simulations should be named lgm_vm_dn, with m = 1 or 2 and n defined according to the definitions above.
Experiments made with a different version or resolution of model from the DECK and historical simulations will also be considered as PMIP4 sensitivity simulations. In addition to running the lgm and pre-industrial experiments with this different model resolution or version, it would be extremely useful to run an abrupt4xCO2 experiment so that the LGM-to-pre-industrial change can be compared to the preindustrial-to-"future" climate change (cf. Fig. 6, Kageyama et al., 2016).

Sensitivity to individual forcings (Tier 2 experiments)
A series of three additional experiments have been designed to disentangle the impact of individual changes in boundary conditions, and thus facilitate the interpretation of the LGM Tier 1 experiment. All three experiments will use the LGM land-sea mask and astronomical parameters, but will use different combinations of ice-sheet and GHG forcings. The experiments are M.  LGM experiment and PMIP4 sensitivity experiments the LGM-PI-ghg experiment, in which all boundary conditions and forcings are set to LGM values except for the GHGs, which are the same as in piControl; the LGM-PI-ice experiment, in which all boundary conditions and forcings are set to LGM values except for the ice-sheet extent, and height is the same as in piControl; and the LGM-PI-ghg_ice experiment, in which all boundary conditions and forcings are set to LGM values except for the GHGs and ice-sheet extent and height, which are the same as in piControl.
Comparison of these sensitivity experiments will allow the impacts of the atmospheric GHG decrease and of the icesheet albedo and topography changes to be disentangled. Provided they are each run to equilibrium, they can be directly compared to the full lgm experiment, allowing the relative importance of different aspects of the change in forcing to be quantified (see e.g. Hewitt and Mitchell, 1997).
4 The lgm experiment: implementing the boundary conditions and model spin-up Table 1 summarizes the implementation of the boundary conditions and forcings and gives check points for each of them.

Atmospheric trace gases
The concentrations of the atmospheric trace gases should be set to -190 ppm for CO 2 , -375 ppb for CH 4 , -200 ppb for N 2 O, and -0 for the CFCs.
-Ozone should be set to its piControl value.
These concentrations have been updated from the PMIP3 values for consistency with the deglaciation protocol , which is based on data from Bereiter et al. (2015) for CO 2 , Loulergue et al. (2008)

Insolation
The astronomical parameters should be set to their 21 ky BP values, according to Berger (1978): The resulting insolation at the top of the atmosphere should then be similar to that displayed in Fig. 2, with a decrease at high latitudes during the summer hemisphere reaching over 10 W m −2 and a mild increase (reaching 3 W m −2 ) between October and April at 40 • N, December and June at the Equator, and mid-January to August at 40 • S.

Ice sheets
The ice sheet can be set to one of the following reconstructions ( Fig. 1): GLAC-1D (Tarasov et al., 2012;Briggs et al., 2014;Ivanovic et al., 2016), ICE_6G-C Argus et al., 2014), or PMIP3 (Abe-Ouchi et al., 2015. GLAC-1D and ICE_6G-C are the most recent reconstructions and are compatible with the set-up of the PMIP4 deglaciation simulation . The use of the PMIP3 ice-sheet reconstruction allows direct comparison with the PMIP3 simulations. These ice-sheet reconstructions significantly differ with each other, in particular in terms of altitude, with differences reaching several hundred metres over North America and Fennoscandia ( Fig. 1 and Ivanovic et al., 2016, Fig. 2). This uncertainty in the boundary conditions results from the different approaches used for the reconstructions, which are summarized in Ivanovic et al. (2016). The implementation of the LGM ice sheets will vary from one model to the other. Here, we give the main implementation steps that have been followed for the IPSL climate Step 0 Step 1 Figure 3. Summary of steps to be followed for the definition of the basic surface types for the atmosphere and ocean boundaries. model (Fig. 3). The details of the implementation may differ for other models, but the same steps should be followed and documented.
The PMIP3 reconstruction files include information about the land fraction (sftlf for surface type land fraction), landice fraction (sftgif for surface type glacier fraction), and difference in orography (diff_orog) that needs to be applied to the piControl orography in order to obtain an lgm orography. This information, in particular sftlf, is not directly available in the GLAC-1D reconstruction and is incomplete in the ICE-6G_C reconstruction (e.g. the Caspian Sea, above the present-day sea level, is missing). The variables available for each reconstruction are listed in Table 2. They can be found on the PMIP4 website (http://pmip4.lsce.ipsl.fr), as provided by the authors of the reconstructions. In particular, the variable names and the resolution have not been modified. In the present step 0, we describe how we compute sftlf, sftgif, and orog from the GLAC-1D and ICE-6G_C data. The IPSL model requires orog at 1/6 • resolution for its gravity wave drag parameterization, which is why we compute diff_orog at this high resolution.
The procedure is as follows ("Pre-pare_LGM_BC_files.py" Python script provided on the PMIP4 website): the input variables listed in Table 1 are read in; these include the land fraction for ICE-6G_C but not for GLAC-1D; for GLAC-1D, the land-sea mask for the present and for the LGM are defined as where topography is positive; small holes (usually one to two isolated grid boxes) in the land-ice fraction are filled, using the "bi-nary_fill_holes" function of the python scipy/ndimage package (for ICE-6G_C, 155 points are filled in, to be compared to the total number of land-ice grid points, with is initially 423 610; for GLAC_1D, 62 points are filled in; the total number of points fully covered by land ice is 23 348); the land fraction is updated to include the land-ice fraction; this land fraction includes unrealistic isolated continental points which are well below sea level (we have considered a threshold of −500 m. There are 23 such points in the ICE-6G_C case, 4 in the GLAC_1D case). These points are filled in using the same function as for the land-ice mask. However, several straits must be re-opened so that the function does not fill in the Red Sea, the Black Sea, the Azov Sea, the Sea of Japan, the Mediterranean Sea, and additionally the Persian Gulf, the Baltic and White seas, the Great Lakes, and the Canadian Archipelago for the present day. "bi-nary_fill_holes" is applied with the appropriate straits opened; then, these are closed again. sftlf is computed following this method for both the present and the LGM; the topography of the points that have been filled in is corrected by averaging the topography of the surrounding points, after removing points well below sea level; and the topography on the continents can be defined for the present and the LGM, and the difference in orography diff_orog can be computed. Similarly, differences in bathymetry can also be computed.
This preliminary step provides the three variables that are necessary to modify the boundary conditions for the atmosphere and the ocean: the land-sea mask, the land-ice mask, and the difference in topography and bathymetry. For the IPSL model, we keep the LGM orography computed at this step for further use.
Step 1. Defining the land-sea mask and the land-ice mask within the climate model In the IPSL model, the coastlines are defined first for the ocean model and then they are used to compute the fraction of land and ocean on the atmospheric grid. We will therefore follow this order here. The procedure is summarized in Fig. 3.
The land-sea mask obtained at the end of step 0 is interpolated on the ocean grid. A threshold of 0.5 is chosen to determine the coastlines. After this first interpolation, the basic features of the LGM coastlines can be checked: presence Table 2. Variables provided with the ice-sheet reconstructions considered for PMIP4.

GLAC-1D ICE-6G_C PMIP3
-HDC: on continents (including ice sheets) and ice shelves: surface altitude (including ice sheets/shelves) on ice-free ocean: bathymetry -HDCB: on continents (including ice sheets) and ice shelves: surface altitude (including ice sheets) of land at locations of the main ice sheets, especially over areas that were glaciated at the LGM but that are covered by oceans today (such as Hudson Bay and the Barents-Kara seas); closure of the Bering Strait, of the straits between the Mediterranean and Black seas, and of the Sahul and Sunda shelves. At this stage, we re-introduce the Caspian Sea in the land-sea mask, using the present-day Caspian Sea. The Caspian Sea is absent from the land-sea masks computed from step 0 because it is higher than global sea level at the LGM. These basic coastlines need polishing, as a function of the ocean model, in order for ocean transport to occur in narrow straits. In particular, the connection from the Red Sea to the Arabian Sea should be checked, as well as of the Sea of Japan to the Pacific Ocean and narrow passages between the Sunda and Sahul shelves. This is detailed for the NEMO ocean in Program 2 given in the Supplement. Once the ocean boundaries are set up, these can be interpolated over the atmospheric grid. The weights required to pass from one grid to the other are computed at the same time.
The land-ice cover is interpolated directly on the atmospheric grid and multiplied by the land-sea mask so that no land ice is defined over the ocean. This might differ for models including a representation of ice shelves.
At the end of step 1, the coastlines are defined for the ocean model, and the land-ice and land-sea masks are defined for the atmospheric model.
Step 2. Implementing the LGM orography The LGM orography is implemented by adding the LGMpresent anomaly in orography computed in step 0 to the pi-Control orography. This is straightforward for models that only require the average orography for each grid point. Additional steps are required for models requiring secondorder moments/minimum/maximum values/slope characteristics for each grid point (e.g. the parameterization proposed by Lott and Miller, 1997). These moments must be computed from a high-resolution orography data set and the anomaly method should be applied for this high-resolution data set before computation of the parameters depending on fine-scale orography. The ice-sheet orography needs to be smoothed before this computation is made, to prevent unrealistic parameters due to the present-day orography (Fig. 4 illustrates the impacts of smoothing the topography for the north-western part of North America). These steps are detailed in program "Prepare_LGM_BC_files.py" (at step 6) given in the Supplement for the LMDZ model. The smoothing is performed with the Gaussian filter provided in the ndimage package, with sigma = 3.
Step 3. Implementing the LGM bathymetry There are two options for implementing the changes in bathymetry. The first option is to use the bathymetry anomalies obtained from step 0 directly and add them to the bathymetry used for the piControl simulations. However, given that the resolution of the ocean models often decreases with depth, this may not be necessary, and a simpler option is to modify the present-day bathymetry by subtracting the mean sea-level drop corresponding to the chosen ice-sheet reconstruction. In this second option, special treatment will be required for straits that are crucial for the ocean circulation and for which the change in bathymetry is significantly different from the mean sea-level drop. The Denmark and Davis straits and the Iceland-Faeroe Rise, for example, must be treated with care, as these are often locations at which the bathymetry for piControl is also adjusted to obtain realistic oceanic currents. The second option is used for the IPSL model, and the corresponding program is provided in the Supplement (program "bathy_lgm.py"). The results are shown in Fig. 5 for the NEMO model used in the IPSL coupled model. Figure 5a shows the changes for global ocean extent, with ocean points disappearing at the location of the LGM ice sheets (e.g. Hudson Bay, Baltic Sea, Barents-Kara seas) and where the present ocean is shallow enough to be sensitive to the LGM sea-level drop (e.g. Bering Strait, Sunda and Sahul shelves, north of Siberia, off the Patagonian eastern coast). Figure 5b shows the global changes in ocean bathymetry, which for this example have been prescribed at the average sea-level drop for the ICE-6G_C reconstruction.  area. We have ensured that the imposed change in bathymetry matches the reconstructed one for the Greenland-Iceland Rise and Iceland-Faeroe Rise.

Freshwater budget: rivers, runoff, and accounting
for positive snow mass balance over the ice sheets The LGM sea-level drop leads to expanded continents and this can mean that prescribed river courses no longer reach the ocean. The North American and European ice sheets also disrupt river courses. At a minimum, the LGM rivers must be set up to ensure they reach the oceans. For instance, the European rivers that today drain into the Nordic and Baltic seas can be redirected to the North Atlantic via the paleo-English Channel (see e.g. Alkama et al., 2006). More realistic river-routing files compatible with the ice-sheet reconstructions will also become available at a later stage. It is highly possible that the snow mass balance over the ice sheets is positive, resulting in a sink of freshwater in the climate model. If this is the case, the average value of the sink (e.g. the average for a 10-year period) should be computed and released to an adjacent ocean, to guarantee closure of the freshwater budget. This should be done following the same procedures as for the DECK experiments or following the procedure advised since PMIP2, which was to compensate for the sink of freshwater by imposing a freshwater flux in broad regions of oceans adjacent to the ice sheets (e.g. the Arctic and North Atlantic north of 40 • N for the North American ice sheet). As this decision might have a large impact on the global ocean overturning circulation, it must be precisely documented (cf. Sect. 4.10).

Vegetation
Models including dynamical natural vegetation should use the corresponding module on all unglaciated continents, in the same way it is used for natural vegetation in other CMIP6 simulations. Modelling groups who do not run with dynamical natural vegetation should use the same vegetation cover as for piControl, extrapolated to the lgm land mask, in their PMIP4-CMIP6 experiment. There is insufficient information to construct a reliable global map of vegetation at the LGM, but one way to take account of LGM vegetation changes in models without dynamic vegetation is to run a biogeography or dynamical vegetation model offline, using climate forcing from the LGM simulation, and to then prescribe the simulated vegetation patterns in the coupled climate model. This ensures that the prescribed vegetation will be consistent with the climate forcing for the given model. These simulations will then be PMIP4 sensitivity experiments (cf. Sect. 3.2.1). A minimum change for models with interactive dust modules will be to remove vegetation from (or only to allow grass in) regions of strong potential dust emissions (cf. Sect. 4.6 below).

Mineral dust
There are several options for implementing dust forcing according to the model's complexity and to the availability of different data sets. Three series of dust data sets are provided. Two of them are based on model simulations (Albani et al., 2014Hopcroft et al., 2015). Both models include a prognostic dust cycle, based on different formulations of the dependency of emissions on wind speed, soil moisture, and vegetation cover arising from the work of Marticorena and Bergametti (1995) and Fécan et al. (1999). In one case (Albani et al., 2014) pre-industrial vegetation is prescribed for physical climate for both PI and LGM climate conditions, but LGM dust emissions at each grid cell are scaled by the non-vegetated fraction, resulting from an offline vegetation reconstruction with BIOME 4 (Kaplan et al., 2003), in equilibrium with LGM climate conditions. In the other case (Hopcroft et al., 2015) a dynamical vegetation model was used to determine the erodible surface. These differences result in different dust emission fields. Furthermore, Albani et al. (2014 further refined their dust emissions by scaling the soil erodibility at the continental scale in order to have a better match to paleodust observations in terms of deposition fluxes. The third data set (Lambert et al., 2015) is a reconstruction of dust deposition, essentially based on geostatistical interpolation of paleodust observations. The three data sets have different specifications in terms of dust size distribution: four size bins spanning 0.1-10 µm diameter (Albani et al., 2014), six size bins spanning 0.0316-31.6 µm radius (Hopcroft et al., 2015), and bulk i.e. integrated over the entire observed size range (Lambert et al., 2015). This has implications for imposing proper constraints on the global dust cycle, e.g. magnitude of emissions (Albani et al., 2014), as well as for dust radiative forcing, when considered in combination with the prescribed dust optical properties (Kok et al., 2017). Therefore modelling groups should carefully account for this aspect when integrating one of these data sets into their model framework.
For models with interactive dust modules but without dynamic vegetation, it is advisable to take into account the more extensive dust sources at LGM. These are described by the "erodibility map" from the Albani et al. (2016) data set and a bare soil map for the Hopcroft et al. (2015) data ( Fig. 6a and b, respectively). For these regions, vegetation must be set to either low vegetation (grasses) or bare soil; otherwise, the source functions should be adapted depending on the precise formulation of the dust emission module in the particular model (e.g. Ginoux et al., 2001) so that dust emissions are allowed. For models that compute the dust radiative forcing from atmospheric dust mass loading, two data sets are available for the LGM: Albani et al. (2014Albani et al. ( , 2016 and Hopcroft et al. (2015). The prescribed LGM mass loading should be implemented as perturbations of the piControl loading, i.e. by either adding an anomaly to these piControl loads or by multiplying them by a ratio, the anomaly, or the ratio being computed from the Albani et al. (2014Albani et al. ( , 2016 or Hopcroft et al. (2015) data sets. Alternatively, modelling groups can compute their own atmospheric dust mass loads offline and use them as prescribed fluxes in their coupled simulations.
Both dust data sets also provide radiative forcing. These should not be used directly because the specified radiative properties of dust vary among models, and using the forcing from the models used to produce the dust fields would be incompatible with other CMIP6 experiments. The dust radiative forcing provided with the data sets is only given with the purpose of broad comparison with the modelling groups' own model output (Fig. 6c and d).
Models including marine biogeochemistry should use LGM dust deposition on the oceans, using the same data set as for the atmospheric forcing ( Fig. 6e and f). If LGM dust atmospheric forcing cannot or is not taken into account, then the Lambert et al. (2015) data set can also be chosen (Fig. 6g).
Modelling groups undertaking the implementation of dust in their models are advised to perform a first trial with an atmosphere-only simulation, as run-away effects involving dust, vegetation, and climate have been experienced by some modelling groups (Hopcroft and Valdes, 2015b). In the latter case, it was the choice of parameters in the dynamic vegetation model which proved to be inadequate.

Other inputs for ocean biogeochemistry models
The global amount of dissolved inorganic carbon, alkalinity, and nutrients should be initially adjusted to account for the change in ocean volume. This can be done by multiplying their initial value by the relative change in global ocean volume. Other features that may need adjustment, given the changes in coastlines and bathymetry, include the amount of nutrients brought by rivers and by boundary exchange at the ocean-sediment interface. Modelling groups must document any such changes in the description of their simulations (cf. Sect. 4.10).

Initialization and spin-up
First, it is suggested to run the atmosphere model separately, using the sea surface temperatures and sea ice from the ocean's initial conditions, in order for the atmosphere to adjust to the topography and surface-type changes. At this stage, it is advised to check that the total atmospheric mass (or globally averaged surface pressure) is the same as for piControl. This run will yield an initial state for the atmospheric component of the model.
The ocean should be initialized with a salinity 1 psu higher than for piControl, which is consistent with the sea-level difference between LGM and piControl (and the volume of freshwater stored in the ice sheets). Similarly, ocean biogeochemistry models should adjust their alkalinity and models including oxygen isotopes should initialize them with a  (Hopcroft et al., 2015). Maps of LGM dust aerosol optical depth (AOD) from the simulations of (c) Albani et al. (2014) and (d) Hopcroft et al. (2015). Maps of LGM dust deposition (g m −2 a −1 ) (e) simulated with the Community Earth System Model (Albani et al., 2014), (f) simulated with the Hadley Centre Global Environment Model 2-Atmosphere (Hopcroft et al., 2015), and (g) reconstructed from a global interpolation of paleodust data (Lambert et al., 2015).
Standard Mean Ocean Water (SMOW) of +1 ‰. The ocean model can be initialized from a piControl experiment or from previous LGM experiments, to minimize spin-up duration.
Practically, the ocean model can be generally initialized from a piControl ocean state with adjusted salinity (and oxygen isotope, if applicable), or from previous LGM experiments (e.g. with well-stratified glacial ocean states), to minimize spin-up duration. Such ocean states, such as described in Werner et al. (2016), which provide 3-D fields of sea temperature, salinity, and associated stable water isotopes on a regular 1 • × 1 • grid, are available on the PMIP4 website and from the PMIP2 and PMIP3 databases.
The model should be spun up until equilibrium. In previous PMIP protocols (in particular http://pmip2.lsce.ipsl.fr) the simulations were considered at equilibrium when the trend in globally averaged SST was less than 0.05 • C/century, the Atlantic Meridional Overturning Circulation (AMOC) was stable and, for models including representations of the carbon cycle or dynamic vegetation, the requirement was that the carbon uptake or release by the biosphere is less than 0.01 Pg C per annum. Recent works give other criteria or recommendations for defining or reaching the equilibrium. For instance, to avoid impacts of potential transient character-istics in the deep ocean on AMOC strength , the equilibrium ocean should ensure that the trend in zonal mean sea salinity in the Southern Ocean (south of the winter sea-ice edge) remains small, especially in the Atlantic sector. Marzocchi and Jansen (2017) show that the AMOC has to be monitored on multi-centennial timescales because variability on the timescales of decades to a century prevents a precise determination of the trends, and hence of whether the model is close to equilibrium or not.
It is required that at least 100 years of data from the equilibrated part of the simulation is stored on the ESGF (Earth System Grid Federation). In order to document the approach to equilibrium, we recommend that the modelling groups monitor the variables listed in Table 3 and save them for model documentation (cf. Sect. 4.10). We recommend these data be saved for at least a few hundred years before the period stored on the ESGF and for the ESGF period, so that the trends in these variables can be better determined. This will help characterize the "ESGF period" within the full simulation and make us aware of possible remnant drifts. The same variables should also be provided for the corresponding pi-Control simulation, for comparison. failure to close the freshwater budget, which can arise from either inadequate compensation for a positive snow mass balance over the ice sheets or from rivers not reaching the ocean, numerical instabilities in the atmosphere, especially near or above the ice sheets, and run-away cooling due to climate-vegetation-dust feedbacks, as reported by Hopcroft and Valdes (2015b). In this case the dynamic vegetation scheme was found to be overly sensitive to temperature, so that grass plant functional types started to die back below 5 • C, resulting in higher albedo, further cooling, and eventual desertification across most of Eurasia in the first LGM simulation with HadGEM2-ES.

Documenting the simulations
The documentation of the simulations should include the model version used, in particular in terms of vegetation and dust representations (interactive, prescribed, or absent), the ice-sheet reconstruction chosen and how it has been implemented, how river routing has been modified and how positive snow mass balance over the ice sheets is dealt with, in particular the regions over which the excess freshwater is applied, the vegetation used in the simulation and how it was obtained and/or implemented, the dust reconstruction used and how it has been implemented, the forcings used (dust, nutrients from rivers and sediments) if ocean biogeochemistry is included in the model, and the spin-up strategy and duration, with documentation of the variables listed in Table 3.
A PMIP4 special issue in GMD and Climate of the Past is open so that groups can publish these documentations. Modelling groups are also encouraged to contribute their simulation and model documentation to the ES-DOC facility.
4.11 "ripf"code for the simulations CMIP6 simulations can be documented through their "ripf" code, these letters standing for "realization", 'initialization", "physics", and "forcing". In practice, each of these letters is followed by a number which indicates after the "r": the simulation number in the ensemble of simulations with the same characteristics; after the "i": the initial method; after the "p": the chosen model's physics; and after the "f": the forcing used for the simulation.
Since there are multiple choices for setting up PMIP4-CMIP6 and PMIP4 LGM experiments, we propose the systematic use of common "f" indices within the CMIP6 "ripf" indices so that the simulations can be distinguished easily from each other. The first digit should describe the ice-sheet reconstruction. It should be set to The third and fourth digits should describe how dust is included in the set-up.
-If no dust forcing can be taken into account, they should be set to 00.
-If dust is prescribed from a PMIP4 data set, they should be set to -11 for the Albani et al. (2014Albani et al. ( , 2016) data set, -12 for the Hopcroft et al. (2015) data set, -13 for the Lambert et al. (2015) data set (for ocean biogeochemistry models only), and -19 for the modelling group's own dust forcing.
-If dust is interactively computed, they should be set to -20 if the surface maps are dynamically simulated using a coupled dynamic vegetation scheme, -21 if the surface maps for emissions are those from Albani et al. (2014Albani et al. ( , 2016, -22 if the surface maps for emissions are those from Hopcroft et al. (2015), and -29 if the surface maps for emissions are produced by the modelling group itself, e.g. by using an offline vegetation model.

Output
The data should be formatted so as to comply with the CMIP6 standards (to be documented in the GMD CMIP6 special issue; cf. Eyring et al., 2016) and PMIP4 data request  so that analyses including other PMIP and CMIP6 simulations can be performed easily. The current list of variables is given in the Supplement but is still subject to potential changes following adjustments of the full CMIP6 list. The PMIP4 data request can be found on the PMIP4 website (https://pmip4.lsce.ipsl.fr/doku.php/ database:pmip4request#the_pmip4_request).

Analyses and outlook
The LGM experiment is a major investment by climate modelling groups, but provides a demanding test of model reliability under extreme and well-documented conditions. Indeed, our experience is that several groups have found model errors while setting up their LGM climate simulations, in particular in the coupling between the atmosphere and the ocean and in the global freshwater budget. The PMIP4-CMIP6 simulations, along with PMIP4 sensitivity experiments and previous PMIP2 and PMIP3 experiments, will create an unprecedented data set about the LGM climate state. With a larger number of simulations, and a better sampling of the forcing uncertainties, we should be able to reach more robust conclusions about, for example, the ability of state-of-the-art climate models to represent a climate very different from the pre-industrial or present climates: benchmarking these simulations will provide a measure of how well models simulate large climate changes, comparable in magnitude to changes expected over the 21st century. Although there are data sets documenting environmental conditions and climate at the LGM, the planned PMIP4-CMIP6 analyses would benefit from the improvement and geographic expansion of these data sets. In addition, there is scope for the creation of new data sets, particularly data sets that can be used to evaluate aspects of the more complex Earth system models that are being run in PMIP4-CMIP6; the relationships between climate or environmental changes at far away locations, or between different features of the climate system: for instance, as alluded to in the introduction, we expect the atmospheric and ocean circulations in the North Atlantic area to be sensitive to the ice-sheet height; the PMIP4-CMIP6 experimental design allows for multi-model studies on this topic; at large scales, the polar amplification and landsea contrasts that have been studied with PMIP2 and PMIP3 experiments could be altered with the PMIP4 more complex simulations including vegetation or/and dust changes; and the potential constraint from the LGM (in particular via the LGM tropical SSTs) on climate sensitivity.
The Tier 2 sensitivity experiments will allow the quantification of the role of individual forcings and feedbacks in climate. This is an essential step in understanding the LGM climate, but also in characterizing and understanding common and/or contrasting features of the most recent past warming (between the LGM and the present) and the predicted future warming.
These are a few examples of possible analyses of the PMIP4-CMIP6 lgm simulations. The analysis of the PMIP4-CMIP6 and PMIP4 sensitivity experiments also relates to other CMIP6 projects and we hope these data will also be analysed by experts from other CMIP6 MIPs. For instance, the understanding of the impacts of the LGM climate forcings and the role of radiative feedbacks is related to CFMIP (Webb et al., 2017) and RFMIP (Pincus et al., 2016). The PMIP4 single forcing experiments can be used in view of the CFMIP experiments testing the impact of uniform lowering of SSTs or CO 2 decrease (in AMIP configuration) and the connection to climate sensitivity for CO 2 increase should be made easier to analyse with these experiments. In terms of diagnostics that can be used to analyse the role of each component of the climate models in setting up the LGM climate, we also expect new studies based on diagnostics developed by the CMIP6 MIPs on the ocean (OMIP, Griffies et al., 2016;Orr et al., 2017), land surface and snow (LS3MIP, van den Hurk et al., 2016), aerosols (AerChemMIP, Collins et al., 2017), sea ice (SIMIP, Notz et al., 2016), and ice sheets (IS-MIP6, Nowicki et al., 2016). It is therefore important to keep the relevant output for these analyses, and the PMIP4 data request has been built based on the lists for these other MIPs.
LGM experiments will also be the starting point for simulations of the last deglaciation, i.e. the transition from the full glacial state to the present interglacial state (21-9 ky BP) and through to the present . Given the large and abrupt changes in AMOC during the glacial period and during the deglaciation, the LGM will also be a reference state for freshwater hosing studies, which will allow further analyses of the relationships between the AMOC state and climate. This is relevant for studying the processes at work during Heinrich events 1 and 2, but also for establishing a coherent view of the LGM physical climate system state throughout all its components. Stated in a different manner, analysing lgm simulations characterized by different AMOCs, obtained through freshwater hosing or not, will help to determine whether all the reconstructions that are available for the different components of the climate system (state of the AMOC, state of the ocean surface, state of the continental surface) are consistent from the point of view of the physics summarized in a climate model, as suggested by studies carried out with one model Klockmann et al., 2016).
This brief outline of possible analyses of the PMIP4-CMIP6 lgm simulations is not meant to be exhaustive, but rather to illustrate how these simulations will contribute to progress on the overarching questions of CMIP6.
Acknowledging CMIP6 according to the instructions given in Eyring et al. (2016), PAGES, and WCRP, which endorse PMIP, as well as the modelling groups which have contributed to the CMIP6 and PMIP4 effort, will be greatly appreciated.