21st century ocean forcing of the Greenland Ice Sheet for modeling of sea level contribution

Changes in the ocean are expected to be an important determinant of the Greenland Ice Sheet’s future sea level contribution. Yet representing these changes in continental-scale ice sheet models remains challenging due to the small scale of the key physics, and limitations in processing understanding. Here we present the ocean forcing strategy for Greenland Ice Sheet models taking part in the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6), the primary community effort to provide 21st century sea level projections for the Intergovernmental Panel on Climate Change 6th Assessment Report. Begin5 ning from global atmosphere-ocean general circulation models, we describe two complementary approaches to provide ocean boundary conditions for Greenland Ice Sheet models, termed the ‘retreat’ and ‘submarine melt’ implementations. The retreat implementation parameterizes glacier retreat as a function of projected submarine melting, is designed to be implementable by all ice sheet models, and results in retreat of around 1 and 15 km by 2100 in RCP2.6 and 8.5 scenarios respectively. The submarine melt implementation provides estimated submarine melting only, leaving the ice sheet model to solve for the resulting 10 calving and glacier retreat, and suggests submarine melt rates will change little under RCP2.6 but will approximately triple by 2100 under RCP8.5. Both implementations have necessarily made use of simplifying assumptions and poorly-constrained parameterisations and as such, further research on submarine melting, calving and fjord-shelf exchange should remain a priority. Nevertheless, the presented framework will allow an ensemble of Greenland Ice Sheet models to be systematically and consistently forced by the ocean for the first time, and should therefore result in a significant improvement in projections of the 15 Greenland ice sheet’s contribution to future sea level change.


Introduction
The rapid response of the Greenland Ice Sheet to climate warming in the past few decades, together with expectations of future climate change, have raised concern that Greenland will contribute significantly to sea level change over the coming decades and centuries (Shepherd et al., 2012;Church et al., 2013;Nick et al., 2013). Greenland contributed ⇠13.7 mm to global mean sea level between 1972 and 2018, with surface mass balance comprising 35-60% of this ice mass loss (van den Mouginot et al., 2019). The remainder derives from discharge from tidewater outlet glaciers, most of which have retreated, accelerated and thinned in recent decades (Rignot and Kanagaratnam, 2006;Khan et al., 2014;Murray et al., 2015). These tidewater glaciers are understood to have responded to climate forcing occurring at their calving fronts, where 5 the ice sheet meets the ocean (Nick et al., 2009;Luckman et al., 2015;Wood et al., 2018). Thus, processes at the ice-ocean boundary and their representations in ice sheet models are a critical component of accurate future sea level projections.
The Ice Sheet Model Intercomparison Project (ISMIP6; Nowicki et al., 2016) is our ::: the community-leading effort projecting future sea level contribution from the Greenland and Antarctic Ice Sheets for the coming 6th Assessment Report of the Intergovernmental Panel on Climate Change (IPCC AR6). ISMIP6 follows a history of similar initiatives, such as SeaRISE 10 (Bindschadler et al., 2013;Nowicki et al., 2013a, b) and ice2sea (Gillet-Chaulet et al., 2012;Goelzer et al., 2013), aimed at bringing together a number of ice sheet models and scientists across disciplines to improve projections of ice sheet mass loss.
Ocean forcing of the Greenland Ice Sheet occurs at around 300 approximately vertical glacier calving fronts around Greenland and at several larger ice shelves and floating ice tongues located in the far north. Ocean forcing is here broadly defined as melting of the ice-ocean boundary (hereafter called submarine melting) and the impact of this melting on calving and glacier retreat. The design of boundary conditions that represent ocean forcing must take into account three sets of processes. First, 25 the transport of ocean heat from the far-field ocean to calving fronts across continental shelves and up long and narrow fjords.
Second, the near-ice circulation which drives heat transfer through the ice-ocean boundary. Third, the impact of submarine melting on iceberg calving and glacier retreat.
Understanding of these key processes has advanced through both observations and models. Considering observations, warm Atlantic-origin water is found on the continental shelf around Greenland either due to transport from the deep ocean to the 30 shelf, often in deep troughs, or due to advection along the shelf by coastal currents (Sutherland et al., 2013;Rykova et al., 2015;Schaffer et al., 2017). The same waters are found adjacent to calving fronts (Straneo et al., 2012) and may enter the fjords by numerous processes including a glacier-driven estuarine-type circulation that may be prevalent in summer (Motyka et al., 2003;Gladish et al., 2015), fjord-shelf exchange driven by winds both inside and outside of fjords (Jackson et al., 2014;Spall et al., 2017) and exchange due to variability in shelf water properties (Mortensen et al., 2014;Carroll et al., 2018). Once warm water reaches the calving front, the transfer of heat across the ice-ocean boundary layer is promoted by the near-ice circulation (Holland and Jenkins, 1999). During summer, the release of ice sheet meltwater into fjords from beneath tidewater glaciers drives localised but vigorous upwelling plumes which are thought to drive rapid submarine melting (Mankoff et al., 2016;Sutherland et al., 2019). These plumes may also fuel a fjord-wide circulation which enhances submarine melting over the full calving front . Submarine melting can shape the calving front, creating regions 5 of undercut and overcut ice (Fried et al., 2019), which may in turn enhance iceberg calving and drive glacier retreat How et al., 2019). Greenland's shelf and fjords, however, remain sparsely observed, especially in winter, and we have very few observations of submarine melt rate. Similarly, significant uncertainty surrounds the dynamic impact of submarine melting on calving due to the difficulty of making the necessary measurements close to dangerous calving fronts.
Many of these processes can be captured by models at the individual fjord or glacier scale. Cowton et al. (2016) and Fraser 10 et al. (2018) have modeled fjord-shelf exchange at Kangerdlugssuaq Fjord in south-east Greenland, while Carroll et al. (2017) modeled fjord water renewal driven by subglacial discharge in an idealised domain. Plumes and the near-ice circulation they generate have been captured by models focused on the part of the fjord within a few kilometers of the calving front Slater et al., 2018). The impact of submarine melting on calving has been studied at high resolution in both idealised and realistic settings Ma and Bassis, 2019;Todd et al., 2019). Yet, the model resolution in these studies, 15 at ⇠500 m for the fjord simulations, ⇠10 m for the plume simulations and ⇠50 m for the calving simulations, is far smaller than the ⇠50 km resolution of AOGCMs (e.g. Watanabe et al., 2010) or ⇠2 km resolution of Greenland Ice Sheet models (Goelzer et al., 2018). Even regional ocean models (e.g. Gillard et al., 2016) do not yet represent fjords and fjord processes.
Thus, climate and ice sheet models do not have sufficient resolution to capture the processes that modulate the effect of the ocean on the Greenland Ice Sheet. 20 At present, therefore, projecting the sea level contribution of the Greenland Ice Sheet requires that we parameterise iceocean processes, but well-validated parameterisations are not readily available. While progress has been made in observing and modeling fjord circulation and fjord-shelf exchange, we as yet lack simple parameterisations or box models that could represent these processes in an efficient fashion (i.e. without resorting to computationally expensive hydrodynamic models).
Conversely, parameterisations exist for the submarine melting induced by plumes (Rignot et al., 2016;Slater et al., 2016) but 25 we still have few observations with which to validate these parameterisations . Lastly, the search for a universal calving law has a long history (e.g. Benn et al., 2007) but as for submarine melting no calving parameterisation has undergone sufficient validation for confident use.
Given the described process uncertainty, the small scale of key processes and the current lack of parameterisations for these processes, projecting ocean-induced ice mass loss from the Greenland Ice Sheet is very challenging. To date, attempts to 30 project future ice discharge from tidewater glaciers have often relied on extrapolation from a few glaciers to the whole ice sheet (Goelzer et al., 2013;Nick et al., 2013;Peano et al., 2017;Beckmann et al., 2019;Morlighem et al., 2019) or have employed ad-hoc methods to mimic the impact of ocean forcing that are not easily relatable to climate warming scenarios (Price et al., 2011;Bindschadler et al., 2013;Fürst et al., 2015). In a single ice sheet model, a significant advance was recently made by Aschwanden et al. (2019), who ran full ice sheet projections that resolve tidewater glaciers and were forced by estimated 35 submarine melt rates but many of the ice sheet models taking part in ISMIP6 do not currently have the resolution of :: or : technical capability for this approach (Goelzer et al., 2018).
In spite of the described difficulties, we present a strategy for simulating the impact of the ocean on the ice sheet that will enable a suite of Greenland Ice Sheet models of diverse capability ::::::::: capabilities : to be systematically forced by future warming scenarios ::::::::::::::::: (Goelzer et al., 2020). We do not aim to solve the problems of process understanding, scale and parameterisation but 5 rather to offer a pragmatic approach based on the current state of knowledge. This approach draws on existing parameterisations for tidewater glacier retreat  and submarine melting (Rignot et al., 2016). The paper proceeds as follows.
An overview of the two-tier strategy for ocean forcing is given, before the subglacial runoff and ocean thermal forcing datasets are described. These time series are combined into projections of glacier retreat and submarine melting. We finally discuss the projected ocean forcing, its temporal evolution and spatial and inter-model variability.

