Interactive comment on “ A fully coupled 3-D icesheet – sea-level model : algorithm and applications ” by B

Abstract. Relative sea-level variations during the late Pleistocene can only be reconstructed with the knowledge of ice-sheet history. On the other hand, the knowledge of regional and global relative sea-level variations is necessary to learn about the changes in ice volume. Overcoming this problem of circularity demands a fully coupled system where ice sheets and sea level vary consistently in space and time and dynamically affect each other. Here we present results for the past 410 000 years (410 kyr) from the coupling of a set of 3-D ice-sheet-shelf models to a global sea-level model, which is based on the solution of the gravitationally self-consistent sea-level equation. The sea-level model incorporates the glacial isostatic adjustment feedbacks for a Maxwell viscoelastic and rotating Earth model with coastal migration. Ice volume is computed with four 3-D ice-sheet-shelf models for North America, Eurasia, Greenland and Antarctica. Using an inverse approach, ice volume and temperature are derived from a benthic δ18O stacked record. The derived surface-air temperature anomaly is added to the present-day climatology to simulate glacial–interglacial changes in temperature and hence ice volume. The ice-sheet thickness variations are then forwarded to the sea-level model to compute the bedrock deformation, the change in sea-surface height and thus the relative sea-level change. The latter is then forwarded to the ice-sheet models. To quantify the impact of relative sea-level variations on ice-volume evolution, we have performed coupled and uncoupled simulations. The largest differences of ice-sheet thickness change occur at the edges of the ice sheets, where relative sea-level change significantly departs from the ocean-averaged sea-level variations.


