Using results from the PlioMIP ensemble to investigate the Greenland Ice Sheet during the mid-Pliocene Warm Period. of the Past

. During an interval of the Late Pliocene, referred to here as the mid-Pliocene Warm Period (mPWP; 3.264 to 3.025 million years ago), global mean temperature was similar to that predicted for the end of this century, and atmospheric carbon dioxide concentrations were higher than pre-industrial levels. Sea level was also higher than today, im-plying a signiﬁcant reduction in the extent of the ice sheets. Thus, the mPWP provides a natural laboratory in which to investigate the long-term response of the Earth’s ice sheets and sea level in a warmer-than-present-day world. At present, our understanding of the Greenland ice sheet during the mPWP is generally based upon predictions using single climate and ice sheet models. Therefore, it is essential that the model dependency of these results is assessed. The Pliocene Model Intercomparison Project (PlioMIP) has brought together nine international modelling groups to simulate the warm climate of the Pliocene. Here we use the climatological ﬁelds derived from the results of the 15 PlioMIP climate models to force an ofﬂine ice sheet model. We show that mPWP ice sheet reconstructions are highly dependent upon the forcing climatology used, with Greenland reconstructions ranging from an ice-free state to a near-modern ice sheet. An analysis of the surface albedo variability between the climate models over Greenland offers insights into the drivers of inter-model differences. As we demonstrate that the climate model dependency of our results is high, we highlight the necessity of data-based constraints of ice extent in developing our understanding of the mPWP Greenland ice sheet.