Overview
We develop two possible implementations for ocean forcing of Greenland Ice Sheet models, referred to as the retreat implementation and the submarine melt implementation (Fig. 1). The retreat implementation is designed to be implementable by all of the ice sheet models taking part in ISMIP6 regardless of resolution, model physics or spin-up procedure. In this imple- 15 mentation, retreat of the ice-ocean boundary is estimated as a linear function of parameterised submarine melting  and is imposed on an ice sheet model through a time-variable ice mask, an approach first suggested by Cowton et al. (2018). The submarine melt implementation instead provides ice sheet modeling groups with fields of subglacial runoff and ocean properties together with a suggested parameterisation for estimating submarine melt from these quantities. Since glacier retreat is given by a competition between frontal ice velocity, calving and submarine melting, the retreat implementation heav-20 ily parameterises ocean forcing by implicitly assuming that all quantities are proportional to submarine melt rate . The submarine melt implementation allows ice sheet models to resolve the competition between velocity, calving and melting, perhaps by implementing a calving law that depends on submarine melt rate.
Both implementations require a parameterisation for submarine melting. Theoretical considerations suggest that melt rates are controlled primarily by local ocean velocity and ocean thermal forcing, the latter defined as the difference between the 25 in-situ temperature and in-situ freezing point (Gade, 1979;Holland and Jenkins, 1999). Near-ice ocean velocities are thought to be highest inside vigorous plumes resulting from the emergence of buoyant subglacial runoff from the grounding line of the glacier (Mankoff et al., 2016). Submarine melt rate parameterisations (Jenkins, 2011;Xu et al., 2013;Slater et al., 2016), therefore, typically include the basic ingredients of subglacial runoff (Q) and ocean thermal forcing (TF). In the retreat implementation we follow Slater et al. (2019) in assuming that submarine melting is proportional to Q 0.4 TF and retreat ( L, 30 in km) is proportional to submarine melting, so that retreat may be estimated as where Q is the mean summer (June-July-August) subglacial runoff (in m 3 s 1 ) and TF is the ocean thermal forcing (in C).
where h is grounding line depth (in m), TF is the ocean thermal forcing (in C) and q is the annual mean subglacial runoff normalised by calving front area (in m d 1 ). We acknowledge the inconsistency of using summer runoff for the retreat implementation and annual runoff for the submarine melt implementation but we emphasise that this makes no practical difference since annual and summer runoff are very closely related, even in the future projections when the melt season becomes longer 10 . The parameterisation for submarine melting is slightly more complex than that for retreat but is functionally very similar.
The chosen parameterisations require the two basic inputs of future subglacial runoff and ocean thermal forcing, which are estimated from CMIP AOGCMs. While it is hoped that some of the new generation of climate models (CMIP6) will be used in ISMIP6, very few CMIP6 simulations were available at the time of writing and given the time constraints of the ISMIP6 15 project it was decided to focus largely on CMIP5, for which the full ensemble is already available. We consider 6 CMIP5 AOGCMs ( Table 1) that represent a subset of the full CMIP5 ensemble but emphasise that the process would be identical for CMIP6 inputs. The 6 CMIP5 AOGCMs have been chosen by selecting AOGCMs with minimal biases in the present day and with the aim of sampling the diversity of projected climate change, as described in Barthel et al. (2019). The focus is on the RCP8.5 scenario, a high greenhouse gas emissions pathway in which radiative forcing reaches 8.5 W m 2 in 2100 (Riahi et al., 20 2011;Nowicki et al., 2016). We also consider a single RCP2.6 simulation (radiative forcing of 2.6 W m 2 in 2100). Each of the CMIP5 AOGCM simulations covers the period 1850-2100, with 1850-2005 considered the historical spin-up period and the emissions forcing applied from 2006-2100 (Taylor et al., 2012). Ice sheet model ocean forcing is delivered for the time period from 1950-2100. The remainder of this methods section describes the calculation of subglacial runoff and ocean thermal forcing from AOGCM output and the combination of these datasets into ice sheet model ocean forcing in the retreat 25 and submarine melt implementations (Fig. 1).

30
(2013)). The most recent version of the model, MAR 3.9.6, is run at 15 km resolution with surface mass balance components (including runoff) statistically downscaled afterwards to 1 km (Franco et al., 2012) to better account for sub-grid topography ( Fig. 2a). Each simulation is forced at its boundaries by 6-hourly output from a CMIP5 AOGCM (Table 1) over the period 1950-2100.