Introduction
Periodical expansion and retreat of continental ice sheets has been the main driver of global sea-level fluctuations during the Pleistocene (Fairbridge, 1961).Similarly, deep-sea benthic δ 18 O records, a proxy for deep-water temperature and ice volume, indicate that the volume of the oceans oscillated throughout the Pleistocene in response to global climate changes (Chappell and Shackleton, 1986;Lisiecki and Raymo, 2005).The separation of the benthic δ 18 O signal into deep-water temperature and ice volume can be deduced by using a combination of ice-sheet models and an air-to-ocean temperature coupling function (de Boer et al., 2013).However, the exact contribution of the different ice sheets to the spatially varying relative sea level (RSL), i.e. the change of the sea surface relative to the solid Earth, is unknown.
One of the best studied intervals in the past is the last glacial maximum (LGM, ∼21.0 kyr ago), for which a wealth of observational data has been collected, for example RSL and extent of the ice sheet.The LGM was a glacial event during which continental ice sheets covered large portions of North America and Eurasia, and when the ice sheets on Antarctica and Greenland extended towards the continental shelf edge (e.g.Ehlers and Gibbard, 2007;Denton et al., 2010).Several well-dated surface geological features of depositional and/or erosive origin constrain the maximum extent of these LGM ice sheets (Ehlers and Gibbard, 2007).The estimated total volume of ice inferred from the RSL data correlates well with the ice-sheet volume increase inferred from the benthic δ 18 O data (Lisiecki and Raymo, 2005).Both B. de Boer et al.: Global coupled 3-D ice-sheet-sea-level model benthic δ 18 O records and surface glacial geological features show the −120 to −130 m relative sea-level low stand that was recorded by submerged fossil coral terraces at Barbados (Peltier and Fairbanks, 2006;Austermann et al., 2013), Tahiti (Bard et al., 1996(Bard et al., , 2010;;Deschamps et al., 2012) and Bonaparte Gulf (Yokoyama et al., 2000).The sea-level rise recorded at these far-field sites during the last ∼ 19 kyr following the melting of the LGM ice sheets is consistent with a decrease of benthic δ 18 O.This marks the transition to the current warmer interglacial (Fairbridge, 1961).
However, several coeval post LGM paleo-sea-level indicators from different regions are found at present at very different elevations above and below the current mean sea level (Pirazzoli, 1991).In particular, a long-term sea-level fall is observed in the proximity of the former LGM ice sheets (Lambeck et al., 1990).Moving slightly away from the formerly glaciated areas, the sea-level trend first switches towards a steep rise (Engelhart et al., 2011), reaching a mid-Holocene high stand (Basset et al., 2005) and then smoothly changes towards a eustatic-like sea-level fall that is observed at the far-field sites like Barbados and Tahiti (Fairbanks, 1989;Peltier and Fairbanks, 2006;Bard et al., 1996Bard et al., , 2010;;Deschamps et al., 2012).This illustrates that the regionally varying sea-level changes resulting from the melting of the LGM ice sheets shows that the spatial variability of sea-level change strongly depends on the distance from the former ice sheets and also on the shape of the ocean basins (Pirazzoli, 1991;Milne and Mitrovica, 2008).
Following the deglaciation of an ice sheet, the solid Earth rebounds upwards beneath the former glaciated area, while the far-field ocean basins experience subsidence as a consequence of the increasing ocean water loading.Therefore, if the ice-sheet thickness variation is the forcing function for the sea-level change, the solid Earth response plays an important role as a modulator of sea-level change.Global sealevel changes during the Holocene and in particular the last 6000 years (Pirazzoli, 1991) show that the solid Earth continued to deform after the North American ice sheet (NaIS) and the Eurasian ice sheet (EuIS) were completely melted.This delayed response implies that the solid Earth behaves like a highly viscous fluid on geological timescales (Ranalli, 1985;Turcotte and Schubert, 2002).Additionally, the current vertical land uplift shown by GPS observations over the formerly glaciated areas of Scandinavia and Hudson Bay also implies that the solid Earth is not in isostatic equilibrium (Milne et al., 2001).
The mean sea surface is also affected by the gravitational pull exerted by the continental ice sheets on the ocean water.During the melting of an ice sheet, the ocean volume and thus the hypothetical eustatic sea level, i.e. the global mean change in sea level, are increasing.However, due to the smaller ice mass there is a reduction in the gravitational pull exerted on the ocean water, which causes a sea-level fall at the ice-sheet margins and a rise, more than it would do eustatically, far away from the ice sheet.This effect is known as self-gravitation (Woodward, 1888), and combined with the solid Earth deformation it attributes a large proportion of the spatial variability of the sea-level change (Farrell and Clark, 1976).Furthermore, due to the rotation of the Earth around its axis, any surface mass displacement together with the solid Earth and geoidal deformations triggers a perturbation of the polar motion that in turn affects the redistribution of melt water in the oceans and ,hence, the mean sea-surface height (Milne and Mitrovica, 1996;Kendall et al., 2005).
All the feedbacks described above make up the complex process known as glacial-hydro isostatic adjustment (GIA), which includes deformation of the Earth and changes of the geoid, and describe any sea-level change that is dictated by ice-sheet fluctuations (Farrell and Clark, 1976).According to the GIA theory, the sea-level change recorded at any location represents the combined response of the solid Earth and of the geoid to the ice-sheet fluctuations; it cannot be directly used as representative of the eustatic sea-level change.GIA feedbacks produce mutual motion of the solid Earth and of the geoid, and hence any land-based sea-level indicator is essentially a RSL indicator as it records the local variation in the vertical distance between the geoid and the bottom of the ocean.
The GIA feedbacks are usually accounted for by solving the gravitationally self-consistent sea-level equation (SLE), which was initially developed by Farrell and Clark (1976), and subsequently updated to include all the GIA feedbacks (Mitrovica and Peltier, 1991).The SLE describes the global RSL change for a prescribed ice-sheet chronology and solid Earth rheological model (Spada and Stocchi, 2007).The SLE has been widely employed to improve and refine our knowledge of the LGM ice-sheet's volume, thickness and extent and their subsequent deglaciation until the present day (e.g.Peltier, 2004;Spada and Stocchi, 2007;Whitehouse et al., 2012a;Stocchi et al., 2013;Gomez et al., 2013).
To explain the observed RSL changes over the past glacial cycles, a global ice-sheet model is needed to calibrate the corresponding ice volume.At the same time, the observed RSL changes are needed to verify the simulated ice volume with the global ice-sheet models.This problem of circularity follows from the fact that the evolution of the ice sheets is coupled to the RSL changes.Most importantly the ice-sheetinduced RSL changes affect the growth and retreat of marine ice sheets, which are in direct contact with the ocean.For example, due to the self-gravitational pull of the ice sheet, the RSL close to the ice sheets actually rises when an ice sheet grows, this will then counteract the advance of the (marine) ice sheet into the ocean.
Thus far most transient simulations of ice sheets have been carried out using a global average sea level (Huybrechts, 2002;Zweck and Huybrechts, 2005;Bintanja and Van de Wal, 2008;Pollard andDeConto, 2009, 2012;de Boer et al., 2013).There have been no studies that simulate a mutually consistent solution of ice volume and regional sea level that include multiple ice-sheet models over longer timescales (Clark et al., 2002;Weber et al., 2011;Raymo and Mitrovica, 2012).During the mid-1970s the importance of including the effect of relative sea-level change on the instability of marine terminating ice sheet was recognised (Weertman, 1974;Farrell and Clark, 1976).This is important because sea-level change has a strong influence on the dynamical behaviour of marine ice sheets (e.g.Gomez et al., 2010aGomez et al., , 2013) ) such as the West Antarctic ice sheet (Pollard and DeConto, 2009).Only recently, Gomez et al. (2013) have succeeded in coupling a single 3-D ice sheet with a sea-level model for simulating the Antarctic ice sheet from the LGM to the present.Henceforth, it is of vital importance to incorporate regional sea-level variations when modelling ice sheets over 10 3 -10 6 years.
In this paper, we present a system of four regional 3-D ice-sheet-shelf models (de Boer et al., 2013) that is fully dynamically coupled to a global GIA model based on the SLE.Where Gomez et al. (2013) employed a similar system for the Antarctic ice sheet, our algorithm represents a method for modelling ice-sheet fluctuations and the related GIA-induced relative sea-level changes on a global scale.In addition, it is dynamically coupled to a deformable Earth model where crustal and geoidal deformations account for self-gravitation, Earth rotation and an adequate treatment of the migration of coastlines.Our new model offers the opportunity to model any ice-sheet and sea-level fluctuation, from the past to the present day as well as into the future.Here, we include a temporal discretisation of past ice-sheet fluctuations with a moving time window that allows us to calculate RSL as a function of the total ice-sheet volume change over the globe over four glacial cycles, starting 410 kyr ago.This allows a comparison with any local record of RSL during this period.

Methods
In this study we present a new system that is based on the dynamical coupling between (i) ANICE, a fully coupled system of four 3-D regional ice-sheet-shelf models (de Boer et al., 2013) and (ii) SELEN, a global scale SLE model that accounts for all the GIA feedbacks (Spada and Stocchi, 2007).In the following, we first describe separately the ANICE and SELEN sub-systems and subsequently introduce the coupling method/algorithm with particular emphasis on spatial (see Appendix A) and temporal discretisation.

