FAMOUS version xotzb (FAMOUS-ice): a GCM capable of energy- and water- conserving coupling to an ice sheet model

The physical interactions between ice sheets and their surroundings are major factors in determining the state of the climate system, yet many current Earth System models omit them entirely or approximate them in a heavily parameterised manner. In this work we have improved the snow and ice sheet surface physics in the FAMOUS climate model, with the aim of improving the representation of polar climate and implementing a bidirectional coupling to the Glimmer dynamic ice sheet model using the 5 water and energy fluxes calculated by FAMOUS. FAMOUS and Glimmer are both low resolution, computationally affordable models used for multi-millennial simulations. Glaciated surfaces in the new FAMOUS-ice are modelled using a multi-layer snow scheme capable of simulating compaction of firn and the percolation and refreezing of surface melt. The low horizontal resolution of FAMOUS compared to Glimmer is mitigated by implementing this snow model on sub-gridscale tiles that represent different elevations on the ice sheet within each FAMOUS grid-box. We show that with this approach FAMOUS-ice 10 can simulate relevant physical processes on the surface of the modern Greenland ice sheet well compared to higher resolution climate models, and that the ice sheet state in the coupled FAMOUS-ice-Glimmer system does not drift unacceptably. FAMOUS-ice coupled to Glimmer is thus a useful tool for modelling the physics and co-evolution of climate and grounded ice sheets on centennial and millennial timescales, with applications to scientific questions relevant to both paleoclimate and future sea level rise. 15 1 https://doi.org/10.5194/gmd-2020-207 Preprint. Discussion started: 23 July 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Continental-scale ice sheets are one of the primary components in the Earth's climate system. The climatic influences that result as they grow and shrink are key features of the global-scale glacial cycles of the last million years. Their waxing and waning volume causes sea-level variations of over 100m in amplitude over these cycles (Spratt and Lisiecki, 2016) through their barystatic contribution (the global-mean effect of changing the mass of the ocean), by altering the gravitational field and 5 rotation of the Earth and by deforming its solid surface (sea level concepts and terminology are reviewed by Gregory et al. (2019)). Ice sheet mass loss accounts for around a third of the present rate of global mean sea-level rise; their contribution is expected to increase and likely eventually dominate sea-level change in the coming decades and centuries (Oppenheimer et al., 2019).
The role of ice sheets in the climate system involves strong bidirectional interactions with both the atmosphere and ocean. 10 As ice sheets evolve they substantially modify the surface radiation balance and temperature by changing the albedo and altitude of the Earth's surface, the supply of fresh water to the ocean, and the circulations of both the atmosphere and the ocean (e.g. Golledge et al. (2019)). The diversity and strength of these feedbacks means that ice sheets should properly be modelled as physically coupled, interactive components in climate models when addressing many scientific questions related to the evolution of the Earth system. Such questions include projections of the impacts of future climate change on multi- 15 decadal timescales and longer due to the increasing influence of ice sheets on the global mean sea-level budget (e.g. Edwards et al. (2014); Nowicki et al. (2016)). However, current understanding of the co-evolving physics of our climate and ice sheets is limited in detail. This hinders our ability to explain major features of climate change in the past and make projections about what may happen in the future.
There are numerous challenges in successfully modelling the physics of ice sheet-climate interactions in the atmosphere-20 ocean general circulation models (AOGCM) that are commonly used for comprehensive climate studies. Many ice sheetrelevant processes have characteristic length scales of kilometres or less (e.g. gradients of precipitation and surface melt on the sloping margins of Greenland, or ocean melt at the grounding line of Antarctica's floating ice shelves). This is much less than the grid-box size of tens or hundreds of kilometres in contemporary AOGCMs (e.g. UKESM1 (Sellar et al., 2019) or CESM2 (Danabasoglu et al., 2020)). Many of the ice sheet features that depend on these small scale processes can, however, evolve 25 on millennial timescales, implying climate simulation lengths that are computationally unaffordable with complex AOGCMs.
The technical structure of most AOGCMs also makes them ill-suited to changing the boundaries of the land or ocean domain as they run. This issue is a particular problem where the evolution of the ice involves the collapse of an ice shelf or a change in global sea-levels. Furthermore, many of the physical processes that are key to modelling the surface physics of polar regions are not well captured by models that are intended for global use, such as stable polar boundary layers, katabatic winds or 30 multi-year snow packs (Connolley and Bracegirdle, 2007).
For all these reasons, the majority of AOGCMs do not include interactive ice sheets, and omit much of the physics required to directly model boundary conditions for ice sheet models. For their part, ice sheet modellers have developed empirical parameterisations to translate the climate fields that AOGCMs do provide into usable boundary conditions (e.g. Reeh (1991); developed to include more sophisticated surface schemes for ice sheet regions, and some have coupled these schemes to ice sheet models (ISM) (e.g. Sellar et al. (2019); Danabasoglu et al. (2020)).
The mismatch in spatial resolution between the grids used for key atmosphere and ice sheet processes has been addressed either through explicitly modelling the local atmosphere at much higher resolution in a regional or nested model (van , or by calculating sub-gridscale surface fluxes and mass balance terms for portions of atmospheric grid-boxes 15 (Ganopolski et al., 2010;Vizcaíno et al., 2013;Ziemen et al., 2014;Sellar et al., 2019). Using sub-gridscale methods is often more computationally affordable than increasing the explicit resolution of the atmosphere model, so allows much longer coupled climate-ice sheet simulations to be carried out. Many of the quantities required for these sub-gridscale calculations are well correlated with surface temperature, and in regions with significant topographic gradients-like the margins of an ice sheet-are also a strong function of altitude. Several models now exist that treat each land surface grid-box in ice sheet 20 areas as a collection of sub-gridscale tiles, each representing a different elevation range within the grid-box (e.g. Vizcaíno et al. (2013); Sellar et al. (2019)). The tiles are aggregated into a grid-box average for communication with the atmosphere, but used individually when required to provide finer-grained information to an ice sheet model. FAMOUS is a relatively coarse resolution AOGCM that has been used successfully in a wide range of climate studies (e.g. for multi-millennial simulations. In previous coupled climate-ice sheet studies with FAMOUS (Gregory et al., 2012;Roberts et al., 2014), the coupling was achieved through a simple positive-degree-day parameterisation and it was later found that the long-term response of the ice was very sensitive to ill-constrained empirical parameters within this scheme. Building on FAMOUS's climate simulation and a framework for sub-gridscale modelling within its land surface scheme, we decided to improve FAMOUS's representation of ice sheet surface physics and to directly couple water and energy fluxes between 30 FAMOUS and the ice sheet model.
In this paper we describe modifications that we have made to the physics and tiling in the land surface of FAMOUS to produce a configuration called FAMOUS-ice, and we show that the results are scientifically useful within the context of coupling to a model of the Greenland ice sheet (GrIS). We do not attempt to model any ocean-ice sheet interaction in this work, as the fjords and marine-terminating glaciers of Greenland are not resolved on the FAMOUS grid. We focus primarily 35 on technical aspects of our work, including the scheme by which FAMOUS-ice has been directly coupled to the Glimmer ice sheet model. A scientific application of FAMOUS-ice, looking at the future stability of the Greenland ice sheet, is being published separately (Gregory et al, 2020, submitted manuscript).

Standard FAMOUS
FAMOUS (Smith et al., 2008;Williams et al., 2013) is a configuration of version 4.5 of the UK Met Office Unified Model 5 (UM4.5 (Gordon et al., 2000)). The atmosphere in FAMOUS has a horizontal resolution of 7.5 • x5 • , with 11 vertical levels. The development described here is based on the most recently released configuration of FAMOUS, xfhcu (Williams et al., 2013).
In the present work, the ocean model is replaced by prescribed sea surface boundary conditions taken from a set of higher resolution AOGCM simulations. This allows us to evaluate the performance of the ice sheet surface mass balance scheme with minimal bias in the wider simulated climate. 10 FAMOUS xfhcu uses the MOSES2.2 land surface model (Essery et al., 2003), an early version of the JULES land surface scheme (Best et al., 2011) used in current UM configurations (e.g. HadGEM3 and the UK Earth System model (UKESM) (Sellar et al., 2019;Williams et al., 2018)). MOSES2.2 shares the same grid resolution as the overlying atmosphere model, and can represent heterogeneity in land surface characteristics by using a set of sub-gridscale tiles in each grid-box, each of which simulates a different type of surface. Each tile covers a particular fraction of the area of its grid-box, but is not explicitly 15 geographically delimited. Effects at the grid-box-scale are given by an average across the set of tiles in that box, weighted by the fraction of the area of the grid-box that each tile covers. Tiles in MOSES are generally used to represent different types of surface vegetation. FAMOUS xfhcu, like most UM configurations with tiles to date, uses 9 different surface types, one of which represents permanently glaciated surfaces. The underlying soil in xfhcu is only represented by entire grid-boxes; all surface tiles in a grid-box share a common sub-surface. Because the glaciated ice surface tile requires its own distinct ice sub-surface, 20 land ice cannot share a grid-box with any other surface tile, and entire grid-boxes must be either entirely glaciated or ice-free.
This means that the spatial representation of the edge of an ice sheet in MOSES2.2 is usually very blocky in a low resolution UM configuration like FAMOUS. MOSES2.2 as described in Essery et al. (2003) has a "zero-layer" snow model, which represents only the bulk properties of the snow. This model assumes that snow has a fixed density (250 kg/m 3 ) and, if deep enough, can insulate the surface 25 from the atmosphere. Snow surface albedo is calculated in four radiation bands (direct and diffuse radiation in visible and infrared bands) from the grain size of the snow, which may evolve with time and surface temperature. The underlying ice tile has a prescribed albedo of 0.75 in all bands. The performance of the MOSES albedo scheme and the lack of internal snow pack refreezing of surface melt were identified as particular shortcomings for the calculation of ice sheet surface mass balance (SMB) in previous work with the UM (Rae et al., 2012). To improve the representation of SMB in glaciated areas in 30 FAMOUS-ice we have significantly modified the capabilities and behaviour of the tiles, snow pack model and snow and ice surface albedos. The body of the snow scheme we describe in the following has largely been backported from version 3 of JULES, with further modifications to both the snow and the tiling system that have been done largely in parallel with work in more recent versions of JULES and UKESM.
3 Tiles for fractional ice extent and sub-gridscale elevation In FAMOUS-ice, we have modified the MOSES2.2 tiling scheme to allow each of the 9 basic surface types to exist at multiple fixed-height elevations within a grid-box. Thus, calculating surface exchange for each tile separately, we simultaneously 5 simulate the different conditions that would be seen on ice surfaces at, for example, sea-level, 100m and 500m, with the same large-scale atmospheric column above the boundary layer for each. In this paper we have configured these height-dependent tiles only on Greenland to minimise computational cost, although they could be used globally in principle. The resulting threedimensional (longitude, latitude, height) arrays allow the climate model output for GrIS to be interpolated in both horizontal and vertical dimensions onto a finer resolution ice sheet model grid, while the atmosphere uses the grid-box mean of the tiles, 10 weighted by their area fractions. All the fine-grid climate information seen by the ice sheet model is thus consistent with the large-scale average in the climate model as it evolves.
In FAMOUS-ice, the ice and soil sub-surface models are allowed to co-exist, so that each grid-box may contain any combination of ice and non-ice surface tiles at the designated altitudes. This allows the atmosphere to see a sub-gridscale, fractional representation of area at the edge of the ice sheet, rather than binary blocks of entire ice/not-ice grid-boxes. The soil and ice 15 sub-surface models still operate at the grid-box-scale and are not independent for each surface tile. This means that tiles at all elevations are coupled to a common sub-surface layer. We will describe later how tiles at different elevations communicate with their sub-surface in such as way as to minimise the spurious flow of information between them. Figure 1 compares the frequency distributions of height surfaces for a sample GrIS topography in FAMOUS-ice, using the standard FAMOUS grid-box mean orography and two choices of distributions of elevation tiles. The finer-scale detail generated 20 for the ice sheet model comes from an increase in the total number of surface tiles simulated and a closer approximation to the true frequency distribution of surfaces with height as the number of elevation classes used within each grid-box increases.
Increasing the number of tiles used implies both computational overhead and an increase in the memory required for I/O, which can be slow in FAMOUS, so it is wise to not to use more tiles than really needed for a given application. In practice, we find the use of 10 tiles, with the same vertical distribution as is used in CESM1 (Vizcaíno et al., 2013) with boundaries at [0, 200, 25 400, 700, 1000, 1300, 1600, 2000, 2500, 3000, 10000] metres, is a good compromise. This increases the number of surfaces calculated on GrIS in FAMOUS by a factor of 10 compared to using the grid-box mean, and adds many more surfaces at low elevations which is important for simulating higher ablation areas.

Downscaling climate to the tiles
In order to reproduce differences in snow pack and ice evolution on different elevation tiles in a grid-box-for example snow/ice 30 mass loss (ablation) at lower, warmer altitudes and net accumulation above the equilibrium line-each tile must experience atmosphere variables adjusted to the particular altitude that it represents. This requires a downscaling procedure that is computationally fast enough to be applied timestep-by-timestep in the atmosphere model, and which conserves the grid-box-mean of the tiles and the physical relationships between the variables being scaled. We require the scheme to provide good results for local altitude adjustments only over a specific ice sheet; the downscaling does not need to work globally. We focus further on downscaling only those fields which are required to produce significant differences in surface mass balance at different 5 altitudes in a single grid-box.
Global, fixed lapse rates are often used to downscale climate model variables onto higher resolution ice sheet surfaces, especially temperature for use in degree-day or other temperature-index schemes (e.g. Roche et al. (2014); Vizcaíno et al. (2013)). Although in preliminary work we tried more sophisticated methods, sampling the local atmospheric lapse rates for each grid-box and attempting to maintain net zero averages in the sum of the adjustments for the tiles, we found the best 10 results for downscaling near-surface air temperature and downwelling longwave radiation onto the elevation tiles came from prescribing spatially constant lapse rates. Optimal values for the lapse rates of 6 K/km (for temperature) and 3.6 W/m 2 /K for longwave radiation were found by comparing lower-atmospheric conditions in FAMOUS-ice downscaled to elevation tiles over GrIS with output from the MAR regional model (Fettweis et al., 2013). This empirically tuned approach is open to question if the ice sheet (in a coupled model) is allowed to evolve to a significantly different state from that the lapse rates were calibrated 15 for. However, we gain confidence from the results of our model in a climate-change experiment in which the Greenland ice sheet is reduced to a small fraction of its present size; we find that the ice sheet area-mean summer-mean surface air temperature change as a function of change in surface altitude is close to our chosen lapse rate (Gregory et al., 2020, submitted manuscript).
Specific humidity is downscaled using a lapse rate derived from the local atmospheric lapse rates for each grid-box. Precipitation is not adjusted for the tiles, either in terms of the magnitude or partitioning between snow and rain, because of complications in robustly maintaining consistent vertical energy budgets in the atmospheric column when changing the amount or phase of moisture. Each tile in a grid-box receives the same amount and phase of precipitation, calculated with reference to the mean grid-box orographic height. Downwelling shortwave radiation and surface momentum fluxes are also not downscaled because 5 in FAMOUS-ice they are found to vary negligibly across the elevation range used for the tiles.
All tiles within a grid-box share the same soil or ice sub-surface, regardless of their elevation, so it is also necessary to adjust the bottom boundary layer temperature of each elevation tile to limit the potential for unrealistic heat flows between the tiles via their shared sub-surface. For example, a high-latitude, cold ice tile with a temperature of 250K at the bottom of the snow pack should not see the same sub-surface boundary temperature as a tile at sea level in the same grid-box where conditions 10 may be causing the snow to melt. A sub-surface temperature lapse rate of -1 K/km was used. We have not attempted to adjust moisture availability to the elevated tile from the soil sub-surface.
Example results show how surface exchange fluxes vary with height on the ice sheet in FAMOUS-ice compared to those from regional climate model (RCM) simulations of GrIS in MAR (Fettweis et al., 2013) and RACMO (Noël et al., 2018) (figure 2). Average vertical gradients of all quantities are of similar magnitude and sign as in the RCMs, improving on the 15 elevation class implementation analysed by Sellevold et al. (2019). Summer surface albedo in FAMOUS is clearly smaller than both RCMs at most heights on the ice sheet -this will be discussed in section 6.

Multi-layer snow-ice scheme
In FAMOUS-ice, we have replaced the MOSES zero-layer snow scheme with the multi-layer snow scheme now available in JULES. This model simulates vertical gradients of temperature and density within the snow pack, mechanical compaction, and 20 the internal percolation, retention and refreezing of melt. The scheme was originally devised to simulate seasonal snow packs, so we have modified it for perennial snow, firn and solid ice as described below. We use it in FAMOUS-ice to model the whole upper portion of the ice sheet that is relevant for SMB and climate.
The pore space available for melt-water retention in the JULES multi-layer snow scheme scales with the depth of the snow pack. Surface melt percolates downward through the snow, being retained where there is unfilled pore space in a layer (and 25 refrozen if the layer is cold enough) or passed on to the layer beneath if there is not. Where snow accumulates year on year and has become many metres deep, scaling the available pore space solely with layer depth leads to the unrealistic retention of all the surface melt. To counter this effect we additionally make pore space availability a function of the density of snow in each layer. Snow layers in our scheme compact under pressure from the overlying snow. Once the layer density increases above a threshold, the potential pore space amount is reduced linearly, reaching zero when the density reaches that of solid ice.
where: S is a liquid-retention parameter = 0.05, dz k is the thickness of snow layer k, ρ water is the density of water = 1000 kg/m 3 and ρ snow,k is the density of snow layer k.
The most recently-released UK Met Office AOGCM, HadGEM3-GC3 (Williams et al., 2018), is their first to use the JULES 5 multi-layer snow scheme. HadGEM3-GC3 is configured with 3 snow layers, with interfaces at 0.04m and 0.16m. This configuration has been shown to simulate seasonal snow well. In the JULES multi-layer scheme, the bottom layer is allowed to thicken as much as is required to allow the total snow mass to be accounted for, but the pack gains no additional vertical resolution. In this 3 layer configuration, the mass of snow in the two thin upper layers can never be enough to compact the lowest layer significantly. The density of the lowest layer thus remains permanently near its initial low value, however much 10 snow accumulates and however long it stays there. A perennial, deep snow pack in this configuration thus does not become compacted and retains all surface melt in its pore space, producing no runoff.
For FAMOUS-ice we retained three thin top layers, and added additional layers beneath to be used for perennial snow in metres. On glaciated tiles, we initialise the snow pack with 100 metres of ice-density snow, which is unable to hold melt-water. 15 The upper layers contain sufficient mass to cause compaction to occur in lower layers. The seasonal signal of temperature variation is resolved in the uppermost layers, while the lowest layers are insulated from variations on sub-annual timescales.
This configuration of the multi-layer snow scheme provides a representation of snow physics both for glaciated areas and for seasonal snow in other areas.
6 Snow albedo 20 The basic FAMOUS radiation scheme divides the incoming shortwave into two spectral bands, one for the visible and one for near-infrared, with distinct albedos. For snow, the albedos are based on the surface grain size. Where the snow pack has completely melted away and the underlying sub-surface is exposed to the atmosphere, albedos for each band are instead prescribed according to sub-surface type; for the glaciated portion of the grid-box, this is bare ice, which has lower albedo than snow. The albedo over GrIS in the original MOSES2.2 FAMOUS xfhcu (Williams et al., 2013) is close to that of fresh 25 snow everywhere, with little spatial or temporal variation. This is attributable to the lack of significant grain-size evolution in the zero-layer snow model and to the chosen albedo being unrealistically high. In both respects, the HadRM3 configuration of UM4.5 is similar to FAMOUS xfhcu, and they contribute to the poor simulation of GrIS SMB noted by Rae et al. (2012).
There are 4 areas where the snow and ice albedos in FAMOUS-ice have been modified: 1. The multi-layer snow scheme itself provides for a more effective coarsening of the snow pack grain-size compared to the old zero-layer snow model, and this improves the seasonal variation of albedo over snow surfaces globally. The increase in sensitivity of the grain-size may be because the thin surface layers are allowed to evolve independently with temperature, rather than having to change the full bulk of the snow pack.
2. For the grain-size-dependent calculation of snow albedo, parameters in the Marshall (1989) implementation of Wiscombe 5 and Warren (1980) have been tuned to make the albedo lower for larger grains.
3. The albedo of the bare ice sub-surface has been tuned to lower values, and includes a dependency on air temperature as the surface rises towards 0 • C.
4. If there is a snow pack on an ice tile but the surface layers of the snow have a density approaching solid ice, then the bare ice albedo is used, rather than the albedo derived from the grain-size calculation, which is not appropriate for dense firn.

10
In JULES the albedo of seasonal snow, a snow , is formulated in four discrete bands (for direct and diffuse components of both visible and near-infrared wavelengths) dependent on the snow grain size, using the Marshall (1989) parameterisation. In FAMOUS-ice, for the visible wavelengths we have tuned the maximum snow albedo to 0.95 and increased the sensitivity to changes in surface grain size, ∆a snow,visible so that our albedo has a grain-size dependency closer to the results of Roesch et al. The surface albedo on ice sheets, a icesheet , in FAMOUS-ice is formulated as where a bare ice is the albedo of a bare ice surface, a max is the maximum albedo of a bare ice surface averaged across the visible and near-infrared bands (set to 0.55), a min is the minimum albedo of a bare ice surface averaged across the visible and near-infrared bands (set to 0.2), ∆a ice is the sensitivity of bare ice albedo to surface temperatures once the surface is expected 25 to melt (set to -0.35 K −1 ), T surface is the surface temperature, T threshold is the threshold grid-box temperature above which melting is expected somewhere at the surface (set to 272K), ρ surface is the density of the snow pack at the surface and m snow is the mass of the snow pack on the ice tile. This formulation allows fresh snow on the ice sheet to have the same grain-dependent albedo parameterisation as for snow on non-ice sheet surfaces, but to transition to a lower albedos characteristic of bare ice surfaces when denser layers of the perennial snow pack are exposed at the surface. Bare ice albedo values have been tuned to match the range of ice albedos present in a MAR simulation of present-day GrIS. The temperature dependency of the ice albedo for temperatures above T threshold is a simple parametersation that mimics the presence of pooled melt-water or biological activity, which have been 5 shown to significantly lower surface albedo in certain regions of the ice sheet (Greuell, 2000;Williamson et al., 2020); these processes are not explicitly included in our surface scheme.
The resulting albedo for the GrIS in FAMOUS-ice does indeed reproduce lower values than in FAMOUS xfhcu at lower elevations and for melting surfaces. However, the FAMOUS-ice simulation now features a widespread low bias higher up on the ice sheet during summer months when compared to MAR. This is likely due to overestimating the altitude at which melt 10 can occur, increasing the snow grain size and triggering the temperature dependence that lowers albedo for higher surface temperatures. FAMOUS also has a low bias in downwelling shortwave over GrIS (fig:4), likely due to excessive cloud over the ice sheet. These two biases act in opposite directions, and the amount of solar radiation absorbed by the ice sheet is itself in line with RCM results. It should be noted that the configurations that use the end members of our range of tunings for ∆a snow,visible both also have a low-biassed albedo in their present-day simulations; this low albedo bias seems a robust feature 15 of FAMOUS-ice simulations that have a realistic ice sheet shape in the present day, but it does not necessarily imply a simple over-or under-estimate of the sensitivity of the albedo to future changes in climate. FAMOUS-ice tends to produce an insufficiently positive SMB between 1000m and 2000m. The seasonal variation in these averages is of a similar range to the RCMs for positive mass balances, but overestimates the monthly variation in negative SMB, especially at low elevations. SMB is primarily a balance between accumulation and liquid runoff from the snow pack, and it is apparent that the difference in the SMB profile between FAMOUS-ice and the RCMs can be attributed to the surface melt (and thus runoff) component. In both FAMOUS-ice and MAR there is relatively little variation in precipitation with 5 height, and the magnitudes match well. In contrast, RACMO simulates more precipitation at low elevations and shows a steady reduction with height above around 500m, which largely accounts for the differences between the SMB-height profiles for RACMO and MAR. This may be because of the higher native resolution of the RACMO simulation in this case resolving more orographic precipitation near the coast of Greenland.
The spatial pattern of the annual average SMB shows a clear correlation with the height of the ice sheet, as expected 10 ( figure 5). In general, the pattern diagnosed from FAMOUS-ice has smaller magnitudes than MAR in both accumulation and ablation zones. The relatively low spatial resolution of the FAMOUS atmosphere cannot simulate the same degree of intense precipitation near the coasts, especially in the south east. The equilibrium line altitude, demarking the transition between areas of net accumulation and ablation, is generally a little higher and further inland in FAMOUS-ice. This reflects the tendency of FAMOUS to both melt too far up the ice sheet and to distribute precipitation across a broad area rather than concentrate it on 15 the slopes at the margins of the ice. In the south east and north west the FAMOUS equilibrium line is close to that of MAR but it differs significantly in other areas.
It is further encouraging that FAMOUS-ice can simulate the magnitude and the basic spatial distribution of internal refreeze of melt within the snow pack ( figure 6). Anomalies in this field reflect the biases seen in the SMB field, with refreezing too far   Fettweis et al. (2013) (using an earlier version of RACMO than that plotted in figures 2 and 4). SMB is the total change in snow pack mass over a year. Ablation is the difference between snowfall and SMB, so includes all processes that remove mass from the snow pack such as sublimation. Since FAMOUS-ice does not route rainfall through the snow, refreeze has been defined under the assumption that all of the rain falling on the snow has run off. The ± uncertainty shown for SMB is the standard error of the time-mean, estimated by assuming annual values to be independent. inland in both the south and the north of GrIS, reflecting the fact that FAMOUS simulates summer melt over too much of the ice sheet in the present day, rather than confining melt to the coasts.
Integrated over the ice sheet as a whole, SMB and its components simulated for the present day are in reasonable agreement with RCMs (table 1), although there is a high bias in melt that produces a lower overall SMB than in the RCMs. FAMOUS-ice reports smaller annual variability in SMB than the RCMs as the climatological SSTs used in FAMOUS-ice exclude interannual Our primary motivation for the model development work we have described is to enable FAMOUS-ice to carry out centennial coupled climate-ice sheet simulations. As in previous studies with FAMOUS (Gregory et al. (2012); Roberts et al. (2014), we use the Glimmer ISM (Rutt et al., 2009). More recent versions of Glimmer (Glimmer-CISM, ), are part of the ice sheet model within the CESM (Danabasoglu et al., 2020). For computational efficiency we use Glimmer in its 5 original shallow-ice formulation without sliding and a fixed internal temperature profile, although the climate model coupling itself does not preclude the use of more sophisticated ice sheet or solid earth physics. In this section we concentrate on technical aspects of coupling FAMOUS-ice and Glimmer directly. The results of a suite of simulations exploring the future evolution of the climate and GrIS in this coupled system are described in Gregory et al. (2020, submitted manuscript).
For coupling, we use our modified JULES multi-layer snow scheme to represent the surface and upper levels of the ice sheet, 10 down to a depth beyond which seasonal climate variations are not important. Thus the top of the ice sheet is considered to be within the domain of FAMOUS, and the dynamic ice underneath is handled by Glimmer. Transfer of information between FAMOUS and Glimmer consists of annual average fields, passed once a year between the bottom of the FAMOUS snow pack at the interface with Glimmer. This proceeds in a series of steps, as follows.

15
SMB at each elevation is diagnosed as the change in snow pack mass in FAMOUS-ice over the course of a year. The absolute mass of the snow pack is irrelevant in this context, as long as it is deep enough to not become exhausted over the year.
At each annual coupling step, once the SMB has been calculated for Glimmer, the FAMOUS snow pack mass is reset to its initial value. This procedure in effect passes mass from FAMOUS to Glimmer or vice versa. An accumulation of snow pack mass in FAMOUS indicates a positive mass tendency for the ice sheet, a loss of snow pack mass (i.e. ablation) a negative 20 tendency. The situation where there is not enough mass available in Glimmer to refill the FAMOUS snow pack (i.e the ice sheet has retreated from a location) will be addressed later. In this way, mass is transferred between the models, and the tiles representing the ice sheet in FAMOUS-ice will always have enough snow present to diagnose an SMB term for Glimmer during the next year of simulation. The amount chosen for the ice sheet snow pack is 100 m of liquid water equivalent (LWE). In the coupled model, we initialise the snow pack on ice sheet tiles as ice-density snow. During the model spin-up, fresh snow builds 25 up on the surface in the accumulation zone and mass is removed by the coupling at the base of the snow pack, while in the ablation zone ice melts at the surface and the coupling adds mass from Glimmer at the base of the snow pack.

Interpolation of SMB to the Glimmer grid
The SMB field on the elevation classes received by Glimmer is transferred onto the Glimmer topography using the same trilinear interpolation as in Vizcaíno et al. (2013). The full SMB field for each elevation class is first bilinearly interpolated 30 onto to horizontal locations of the Glimmer grid. Each point on the Glimmer grid then interpolates vertically between the grids of the elevation classes above and below its topographical height. If the Glimmer topography is below the mid-height of the lowest FAMOUS elevation class, or above the mid-height of the highest elevation class, then no vertical interpolation can be done and the value from the relevant elevation class is used unaltered. The same interpolation procedure is applied to the temperature field at the bottom of the FAMOUS snow packs, to provide a surface temperature boundary condition for Glimmer.
Linear interpolation is not generally conservative, and the interpolated output can be globally scaled to ensure that the total mass given to Glimmer matches the mass change in FAMOUS-ice, although this has not been done in the illustrative 5 simulations here. Non-conservation may also arise from the area of the ice sheet being differently represented in each model e.g. due to mismatches of coastline, in which case imposing a numerical scaling on the SMB seen by the ice sheet model to make it conform to the SMB in FAMOUS-ice may not provide the best estimate of the surface forcing that the climate would really produce for an ice sheet of that precise shape. A local adjustment to conserve the mass change of the Glimmer points located in each FAMOUS-ice grid-box was also tested, but was found to heavily imprint the outline of the 5 • x7.5 • FAMOUS-10 ice grid on the 20km Glimmer SMB field. We thus do not use this correction in FAMOUS-ice, but more acceptable results would be expected if such a technique were used with a less-coarse grid climate model as a source. In practice, in the example simulations shown above, the non-conservation due to regridding between the models was less than 2% of integrated GrIS SMB taken over a 10 year mean.
Using these boundary conditions from FAMOUS-ice, Glimmer simulates the evolution of the ice sheet state. For synchronous 15 coupling with the climate, Glimmer simulates a year of ice flow before transferring information back to FAMOUS. However, under the assumptions that the change in ice sheet geometry over N years, where N is a small number, is too small to have a significant effect on the climate and SMB calculated by FAMOUS-ice, and that the global climate is either constant or changes negligibly within N years, it is also possible to use an asynchronous coupling scheme where the ice is allowed to evolve for N years under the same climate boundary conditions. In our experience, running N =10 ice sheet years for every year of climate 20 simulated does not significantly affect the evolution of the coupled climate-ice system, but running N=100:1 does result in differences compared to a synchronously coupled run.

Update of FAMOUS following Glimmer
After timestepping the ice sheet for the required number of years, Glimmer passes several fields to FAMOUS: the mean orographic height for each FAMOUS-ice grid-box; the fractions of ice and non-ice covered areas in each elevation class to 25 define the tile fractions in FAMOUS-ice, and some fields describing the sub-gridscale distribution of the orography within each grid-box for use by parameterisations of atmospheric gravity-wave and boundary layer drag. All these are calculated from the updated topography on the Glimmer grid.
Glimmer not only changes the height of the ice sheet, but can also produce ice in formerly unglaciated areas or make areas ice-free if the SMB provided is sufficiently negative. All of these effects are represented in FAMOUS-ice through changes in 30 the area fractions of the sub-gridscale tiles, which in turn affect the grid-box mean surface properties. For structural reasons within FAMOUS, it is not possible to change the land-sea mask of the climate model as it runs. Glimmer may thus only move the ice within the boundaries of FAMOUS-ice's pre-defined Greenland land area, and our coupled system is unable to simulate the growth of floating shelves, displacement of ocean by advancing grounded ice or changes in sea-level. Any ice that Glimmer moves beyond the coastline defined by FAMOUS-ice is calved to the ocean.

Creating or destroying glaciated grid-boxes
An unglaciated grid-box in Glimmer may be occupied by ice in one of two ways. Glimmer may move ice into it, or it may accumulate a sufficient depth of meteoric snow to pass into the regime where we would want Glimmer to be able to treat it 5 dynamically. To allow this latter process, information about the snow mass on unglaciated tiles in FAMOUS is also passed to Glimmer and interpolated in the same trilinear fashion as the SMB. Where the resultant field has a snow depth greater than a given threshold on a grid-box that is not already occupied by ice in Glimmer then a new ice point is initialised. To conserve mass, the unglaciated tile in FAMOUS-ice that sourced this snow has an equivalent mass subtracted from it.
Conversely, as the ice in Glimmer thins then we need a method for deglaciating grid-boxes. When the ice becomes so thin 10 it is no longer dynamically active and potentially unable to fully resupply the FAMOUS snow pack at the start of a coupling period we remove all ice from that grid-box in Glimmer and increase the area fraction of the corresponding unglaciated tile in FAMOUS-ice accordingly. The mass of ice removed from Glimmer is also added to the snow pack in the unglaciated tile in FAMOUS-ice. The threshold for creation or removal of ice from Glimmer is set at 50m liquid water equivalent in the example simulation used here. 15 Since snow in FAMOUS is held as a mass-per-area term, a change to the area fraction of a tile will result in a change in the absolute mass of snow on that tile, regardless of any additional adjustments that may also be required by the processes described above. Full conservation of the multi-layered, multi-prognostic snow pack properties as the ice sheet area waxes or wanes is a complex operation. There are many science questions that do not require a model to be able to initiate dynamically active ice in new areas purely from meteoric snow, nor which require water mass to be conserved exactly across the system 20 components. In such cases we allow ice to spread to new areas only by Glimmer ice sheet dynamics, and diagnose the degree of non-conservation implied by not adjusting the snow packs to account for changes in area fraction.

A simulation of the present-day Greenland Ice Sheet in FAMOUS-ice
It has been shown that the SMB calculation of FAMOUS-ice for GrIS in a modern day climate compares reasonably well to other models. Once coupled to a ice sheet model however, biases in the simulated SMB may result in an unacceptable 25 simulation of the dynamics and geometry of the ice.
In a run of 5000 years of ice sheet evolution (500 years of climate simulation and an asynchronous time factor of N=10) of FAMOUS-ice coupled to Glimmer forced by the same modern greenhouse gases and sea surface conditions described in section 7, GrIS volume and area grow to reach a steady state about 10% larger than observed (figures 7, 8). The resolution and shallow-ice numerics of Glimmer mean that individual outlet streams and Greenland fjords are not well reproduced, and 30 the requirement that ice only calves when it reaches the coast of the climate model means that the ice sheet expands toward the coast. As it does so it slumps in the centre and builds up higher shoulders at the margins to allow enough ice to calve to maintain balance with the climate model SMB ( figure 7, red lines, figure 8a). When calving is instead enforced at the current Figure 7. Spinup of the Greenland ice sheet in the coupled FAMOUS ice. ice sheet volume (x10 x km 3 LWE, solid lines, left axis) and area (x10 6 km 2 , dashed lines, right axis). Black: calving imposed near the observed ice edge; red: ice sheet allowed to calve at the coast.
Asynchronous coupling was used: the time-axis shows climate model years, the ice sheet evolves 10 years for every 1 year of climate. ice-edge, making up for the lack of explicit fjord calving at this resolution of Glimmer, the ice sheet does not spread in this way and the volume is also more realistic ( figure 7, black lines, figure 8b).
More details of the simulation of GrIS in coupled FAMOUS-ice and its sensitivity to climate change can be found in Gregory et al (2020, submitted ms).
We have described modifications to the FAMOUS climate model that enable it to calculate a realistic surface mass balance for ice sheets, and its coupling to an model of dynamic ice sheet flow. Together, these developments allow coupled simulations of the multi-millennial evolution of climate and grounded ice sheets to be carried out with more fidelity to the physical processes in reality. 5 FAMOUS-ice has been used for studies of the future evolution of GrIS, including an analysis of the ice sheet's stability and the presence of tipping points (Gregory et al. (2020, submitted manuscript). This same framework will be used for paleoclimate studies, where the co-evolution of and feedbacks between the climate and ice sheets are essential to understanding the glacial cycles of the last millennia. A necessary future development is the ability to simulate sea-level rise and the spread of ice onto coastal shelves, which is currently not possible as the land-sea mask in FAMOUS is fixed. 10 The configuration of FAMOUS-ice described in this paper is limited to simulating ice sheets that are predominantly grounded and not subject to direct influence from contact with the ocean. The shallow ice approximation used here in Glimmer is also not suited to simulating streaming ice. Techniques for coupling models of floating ice shelves and marine-terminating glaciers to the ocean is a subject of current research (e.g. Asay-Davis et al. (2016)), and it is not yet clear how this should best be done in any modelling system, let alone at coarse spatial resolution and over centennial or millennial timescales. Therefore, whilst 15 FAMOUS-ice can justifiably be applied to Greenland and the Laurentide ice sheet of the Quaternary glacial cycles, it cannot simply be used to simulate Antarctica or other paleo ice sheets which are predominantly affected by interactions with the ocean. It would be possible to use a parameterisation of the influence of ocean temperature on neighbouring ice sheets, such as in Favier et al. (2019), and this will be pursued in future work. Elevation tiles are also not as effective for downscaling a coarseresolution climate simulation on Antarctica, where surface conditions are less tightly dependent on altitude than on Greenland, 20 and more determined by small-scale dynamic processes such as katabatic winds, blowing snow and synoptic variability. In this case accumulation needs to be redistributed both vertically between tiles and horizontally between grid-boxes.
The code for the coupled system we have developed for FAMOUS-ice may be used with some other models. HadCM3 (Gordon et al., 2000) and the HadRM3 and PRECIS regional models (Jones et al., 2004) rely on the same UM4.5 code-base as FAMOUS, and it would be possible to directly transfer the adaptations to the tiles and the multi-layer snow scheme to 25 these models. HadCM3 has a better climate simulation than FAMOUS and is still used extensively in the UK climate research community, although it is more computationally expensive. On the ice side, the FAMOUS-ice coupling to Glimmer could alternatively accommodate the BISICLES ice sheet model (Cornford et al., 2013). BISICLES is a variable resolution, adaptive mesh ISM whose numerics allow it to simulate ice streams, floating ice shelves and grounding line retreat in a computationally affordable manner. Using BISICLES would remove the restrictions inherent in Glimmer's shallow-ice approximation. in a manner very similar to that described here. Although the climate simulation of UKESM1 contains considerably more physical detail than FAMOUS-ice, UKESM1's computational cost makes it impossible to conduct the multi-millennial climate simulations that are needed to resolve many questions of the co-evolution of climate and ice sheets. Simpler models such as FAMOUS are still needed for many of the scientific questions we urgently need to answer about how the Earth system functions. The mods required to configure the PUM as FAMOUS-ice can be downloaded from https://puma.nerc.ac.uk/svn/UniCiCles_svn/ UniCiCles/tags/FAMOUS-ice_SEG, along with a fork of version 1.9 of Glimmer-CISM it can couple with. Namelists and configuration details that specify the FAMOUS-ice simulation illustrated in this work are available via job xotzb in the central UMUI service provided by NCAS CMS at puma.nerc.ac.uk. Additional support for obtaining and running the PUM and 15 FAMOUS-ice simulations in practice is available via the NCAS CMS helpdesk at http://cms.ncas.ac.uk/wiki/CmsHelpdesk or via the authors.

Data availability
Simulation output from the FAMOUS-ice run used in the illustrative plots can be downloaded from http://gws-access.jasmin.ac.uk/ public/ncas_climate/rssmith/sgfjb 20

Author contributions
RS and JG designed the coupling schemes, which were coded in FAMOUS and Glimmer by RS and SG. SG improved the model and conducted the simulations. All authors contributed to the simulation analysis. The manuscript was written by RS.