Hydrological drainage basins
Both the retreat and submarine melt implementations use an estimate of subglacial runoff per tidewater glacier, which requires a hydrological drainage basin for each glacier (Fig. 1). These basins are delineated based on the hydrological potential (Shreve,5 1972): where ⇢ w = 1000 kg m 3 and ⇢ i = 910 kg m 3 are the densities of freshwater and ice respectively and g = 9.81 m 2 /s is the gravitational acceleration. Bed topography, b (m) and ice thickness, h (m), come from BedMachinev3 (Morlighem et al., 2017).
The variable f represents the ratio of subglacial water pressure to ice overburden pressure. Based on limited borehole pressure 10 records we set f = 1 (Meierbachtol et al., 2013;Andrews et al., 2014;Doyle et al., 2018) but acknowledge that different values of f can alter drainage pathways (e.g. Chu et al., 2016;Moyer et al., 2019). By performing flow routing on (Schwanghart and Scherler, 2014), we identify the area of the ice sheet that drains subglacial water to a given tidewater glacier calving front, defining hydrological drainage basins for each tidewater glacier around the ice sheet (Fig. 2b). For simplicity the hydrological drainage basins are assumed to be constant in time. Given the high density of moulins observed around the margins of the ice 15 sheet during summer (e.g. Yang and Smith, 2016), we assume that all surface meltwater drains to the ice sheet bed close to where it melts. The subglacial runoff for each glacier can then be estimated by summing the surface runoff from MAR over the hydrological drainage basin for each glacier (Fig. 2b). Studies that have assessed subglacial runoff from fjord observations find agreement between their oceanographic estimates and the method described here (Jackson and Straneo, 2016;Mankoff et al., 2016;Jackson et al., 2017).

Present-day bias correction
The :::: Many : CMIP5 AOGCMs may deviate considerably from the observed present-day climate in both the atmosphere and ocean. For example, Menary et al. (2015) show CMIP5 ocean temperature biases can exceed 2 C in the Labrador Sea. If the AOGCM-simulated atmosphere is substantially colder than observations, runoff will be underestimated in MAR when forced by the AOGCM in question . Since in the ISMIP6 exercise we wish to sample uncertainty in future 25 projections rather than the representation of the present day, we perform a bias correction of the projected subglacial runoff at each glacier to ensure it agrees with our best estimate of present-day runoff (Fig. 1). This bias correction furthermore ensures a continuous transition from present to future forcing, which is desirable as the ice sheet models have been initialised to the present-day forcing (Goelzer et al., 2018).
Present day is defined as the time period 1995-2014. For our best estimate of runoff in the present day we use a 5.5 km 30 resolution regional climate simulation using RACMO2.3p2, forced at its boundaries by ERA-Interim atmospheric reanalysis . We ensure that the projected runoff (Q P ROJ ) agrees with the RACMO runoff (Q RACM O ) in the present day by bias-correcting the projected runoff for each glacier (j) as follows: where the 1995-2014 in parentheses indicates the mean value between 1995 and 2014. We assume that the bias remains constant in time. An example of this procedure for Helheim Glacier in SE Greenland under MIROC5 in an RCP8.5 scenario is shown in Fig. 2c. In this case the JJA runoff estimated from MAR forced by MIROC5 is decreased by 55 m 3 s 1 to 316 m 3 s 1 to bring it 5 into agreement with the temporally averaged RACMO2.3p2 output over the period 1995-2014. Note that we do not expect the interannual runoff variability in MAR forced by MIROC5 to agree with RACMO2.3p2 forced by ERA-Interim (Fig. 2, inset) because MIROC5 is a free-running climate model whereas ERA-Interim is an atmospheric reanalysis.
Over all glaciers and all CMIP5 AOGCMs considered (Table 1), the mean bias correction is +2 m 3 s 1 with a standard deviation of 56 m 3 s 1 and a minimum and maximum correction of -527 and +519 m 3 s 1 , respectively (Fig. S1). As a fraction 10 of the present-day runoff the mean bias correction is +0.13 with a standard deviation of 0.47. Bias corrections for the largest glacier by ice flux in each sector and for all models are shown in Fig. S1.
It would be better to use MAR forced by ERA-Interim for our best estimate of present day, because it is MAR that is used for the forward projections. If we define the interannual runoff variability as the standard deviation of the detrended projections, we find a mean interannual variability across all glaciers and AOGCMs of 74 m 3 s 1 . Given that the bias correction (i.e. the 15 difference between RACMO and MAR in the present day) is typically smaller than the interannual variability of the projections, the use of RACMO for the present day does not cause any inconsistency in methodology :::::: practice.

Defining ocean thermal forcing
Due to a lack of parameterisations that can capture fjord-shelf exchange and fjord circulation without resorting to full hydro-20 dynamic models, we take a simplified approach to estimating ocean thermal forcing in which the forcing experienced by the glacier is directly related to far-field ocean properties. As such, we are hard-wiring tidewater glaciers to respond to large-scale ocean changes at the expense of most of the local details that we cannot currently account for. Specifically, we spatially average ocean properties over predefined ocean regions and use these properties to force all tidewater glaciers in the same region ( Fig. 1). For the retreat implementation, the far-field ocean properties are furthermore depth-averaged (section 2.4.1) while for 25 the submarine melt implementation, the far-field ocean properties are extrapolated into fjords taking account of bathymetry (section 2.5.1).

Choice of ice-ocean sectors and spatial averaging
The ice sheet and surrounding ocean were divided into 7 ice-ocean sectors ( Fig. 3a) over which ocean properties were spatially averaged ( Fig. 1). Each sector is hereafter referred to by its acronym (Fig. 3a), where SW is south-west Greenland, CW is 30 central-west Greenland, NW is north-west Greenland, NO is northern Greenland and similarly for the eastern side of the ice sheet. The sectors, identical to those considered in Slater et al. (2019), were chosen as regions with similar ocean properties largely defined by ocean bathymetry (e.g. Denmark, Fram, Nares Straits) and consistent with the boundaries of commonly used ice sheet drainage basins (e.g. Mouginot et al., 2019) once extended into the ice sheet (see Slater et al. (2019) for a more in depth description). The small region in CE Greenland is a transition zone between the warm Atlantic waters in the Irminger basin to the south and cool Arctic waters in the Nordic Seas to the north and, as such, was split from the SE and NE 5 Greenland sectors. Each ice-ocean sector extends to the centre of the offshore ocean basin or strait, except for in the Arctic Ocean, Greenland Sea and Labrador Sea where the ocean basin is very large and the sector boundary is located approximately 150 km beyond the shelf break (Fig. 3). With these choices we sample the water masses that interact with the ice sheet but not those that are recirculating (e.g. in western Baffin Bay). Extending the sectors beyond the shelf break also allows us to access many more ocean observations (Fig. S2), which provides greater confidence in the calibration of the retreat parameterisation To obtain sector ocean properties, monthly CMIP5 AOGCM outputs of modeled ocean potential temperature (T ) and practi- 15 cal salinity (S) are first temporally averaged to annual means ( Fig. 3a). Temperature and salinity are then linearly interpolated onto a regular grid with 50 km spatial and 50 m depth resolution (Fig. 3b). Sector ocean properties are finally obtained by taking a simple spatial average over all regular grid points inside a given sector to give a single temperature and salinity profile for each ice-ocean sector for each year (e.g. Fig. 3c).