The ANICE regional ice-sheet-shelf model
ANICE is a 3-D coupled ice-sheet-shelf model (Bintanja and Van de Wal, 2008;de Boer et al., 2013de Boer et al., , 2014)).It is a shallowice model, for which we use approximate equations for sheet and shelf flow.These approximations are based on the shallowness of a large ice body, with horizontal scales far larger than the thickness of the ice.In ANICE, we apply two approximations, the shallow-ice approximation (SIA) (Hutter, 1983) that is used as the basis for land based ice flow, and the shallow-shelf approximation (SSA) (Morland, 1987) that is used for floating ice and sliding velocities (de Boer et al., 2013).The latter is computed for both grounded and floating ice; thus, incorporating the transition zone from sheet to shelf (Winkelmann et al., 2011;de Boer et al., 2013).Within this framework we incorporate four separate icesheet models for the regions with major ice sheets during the Pleistocene: the Antarctic ice sheet (AIS), the Greenland ice sheet (GrIS), the NaIS and the EuIS.The models are solved separately on a rectangular x-y grid (see Table 1).Ice temperatures and velocities are solved in three dimensions with 15 grid points in the vertical, which are scaled with ice thickness and have a higher resolution at the bottom, starting with 1 % and increasing to 10 % at the top.
We adopt the initial basal and surface topographies for Antarctica from the ALBMAP data set (Le Brocq et al., 2010) and for Greenland from Bamber et al. (2001).The initial topography for Eurasia and North America is based on a high-resolution present-day (PD) topography data set (SRTM30_PLUS; Becker et al., 2009).For the initial climate forcing, we use the PD meteorological conditions from the ERA-40 Re-analysis data set (Uppala et al., 2005).We calculate monthly averages from 1971 to 2000 for precipitation (in m w.e.yr −1 ), 2 m surface-air temperature ( • C), and 850 hPa wind fields (in m s −1 ).The surface topography for the EuIS and NaIS and the ERA-40 climate fields are interpolated on the rectangular ANICE grids with an oblique stereographic projection using the OBLIMAP programme (Reerink et al., 2010).The AIS, EuIS and NaIS models incorporate grounded and floating ice and a sub-shelf melting parameterisation, whereas for the GrIS we only consider ice on land (de Boer et al., 2013).All four models are solved with their own internal time step varying between 1 and 5 years, depending on the stability criterion, i.e. the ice velocity cannot exceed the grid scale (Table 1) divided by the time step.
The uncoupled ice-sheet model ANICE accounts for the regional bedrock deformation that is calculated from variations in ice thickness and ocean water by means of a two layer flexural Earth model, a flat elastic lithosphere (EL) resting over a viscous relaxed asthenosphere (RA), i.e. the ELRA model (Le Meur and Huybrechts, 1996).The upper layer mimics the elastic lithosphere and therefore accounts for the shape of the deformation.The time response of the bedrock deformation is controlled by the lower viscous asthenosphere, with a constant response time of 3 kyr.The rate of the vertical bedrock movement is proportional to the B. de Boer et al.: Global coupled 3-D ice-sheet-sea-level model deviation of the profile from the initial equilibrium state and inversely proportional to the relaxation time (Le Meur and Huybrechts, 1996;de Boer et al., 2013).Within the uncoupled ANICE system, eustatic sea-level change is internally calculated from ice-volume changes relative to PD (de Boer et al., 2013(de Boer et al., , 2014)).

Model forcing
To simulate the evolution of the ice volume through time, we use benthic oxygen isotope δ 18 O data as an input, which is a proxy for changes in ice volume and deep-water temperature (Chappell and Shackleton, 1986).Here, we used the LR04 benthic δ 18 O stack of 57 deep-sea sediment records (Lisiecki and Raymo, 2005) with an inverse procedure to separate the benthic δ 18 O data into an ice volume and deep-water temperature component (de Boer et al., 2013).Since this data set uses globally distributed records of benthic δ 18 O data, we assume the record represents a global average climate signal (de Boer et al., 2013).From the benthic δ 18 O data, a surface-air temperature anomaly relative to PD is derived using an inverse procedure (Bintanja and Van de Wal, 2008;de Boer et al., 2013de Boer et al., , 2014)).The method is based on the assumption that both ice volume and deep-water temperature are strongly related to the mid-latitude-to-subpolar Northern Hemisphere (NH) surface-air temperature.This continental mean (40 to 80 • N) temperature anomaly (hereafter NH temperature anomaly) controls the waxing and waning of the EuIS and NaIS (Bintanja et al., 2005).The procedure linearly relates the NH temperature anomaly to the difference between the modelled and observed benthic δ 18 O 100 years later given by Here, 100 years is the time resolution of the δ 18 O forcing record, T NH is the mean NH temperature anomaly over the preceding 2000 years (2 kyr) and the second term on the right-hand side represents the temperature response to changes in the δ 18 O record.The modelled benthic δ 18 O is calculated using ice volume, ice-sheet δ 18 O and deep-water temperatures relative to PD for every 100 years.The length of the mean window of 2 kyr and the scaling parameter of 20 were optimised by minimising the difference between modelled and observed δ 18 O, i.e. the observed δ 18 O record must be accurately followed (de Boer et al., 2013).The calculated NH temperature anomaly is forwarded to the ice-sheet models and uniformly added to the surface temperature field.Within each model the surface temperatures are also corrected for surface height changes with a temperature lapse rate.As a result, the model computes ice volume and temperature consistent with the benthic δ 18 O forcing (Fig. 1).For a full description of ANICE, see de Boer et al. (2013).