404
A. M. Dolan: Investigation of the Greenland Ice Sheet via results from the PlioMIP ensemble mate and sea level (Church et al., 2013;Masson-Delmotte et al., 2013;Vaughan et al., 2013). Therefore, there is increasing interest in understanding the nature and behaviour of the major ice sheets during warm intervals in Earth history. The Pliocene Epoch is a particularly well documented pre-Quaternary environment which has become the focus for intense study within the Pliocene Model Intercomparison Project (PlioMIP; Haywood et al., 2010Haywood et al., , 2011. A warm period within the Late Pliocene (referred to as the mPWP; 3.264 to 3.025 million years ago; Dowsett et al., 2010) is predicted to have been between 2 and 3 • C warmer than preindustrial Lunt et al., 2010) and estimates of atmospheric carbon dioxide (CO 2 ) concentrations suggest levels of up to 450 ppmv (Pagani et al., 2010;Seki et al., 2010). The IPCC Fifth Assessment Report states with high confidence that global mean sea level was above present (up to 20 m) during warm intervals of the mPWP (Masson-Delmotte et al., 2013) and individual records of sea level high-stands (∼ 20 m) support the reduction in the extent of the ice sheets at this time (e.g. Miller et al., 2012;Rovere et al., 2014;Rohling et al., 2014). Although in recent literature the terms mid-Piacenzian or Late Pliocene warm event have been used to describe this interval, here we retain consistency with the original PlioMIP naming convention when discussing model simulations and use the mPWP.
Proxy records of palaeotemperature derived from ice cores (Dahl-Jensen et al., 1998;Cuffey and Marshall, 2000;Johnsen et al., 2001;Rasmussen et al., 2006) and numerical modelling Overpeck et al., 2006;van de Berg et al., 2011;Born et al., 2012;Quiquet et al., 2013;Stone et al., 2013) of more recent interglacials demonstrate that the Greenland ice sheet (GrIS) has a large sensitivity to high-latitude warming. However, there is little proximal evidence to indicate the volume or extent of the GrIS during the warmest intervals of the mPWP. Evidence of long-lasting subaerial soil formation at the base of the central Greenland Ice Sheet have been suggested as evidence for persistent reduction in Pliocene ice, but the soils have not been positively dated as relating to this period (Bierman et al., 2014). The presence of forest fragments in the Kap København Formation in the far North of Greenland up until 2.4 Ma (Funder et al., 2001) suggests that this area may have been ice-free through intervals of the Late Pliocene. Fragments of evergreen taiga forest in Pliocene sediments at Ile de France (Bennike et al., 2002) also suggest that ice marginal regions were much warmer during the Pliocene. Records from the central Labrador Sea suggest that landmasses adjacent to Greenland, such as Ellesmere and Baffin Island, show a predominance of evergreen forest during intervals of the Pliocene (De Vernal and Mudie, 1989; Thompson and Flemming, 1996;Ballantyne et al., 2006;Csank et al., 2011). Additionally, temperature estimates from peat deposits in the Canadian High Arctic (Beaver Pond) suggest elevated Pliocene Arctic temperatures (Ballantyne et al., 2010). Although there are no mPWP temperature records from Green-land against which to test the climate model simulations, there is generally a cool bias in the models compared to the available data in the northern high latitudes Salzmann et al., 2013).
While useful, proxy evidence is too sparse and uncertain to enable a detailed reconstruction of the extent and location of mPWP ice sheets. Therefore, a variety of modelling frameworks have been adopted in order to simulate the mass balance of the GrIS and reconstruct potential ice sheet configurations during the mPWP (Lunt et al., 2008bHill, 2009;Hill et al., 2010;Dolan et al., 2011;Koenig et al., 2011Koenig et al., , 2014a. These modelling frameworks have generally included the offline coupling of an ice sheet model (ISM) to a climate model, and have been limited to the use of three climate models; the UKMO UM (UK Met Office Unified Model; e.g. Hill et al., 2010;Dolan et al., 2011), GENE-SIS (e.g. Koenig et al., 2011) and multiple versions of CAM (Community Atmosphere Model; e.g. Yan et al., 2014). Although all available simulations suggest that the GrIS was reduced in size during the mPWP, the model dependency of the results is yet to be robustly assessed. The extent to which ice sheet reconstructions are dependent on the ISM employed is addressed through a sub-project of PlioMIP, entitled the Pliocene Ice Sheet Modelling Intercomparison Project (PLISMIP; Dolan et al., 2012). Results from Koenig et al. (2014b) suggest that ISM dependency is low. Here, we will address the question of climate model dependency utilising climate model outputs from PlioMIP (Chan et al., 2011;Bragg et al., 2012;Contoux et al., 2012;Stepanek and Lohmann, 2012;Yan et al., 2012;Kamae and Ueda, 2012;Zhang and Yan, 2012;Zhang et al., 2012;Chandler et al. 2013;Rosenbloom et al., 2013) to force the British Antarctic Survey ISM (BASISM). Results from PlioMIP present a unique opportunity to sample differences in model predictions of climate and how this impacts on our reconstruction of the GrIS.
Initially a summary of the PlioMIP experimental design will be provided, followed by a description of the offline coupling method adopted for the ISM simulations in this study, which will include details of the climate differences over Greenland as derived from the PlioMIP ensemble. A discussion of the differences between equilibrium-state ice sheet simulations using the climatological forcing from the 15 different climate model experiments in the PlioMIP ensemble will follow and we will conclude with an assessment of the potential causes of any discrepancies and suggestions for future modelling strategies of the mPWP GrIS.
The aims of this paper can be summarised as -to assess the extent to which GrIS reconstructions for the mPWP are dependent upon the climate model used to force the ISM; -to understand the potential reasons for any differences between the simulated GrISs by considering factors  . Regarding the LSM, "preferred" refers to a LSM that has been entirely altered to meet the PlioMIP boundary conditions, whereas "alternate" is where modelling groups have had to use more similar to modern LSM. More comprehensive details of each model, and their implementation of the LSM, can be found in  and the individual references listed in this In order to systematically examine uncertainties in numerical model predictions of the mPWP, PlioMIP  was initiated as a component of PMIP (Palaeoclimate Model Intercomparison Project). PMIP's aim is to provide a means for coordinating palaeoclimate modelling and model-evaluation activities in order to understand the mechanisms of climate change and the role of climate feedbacks under past climate conditions (Braconnot et al., 2012). Previous comparisons of mPWP simulations had been limited to at most three different climate models and had incorporated different approaches to implementing the mPWP boundary conditions (e.g. Haywood et al., 2000Haywood et al., , 2009. PlioMIP established the design for two initial experiments. Experiment 1 used atmosphere-only climate models (AGCMs) and is detailed fully in Haywood et al. (2010). Experiment 2 utilised coupled atmosphere-ocean climate models (AOGCMs) and is described in Haywood et al. (2011). Here the atmospheric and topographic fields from both the AGCMs and the AOGCMs in PlioMIP (Table 1) will be used to force an offline shallow ice approximation ISM (BASISM; see Sect. 2.2).
The boundary conditions applied to all climate models of PlioMIP are described specifically in Haywood et al. (2010Haywood et al. ( , 2011 respectively. In brief, both experiments utilised the US Geological Survey PRISM3 boundary condition data set . PRISM3 is an improved data set in terms of data coverage compared to its predecessor (PRISM2; Dowsett et al., 1999) and includes information on monthly sea surface temperatures (SSTs) and sea-ice distributions, vegetation cover, sea level, ice sheet extent and topography. Vegetation cover is based on the palaeobotanical reconstruction of Salzmann et al. (2008) and topography is derived from the Sohl et al. (2009) palaeogeographic reconstruction. The PRISM3 ice sheets applied in the climate models were derived from offline ISM experiments forced with climatological fields from the Hadley Centre Atmosphereonly climate model ( Fig. 1; HadAM3; Hill, 2009), and represent an ice sheet that is consistent with the rest of the PRISM3 reconstruction. For the AGCMs the SST and sea ice distribution was fixed according to PRISM3, whereas the AOGCMs predicted their own mPWP sea surface conditions.
In all of the PlioMIP experiments, the atmospheric concentration of CO 2 was set to 405 ppmv (Haywood et al., (Hill, 2009;Dowsett et al., 2010). The forcing climatology for this ice sheet reconstruction is a HadAM3 simulation with PRISM2 boundary conditions (as described in Salzmann et al., 2008). 2010, 2011a). This is slightly higher than the previous standard PRISM2 level (400 ppmv), but still falls well within the uncertainty limits of current CO 2 proxy records (e.g. Pagani et al., 2010;Seki et al., 2010;Bartoli et al., 2011). All other trace gases were specified at a pre-industrial concentration and the selected orbital configuration was unchanged from modern .
Each of the PlioMIP models were set up with PRISM3 boundary conditions as described above and then run for a minimum integration length of 50 years for the AGCMs and 500 years for the AOGCMs. Average climatological forcing fields were derived from the final 30 years of the simulation. Each modelling group's standard pre-industrial simulation was used as a control run.
Details of participating groups and climate models can be found in Table 1. Simulations from seven AGCMs and eight AOGCMs were completed and results submitted to PlioMIP. The AGCMs and AOGCMs sample differing levels of complexity and resolution, from higher-resolution IPCC AR5class models to intermediate-resolution models . . mPWP minus pre-industrial mean annual surface air temperature ( • C) over Greenland for the PlioMIP ensemble using atmosphere-only (AGCMs) and coupled atmosphere-ocean climate models (AOGCMs). Temperature plotted on the original climate model resolution.

Climatological forcing over Greenland
Greenland mean annual temperature and precipitation, and summer temperature anomalies between the mPWP and the pre-industrial for each of the PlioMIP AGCMs and AOGCMs are shown in Figs. 2-4. Over Greenland simulated mPWP climates from the AGCMs show an increase in mean annual temperature of between 11.9 and 14.1 • C, whereas the range predicted from the AOGCMs is much greater (5.3 to 12.8 • C; Table 3). For the AGCMs, mPWP mean annual precipitation levels over the Greenland region (Table 3) increase compared to pre-industrial in models. All AOGCMs show an increase in mPWP precipitation of between 0.2 and 0.8 mm day −1 . Simulated mPWP summer temperatures were on average 12.6 • C warmer over Greenland for the AGCMs Figure 3. mPWP minus pre-industrial mean annual precipitation (mm day −1 ) over Greenland for the PlioMIP ensemble using atmosphere-only (AGCMs) and coupled atmosphere-ocean climate models (AOGCMs). Precipitation plotted on the original climate model resolution. and 12.3 • C warmer in the AOGCMs. However, the average for the AOGCMs is lowered due to MRI-AOGCM simulating a warming of 4 • C, whereas all other AOGCMs fall between 11.6 and 15 • C of warming in the mPWP relative to the pre-industrial control simulation.

Ice sheet modelling framework
In this study we used BASISM, which has previously been applied to study mPWP ice sheets (Hill et al., 2007Hill, 2009;Dolan et al., 2011). BASISM is a finite difference, thermomechanical, shallow ice approximation (SIA) ISM, utilising an unconditionally stable, implicit numerical solution of the non-linear simultaneous equations of ice flow. BASISM is similar to other SIA models described by Huybrechts (1990), Ritz et al. (2001) and Rutt et al. (2009) and . mPWP minus pre-industrial mean July surface air temperature ( • C) over Greenland for the PlioMIP ensemble using atmosphere-only (AGCMs) and coupled atmosphere-ocean climate models (AOGCMs). Temperature plotted on the original climate model resolution.
a more detailed discussion of the numerical formulations behind BASISM can be found in Hindmarsh (1993Hindmarsh ( , 1996Hindmarsh ( , 1999Hindmarsh ( , 2001. As well as the internal glaciological dynamics, interactions with the bedrock are simulated with a simple model of elastic rebound, with a rebound timescale of 3000 years (Le Meur and Huybrechts, 1996). The bedrock height for all initial conditions is recalculated using this model, on the assumption of isostatic equilibrium and then the bedrock is allowed to dynamically evolve and adjust during subsequent ice sheet changes.
For this study, BASISM was run on a 20 km × 20 km grid, with 21 vertical layers, in a domain covering the modern grounded GrIS. The ISM is forced using climatological fields of mean annual temperature (Fig. 2) and precipitation (Fig. 3) and warmest mean monthly temperature (July;   Hill (2009). An exponential function is used to convert temperatures into the number of positive degree days (PDD; Reeh, 1991), which shows a high level of correlation between warmest month temperatures and observations of present-day melt (Hill, 2009). Bilinear interpolation was used to downscale the meteorological fields from the original climate model grid onto the higher-resolution ISM grid. Downscaling is problematic in that the coarse horizontal resolution of the climate model is inadequate to resolve the steep topographic slopes around the edges of Greenland (Thompson and Pollard, 1997;Ridley et al., 2005). This is partly addressed by applying a uniform and constant lapse rate correction to resolve for the difference in climate model and ISM topography, both in the initial conditions and as the ice sheet surface evolves during the simulation. The standard lapse rate used within BASISM is −6.0 • C km −1 , which lies within modern observations of lapse rates on Greenland (Steffen and Box, 2001;Hanna et al., 2005). Currently, there is no known similar simple relationship between precipitation and altitude. Precipitation over the GrIS is highly non-linear, with synoptic patterns of atmospheric circulation tending to drive patterns of accumulation (Schuenemann and Cassano, 2009). Combined model simulations and tree-ring isotopes have shown that the dominant patterns of Pliocene North Atlantic atmospheric circulation are likely to have remained similar to today . Although some direct effects of altitude and temperature will occur as the ice sheet evolves, one of the key changes will be feedbacks on atmospheric circulation, which can only be modelled in a coupled ice sheet-climate model (Mayewski et al., 1994). Where downscaling methods do exist (e.g. Ritz et al., 1997), the ratio of precipitation change with temperature change is poorly constrained (Charbit et al., 2002). Therefore, no correction for precipitation has been made within the ice sheet modelling experiments presented here.
The PDD method was employed to convert the climate fields into a melt rate (Reeh, 1991;Braithwaite, 1995), and is well established in coupled atmosphere-ice sheet palaeoclimate modelling studies (e.g. DeConto and Pollard, 2003;Lunt et al., 2008aLunt et al., , b, 2009). The PDD technique assumes that the melting of the ice sheet surface can be fully described by three physical constants (melt rate or PDD factor of ice and snow and the maximum fractional refreezing rate; W max ) and the temperature record. Although many other factors could contribute this method has been shown to have some physical justification (Ohmura, 2001). Standard PDD parameters for ice (∝ i ) and snow (∝ s ) are set to ∝ i = 8 mm day −1 • C −1 and ∝ s = 3 mm day −1 • C −1 respectively, which is within observations of different modern-day climates (Braithwaite, 1995). Further developments of the PDD method have been used in previous studies, but they rely on additional glaciological parameters that may not remain constant in palaeoclimate scenarios, thus it is unclear Table 2. The three glaciological parameters and the values considered in the ice sheet modelling simulations. By varying each glaciological parameter independently, while holding the others constant, there are a total of 48 sensitivity experiments performed for each ice sheet model simulation.
Lapse rate PDD factor snow PDD factor ice how to assign them for the mPWP (Janssens and Huybrechts, 2000;Tarasov and Peltier, 2002). The aforementioned "standard" glaciological parameters (i.e. lapse rate, and the PDD factors of ice and snow) used in BASISM were originally tuned for a HadAM3 experiment (Hill, 2009), so that the best representation of the modern GrIS and East Antarctic Ice Sheet (EAIS) were simulated. However, these parameter values are still poorly constrained, resulting in highly variable ice sheet volumes and extents depending on the exact values prescribed (Ritz et al., 1997;Lunt et al., 2008b;Stone et al., 2010). Stone et al. (2010) demonstrated that the ice sheet extent is predominantly dependent on the PDD factors and the lapse rate, and therefore we have chosen to vary these parameters in order to obtain an additional estimate of uncertainty on our ISM-based GrIS reconstructions.
The typical annual lapse rate used for a variety of studies on Greenland (e.g. Ridley et al., 2005;Huybrechts and de Wolde, 1999;Vizcaíno et al., 2008) ranges from −6.0 to −8.0 • C km −1 and therefore here we will test values within this range ( Table 2). The PDD parameter values for ice and snow vary much more within the literature and previous modelling studies. The standard value for ice used by many modellers is 8 mm day −1 • C −1 (e.g. Huybrechts and de Wolde, 1999;Ritz et al., 1997), although Braithwaite (1995) suggested that the value could be as much as 20 mm day −1 • C −1 . Modelling studies for the Late Pliocene Greenland (e.g. Lunt et al., 2008b) have tested a range of PDD parameters from low PDD factors (∝ i = 8 mm day −1 • C −1 and ∝ s = 3 mm day −1 • C −1 ; the same as BASISM standard) to very high PDD factors (∝ i = 64 mm day −1 • C −1 and ∝ s = 24 mm day −1 • C −1 ) and have shown that the higher end of these ranges do not lead to a good simulation of the modern GrIS.
Although it is possible to use statistical methods such as Latin Hypercube Sampling (LHS) to define random plausible parameter sets within a given range (e.g. Stone et al., 2010), here we simply choose to co-vary parameters. Table 2 shows the parameter values tested here, equating to 48 parameter permutations for each simulation based on the forcing from a specific climate model. In every ISM simulation, absolute temperatures and precipitation values were used to force the ISM, and no correction was made to account for temperature biases in each model's simulation of the pre-industrial (see Lunt et al., 2009). BASISM was run for 50 000 years, which is enough time for the simulated ice sheet to come into geometric and thermal equilibrium with the forcing climate.
Prior to simulating the mPWP GrIS, control cases were run in order to enable an assessment of the modelling framework for the pre-industrial. For the pre-industrial simulations, BA-SISM was initialised from a modern ice configuration. Ini-tially it is useful to determine whether the pre-industrial control climate from each model produces a sensible reconstruction of the present Greenland ice sheet using BASISM with the range of glaciological parameters that are identified in Table 2. In order to analyse the ice sheet geometries from the 48 experiments undertaken for each of the PlioMIP climate models, we have chosen two performance metrics to investigate for each model simulation. Following the methods of Stone et al. (2010), the difference in total ice volume compared to estimated modern volume will be used as an overall diagnostic of how well each simulation reconstructs the observations of ice thickness. The second performance metric will be the root mean square error (RMSE), which is a measure of the spatial fit of the ice sheet thickness reconstruction  Table 2). The coloured circles are superimposed on the grey circles, so when the GrIS volume is similar, the grey circles (or individual colours) will not be visible.
over the Greenland domain. The RMSE describes the magnitude of the differences between two fields (e.g. observed ice thickness and simulated ice thickness). In both cases, lower values describe a better match between the modelled and the observed GrIS. We use the digital elevation model (DEM) of Bamber et al. (2001), interpolated on to the ISM grid (20 km resolution), to calculate observed ice sheet volume and thickness. This technique will also allow the definition of optimal parameter sets (within the envelope of parameter values tested) which gives each forcing climate model the "best" estimate of the present GrIS. These parameter sets were than used with each of the climate forcings from the PlioMIP ensemble.

Pre-industrial control Greenland ice sheets
For the pre-industrial control experiments, BASISM was initialised from the modern GrIS. Figure 5a (AGCMs) and 5b (AOGCMs) summarise the sensitivity of modelled GrIS volume to the three tuneable glaciological parameters (Table 2). For most of the PlioMIP climate model inputs, the choice of parameter values for atmospheric lapse rate and the PDD factors of ice and snow have little impact on the resulting GrIS volume (with the exception of HadAM3 and the fully coupled version of MIROC, where the final volume depends on the choice of the parameter set; Fig. 5b). This is due to the modern ablation zone being constrained to the steep slopes on the periphery of the ice sheet and the constraints applied at the ice sheet grounding line. The parameter set for each PlioMIP model that gives the optimal ice sheet in terms of total ice volume or RMSE of ice thickness for steady-state conditions in comparison to modern observations, is also shown in Fig. 5a and b. Based on the diagnostics chosen here, the optimal parameter sets are never equal to the standard parameter values used within BASISM, although the impact of this on the pre-industrial GrIS is minimal.
For ease of comparison, if we consider using the standard BASISM parameters, all forcing climatologies produce a GrIS which is similar to modern observations. However, the ISM consistently overestimates ice volume by between +3 and +17 %. Comparing the spatial differences between Bamber et al. (2001) and the PlioMIP-based ISM simulations, there are similar biases in elevation (Fig. 6) between the different climate forcings. Over central Greenland, some BASISM simulations produce ice sheets that are too low (∼ 200 to 400 m) in comparison to observations (Fig. 6) although others -notably CAM3.1, COSMOS (AGCM and AOGCM), NorESM (AGCM and AOGCM) -are very close to observations in these regions. Consistent with other ISMs (e.g. Koenig et al, 2014b), all BASISM simulations produce ice sheets that are too high (up to ∼ 800 m) at the ice sheet margins (Fig. 6). These largest deviations from observations occur in the regions of fast ice sheet flow around the ice sheet margins. At these locations there are inherent problems with both the relatively coarse resolution climate model and the ISM when simulating areas of steep topography and complex dynamics. Additionally, as a large proportion (∼ 40 %) of the ice loss in Greenland occurs through iceberg calving (Huybrechts et al., 1991) and such grounding line physics are omitted from this SIA ISM, it is expected that ice loss at the margin would be underestimated (Fig. 6). For smaller simulated ice sheets where ice terminates on land (such as those in the mPWP e.g. PRISM3; Fig. 1), problems associated with ice dynamics such as calving are anticipated to have less of an influence on the reconstruction. The ranking between the simulations depends upon the choice of metric (volumetric or spatial) and thus no one climatological forcing stands out as giving the best representation of the present GrIS. Therefore these metrics will be considered separately in the analysis of mPWP results. RMSE values for each PlioMIP model based on the optimal parameter sets range from 250 to 305 m, and there is no discernible difference in skill at reproducing the modern GrIS between the AGCMs (Fig. 5a) and the AOGCMs (Fig. 5b). Considering both the AGCMs and AOGCMs, the parameter set for each model which gives the smallest RMSE, simulates a difference in volume between the models of 3.01×10 6 km 3 and 3.47×10 6 km 3 . Using the standard parameter set used in BA-SISM, the volume difference for the pre-industrial is similar (3.02 × 10 6 km 3 and 3.44 × 10 6 km 3 ). In summary, none of the simulated ice sheets show any significant biases beyond those inherent when using a SIA ISM (see also Ritz et al. 1997;Saito and Abe-Ouchi, 2005). This provides confidence in the results of the mPWP ISM simulations using the same modelling framework.

mPWP Greenland ice sheets
For the mPWP simulations, BASISM was initialised from the PRISM3 ice configuration  Fig. 1), con-  (Bamber et al., 2001) is shown for comparison. Redfilled circles show the standard parameter set used within BASISM (∝ i = 8 mm day −1 • C −1 and ∝ s = 3 mm day −1 • C −1 , lapse rate = −6 • C km −1 ) and blue-filled circles show the parameter set that gives a volumetric reconstruction closest to observed. Yellow-filled circles show the parameter set that gives the smallest RMSE in terms of simulated ice sheet thickness. Grey circles show the sensitivity of the ice sheet volume to different values of lapse rate and the PDD factors for ice and snow (see Table 2). The coloured circles are superimposed on the grey circles, so when the GrIS volume is similar, the grey circles (or individual colours) will not be visible. sistent with the climate model forcing. Figure 7 shows the simulated GrIS volume for each of the PlioMIP ensemble members using the different glaciological parameters listed in Table 2. In contrast to the pre-industrial ice sheets, mPWP simulations are much more sensitive to the chosen parameter values within the ISM. This is consistent with results presented by Robinson et al. (2011) using a different modelling framework, which show that the modern GrIS is less sensitive to changes in melt parameters than ice sheet reconstructions for the warmer-than-modern Eemian Interglacial (ca. BP 130-115 ka). In all cases, the use of the standard, the volumetrically optimal or the spatially optimal paramewww.clim-past.net/11/403/2015/ ters within BASISM has a significant impact on the resulting mPWP GrIS reconstruction (Fig. 7). Figure 8 shows the surface mass balance (SMB) calculated by BASISM for the PlioMIP climatologies from the initial ISM time step. BASISM simulates a positive SMB over the PRISM3 ice sheet region for the majority of PlioMIP climate forcings and over the southern and western parts of Greenland, net ablation of up to 10 m yr −1 is predicted. In MRI-CGCM2.3 (AOGCM), the cold summer mPWP temperatures ( Fig. 4; Table 3) lead to accumulation over most of the landmass of Greenland (Fig. 8). Conversely, the high summer temperatures exhibited in the NorESM-L models means that the GrIS area experiences only ablation, even over the centre of the PRISM3 GrIS. Figure 9 shows the spatial distribution of the GrIS when BASISM (standard parameter set; red dots in Fig. 7) is forced with atmospheric input fields from each of the PlioMIP models. These results show large differences in both the ice thickness and extent from one simulation to another. Using the AGCMs, ice cover ranges from no ice (NorESM-L) to modern extent (COSMOS, MIROC4m and MRI-CGCM2.3). The absence of ice in the NorESM-L reconstruction is due to the Figure 9. BASISM reconstructions of the mPWP GrIS for individual AGCM and AOGCM forcings. All BASISM simulations were forced with climate model fields (i.e. temperature and precipitation) that were downscaled by a bilinear interpolation method from the original model grid to 20 km × 20 km resolution. GCM specific topography was also used and the ISM simulations were initialised from the PRISM3 ice sheet configuration (Fig. 1). The ice sheet configurations relate to the volumes (red-filled circles) shown in Fig. 7 which use standard glaciological parameters. fact that summer temperatures remain above freezing even when a lapse rate correction has been applied (that accounts for the differences in altitude between the GCM and the ISM grid). Therefore, no ice is able to survive the melt season in this simulation (Fig. 9). The ice sheet reconstructions using CAM3.1 (0.77 × 10 6 km 3 ) and LMDZ5A (1.47 × 10 6 km 3 ) provide ice sheets that are comparable in terms of volume to Table 4. GrIS diagnostics for the PlioMIP simulations, including volume, sea level equivalent (SLE) and ice area using the standard BASISM parameters. Values are given as a difference from the simulated pre-industrial GrIS, when the same GCM pre-industrial forcing climatology is used. For example, negative volume or area means that the GrIS reduces in size compared to the GCM preindustrial control. All simulated volumes (for each parameter set) can be found in the Supplement (Table S1) the PRISM3 GrIS (1.07×10 6 km 3 ), although the distribution of ice is most similar in LMDZ5A (Fig. 9). All AOGCMs produce some ice over Greenland during the mPWP (Fig. 9) and seven of the eight reconstructions show a reduction in volume in comparison to the GCM-specific preindustrial counterpart (Table 4). Ice is distributed in these seven reconstructions as two ice caps, one in the South of Greenland and one spreading out from the mountains of East Greenland. The simulation performed using MRI-CGCM2.3 (AOGCM) produces a GrIS of modern extent with an overall increase in modelled volume relative to the pre-industrial control (+6.3 %; Table 4). This is consistent with the MRI-CGCM2.3 (AOGCM) simulated mPWP temperature over Greenland, which is on average 9 • C warmer than the MRI-CGCM2.3 pre-industrial. Nevertheless, the absolute mPWP temperatures remain much colder than those simulated within the rest of the ensemble, and are actually more akin to the range of pre-industrial temperatures simulated by the other models (Table 3). At the other extreme, NorESM-L produces a GrIS which is reduced in areal extent by 1.41 × 10 6 km 2 (equivalent to a simulated sea level increase of > 7 m). GISS ModelE2-R, HadCM3 and IPSLCM5A produce relatively similar ice sheet configurations over Greenland with the Northern ice cap not extending across to West Greenland. However, the ice sheets reconstructed by CCSM4, COSMOS and MIROC4m either reach or stretch to within ∼ 60 km of the Baffin Bay coastline (Fig. 9). In terms of areal extent and volume, the IPSLCM5A and the GISS ModelE2-R ice sheet reconstructions are the closest to the original PRISM3 GrIS.

Discussion
To date, only a few studies (e.g. Charbit et al., 2007;Quiquet et al., 2012;Yan et al., 2013) have explicitly tested the sensitivity of an ISM to atmospheric input fields, with more studies focusing on parametric uncertainty within ice sheet modelling (e.g. Marshall et al., 2002;Tarasov and Peltier, 2004;Hebeler et al., 2008;Stone et al., 2010). In this study we have tested the climate model dependency of ice sheet reconstructions using output from multiple mPWP climate models. The simulated realisations of the mPWP GrIS reveals significant differences from one simulation to the other, with respect to both the simulated ice volume and ice-covered area, and to the shape and spatial distribution of the ice sheet.

Understanding climate model differences
By comparing the ISM output ( Fig. 9) with GCM-predicted mPWP climate forcing (Figs. 2-4) and the calculated SMB fields (Fig. 8), it is clear that some of the major variations are reflected in the differences in temperature and precipitation distribution amongst the model ensemble. This is in agreement with the study of Charbit et al. (2007), who demonstrated that variability in climate forcing through the last glacial-interglacial cycle induced large differences in simulated Northern Hemisphere ice sheets.
In order to better understand the mechanisms that cause inter-climate model differences in temperature, a more indepth analysis is required that investigates how changes in the energy balance lead to a redistribution of global heat (e.g. Heinemann et al., 2009;Lunt et al., 2012). Hill et al. (2014) have performed such an analysis on the AOGCM (Experiment 2) results from PlioMIP and have shown that the dominant control on annual mean temperature changes in the Arctic regions is related to the clear sky albedo in each model. All AOGCM simulations show that the strongest warming signals come from clear sky albedo (α), although the range in the magnitude of this warming is large (3-12 • C; Hill et al., 2014). Clear sky albedo reflects changes on the Earth surface such as vegetation, snow cover and ice (both terrestrial ice and sea ice). also be attributed to the different albedo properties of snow in each of the climate models (Table 5). For example, some climate models have deep-snow albedo values that are dependent on temperature (e.g. HadCM3, MRI-CGCM2.3 and COSMOS) but the range of maximum and minimum albedo values are not always identical (e.g. the minimum albedo in COSMOS is 0.6 whereas in MRI-CGCM2.3 it is 0.64). Moreover, not all climate models account for factors that influence snow albedo such as the aging of snow or the radiative effects of darkening snow. The differences in the snow albedo schemes implemented in the ensemble may help to explain the differences shown in the mPWP experiments especially over the GrIS region (Fig. 11).
In the ice-free regions of Greenland prescribed in PlioMIP, modelling groups were asked to implement the Salzmann et al. (2008) vegetation reconstruction. Due to the challenging nature of this task different implementation methods were used within the modelling groups. The vegetation distribution was given to the groups in terms of the BIOME4 biome or mega-biome types (Salzmann et al., 2008). However, most modelling groups were unable to implement the data set exactly and instead mapped the plant-functional types onto their own biome scheme. In some cases (e.g. with the GISS ModelE2-R) this meant that distinct biome types within BIOME4 became merged into broader categories within an individual model scheme . It is therefore likely that the albedo properties of the altered vegetation types could be quite different between models, which may be an important factor in the clear sky albedo differences shown in Fig. 11.
The impact of differing albedo schemes over Greenland can be seen clearly in the MRI-CGCM2.3 (AOGCM) reconstruction of the mPWP GrIS (Fig. 9). The high albedo values relative to other models are associated with much colder GISS ModelE2-R

Unknown
Yes -following Loth and Graf (1998) Yes -following Wiscombe and Warren (1980) Yes -following Wiscombe and Warren (1980) Yes -following Warren and Wiscombe Schmidt et al.
HadAM3/ HadCM3 Yes -Uses land surface energy scheme MOSES1 (Cox et al., 1999) and albedo of snow is temperature dependent

MIROC4m
Unknown Yes -following Wiscombe and Warren (1980) Yes -following Wiscombe and Warren (1980) Yes -following Wiscombe and Warren (1980) Yes -following Wiscombe and Warren (1980) Numaguti et al.  Table 3) and lead to the reconstruction of a modern-sized mPWP GrIS (Table 4; Fig. 9). High albedo values in the AOGCM version of MRI-CGCM2.3 are also consistent with results from Hill et al. (2014), which show this model as having the least contribution to mPWP warming from clear sky albedo.
In the AOGCMs, it is also useful to consider differences in predicted SSTs and sea-ice around the Greenland region. Hill et al. (2010) and Koenig et al. (2014) have shown minimal and large responses respectively of the GrIS to fixed SSTs within a climate model. Whilst these studies are not directly comparable (due to the use of different modelling frameworks and different initial conditions), Hill et al. (2010) suggest that the GrIS volume is relatively insensitive to changes in SSTs, with alterations in precipitation being the dominant forcing of the small changes (< 20 % of present GrIS volume). However, Koenig et al. (2014) have demonstrated a greater sensitivity of the GrIS to changes in temperature incurred by fixed SST and sea-ice boundary conditions in the climate model. Also, Ballantyne et al. (2013) have shown that Arctic continental temperatures in general (including those over Greenland) are highly sensitive to the prescription of sea-ice conditions within a model. For the AGCMs the mPWP albedo values over the sea-ice region around the coast of Greenland are very similar, reflecting the prescribed sea-ice conditions in these models (including a sea-ice-free summer; Fig. 11; see also Haywood et al., 2010). Minor albedo differences in the AGCMs are attributed to the varying sea-ice albedo schemes used in the models. Conversely, in the AOGCMs, where the models can freely simulate seaice conditions, there are significant differences in albedo values which reflect the variability in sea-ice predictions in this region. Howell et al. (2015) have performed an in-depth analysis of the differences in Arctic sea-ice predictions within the PlioMIP AOGCM ensemble. It is possible to draw correlations between some models' sea-ice and GrIS reconstructions. For example, the higher summer temperature in July in NorESM-L may be partially attributed to the greatly reduced sea-ice and increased SSTs over the sub-polar North Atlantic, whereas using CCSM4, which retains a substantial sea-ice cover in the Arctic during summer, produces one of the largest predicted GrISs ( Fig. 9; F. W. Howell, personal communication, 2014). Whilst the differing conditions in the surrounding oceans offer some explanation as to the different GrIS predictions from the PlioMIP AOGCM ensemble, they do little to shed light upon the reasons for inter-model differences within the AGCMs. Thus it is difficult to promote sea-ice and SSTs as the sole fundamental control on the extent of the GrIS based on the results presented here.
One further potential contributor to the inter-model differences between ice sheet reconstructions could be the differences in resolution within the PlioMIP ensemble, as GCM resolution (within one model) has been shown to impact on the simulated climate (Roeckner et al., 2006). On one Figure 12. mPWP mean annual temperature ( • C) and precipitation (mm day −1 ), and mean July temperature differences simulated between the AOGCMs and AGCMs over Greenland for comparable models from the PlioMIP ensemble (AOGCM climate minus AGCM climate).
hand, there are multiple scenarios presented here where the GCM horizontal resolution is comparable (i.e. COSMOS and NorESM-L, MIROC4m and MRI-CGCM2.3), but the simulated ice sheet is very different (Fig. 9). However, it is also noticeable that the extent of the prescribed GrIS within each of the PlioMIP models is slightly different due to the model resolution (Table 1). This can be seen most clearly when considering the southward and eastward extent of the regions of accumulation (where the model predicts a positive SMB) in Figure 8. In general such regions of positive SMB track the shape of the prescribed PRISM3 ice sheet in the GCM, and the overall area of accumulation will have an influence on the final GrIS volume.
Where possible, it is also interesting to contrast the results obtained from using a fully coupled version of the model to those obtained using the atmospheric component of the same model (Fig. 12). Six model lineages can be considered in this way: COSMOS (AGCM/AOGCM, Stepanek and Lohmann, 2012), HadAM3/HadCM3 , LMDZ5A/IPSLCM5A (Contoux et al., 2012), MIROC4m (AGCM/AOGCM, Chan et al., 2011), MRI-CGCM2.3 (AGCM/AOGCM, Kamae and Ueda, 2012) and NorESM-L (AGCM/AOGCM, Zhang et al., 2012). As the AOGCM experiments incorporate a dynamic ocean, there is no reason to anticipate that the reconstructed ice sheets will necessarily be comparable when only the atmospheric component of the model is employed. Of the six climate models, four simulate a larger GrIS using the AGCM component than the AOGCM (COSMOS, LMDZ5A/IPSLCM5A, HadAM3/HadCM3 and MIROC4m; Fig. 9). Larger ice sheets are generally associated with the decrease in summer temperatures and increase in precipitation levels in the AGCMs (Fig. 12).
In summary, there are substantial differences in the predicted volumes of the GrIS when forced with multiple climate model predictions (performing a standard experiment), which suggests that the climate model dependency of ISM results is high. However, it is difficult to ascertain why the modelled differences occur between the PlioMIP simulations, although we have shown that the clear sky albedo within each model may be an important factor. In contrast, Koenig et al. (2014b) derive a much lower inter-ISM spread when reconstructing the GrIS during the mPWP, which suggests that relative to climate model dependency, ISM dependency is low. This also gives us confidence that the BASISMbased ice sheet predictions presented here would also hold true if repeated with a different ISM (see also Yan et al., 2014).

Understanding the mPWP Greenland ice sheet
Our results show a high climate model dependency of ISM simulations over Greenland, which implies that the PRISM3 ice sheet configuration (Hill, 2009;Dowsett et al., 2010) is likely dependent on the climate model used within the modelling framework (in this case HadAM3). A better estimation of the GrIS during the mPWP might be derived from considering a "mean" modelled ice sheet, rather than a single reconstruction. A number of studies have shown that a multimodel average often out-performs any individual model if compared to observations (Knutti et al., 2010). This has been demonstrated for mean climate (Gleckler et al., 2008;Reichler and Kim, 2008), but also in regional climate model assessments of the mPWP . A similar approach was taken for defining the Last Glacial Maximum (LGM) ice sheet configuration in the Northern Hemisphere. In the PMIP3/CMIP5 LGM experiments a blended product was obtained by averaging three different ice sheet reconstructions, because of the uncertainties associated with each individual reconstruction (PMIP3, 2010). Here we have calculated an un-weighted multi-model mean (MMM), which is the average of simulations in our multi-model ensemble, treating all models equally.
We have calculated MMMs for both the ice sheet configurations derived using the standard BASISM glaciological parameters and the parameter sets that give the best GrIS reconstruction in terms of modern volume. Figure 13a displays the differences between the calculated MMMs for the AGCM and AOGCM simulations. Present-day observations suggest that if the modern GrIS entirely deglaciated, global sea would rise by around 7.36 m (Bamber et al., 2013). The mPWP GrIS MMM volumes are equivalent to a range in global sea level rise of 2.2 to 4.4 m (Fig. 13a). Due to the difficulties in creating a spatially consistent MMM GrIS, possible ice sheet configurations (taken from the BASISM ensemble of predicted ice sheets) that are approximately equal to the largest and smallest MMM volume are shown in Fig. 13b. It is notable that the smallest MMM ice sheet is very similar to the PRISM3 GrIS boundary condition prescribed in the PlioMIP climate models (Fig. 1), with the exception of the ice cap on Southern Greenland.
There are nevertheless a number of problems with this approach that suggest that caution should be applied when interpreting these results. Firstly, given sea level records and proximal estimates of Greenland ice, it is unlikely that a modern-extent GrIS prevailed during the warmest parts of the mPWP. In the case of the AOGCM "best-fit" parameters, the removal of the large MRI-CGCM2.3 ice sheet reconstruction would make the ensemble spread significantly smaller and also impact upon the calculated MMM (the alternative MMM ice sheet reconstruction in this case would be equivalent to a 5.1 m sea level rise rather than 4.4 m).
Secondly, Contoux et al. (2015) highlight the possibility that the use of the PRISM3 GrIS as a climate model boundary condition for the experiments presented here might bias or precondition the subsequent ISM experiments towards a PRISM3-like GrIS. Contoux et al. (2015) show that when an ice-free Greenland is prescribed in the IPSLCM5A cli-mate model, the subsequent ISM reconstruction is smaller than the PRISM3 GrIS and restricted to the East Greenland Mountains and the southern tip of Greenland. This is supported by the inter-model assessment presented in Koenig et al. (2014b). When prescribing an ice-free Greenland during the mPWP in HadAM3, five SIA ISMs reconstruct a mean ice loss equivalent to a ∼ 7 m global sea level rise. However, when using the same set of boundary conditions to this study (i.e. PRISM3 ice in HadAM3), the contribution of the GrIS to sea level rise ranges between 2.2 and 1.6 m as a MMM (see Koenig et al. 2014b for further details). This highlights the impact of the choice of the initial ice configuration in the climate model. However, without a fully coupled ice-sheetclimate model, this is a difficult problem to overcome. Given the modelling framework adopted here, and the likely presence of ice on Greenland, it is essential to prescribe an ice sheet in the climate model, which requires a number of a priori assumptions regarding ice distribution. Not only does this have implications for our understanding of the GrIS during warm interglacials of the Late Pliocene, but an incorrect representation of the ice sheets in general may have a negative impact when assessing global climate model simulations against proxy data from the mPWP/Late Pliocene (e.g. Dowsett et al., 2012Dowsett et al., , 2013Salzmann et al., 2013).
A final caveat to this research is derived from the uncertainty as to whether a good simulation of the modern GrIS (when compared to observations) necessarily implies a realistic representation of the mPWP GrIS.  found that when simulating the Eemian GrIS (where significantly more constraints are available than for the Late Pliocene), the ISM simulation that gave the most realistic modern ice sheet gave an entirely unrealistic ice sheet for the Eemian when compared with proxy data. This highlights the need for further palaeodata constraints regarding the extent and thickness (where possible) of the Late Pliocene GrIS in order to thoroughly assess the results presented here.

Climate model boundary conditions for PlioMIP phase 2
The final aim of this study and the wider PLISMIP project  is to inform decisions regarding the ice sheet boundary conditions to be prescribed in the second phase of PlioMIP . The high climate model dependency of the GrIS shown here now brings into question the suitability of the PRISM3 GrIS in PlioMIP Phase 1, as this was the result of a one climate model/one ISM modelling framework. However, the broad range in the MMM ensemble presented here and the problems associated with a priori assumptions necessary to undertake this modelling framework suggest that the simple use of a MMM GrIS is inappropriate.
It is therefore likely that future GrIS reconstructions will be based on a combination of climate/ice sheet modelling results (e.g. Koenig et al., 2014b;Contoux et al., 2015 and those presented here) and data-based constraints. Evidence of vegetation suggesting ice-free conditions can be found in North Greenland (Funder et al., 2001), at Ile de France (Bennike et al., 2002), on Ellesmere Island and the Canadian Archipelago (De Vernal and Mudie, 1989;Thompson and Flemming, 1996;Ballantyne et al., 2006;Csank et al., 2011), and these offer limited constraints on a mPWP GrIS reconstruction. More recently Bierman et al. (2014) have shown a preservation of a preglacial landscape under the centre of the GrIS at the site of the Greenland Ice Sheet Project 2 core. They suggest that the soils which formed at the base of the core (at the onset of Northern Hemisphere Glaciation around 2.7 Ma) could have been subaerially exposed for between 200 000 and 1 million years, which has been suggested to imply that this region was potentially ice free in the warm Late Pliocene. Additionally, a recent reassessment of pollen derived from ODP Hole 646B off Southwest Greenland (de Vernal and Mudie, 1989) confirms that southern Greenland would have been vegetated (boreal and cool-temperate conditions) during warm parts of the Pliocene (A. de Vernal, personal communication, 2014).
Combined, the proxy-based evidence and the modelling work done to date would suggest that a smaller ice cap (in relation to PRISM3), centred on the eastern Greenland Mountains, is the best available estimation of a warm interglacial mPWP GrIS configuration. Clearly however, there is a critical need for further data pertaining to ice extent (e.g. Bierman et al., 2014) or potentially the Greenland climate (such as vegetation records) in order to more accurately constrain this reconstruction.

Conclusions
The Pliocene Ice Sheet Modelling Intercomparison Project  was initiated in order to ascertain the degree to which ice predictions over Greenland are influenced by the choice of ISM and climate model. Whilst Koenig et al. (2014b) have shown that ISMs are generally relatively consistent in their predictions when forced with the same climatology, here we show that the choice of climate model significantly affects the predicted GrIS. Ice sheet reconstructions using forcing from the PlioMIP AGCMs and AOGCMs range from a larger-than-modern GrIS to an ice-free Greenland. Such a result demonstrates the difficulty in using only one climate model to draw conclusions regarding ice sheet stability in the warm Late Pliocene and highlights the need for an alternative ice sheet reconstruction going forward with PlioMIP Phase 2.
The Supplement related to this article is available online at doi:10.5194/cp-11-403-2015-supplement.