Present-day bias correction 20
As for the subglacial runoff, we bias-correct the ocean properties to ensure consistency with observations in the present day We obtain annual profiles per ice-ocean sector from EN4 in the same fashion as for the CMIP5 AOGCM projected profiles.
While for subglacial runoff we bias-corrected a single value, here we must bias-correct a whole temperature or salinity profile.
Rather than applying a different bias correction at each depth level, we apply a single bias correction to the whole profile based 30 on the observed bias in the 200-500 m depth range. Specifically, we bias-correct ocean temperature ( Fig. 3c) as follows Salinity is bias-corrected in exactly the same fashion. Since the vertical structure of the ocean can vary in time in the CMIP5 5 AOGCMs, we felt a depth-varying bias correction could lead to unphysical profiles and that a single-valued correction, centered over the depth range most relevant to tidewater glacier grounding lines (Morlighem et al., 2017), was preferable. As for the runoff, the bias correction is assumed constant in time. The magnitude of these corrections can be significant. For example, in MIROC5 RCP8.5 the temperature bias correction for SE Greenland is 1.4 C (Fig. 3c). Over all sectors and CMIP5 AOGCMs considered, the mean temperature bias correction is +0.1 C with a standard deviation of 1.5 C and a minimum and maximum 10 correction of -3.1 and +3.2 C, respectively (Fig. S3).

Calculation of ocean thermal forcing
To calculate the thermal forcing that enters the retreat parameterisation in Eq.
(1), profiles of ocean temperature and salinity (e.g. Fig. 3c) are first converted to profiles of thermal forcing (Fig. 1). The thermal forcing (TF) is for the retreat parameterisation 15 defined as the elevation of the potential ocean temperature T above its local freezing point T f where in the second equality we have employed a linearised expression for the local freezing point in terms of the practical salinity S and depth z and the constants take values 1 = 5.73 ⇥ 10 2 C psu 1 , 2 = 8.32 ⇥ 10 2 C and 3 = 7.61 ⇥ 10 4 C m 1 (Jenkins, 2011). As before, i indexes the ice-ocean sector. 20 In-keeping :: In :::::: keeping : with the simple philosophy of the retreat parameterisation, the profiles of thermal forcing TF i (z, t) are finally depth-averaged between 200 and 500 m depth, this being the depth range most relevant to tidewater glacier grounding lines in Greenland (Morlighem et al., 2017). The final thermal forcing entering Eq. (1) in the retreat implementation is a single value per ice-ocean sector per year, for each CMIP5 model considered (Table 1).

Glacier-by-glacier projection of retreat
is the ice-ocean sector i from 1 to 7 in which the glacier j is situated (Fig. 4a). Since this time series has high interannual variability and since for ISMIP6 we are most interested in the multi-decadal sea level contribution, the time series is smoothed using a 30 20-year centered moving average (Fig. 4a). Lastly, in the CMIP6 and ISMIP6 frameworks (Nowicki et al., 2016;Eyring et al., 9 2016) the projections begin in 2015 and we project retreat relative to 2014. Thus for each glacier j, projected retreat L j (t) is given by where both terms on the right-hand side refer to the smoothed time series. We generate 10 4 possible future retreat trajectories for each glacier (Fig. 4b) by sampling 10 4 values of  from its distribution obtained from observations .

Averaging retreat per ice-ocean sector
Due to limitations of the retreat parameterisation -principally its lack of ability to capture individual glacier effects related to bed topography -it is most appropriate to apply retreat averaged over a population of glaciers rather than on an individual glacier basis . From the ice sheet model perspective, this is also preferable because the state of the ice sheet may differ significantly from the observed ice sheet (Goelzer et al., 2018). Thus, identifying individual glaciers in a given ice 10 sheet model is not trivial so that applying retreat to individual glaciers is also difficult. An obvious solution is to impose a given retreat over a predefined geographical region (or ice-ocean sector), which means averaging retreat over a population of glaciers.
A potential issue is that under the retreat parameterisation (Eq. (1)), glaciers with large hydrological catchments (typically glaciers such as Jakobshavn Isbrae or Helheim) undergo large changes in subglacial runoff and have large projected retreat 15 relative to smaller glaciers. This is considered an important feature of the retreat parameterisation . Each ice-ocean sector (Fig. 3a) typically has a small number of large glaciers and a large number of small glaciers, such that taking a simple mean of the projected retreat over the glaciers in a sector will result in a trajectory that is much closer to that of the small glaciers than the large glaciers. This is problematic because the primary objective of ISMIP6 is sea level contribution and for Greenland this is dominated by the largest glaciers (Enderlin et al., 2014). To address this problem, we take an ice 20 flux-weighted mean over glaciers in a sector (Fig. 1). Specifically, we define the retreat for each sector i as where f j is the 2000-2010 mean observed ice flux (Enderlin et al., 2014;King et al., 2018) and the sum runs over all glaciers j in ice-ocean sector i. This ensures that the largest glaciers are treated as the most important when generating a retreat projection per sector. Since we have 10 4 retreat trajectories for each glacier (Fig. 4b), this procedure produces an ensemble of 10 4 ice 25 flux-weighted retreat trajectories for each ice-ocean sector. As expected, the median retreat of this ice flux-weighted ensemble is larger than the median retreat that would have been obtained by taking a simple mean over glaciers in a sector (Fig. 4c).

Low, medium and high scenarios
Given the large uncertainty associated with tidewater glacier response to climate forcing and the need to quantify uncertainties on future sea level contributions, it is desirable to provide a range of projected retreat that brackets the uncertainty associated 30 with the retreat implementation. For each CMIP5 AOGCM we identify a low, medium and high retreat scenario (Fig. 1). From the ensemble of 10 4 ice flux-weighted retreat trajectories for each ice-ocean sector, we define the medium retreat scenario as the trajectory with the median retreat at 2100 and the low and high retreat scenarios as the trajectories with the 25th and 75th percentile retreats at 2100 (Fig. 4c).

Extrapolation of ocean properties into fjords
In the submarine melt implementation, we account for the effects of fjord bathymetry and grounding line depth on the thermal forcing experienced by the glacier (Fig. 1). This is achieved by extrapolating the ocean property profiles (e.g. Fig. 3c) into fjords and below the present-day ice sheet by taking into account ocean bathymetry and subglacial topography in the same manner as This procedure is repeated for all fjords around the ice sheet, including below the present-day ice sheet, so that ocean conditions at calving fronts will be available to ice sheet models after calving fronts have retreated.