The SELEN global SLE model
SELEN solves the SLE by using the pseudo-spectral method (Mitrovica and Peltier, 1991;Spada and Stocchi, 2007;Stocchi et al., 2013) and calculates the deformations of the solid Earth (U ), the geoid (N; mean sea surface at rest) and relative sea level (S) as a function of time on a global scale: The two fundamental inputs for the SLE are (i) the icesheet thickness chronology that represents the forcing function for sea-level change, and (ii) the solid Earth rheological model that describes the response of the solid Earth and the geoid to the melt-water redistribution.The Earth is assumed to be spherically symmetric, self-gravitating and radially stratified.
The rheology model only accounts for radial variability.The outer shell is assumed to be perfectly elastic and mimics the lithosphere.The mantle is discretised into n Maxwell viscoelastic layers (linear rheology) while the inner core is assumed to be inviscid.Our default settings for the coupled ANICE-SELEN system are an elastic lithosphere thickness of 100 km and a n = 3 layer Earth model with a viscosity for the shallow upper layer of 3 × 10 20 Pa s, a transition zone of 6 × 10 20 Pa s and a lower layer of 3 × 10 21 Pa s.We adopt the normal mode technique to generate the response of the Earth to variations of land ice and water loading (Peltier, 1974).
The spatio-temporal variations of ice-sheet thickness represent the a priori forcing function that drives the corresponding self-consistent RSL changes.At the core of the SLE is the concept that any local RSL change depends upon all surface mass displacements (both ice and melt water) that have occurred since the beginning of the ice-sheet  (de Boer et al., 2013).Every 1000 years ANICE forwards grounded ice thickness, the Ice load given in the input array (IA), to SELEN, which computes the gravitationally self-consistent sea level and bedrock topography adjustment that are coupled back to ANICE; in terms of the RSL, S is given in the output array (OA).chronology anywhere on the Earth.Recent improvements account for the dynamical feedback from the solid Earth rotation and the lateral migration of coastlines, also known as the time-dependent ocean function (Milne and Mitrovica, 1996;Mitrovica and Milne, 2003;Kendall et al., 2005).We solve the SLE with a pseudo-spectral numerical scheme (Spada and Stocchi, 2007) that we truncate at a spherical harmonic degree of order 128 to save computation time.Moreover, the SLE is solved by means of an iterative procedure where, at the first iteration, the RSL change S is assumed to be eustatic.After 3 iterations, the solution has converged and S is regionally varying (non-eustatic, non-globally uniform) according to GIA feedback (Farrell and Clark, 1976;Mitrovica and Peltier, 1991;Spada and Stocchi, 2007),

The fully coupled system of ANICE-SELEN
In the following, we describe the dynamical interaction between the four regional 3-D ice-sheet-shelf models, which define the ANICE sub-system (see Sect. 2.1), and the gravitationally self-consistent SLE, which is solved by means of SELEN (see Sect. 2.3) as illustrated in Fig. 2. In the coupled ANICE-SELEN system, the RSL change that is provided to ANICE includes bedrock deformation and changes in the sea surface and thus replaces the regional ELRA model that is used for the stand-alone ANICE simulations.According to the SLE, solid Earth and geoid deformation at each point in space and time linearly depend on all the ice-sheet thickness variations and on the corresponding changes in the ocean loading that have occurred until that time.Having ANICE-SELEN fully and dynamically coupled implies that information is exchanged between the two sub-systems through time.ANICE provides SELEN with ice-sheet thickness variation in space and time, while SELEN returns the corresponding RSL change (representing both variations in U and N , see Eq. 2) to ANICE.The two means of communication between ANICE and SELEN are the input array IA(λ, θ ), which carries information about ice-sheet thickness variation in space, and the output array OA(λ, θ ), which retrieves the RSL change at each element of the four ANICE sub-domains.Both are a function of latitude (λ) and longitude (θ).The output array, containing RSL change, is used within ANICE to update the topography for the next time step.This procedure is repeated with a coupling interval t C = 1 kyr (Fig. 2 and Table 2).Before the coupling starts, AN-ICE is spun up for 1 glacial cycle in the uncoupled mode without SELEN.In the uncoupled ANICE sub-system, each regional ice-sheet model deforms its own regional topography independently from the other three ice sheets.Together, the four regional ice-sheet models contribute and respond to the global eustatic sea-level change.The latter is internally calculated from the changes in ice volume and is the only means of connection among the four ice sheets.When the coupling starts at 410 kyr ago, the ELRA model is switched off and all four regions use the spatially varying RSL as provided by SELEN, which implicitly includes the deformation of the Earth.

Spatial discretisation
The execution of the algorithm starts with the discretisation of the Earth surface into almost equal-area hexagonal elements.The number of hexagons, i.e. the spatial resolution of the global mesh used within SELEN, depends on the parameter RES (Spada and Stocchi, 2007) (see Appendix A).The time steps of the moving time window 10 × 1, 2 × 5, 3 × 20 kyr A RES value of 60 was adopted, which results in 141 612 hexagonal elements, and each element approximately corresponds to a disc with a half amplitude of α = 0.304 angular degrees (see Appendix A).We employ the surface interpolation routine grdtrack from GMT (Wessel and Smith, 1991) to project ETOPO1 topography on the global mesh (Amante and Eakins, 2009).For each element the grdtrack routine provides a value for the bedrock topography as well as a value for the ice elevation that is non-zero wherever ice is currently present.Wherever the bedrock height is negative and the ice elevation is non-zero, we evaluate whether the ice is grounded or floating.This is essential for defining the ocean function (OF), which describes if an element belongs to the ocean (OF = 1) or to the land (OF = 0) (Milne et al., 2002).
Once the initial global topography file is generated, we update this field by projecting the four initial ANICE topographies and ice-sheet thickness on the SELEN grid.However, the three Northern Hemisphere regional ice-sheets models (NaIS, EuIS and the GrIS) share overlapping regions.We therefore define a hierarchical procedure where the topography and ice-sheet thickness values to be interpolated on the elements of the global mesh are, firstly, those from the NaIS, secondly, the EuIS and, finally, the GrIS and AIS (see Supplement).The ANICE grid points and SELEN ice elements are shown in Fig. 3 and the specific number of x and y grid points and SELEN elements of each ice-sheet model grid are provided in Table 1.
The geographical coordinates of the elements that are initially updated with the four separate ANICE topographies and ice-sheet thickness are stored in the input array.These are the elements that could potentially be affected by ice-thickness variation through time, and consequently are recognised by SELEN as ice-sheet elements (Fig. 3b, d).The ice-sheet thickness is initially zero in ice free areas and nonzero wherever there is currently grounded ice, i.e. on Greenland and Antarctica.This initial array is the projection of topography and ice thickness of the four ANICE sub-domains on the global hexagonal mesh and represents an interglacial stage from which all of our simulations start.At each coupling time step t C of the simulation, the array is updated with new ice-sheet thickness values according to ANICE, and the information is passed to SELEN for the computation of the GIA-induced RSL changes.The latter are returned to ANICE by means of the output array that stores the  1. geographical coordinates of the centroids of the equal-area elements of the four ice-sheet regions.