25
In line with the more complex nature of the submarine melt implementation relative to the retreat implementation we use full, non-linear TEOS-10 routines (McDougall and Barker, 2011) to convert ocean property profiles to ocean thermal forcing profiles ( Figs. 1 and 5d). Specifically, the CMIP5 quantities of depth, practical salinity and potential temperature are converted to pressure, absolute salinity and in-situ temperature using the 'gsw_p_from_z', 'gsw_SA_from_SP' and 'gsw_t_from_pt0' routines, respectively. A full three-dimensional, time-varying thermal forcing field TF(x, y, z, t) is obtained as 30 TF(x, y, z, t) = T (x, y, z, t) T f (x, y, z, t) where T is the in-situ temperature and T f is the in-situ freezing point that depends on pressure and absolute salinity as defined by the 'gsw_t_freezing' routine. Lastly, we collapse the three-dimensional thermal forcing field to two-dimensions by considering only the value at the ocean bottom, so that the final thermal forcing field (TF) is defined at annual resolution on a 1 km x-y grid covering Greenland (Fig. 6a). The motivation for using the ocean bottom value is that this is the thermal forcing experienced by the grounding line of a glacier if its calving front was located in the grid cell in question. Furthermore, plumes 5 upwell deep waters towards the fjord surface so that the temperature profile within the plume is well approximated by the value at the ocean bottom (Mankoff et al., 2016). Lastly, ::: We :::: note ::: that : the submarine melt parameterisation is non-linear in TF (Eq. 2) so that annual mean melt is not equal to melt calculated from annual mean TF. The difference is, however, less than 1% and it is, therefore, justified to use annual mean TF.

Assignation of runoff to drainage basins 10
The treatment of subglacial runoff is initially the same as for the retreat parameterisation. Once the time series of bias-corrected subglacial runoff has been obtained for each marine-terminating glacier (section 2.2), this runoff is distributed onto a 1 km x-y grid by assigning the total runoff for each hydrological basin (Fig. 2b) to every grid point lying inside the basin (Figs. 1 and   6b). In this way, as a calving front retreats over the x-y grid, the calving front submarine melt rate may be obtained by sampling the ocean thermal forcing and subglacial runoff from the grid point at which the calving front is currently located. We assume 15 that the hydrological drainage basins remain fixed in time at their present-day extent. Extending the runoff field beyond the present-day ice sheet is desirable to allow for potential calving front advance in the simulations, or to accommodate models whose initial ice extent is larger than observations. We choose to extrapolate subglacial runoff values beyond the present-day ice sheet by three 1 km grid cells using an iterative buffering approach. First, we sort the drainage basins by area from largest to smallest. For each iteration, we buffer runoff values by one 1 km grid cell around each basin, starting with the largest basin 20 and ending with the smallest basin. We fill only empty grid cells such that if a grid cell has already been populated by a runoff value from a larger basin, we do not overwrite that value. In this way, grid cells that are adjacent to two drainage basins are filled with runoff values from the larger basin. After the third iteration, we are left with a field of annual cumulative basin runoff values that have been extrapolated by three 1 km grid cells beyond the present-day ice sheet extent.
The submarine melt parameterisation Eq.
(2) takes as input the subglacial runoff normalised by the submerged area of the  to poor bed topography in some regions, which typically means unrealistically shallow topography in the region of a calving front, we impose a minimum submerged surface area of 0.2 km 2 , equivalent to a glacier of width 2 km and grounding line depth 100 m.

Application of submarine melt parameterisation
Armed with both ocean thermal forcing and subglacial runoff fields defined at annual resolution on 1 km grids and with the submarine melt rate parameterisation Eq.
(2), submarine melt rates may be estimated for the time period 1950-2100 and for each CMIP5 model ( Fig. 1 and Table 1). While this defines a submarine melt rate on every grid cell where both ocean thermal forcing and runoff are defined (Fig. 6c), the intention is that the ice sheet model applies this submarine melt rate only when the 5 model has a calving front within this grid cell. In this way, the ice sheet models may apply a time-varying submarine melt rate to calving fronts around the ice sheet as these calving fronts retreat over the coming century.

Results
We here present the Greenland Ice Sheet ocean forcing arising from the choices and steps made in section 2. The intention is to highlight temporal evolution of the forcing, together with spatial and model-to-model variability, as these factors will drive 10 variability in sea level projections once implemented in an ice sheet model. The results are discussed with the same structure as section 2 and Fig. 1.