Temporal discretisation
In SELEN, the temporal discretisation is performed assuming that the variables vary stepwise in time (Spada and Stocchi, 2007).Usually, the late Pleistocene ice-sheet time histories that are available from literature (e.g.Peltier, 2004) are discretised into time steps of 500 or 1000 years.Provided that the solid Earth behaves like a Maxwell viscoelastic body (see Sect. 2.3), the RSL change induced by the ice thickness variation between two consecutive times accounts for (i) an immediate elastic part that occurs as soon as the second icesheet thickness is loaded and (ii) a viscous part that depends on the mantle viscosity profile and on the length of the time step t S , the time step at which the viscous response is discretised.
When using SELEN for a prescribed a priori ice-sheet chronology, the spatio-temporal discretisation of ice-sheet thickness is assimilated at once (Spada and Stocchi, 2007).Consequently, given the time step t S , the total number of time steps and the load Love numbers (Peltier, 1974), the RSL change is computed by means of spatio-temporal convolutions over the surface of the Earth.Accordingly, the change in RSL at any location on the Earth and at any time since the beginning of the ice-sheet chronology is determined by all the ice and ocean load variations that have occurred until that time step (see Sect. 2.3).This implies that, by assuming a predefined mantle viscosity profile, it is possible to compute RSL changes at any time t after the end of the ice-sheet chronology as a consequence of the mantle viscous relaxation, which is an exponentially decaying function of time (Peltier, 1974).
When coupling ANICE-SELEN, a problem arises because the ice thickness variation through time is not known a priori.The ice-sheet thickness variation is only known until the time SELEN is called by ANICE, which is done with an interval of t C = 1 kyr.This implies that any time ANICE calls SE-LEN to compute the bedrock deformation and the sea-surface variation for a specific time t > −410 kyr (the first time that ANICE calls SELEN), all the deformations triggered by the previous time steps are required.Hence, any time t SE-LEN is called, the SLE must be solved starting again from t = −410 kyr (the first ice thickness change).As a consequence, the arrays carrying the SLE results grow throughout the simulations.This is not a big problem when simulating short ice-sheet fluctuations like the post LGM melting, but it is definitely a limitation when simulating multiple glacial cycles.
To avoid this problem, we take advantage of the linearity of both the SLE and of the rheological model.In particular, we use the fact that the viscous response of the bedrock deformation exponentially decays with time and ceases once isostatic equilibrium has been reached.At any time t when ANICE calls SELEN, the bedrock deformation U (t) and the geoid change N (t), due to the ice-thickness change I (t) = H (t) − H (t − t S ), are computed between t = t + t S and a predefined t = L, where L is the total length (in kyr) of a moving time window (see Table 2).Here, H (t) is the ice thickness at time t, and I (t) is the change in ice thickness relative to the previous time step.
We call this temporal discretisation scheme the "moving time window".The length of the moving time window L, i.e. how far into the future SELEN solves for the RSL change, is a free parameter.The longer the moving time window, the more accurate the results will be, because more information from the past is taken into account.In order to maintain a long enough moving time window and to save CPU time, it is important to consider how many time steps NT (Number of time steps of the moving time window) of t S are used to define the moving window.If the length of the moving window allows for a longer memory, the number of time steps allow for an accurate discretisation of the RSL change.
Figure 4 illustrates this process for a 20 m thick ice sheet that is added at time t = 0 (inset of Fig. 4a; using a schematic set-up as shown in Fig. 5).SELEN computes the bedrock deformation from t = 0 to t = L, the length of the moving time window that is set to L = 80 kyr.The bedrock deformation is computed at NT = 15 time points in the future, with NT the number of time steps t S of the moving time window (see Table 2).The time steps are heterogeneous, i.e. 10 steps of 1 kyr, 2 steps of 5 kyr and 3 steps of 20 kyr.The discretisation time step t S is thus an array of length NT = 15.The black squares show the predicted bedrock deformation at each time step.Then, the bedrock deformation is interpolated within the total window of 80 kyr to have a discretised solution at the resolution of the coupling interval t C = 1 kyr (Fig. 4b).At the following time t = t + t C , another 20 m of ice is added above the initial ice layer, and the bedrock deformation due to this extra mass is computed again in the same way.The new array is summed to the previous one to incorporate the viscous deformations of the initial ice-thickness variation (Fig. 4c).This process is carried on throughout the whole simulation so that the memory of previous ice thickness variations is maintained.
Two auxiliary arrays, (auxiliary sea level) AS (λ, θ, t) and (auxiliary ocean function) AOF(λ, θ, t), are generated to store the following L kyr of RSL changes and ocean function variations with a temporary resolution of t C , respectively.The auxiliary arrays are generated using the ice and water loading at time t and are both discretised into NT time steps.At the end of each iterative step of the SLE, the ocean functions are updated using the current computed RSL changes, S(λ, θ, t), and the predicted RSL change as stored in the AS(λ, θ, t) array, that includes all past variations of S. After the SLE is solved, the newly obtained ocean function is then stored in the AOF and the calculated RSL change of the current time is added to the AS (Fig. 4c).This is necessary to account for the variation of coastlines.The output array OA(λ, θ ) that is sent back to ANICE only stores the RSL change for the current time including the past variations, AS(λ, θ, 0).
Throughout a full ANICE-SELEN simulation, the role of the auxiliary arrays, AS(λ, θ, t) and AOF(λ, θ, t), is to account for the response to past ice-sheet thickness variations.This avoids the computationally expensive problem of performing, at any call from ANICE, a full-temporal convolution, since the first-time SELEN is called t = −410 kyr.The auxiliary arrays are consequently updated at any call from ANICE to SELEN, t C = 1 kyr, to store the contributions of each ice-sheet thickness variation simulated with ANICE over a period of L kyr, the length of the moving time window.