Future subglacial runoff
For both implementations, projected subglacial runoff is prescribed for each tidewater glacier using its hydrological drainage basin. We visualise variability in runoff by considering runoff for the largest glacier by ice flux in each sector (Table S1; Lastly, we consider model-to-model variability in projected runoff by averaging over the largest glaciers by sector (Fig. 7c).
The only RCP2.6 scenario considered shows a moderate increase in runoff until 2050 before a return to present-day values by 2100 ( Fig. 7c). All RCP8.5 simulations exhibit a similar temporal evolution and show a significant increase in runoff during the 30 coming century. HadGEM2-ES has the highest runoff at ⇠2000 m 3 s 1 in 2100, with IPSL-CM5A and MIROC5 giving similar results. NorESM1-M and ACCESS1-3 have medium runoff and CSIRO-Mk3-6-0 has the lowest runoff at ⇠1150 m 3 s 1 by 2100. The multi-model spread in runoff anomaly at 2100 is ⇠850 m 3 s 1 , around 50% of the multi-model mean of ⇠1650  There is significant regional variability in projected ocean warming in the MIROC5 RCP8.5 simulation (Fig. 8a). The NE sector stands out with a thermal forcing increase of nearly 5 C, while all other sectors exhibit an increase of between 1 and 3 C. Ocean warming in the NE sector amounts to an increase of 150% in thermal forcing relative to the 1995-2014 baseline period (Fig. 8b). The SE and SW sectors see the smallest relative increase amounting to only ⇠20%. We do note, however, that regional ocean warming differs substantially across CMIP5 AOGCMs (Yin et al., 2011;Barthel et al., 2019, Figs. S6 and 15 S7). The NE sector sees the most warming in MIROC5, HadGEM2-ES and IPSL-CM5A-MR, but the CW and NW regions see equivalent or greater warming in the other three RCP8.5 models. It is also interesting to note that the relative increase in runoff (Fig. 7b) is much larger than the relative increase in ocean thermal forcing (Fig. 8b).
We consider ocean warming at the ice sheet scale by taking a mean over the 7 sectors for each CMIP5 AOGCM (Fig. 8c).
For MIROC5 RCP2.6, there is moderate warming of nearly half a degree which persists until the end of the century. This is 20 mostly driven by significant warming in the CW and NW sectors (Fig. S6a) that exceeds warming in these sectors in some of the RCP8.5 simulations (Figs. S6d and S6f). Given the large inter-model variability in ocean warming, this warming feature is likely to be specific to MIROC5 rather than being more broadly representative of RCP2.6 simulations. Among the RCP8.5 simulations, CSIRO-Mk3-6-0 shows the most warming by 2100, reaching 2.8 C above the present-day value. HadGEM2-ES shows the least warming, reaching 1.9 C by 2100. The multi-model spread in thermal forcing anomaly by 2100 is 0.9 C, around 25 35% of the multi-model mean of ⇠2.4 C. Relative to the 1995-2014 baseline value of 4.6 C, thermal forcing is projected to increase by a factor 0.4-0.6 this century under RCP8.5 (Fig. 8c).

Retreat implementation forcing
Projected sector retreat combines the runoff anomaly per glacier (section 3.1), the thermal forcing anomaly per sector (section 3.2) and the ice flux of all glaciers in the sector (section 2.4.3). Thus, sector-to-sector variability in projected retreat arises 30 due to both variability in regional climate and differences in the population of glaciers in each sector.
For the MIROC5 RCP8.5 simulation, the SW sector has the largest retreat (Fig. 9a) because it has a small number of glaciers (Table S1) each experiencing a large increase in subglacial runoff (Figs. 7a-b). The projected retreat for the CW sector is also high (Fig. 9a), partly due to large projected retreat for Jakobshavn, which dominates the sector-average retreat because it alone accounts for around half of the present-day ice flux in the CW sector (Table S1). Projected retreat is smallest for the NW and NO sectors (Fig. 9a) because these sectors comprise a large population of smaller glaciers (Table S1) and experience the least absolute increase in subglacial runoff (Fig. 7a). Fig. S8 shows equivalent plots to Fig. 9a for all other CMIP5 AOGCMs, in which the spatial patterns of retreat are similar in almost all models with large projected retreat for SW and CW and smaller 5 retreat for NW and NO. Note that Fig. 9a shows only the medium retreat case for each sector; low and high projections are plotted in Fig. S9.
To provide an ice sheet-wide view of retreat per CMIP5 AOGCM, we combine the sector-by-sector projections (e.g. Fig. 9a) into an ice sheet projection by weighting according to the present-day ice flux (Table S1). The resulting projections (Fig. 9b) are not used to force the ice sheet models (the ice sheet models are forced by the sector-by-sector projections) but they do 10 illustrate multi-model variability in projected retreat. The RCP2.6 simulation considered shows moderate retreat of ⇠2 km until 2050 and then a stabilisation of terminus positions (Fig. 9b). The retreat is largely driven by significant ocean warming in the CW and NW sectors (Figs. S6a and S8a).
The RCP8.5 projections show ⇠15 km of retreat by 2100. The retreat rate generally increases throughout the century, so that ⇠4 km of retreat occurs before 2050 and ⇠11 km between 2050 and 2100. The multi-model spread in retreat by 2100 is only 2 15 km, or 15% of the multi-model mean. The largest retreat is projected using CSIRO-Mk3-6-0 and the least using HadGEM2-ES, although all models are similar. In contrast, the spread in projections resulting from the low and high retreat cases for a given model is generally large. For the MIROC5 RCP8.5 projections, the difference between the low and high retreat cases at 2100 is 14 km, much larger than the multi-model spread (Fig. 9b). The same is true for the low and high cases in all other RCP8.5 models (not shown).

Submarine melt implementation forcing
Projections of submarine melt rates are obtained by combining ocean thermal forcing, runoff accumulated over each glacier's subglacial drainage basin and a calving front submerged area (Eq. (2)). To illustrate the results, we show melt rates for the glacier with the largest ice flux in each region ( Fig. 10; Table S1). These projections do not take into account the motion of glacier termini and, thus, isolate the change in melt rates due solely to changes in future atmospheric and ocean forcing.  (Figs. 10g-h). There is significant spread in projected melt rates in 2100 for the RCP8.5 scenarios, typically amounting to 25-50% of the multi-model mean but substantially more for Humboldt Glacier. When considering a mean over the 7 glaciers, the multi-model spread under RCP8.5 is much smaller than at individual glaciers, with the mean melt rate increasing from ⇠2 md 1 in the present day to ⇠6 md 1 in 2100 (Fig. 10a). Under the RCP2.6 scenario, melt rates show only moderate increases until around 2050, followed by stabilisation or decrease (Fig. 10). Projected RCP2.6 melt rates in 2100 are lower than the present day for Kangiata Nunata . In general, RCP2.6 melt rates do not depart significantly from RCP8.5 melt rates until around 2050.
A similar picture emerges when a larger population of 125 glaciers is considered. Fig. 11 shows histograms of the relative 5 increase in submarine melt rate between a twenty year period the end of the century (2081-2100) and the present day (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) under all 6 of the RCP8.5 models considered. For example, since we consider 58 glaciers in NW Greenland, Fig. 11a has a total count of 58 x 6 = 348. In SE and SW Greenland, melt rates increase by at most 170%. These regions already experience a warm ocean and atmosphere in the present day and so large increases in absolute melt rate (Figs. 10g-h) appear as smaller relative increases in submarine melting. Moving north, CE, CW and NW Greenland experience increases up to ⇠400% while 10 the NE and NO sectors have the largest relative increases in melting, reaching over 1000%. These northerly regions have a particularly cold ocean in the present day and currently experience very little submarine melting (e.g. Fig. 10b). Thus any increase in absolute melt rate can constitute a very large relative increase.
The spread in relative melt rate increase within regions (Fig. 11) arises from a number of factors. The glaciers in each region have diverse grounding line depths, submerged in fjords with differing sill depths. Thus, glaciers with deep grounding lines 15 that are directly exposed to the ocean are responding to different water masses than glaciers that are grounded in shallow water or protected from the ocean by shallow sills. If these water masses evolve differently over the coming century then adjacent glaciers may experience very different ocean forcing even within the same CMIP5 AOGCM. A second source of variability is that from the 6 CMIP5 AOGCMs themselves, which can differ substantially on the evolution of ocean temperature within a given sector (Fig. S7).

Retreat and submarine melt implementations
The two implementations have distinct advantages and disadvantages. The retreat implementation has the advantage of being accessible to all ISMIP6 ice sheet models :::::::::::::::::: (Goelzer et al., 2020) and has been empirically validated by tuning to match observed glacier retreat over the past 60 years . In addition, it replaces the need for a representation of calving, the pa-25 rameterisation of which remains a large source of uncertainty (Benn et al., 2017). On the other hand, the retreat implementation does parameterise terminus position in a very constraining manner: it does not allow for modeled ice dynamics to influence the terminus position and it takes no account of bed topography, which is known to be an important factor in determining the response of an individual glacier to an ocean perturbation (e.g. Catania et al., 2018). These issues motivate the second proposed implementation.

30
The submarine melt implementation places less constraints on the interaction between the ocean and ice sheet by specifying only the submarine melt rate (or more precisely, the subglacial runoff, ocean temperature and a parameterisation to combine these quantities to estimate submarine melt rate). The representation of calving and its possible coupling to submarine melting is left to the ice sheet model. This implementation has the advantage that the important interactions between submarine melting, calving, ice dynamics and bed topography can be resolved by the model (e.g. Aschwanden et al., 2019;Morlighem et al., 2019).
The disadvantages are that there is large uncertainty in the submarine melt rates obtained from the parameterisation and we still lack a good understanding of -and parameterisation for -calving. Furthermore, the submarine melt implementation may be considerably more computationally expensive and technically challenging to implement than the retreat implementation.

Variability in projections
The projected relative increase in subglacial runoff (factor 2.5-4.5, Fig. 7) is much higher than for ocean thermal forcing (factor 0.4-0.6, Fig. 8) for all models under an RCP8.5 scenario. Yet both forcings contribute significantly to the retreat and submarine melt rate projections due to the form of the retreat and submarine melt parameterisations (Eqs. (1) and (2)). The subglacial runoff Q appears sub-linearly in these parameterisations, while the thermal forcing (TF) appears approximately linearly, so 10 that the impact of increasing thermal forcing on projected retreat and submarine melting is larger than the impact of increasing runoff by an equivalent relative amount.
There also appears to be some compensation occurring between atmosphere and ocean in the 6 AOGCMs we have considered. The multi-model spread by 2100 in projected subglacial runoff is ⇠50% and in thermal forcing is ⇠35%, but the spread in projected retreat and submarine melting is only ⇠15% (Figs. 7c, 8c, 9b and 10a). The model that has the most ocean warming 15 (CSIRO-Mk3-6-0) has the least runoff increase and the model that has the least ocean warming (HadGEM2-ES) has the most runoff increase (Figs. 7c and 8c). Due to the form of the retreat and submarine melt parameterisations (Eqs. (1) and (2)), the atmosphere and ocean projections can compensate each other, reducing the multi-model spread in the retreat and submarine melt projections. Coupled with the large uncertainty in the linear coefficient  appearing in the retreat parameterisation , the spread in projected retreat due to the low and high retreat cases (section 2.4.4) is, therefore, much larger than 20 the spread in projected retreat or submarine melting due to AOGCM selection (Figs. 9b and 10a). It can therefore be expected that the spread in sea level projections arising from the use of the low and high retreat scenarios will be larger than from the medium retreat or submarine melt rate scenarios forced by different CMIP5 AOGCMs. In terms of sampling uncertainty on future sea level within the implementations presented here, it may be more beneficial to prioritise ice sheet simulations sampling uncertainty in coefficients in the parameterisations rather than considering additional AOGCMs. We note that this may not be 25 true for a different ocean forcing implementation and that we have only considered 6 CMIP5 AOGCMs in this study (Table 1) that are a selected subset of the larger CMIP5 ensemble (Barthel et al., 2019). It is possible that use of other CMIP5 AOGCMs would lead to a greater spread in projected retreat and submarine melting.
Examination of the projected submarine melt rates also suggests the possibility for sector-by-sector compensation. For example, CSIRO-Mk3-6-0 has the highest projected melt rates of any RCP8.5 model at Kong Oscar, Jakobshavn and Helheim 30 but is close to the lowest projected melt rates at Humboldt, Daugaard-Jensen and Kangiata Nunata Sermia (Fig. 10). There is no individual CMIP5 AOGCM that gives high melt rates in every single sector or at every single glacier; rather a model that gives high melt rates in a certain sector often gives lower melt rates in another sector. As a result, taking a mean of the projected RCP8.5 melt rates over 7 large glaciers gives trajectories that lie within a narrow envelope (Fig. 10a). Once again, this may act to reduce the spread in the projected sea level from ice sheet models forced by these melt rates. The response of individual sectors or glaciers may differ substantially between CMIP5 AOGCMs, but Fig. 10 suggests that glaciers and sectors may compensate one another, leading to a similar sea level contribution from the full ice sheet under each CMIP5 AOGCM.
The dynamic sea level contribution is not, however, directly related to the magnitude of retreat or submarine melt rate.
For example, although the SW sector has the largest projected retreat, it contains relatively few tidewater glaciers and these 5 glaciers currently account for <4% of Greenland's ice discharge (Table S1). It is therefore unlikely to be a major source of dynamic sea level contribution in the future. In contrast, the NW region has the smallest projected retreat but has a large number of tidewater glaciers that currently account for ⇠20% of Greenland's ice flux (Table S1) and is much more likely to be a significant dynamic contributor to sea level. Within the submarine melt implementation there is also the possibility for nonlinear or threshold response of glaciers to submarine melting, where small changes in forcing may result in large excursions in 10 terminus position and mass loss .

Impact of bias corrections
In order to provide continuous ocean forcing from the present day into the future and to ensure we sample uncertainty in future climate projection rather than representation of the present day, we bias-corrected the subglacial runoff (section 2.2.3) and ocean thermal forcing (section 2.3.3). Due to the non-linearity of the retreat and submarine melt parameterisations, the 15 bias corrections do impact the projected retreat and submarine melt rate. Where there exists uncertainty in the present-day quantities, for example the ocean thermal forcing in CW, NW and NO Greenland, this leads to uncertainty in the projections.
Compared to the situation in which no bias correction is performed, the bias correction can change RCP8.5 projected retreat by 2100 by up to a few kilometers, or around 0 to 20% of the typical retreat by 2100 of 15 km (Figs. 9b and S10). The bias correction is equally likely to increase or decrease the projected retreat (Fig. S10). There are a few instances where the impact of 20 the bias correction is larger, for example in NorESM1-M the bias correction decreases projected retreat by 36% in SE Greenland and increases it by 20% in CW Greenland (Fig. S10). These follow from the large ocean thermal forcing bias corrections applied to this model (Fig. S3). The bias corrections can, therefore, contribute to sector-by-sector differences in retreat projections but do not overall increase or decrease projected retreat. Since the retreat and submarine melt parameterisations have a similar form, the impact of the bias correction on submarine melting will be similar.

Missing processes and priorities for future improvement
Due to the complexity and timescale of the exercise we have had to make a number of simplifications of complex processes in order to deliver the ocean forcing to the ice sheet modeling groups. One key simplification is our treatment of the ocean thermal forcing experienced by tidewater glaciers. Since the CMIP5 ::::: CMIP AOGCMs do not resolve Greenland's fjords, we have had to bridge the gap between the continental shelf and calving fronts. In the retreat parameterisation, the ocean thermal 30 forcing applied to glaciers is a spatially-and depth-averaged value from the continental shelf. Thus, we have neglected spatial gradients in ocean temperatures within the chosen sectors , the processes responsible for transporting and transforming ocean waters between the shelf and calving front (Motyka et al., 2003;Straneo et al., 2010;Mortensen et al., 2011;Jackson et al., 2014;Gladish et al., 2015) and the diverse grounding line and sill depths of glaciers and fjords in Greenland (Morlighem et al., 2017). We do note that the retreat parameterisation Eq. (1) was tuned based on observations from 1960 to present using the same definition of ocean thermal forcing and so, to some extent, all of these processes will have fed into the empirical tuning. This definition of ocean thermal forcing nevertheless neglects much of the individuality of glacier-fjord systems, essentially linking groups of glaciers to large-scale ocean changes only. 5 In the submarine melt implementation, the effect of sills and grounding line depth is taken into account by retaining the depth-variability of ocean conditions and extrapolating these properties into fjords based on the bathymetry. Certainly, the presence of sills is known to modify fjord water properties substantially by blocking access of dense waters to the calving front (Gladish et al., 2015), but this extrapolation remains a simplification due to vertical mixing within fjords (e.g. Inall et al., 2014) and because periodic dense inflows over sills have been observed in Greenland (Mortensen et al., 2011). Therefore, 10 both the retreat and submarine melt implementations would be improved with methods to quantify water mass transformation between the shelf and calving fronts. Such methods might take the form of very high-resolution regional ocean modeling or, perhaps more practically for efforts such as ISMIP6, simple parameterisations or fjord box models. Knowledge of fjord and shelf bathymetry is a prerequisite for these improvements but is currently incomplete and is, therefore, an additional priority.
Both implementations also assume that submarine melting is the primary climate forcing experienced by the calving fronts 15 of tidewater glaciers. This assumption derives from the literature consensus on the important role played by submarine melting in the recent retreat of tidewater glaciers in Greenland (Holland et al., 2008;Straneo and Heimbach, 2013;Fried et al., 2015;Cowton et al., 2018), yet other processes may play a substantial role. In particular, the buttressing provided to glaciers by ice mélange may be sufficient to suppress calving (Amundson et al., 2010;Robel, 2017), has been implicated in rapid glacier retreat (Christoffersen et al., 2012;Moon et al., 2015;Bevan et al., 2019) and is found to be more influential than submarine 20 melt in some models (Krug et al., 2015;Todd et al., 2018). Future ice sheet-ocean forcing efforts might, therefore, look to incorporate the impact of ice mélange buttressing.
Once submarine melting is assumed to be the primary ocean forcing, it must be parameterised, as has been done in Eqs. (1) and (2) for the retreat and submarine melt implementations, respectively. The form of both parameterisations derives from the physics of plumes, which, aside from the submarine melt they induce, are relatively well understood from theory, laboratory 25 and observational work (Morton et al., 1956;Jenkins, 2011;Jackson et al., 2017). Observations of submarine melting with which to constrain key constants in melt parameterisations are, however, severely lacking. Our first direct observations of submarine melting were obtained very recently in Alaska Sutherland et al., 2019) and suggest we may currently be underestimating submarine melt rates, especially outside of plumes. For the retreat implementation, uncertainty in melt parameterisations is less of an issue because the parameterisation assumes proportionality between glacier retreat and 30 submarine melt rate, and since glacier retreat is easily observable, we have good observations to tune the linear coefficient  . This is not the case for the submarine melt implementation, though ice sheet models typically do a spin-up simulation in which they tune their model to try to match present-day ice sheet extent, which may go some way to reducing their sensitivity to uncertainty in the melt parameterisation. It is clear, however, that observations of submarine melting and further work building on Sutherland et al. (2019) would be valuable for reducing uncertainties on sea level contribution in efforts beyond ISMIP6.