Schematic test with the moving time window
As described in Sect.3.2, SELEN is called by ANICE every t C = 1 kyr.The length of the moving window L and the length of the time steps of the moving time window, t S , must be multiples of the coupling interval.We adopt a heterogeneous set of time steps NT to include past variations of GIA.As an example of this algorithm, we use a schematic experiment with an predefined ice load over 480 kyr to demonstrate how the moving time window works (Fig. 6b).We have used an axisymmetric land-ocean configuration that consists of two polar continents, separated by a homogeneous ocean (Fig. 5) using a 2-layer Maxwell viscoelastic Earth model.The coastlines are fixed and an axisymmetric ice load is located on the south pole with a cylindrical shape and a linear varying ice thickness as shown in Fig. 6b.Since the evolution of ice loading is known a priori, we can easily solve the standard SLE solution, for which one complete convolution of the SLE is needed over all time steps of the 480 kyr schematic experiment.These results are then used as a reference solution for the moving-time window experiments.
To test the accuracy of the moving-window technique, we have run a series of simulations that use a linear temporal interpolation between the heterogeneous time steps of the moving window.The moving window here covers the entire length of the simulation, but it consists of NT = 15 time steps.During the first 10 steps, t S is 1 kyr and thereafter five heterogeneous steps are used to complete the 480 kyr window.Figure 6a shows the normalised residual of the RSL change computed with the moving-time window.We have computed this as NormRes = (S mw − S full )/(S full ), where S mw is the RSL change calculated with the moving time window and S full is the RSL change computed with the standard SLE solution.Clearly, the largest differences between the moving window method and the standard SLE solution are located close to the ice sheet and in particular on top of the forebulge area (Fig. 6a and c).Here the GIA signal is more complicated than in the ice-covered area and in the far-field sites because of the lithosphere flexural response.
For the fully coupled experiments, we used an empirically derived window of L = 80 kyr with NT = 15.When using shorter time windows with this schematic set-up, information from past changes in ice-sheet variations is lost, whereas long windows do provide more information from the past but take more computational time.Similarly, we have performed a few small tests with the schematic set-up using shorter coupling intervals of t C = 200 or 500 years.Results of these tests indicate that a shorter coupling interval does not lead to a large improvement of the results.Although the higher time resolution resolves the initial exponential decay of the deformation better, the coarser resolution of the consecutive time steps of the moving time window with NT = 15 results in a larger deviation from the full solution (not shown).The window of 80 kyr and the coupling interval of 1 kyr used in the fully coupled experiments are therefore chosen as a trade-off between including sufficient memory of the deformation of the solid earth and computational time of the full 410 kyr run.

Simulations with coupled system over 410 kyr
Our simulations with the coupled ANICE-SELEN system provide variations of regional sea-level through time (see Supplementary Movie) using the default set-up as described in Sect. 2. As we show in Fig. 7a, RSL varies significantly between different locations and can be quite different from the eustatic curve (Fig. 7b).For the far-field site (Red Sea), the RSL is quite similar to the eustatic curve (black shading), although values are a bit less negative.The largest differences relative to the eustatic curve are found in the Antarctic Peninsula.Here, the change in RSL is always smaller than the eustatic curve due to the isostatic depression in response to the increase in the local ice load, a similar process occurs for western Europe.Differences in RSL can reach up to 100 m, for example between East Coast USA (green) and western Europe (blue), which are both relatively close to the large ice sheets in the NH.In particular the largest deviations from the eustatic curve occur during glacial maxima, directly after the LGM and MIS 6 (the penultimate glacial maximum).This is highlighted by the vertical dashed lines in (Fig. 7b), where the predicted RSL for East Coast USA shows a dip due to a lagged response of the collapse of the forebulge of the NaIS.In comparison, during interglacials local peak values are higher than eustatic, as indicated by the two vertical dashed lines at the Eemian and MIS 9 (310 kyr ago), which supports recent numerical simulation of the GIA correction for the MIS 11 interglacial (Raymo and Mitrovica, 2012).

Comparison with the eustatic solution
The initial set-up of the ANICE model as described in de Boer et al. (2013) calculates the change in sea level from the eustatic contributions of the four ice sheets relative to PD.In Fig. 1b, the four contributions of the ice sheets are shown over the 410 kyr time period.Clearly, the largest contributions arise from the NH ice sheets on Eurasia and North America.When we include the regional sea-level variations, the local evolution of ice thickness will obviously change due to the self-gravitation effect, especially for the marine parts of the ice sheets.
In Fig. 8a-d, we compare the modelled ice volume of the coupled ANICE-SELEN simulation with a simulation that is not coupled to SELEN (ice volume from de Boer et al., 2014).The largest differences occur during the glacial periods, especially for the AIS (Fig. 8a).For Antarctica, these differences are mainly observed in the marine sectors of the ice sheet, i.e.West Antarctica.Here, including the gravitationally self-consistent sea-level change reduces the growth of the ice sheet relative to the non-coupled run.As a result, with only eustatic variations (dashed line in Fig. 8a), the ice sheet grows significantly larger during a glacial period.Thus by including the self-gravitation effects and RSL changes, the growth of the West Antarctic ice sheet results in a local increase of sea level rather than a eustatic fall, which induces a slower advance of the ice sheet and thus a smaller ice volume.This self-stabilisation mechanism has been identified previously in coupled model simulations for Antarctica by Gomez et al. (2013).
The gravitationally self-consistent solution of the SLE provides a much more realistic behaviour of the response of the solid Earth to changes in ice and water loading.The viscoelastic Earth model accounts for the response on multiple timescales and provides a global solution, whereas a single response timescale of 3 kyr is used in the uncoupled solution of ANICE.For all four ice sheets (Fig. 8a-d), our current set-up of SELEN provides a lower response of the bedrock relative to the flexural Earth model used in the uncoupled ANICE simulation (dashed lines).This results in a lower total ice volume for the coupled solution, especially for the NaIS (Fig. 8d).Because the coupled simulation takes into account the change of the coastline over the globe (i.e. the time-dependent ocean function), the area of the total ocean is reduced by about 5 % of the global surface area during glacial maxima (Fig. 8f).Consequently, the total eustatic sealevel change of the two simulations (Fig. 8e) is coincidentally quite similar over the whole 410 kyr period.

Rotational feedback
An important aspect of the gravitationally self-consistent solution of the SLE is the rotational feedback, which is a new feature in SELEN.The changes in the mass distribution of ice and water induce a shift in the position of the rotational axis (polar wander) that has an ellipsoidal form (e.g.Gomez et al., 2010b).The difference as shown in Fig. 9b is described by the spherical harmonics of degree 2 (e.g.Mound and Mitrovica, 1998) (see also the Supplement Movie).As is shown in Fig. 9b, the positive contribution of the degree 2 signal is centred in the North Atlantic ocean and is related to the large increase in ice volume in the NH, which thus adds several metres to the fall in sea level during the LGM.These regional differences result in differences in the local ice thickness (Fig. 9d), but a minimal change in total ice volume.The addition of the rotation feedback, which is a significant contribution to the RSL change reaching up to 5 m or higher (Fig. 9b), is required for the correct interpretation of RSL data.In addition, there is a clear dynamical response of the ice sheets (Fig. 9d) to the differences in RSL, which results in large and significant changes in local RSL values close to the ice sheets.

Discussion and conclusions
In this paper we have presented a fully and dynamically coupled system of four 3-D ice-sheet-shelf models (de Boer et al., 2013(de Boer et al., , 2014) ) and a glacial isostatic adjustment model based on the SLE (Spada and Stocchi, 2007).The two key aspects of the coupling algorithm are the spatial discretisation and related interpolation of the ice volume from the four different regional ice-sheet-shelf models, and the temporal discretisation scheme with the related time interpolation.This system is the first fully coupled global ice-sheet-sea-level model available.Here, we have provided a simulation of the global solution of ice volume and relative sea-level variations over the past four glacial cycles.
The key aspect of our results is the dynamical response of the ice sheets to changes in RSL, which includes both the deformation of the bedrock in response to ice and water loading and the geoidal deformations.When an ice sheet grows, due to the self-gravitational pull of the ice sheet the RSL close to the ice sheets actually rises, whereas the global mean sea level falls.The self-gravitational pull thus acts to stabilise the ice sheets, as has also been shown by Gomez et al. (2013) with a coupled ice-sheet-sea-level model for Antarctica.Henceforth, ice volume is lower during glacial periods.Overall the coupled model results in lower ice volume relative to an uncoupled simulation that uses eustatic sea level derived from ice-volume changes only.We also include a time-dependent ocean function that accounts for the changes in the coastlines over the globe.This leads to a significant reduction in the ocean area during the glacial maxima and hence results in a nearly equal eustatic sea-level change compared to the uncoupled simulations.
The use of the 3-layer Maxwell viscoelastic Earth model gives a lower response in bedrock deformation due to the ice loading relative to the simplified model used in earlier studies (de Boer et al., 2013)  parameters, whereas several other studies clearly show the apparent sensitivity to varying the 1-D structure of the Earth model (e.g.Stocchi et al., 2013;Whitehouse et al., 2012a), which is a simplified version of the complex 3-D Earth structure in itself (van der Wal et al., 2013).Additionally, the adopted ice-sheet model parameters can also be investigated.For example, the mass balance parameters we use in ANICE (see de Boer et al., 2013) can be tested within a certain range of a physical parameter space (e.g.Fitzgerald et al., 2012).Similarly, ice flow and basal sliding can be varied (e.g.Maris et al., 2014).In a future study we will investigate the parameter space for both the Earth and ice-sheet models, and we will compare with observational data (see for example Whitehouse et al., 2012b;Briggs et al., 2013) on a global scale (Tushingham and Peltier, 1992).
Our results presented here provide a first overview of what can be achieved with our coupled ice-sheet-sea-level model.Similar to the choice of model parameters as mentioned above, our results are also sensitive to several other assumptions that are naturally necessary within a modelling framework.For example, our current set-up of a forward model starting from 410 kyr ago results in a final (t = 0 kyr) topography which is not in coherence with the actually PD topography.In future work we will include an additional correction for the difference between the final topography of our coupled experiment and the actual present-day topography, as has been suggested by Kendall et al. (2005).Second, to capture the full glacial contribution of the GrIS and its connection to the NaIS, we aim to include both ice sheets within the same ANICE model domain and thus also incorporate ice shelves for the GrIS and the possibility of merging of the ice sheets in North America and Greenland.
Lastly, to address the sensitivity tests raised in the previous paragraph a future study will include a thorough comparison with observational data such as near field RSL data (e.g.Whitehouse et al., 2012b), RSL over the glacial cycle (e.g.Deschamps et al., 2012;Grant et al., 2012;Austermann et al., 2013) and ice extent (e.g.Hughes et al., 2013).As shown in de Boer et al. (2013), the ANICE uncoupled model already compares reasonably well with other ice-sheet models and observations of sea level.
We presented here a complete dynamic system of four regional ice-sheet models and a global solution of the gravitational self-consistent sea level over time.Within this system, ice volume and global RSL changes are dynamically coupled.As a result, both the influence of the RSL on ice-sheet growth or retreat, and the change in RSL from changes in ice volume are taken into account within a consistent framework.We have developed a moving time window algorithm to account for past ice-sheet fluctuations.This allows us to calculate RSL and ice volume over the globe over four glacial cycles, starting 410 kyr ago.Our simulations show that especially during periods of rapid changes of sea level relative to PD, differences between regions can be very large; thus, showing the importance of this coupled system for modeldata comparison on a regional scale.