Summary
The Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6) constitutes the primary community effort to produce ice sheet sea level projections for the next Intergovernmental Panel on Climate Change Assessment Report (IPCC AR6). ISMIP6 is 5 the first effort to develop a multi-model ensemble of Greenland Ice Sheet models forced by ocean boundary conditions derived from CMIP AOGCMs. Such a strategy is demanding to design due to the evolving nature of our process understanding and ice sheet model technical capabilities. With these challenges in mind, we have proposed two ocean forcing strategies, called the retreat implementation and the submarine melt implementation. By combining these strategies with projected climate from selected CMIP AOGCMs, we have derived ocean boundary conditions for Greenland Ice Sheet models to run 21st century 10 projections ::::::::::::::::: (Goelzer et al., 2020).
In the retreat implementation, retreat is projected using a process-motivated but empirically-calibrated parameterisation that combines subglacial runoff and ocean thermal forcing to estimate tidewater glacier retreat . Retreat is projected for each individual tidewater glacier but for simplicity is applied to the ice sheet homogeneously within each of 7 sectors. Under a high greenhouse gas emissions RCP8.5 scenario, projected retreat that will be applied to the ice sheet 15 models amounts to around 15 km by 2100 with a range of 10-25 km in low and high scenarios. Under a low emissions RCP2.6 scenario, retreat of only ⇠1 km will be prescribed. In the submarine melt implementation, fields of subglacial discharge and ocean thermal forcing covering Greenland are provided, together with a recommended parameterisation that may be used to estimate submarine melt rate wherever a calving front is located. Under RCP8.5, projected melt rates in 2100 are a factor ⇠3 higher than the present day but remain relatively constant under RCP2.6. The sea level contributions resulting from these two 20 implementations will be determined by the modeled dynamic response to these forcings.
The proposed implementations are driven by process understanding but are also pragmatic and have necessarily neglected certain processes or made use of poorly-constrained parameterisations. Foremost amongst these are fjord processes and the transformation of ocean waters between the continental shelf and glacier calving front and the parameterisation of submarine melting. These issues are to some extent ameliorated through tuning, both in the described implementation and at the level of      (2). Note that the thermal forcing and melt rate values in the ice sheet interior are included only to show that the submarine melt implementation defines melt rate everywhere that is below sea level and connected to the ocean. An ice sheet model would only apply these melt rates if the ice sheet margin retreats into the interior, which is unlikely by 2100. Also note that runoff values are only plotted for marine-terminating hydrological basins.  are combined according to their present-day relative ice flux (Table S1). Also shown in the shading is the low and high retreat projections for MIROC5 RCP8.5. Note that the ice sheet models are forced on a sector-by-sector basis, so the projections in (b) are not used to force any models but are included to give a sense of the multi-model variability. See Fig. S8 for full plots of all projections.  Data availability. The bed topography and bathymetry used in this work may be downloaded from https://nsidc.org/data/IDBMG4 (last access September 2019). Information on the RACMO2.3p2 SMB data can be found at http://www.projects.science.uu.nl/iceclimate/models/greenland.php (last access April 2019). EN4.2.1 oceanographic data is available at https://www.metoffice.gov.uk/hadobs/en4/download.html (last access April 2019). CMIP5 model output is available at https://esgf-node.llnl.gov/projects/esgf-llnl/ (last access April 2019). The MAR based future subglacial discharge projections are available on ftp://ftp.climato.be/fettweis/MARv3.9/ISMIP6/GrIS/ (last access April 2019). TEOS-10 5 routines may be found at http://www.teos-10.org (last access September 2019). Further information on the ISMIP6 project may be found at http://www.climate-cryosphere.org/activities/targeted/ismip6 (last access April 2019). We intend to make all code used to create these projections freely available on the ISMIP6 GitHub page at https://github.com/ismip (last access September 2019). All of the projection datasets described in this paper are freely available from the ISMIP6 ftp server; access can be obtained by emailing ismip6@gmail.com.
Author contributions. DS undertook the majority of the analysis, processing, writing and creation of figures. DF processed past and future 10 runoff and contributed to writing and figures. FS led the ISMIP6 ocean forcing and provided oversight at all stages of the process. HG provided invaluable guidance on the implementation of the described ocean forcing in ice sheet models. CML provided CMIP5 model output and expertise. MM performed the extrapolation of ocean properties into fjords in the submarine melt implementation. XF ran MAR simulations forced by the selected CMIP5 models. SN coordinated the ISMIP6 effort. All authors took part in extensive discussion of the methodology and edited the manuscript.