Figure 1 .
Figure 1.The uncoupled ANICE simulations using four ice-sheetshelf models.(a) in black the LR04 (Lisiecki and Raymo, 2005) benthic δ 18 O stack with the two separate contributions of ice volume (blue) and temperature (green).(b)The global eustatic sea level from ice volume in black with the four separate ice-sheet contributions of Eurasia (red), North America (blue), Antarctica (orange) and Greenland (green).Results are the same as shown inde Boer et al. (2014).

Figure 2 .
Figure 2. Scheme of the modelling framework.A coupling interval of 100 years is indicated by the black arrows, red arrows indicate a coupling every t C = 1000 years.The model is forced with benthic δ 18 O data, from which a NH temperature anomaly T NH is computed and forwarded to ANICE and the deep-water temperature module.ANICE computes the separate contributions of ice volume and deep-water temperature to benthic δ 18 O, which are sent back to the inverse routine every 100 years(de Boer et al., 2013).Every 1000 years ANICE forwards grounded ice thickness, the Ice load given in the input array (IA), to SELEN, which computes the gravitationally self-consistent sea level and bedrock topography adjustment that are coupled back to ANICE; in terms of the RSL, S is given in the output array (OA).

Figure 3 .
Figure 3.The four separate ANICE rectangular grid points for (a) the NH and (c) for Antarctica.The corresponding SELEN hexagonal elements for (b) the NH and (d) Antarctica.The colours correspond to each ice sheet: blue -NaIS, red -EuIS, green -GrIS and orange -AIS.The numbers of ANICE grid points and SELEN elements are shown in Table1.

Figure 4 .
Figure 4. Bedrock deformation according to a sequential increase of ice thickness on the south pole (Fig. 5) at a colatitude of 18 • .Every 1 kyr the ice thickness is increased with 20 m, any 20 m increase of ice thickness contributes to 80 kyr of viscoelastic crustal deformation.(a) At t = 1 kyr the predicted bedrock deformation at the 15 time steps of the moving time window.(b) Light grey markers indicate the fully discretised solution that is stored at t C =1 kyr resolution.(c) The predicted deformation for five consecutive time steps.The total solution, including past deformations and the elastic response is shown in red.Insets for panels (a)-(c) show the implied ice thickness variations, steps of 20 m per 1 kyr.

Figure 5 .
Figure 5.A slice of the schematic Earth with 2 polar continents as used in the moving time window experiments.LT: lithosphere of 100 km; UM: upper mantle, a viscosity of 10 21 Pa s; LM: lower mantle, 2 × 10 21 Pa s; CO: inviscid core.

Figure 6 .
Figure 6.(a) Normalised residual of RSL change with the moving-time window relative to the standard SLE solution (Eq.3).Simulations are performed with a schematic Earth with two polar continents (Fig. 5), and an ice sheet on the south pole.The y axis shows the time, the x axis the colatitude ( • ) relative to the south pole.(b) The ice-thickness variations that are applied as a cylindrical shaped ice sheet up until 18 • colatitude (vertical red dashed line in panel (a).(c) The RSL at a colatitude 20 • , in black the full standard SLE solution, and in red the solution with the moving time window.
Figure 7. (a) RSL at four different locations as predicted by the coupled model.Red sea (red), Antarctic Peninsula (orange), western European coast (blue) and East Coast of the USA (green).The black line indicates the eustatic RSL change, calculated as the global mean change of RSL for the entire ocean.(b) RSL minus eustatic for each of the four locations.The four locations are indicated with coloured dots in Fig. 9. Vertical dashed lines indicate key periods, from left to right: the LGM (18 kyr ago), Eemian (122 kyr ago), marine isotope stage (MIS) 6 (138 kyr ago) and MIS 9 (310 kyr ago).

Figure 8 .
Figure 8.A comparison of the coupled ANICE-SELEN solution with model runs using the eustatic sea level.(a) The ice volume of the AIS, (b) ice volume of the EuIS, (c) ice volume of the GrIS and (d) ice volume of the NaIS.For all figures the coupled solution is shown by the solid line and the dashed line represents the model runs using the eustatic sea level (derived from ice volume as in de Boer et al., 2014).(e) The eustatic RSL change from the coupled run with SELEN is in red and the eustatic sea level from de Boer et al. (2014) is in green.(f) The evolution of the time-dependent ocean function, shown on the left y axis is the total ocean area, the right y axis shows the percentage of ocean covered grid points over the globe.

Figure 9 .
Figure 9. Results of a coupled ANICE-SELEN run at the last glacial maximum (here 18 kyr ago).(a) RSL change with respect to the eustatic (= 111.3 m below PD) including rotational feedback.(b) The difference in RSL of a run using rotational feedback (as in a) with a run without rotational feedback.(c) The total ice loading from ANICE (= 112.8 m s.e.) including rotational feedback.(d) The difference in ice loading of a run using rotational feedback (as in c) with a run without rotational feedback.In panel a the coloured dots indicate the locations illustrated in Fig. 7.A full time evolution of the 410 kyr long simulation of these maps is shown in the Supplement Movie.

Table 1 .
Separate model parameters for the four ice-sheet models.

Table 2 .
Model parameters for time discretisation.