Capturing the Interactions Between Ice Sheets , Sea Level and the Solid Earth on a Range of Timescales : A new “ time window ” algorithm

Retreat and advance of ice sheets perturb the gravitational field, solid surface and rotation of the Earth, leading to spatially variable sea-level changes over a range of timescales (~O100-6 years), which in turn feed back onto ice sheet dynamics. Coupled ice-sheet – sea-level models have been developed to capture the interactive processes between ice sheets, sea level and the solid Earth, but it is 10 computationally challenging to capture short-term interactions (~O100-2 years) precisely within longer (~O103-6 years) simulations. The classic coupling algorithm assigns a uniform temporal resolution in the sea-level model, causing a quadratic increase in total CPU time with the total number of input ice history steps, which increases with either the length or temporal resolution of the simulation. In this study, we introduce a new “time window” algorithm for sea-level models that enables users to define 15 the temporal resolution at which the ice loading history is captured during different time intervals before the current simulation time. Utilizing the time window, we assign a fine temporal resolution (~O100-2 years) for the period of ongoing and recent history of surface ice and ocean loading changes and a coarser temporal resolution (~O103-6 years) for earlier periods in the simulation. This reduces the total CPU time and memory required per model time step while maintaining the precision of the model 20 results. We explore the sensitivity of sea-level model results to the model's temporal resolution and show how this sensitivity feeds back onto ice sheet dynamics in coupled modelling. We apply the new algorithm to simulate the sea-level changes in response to global ice-sheet evolution over two glacial cycles and the rapid collapse of marine sectors of the West Antarctic Ice Sheet in the coming centuries, providing appropriate time window profiles for each of these applications. The time window algorithm 25 improves the total CPU time by ~50 % in each of these examples, and this improvement would increase with longer simulations than considered here. Our algorithm also allows coupling time intervals of annual temporal scale for coupled ice-sheet – sea-level modelling of regions such as the West Antarctic that are characterized by rapid solid Earth response to ice changes due to the thin lithosphere and low mantle viscosities. 30 https://doi.org/10.5194/gmd-2021-126 Preprint. Discussion started: 17 May 2021 c © Author(s) 2021. CC BY 4.0 License.

computationally challenging to capture short-term interactions (~O10 0-2 years) precisely within longer (~O10 3-6 years) simulations. The classic coupling algorithm assigns a uniform temporal resolution in the sea-level model, causing a quadratic increase in total CPU time with the total number of input ice history steps, which increases with either the length or temporal resolution of the simulation. In this study, we introduce a new "time window" algorithm for sea-level models that enables users to define 15 the temporal resolution at which the ice loading history is captured during different time intervals before the current simulation time. Utilizing the time window, we assign a fine temporal resolution (~O10 0-2 years) for the period of ongoing and recent history of surface ice and ocean loading changes and a coarser temporal resolution (~O10 3-6 years) for earlier periods in the simulation. This reduces the total CPU time and memory required per model time step while maintaining the precision of the model 20 results. We explore the sensitivity of sea-level model results to the model's temporal resolution and show how this sensitivity feeds back onto ice sheet dynamics in coupled modelling. We apply the new algorithm to simulate the sea-level changes in response to global ice-sheet evolution over two glacial cycles and the rapid collapse of marine sectors of the West Antarctic Ice Sheet in the coming centuries, providing appropriate time window profiles for each of these applications. The time window algorithm 25 improves the total CPU time by ~50 % in each of these examples, and this improvement would increase with longer simulations than considered here. Our algorithm also allows coupling time intervals of annual temporal scale for coupled ice-sheet -sea-level modelling of regions such as the West Antarctic that are characterized by rapid solid Earth response to ice changes due to the thin lithosphere and low mantle viscosities.

Introduction
It is well established that sea-level changes in response to ice-sheet changes feed back onto the evolution of ice sheets (e.g., Gomez et al., 2012;de Boer et al., 2014;Konrad et al., 2015;Larour et al., 2019).
Changes in grounded ice cover perturb the Earth's gravitational field, rotation and viscoelastic solid 35 surface, leading to spatially non-uniform changes in the heights of the sea surface geoid and the solid Earth, i.e., sea-level changes (e.g., Peltier, 1974;Farrell and Clark, 1976;Woodward, 1888;Mitrovica and Milne, 2003). Sea-level changes occur as an instantaneous response to the surface (ice and water) loading changes associated with elastic deformation of the solid Earth and changes in gravity and rotation, followed by a slower response over tens of thousands of years due to the viscous mantle flowing back 40 towards isostatic equilibrium, once again accompanied by gravitational and rotational effects. The spatial and temporal scales of the solid Earth response to ice loading changes depend on the rheological structure of the elastic lithosphere and viscoelastic mantle, which are both radially and laterally heterogeneous (e.g., Dziewonski and Anderson, 1981;Morelli and Danesi, 2004;Nield et al., 2014;An et al., 2015;Lloyd et al., 2020). The contribution from viscous deformation to sea-level changes in regions with 45 thinner lithosphere and lower mantle viscosities such as the West Antarctica occurs on shorter timescales (O ≤ 10 2 years; e.g., Barletta et al., 2018) than they do in regions with thicker lithosphere and higher mantle viscosities such as North America (e.g., Mitrovica and Forte, 2004) and is more localized to the loading changes, calling for higher spatiotemporal resolution in modelling applications in these regions. 50 Spatially variable changes in the sea surface geoid and the solid Earth (i.e., sea level) have different dominant mechanisms in influencing ice sheets in marine and continental settings. The evolution of a marine-based ice sheet is strongly dependent on the slope of bedrock underneath the ice sheet and local ocean depth at the grounding line (e.g., Weertman, 1974;Thomas and Bentley, 1978;Schoof, 2007).
Thus, deformation of the bedrock beneath the ice and sea level changes at the grounding line in response 55 to the marine-based ice sheet's growth and retreat affect the ice flux across the grounding line (Gomez et al., 2010(Gomez et al., , 2012. In the continental setting, solid Earth deformation beneath an evolving land-based ice sheet alters the slope and elevation of the ice surface in the atmosphere. This, in turn, influences the ice sheet's surface mass balance (e.g., Crucifix et al., 2001;Han et al., 2021;van den Berg et al., 2008).
The interactions between ice sheets, sea level and the viscoelastic solid Earth are active over a range of timescales, and several studies have developed coupled ice sheet-sea level models to investigate these interactions (Gomez et al., 2012(Gomez et al., , 2013deBoer et al., 2014;Konrad et al., 2014). Studies have applied coupled modelling to simulate the evolution of the Antarctic Ice Sheet (AIS) during the last deglaciation (Gomez et al., 2013Pollard et al., 2017), the Pliocene  and the future 65 (Gomez et al., 2015;Konrad et al., 2015;Larour et al., 2019), the evolution of the Northern Hemisphere ice sheets over the last glacial cycle (Han et al., 2021) and the global ice sheets over multiple glacial cycles (deBoer et al., 2014). These studies capture the interactions between ice sheets and sea level at a temporal resolution of as short as 50 years for the millennial timescale simulations and 200 years for the glacial timescale simulations, but moving to longer simulations or greater spatiotemporal resolution 70 presents a computational challenge, requiring a trade-off between those two in simulations.
There is a need to overcome this challenge in order to understand ice sheets and sea-level changes over a wider range of timescales and in greater detail, especially as the spatiotemporal resolution and extent of paleo records improves (e. g., Khan et al., 2019;Rovere et al., 2020;Gowan et al., 2021). Motivations 75 include running simulations over longer time periods in the past (e.g., from the warm mid-Pliocene to the modern), or in higher spatiotemporal resolutions in order to accurately capture rapid paleo ice-sheet variability and sea-level rise events observed in geological records (e.g., Ice Rafted Debris events -Weber et al., 2014;Meltwater Pulse 1A event -Fairbanks, 1989;Deschamps et al., 2012;Brendryen et al., 2019).
Furthermore, the present-day WAIS sits atop rapidly responding bedrock (e.g., Barletta et al., 2018;Lloyd 80 et al., 2020;Powell et al., 2020) and is under the threat of catastrophic collapse in a warming climate (e.g., SROCC, 2019). To capture the dynamics of such rapid retreat of ice sheets and the associated sea-level changes, models may need to employ annual-decadal scale resolution (e.g., Larour et al., 2019).
The computational challenge introduced above arises only in a coupled ice sheet -sea level modelling 85 context where (unlike in stand-alone sea level modelling applications where ice cover changes are prescribed, e.g., Peltier 2004;Lambeck 2014) the ice cover changes are unknown at the start of a simulation and predicted by the ice sheet model as the simulation progresses. That is, with the standalone ice age sea-level model algorithm described in Kendall et al. (2005), the model takes in the full history of ice loading at the start of the simulation and computes associated sea-level changes across all 90 time steps and outputs results at once at the end of the simulation. On the other hand, in a coupled icesheet -sea-level simulation, ice cover changes are predicted by a dynamic ice-sheet model and provided to the sea-level model. The sea-level model, in turn, provides updated bedrock elevation and sea surface heights to the ice-sheet model. This exchange of models' outputs happens at every coupling time interval, which necessitates a 'forward modelling' scheme (as described in Gomez et al. 2010) in the sea-level 95 model: The sea-level calculation at every new coupling time interval requires the history of ice loading since the beginning of the coupled simulation as input. The classic forward sea-level modelling algorithm adopted in coupled models employs a uniform temporal resolution throughout a simulation, which leads to a quadratic increase in the amount of surface loading history with the length of a simulation. The sealevel calculation thus becomes more computationally expensive as the simulation progresses and can 100 make very long or very finely temporally resolved simulations computationally infeasible (e.g., Using this framework, coupled simulations in past studies have been limited to 40-125 ky with a temporal resolution of 200 years; Gomez et al., 2013;Han et al., 2021).
To overcome this challenge, de Boer et al. (2014) presented what they called a "moving time window" 105 algorithm in a sea-level model (SELEN, Spada and Stocchi, 2007) and performed coupled ice-sheetsea-level model simulations over four glacial cycles (410 ky). Using the characteristics of exponentially decaying viscous deformation of the mantle, de Boer et al. (2014) interpolated "future" viscous deformation associated with ongoing surface loading changes and added up the interpolated values at later time steps to obtain deformation due to "past" loading changes. They also applied a coupling time 110 interval of 1 ky, while other studies have suggested that simulations over the last deglaciation require coupling intervals of at least 200 years (Gomez et al., 2013) and coupled simulations of future retreat of the Antarctic ice sheet have adopted coupling times of tens of years or less (Gomez et al., 2015;Pollard et al., 2017;Larour et al., 2019) to capture the decadal to centennial-scale interactions between ice sheets, sea level and the solid Earth.
In this study, we develop a new time window algorithm, which takes a different approach from deBoer et al. (2014) to overcome the computational challenges posed in coupled ice-sheet -sea-level modelling.
We modify the classic forward model sea-level algorithm introduced in Gomez et al. (2010) by systematically reducing the temporal resolution of earlier ice history while maintaining high resolution in 120 recent loading. We present the algorithm in Section 2 and perform a suite of simulations with idealized ice sheet evolution and bedrock geometry to show how the temporal resolution of a sea-level model influences the predicted sea level (Section 3.1) and its influence on Northern Hemispheric Ice Sheet dynamics through the last glacial cycle (Section 3.2). Next, we apply our time window algorithm to simulate sea-level changes due to the evolution of the global ice sheets over the last two glacial cycles 125 and due to future Antarctic Ice Sheet evolution in the coming centuries (Section 3.3), presenting appropriate time window parameters for each scenario. We finish with a discussion of the results and concluding remarks in Section 4.

130
We incorporate our time window algorithm into the forward sea-level model presented in Gomez et al. (2010), which draws on the theory and numerical formulations in Milne (2003), Kendall et al. (2005), and Mitrovica et al. (2005). In the forward modelling, at every time step tj, the sea-level model performs a one-step computation between times tj-1 and tj of the global sea-level change associated with ongoing (between tj-1 and tj) and past (between t0 and tj-1) ice loading changes. The numerical form where j is an index for the current time step, represents changes over a single time step (e.g., from tj-1 140 to tj) and ∆ represents the total change since the initial time t0. Thus, !, and ! represent changes in ice thickness (I) and ocean loading (S) between tj-1 and tj, and ∆ ! represents the change in ocean loading before the current time step between t0 and tj-1. and ∆$ ! % represent the geographically non-uniform and uniform components of the globally defined total sea-level change, respectively. * represents an ocean-mask function, defined as 1 where sea level is positive and there is no grounded ice, or zero 145 otherwise. T0 represents initial topography at t0, where the topography is defined as the negative of the globally defined sea level (∆ ! = + ∆$ ! % ). Since the focus of this paper is to modify the traditional implementation of the sea-level equation, we refer readers to Mitrovica and Milne (2003), Kendall et al. (2005) and Gomez et al. (2010) for the detailed derivation of this equation and implementation of its numerical algorithm.
150 Figure 1a represents the classic forward sea-level model algorithm (Gomez et al., 2010) where the time interval dt between each time step of the ice history is uniformly fixed throughout a simulation. By the end of the simulation, the total number of ice history steps (Nj) considered in the calculation across the final time step from tj=f-1 to tj=f is simply the length of the simulation (L_SIM) divided by the prescribed 155 time interval (dt). Thus, the number of time steps in a simulation increases either by performing a longer simulation (i.e., larger L_SIM) or increasing the temporal resolution (i.e., smaller dt) of a simulation, and the CPU time increases quadratically. resolve the prescribed non-uniform time steps (see Fig. 1b-1). It also creates an array of iceload file numbers that, by convention, start at 0 and increment by 1 each time step forward. When the simulation begins (i.e., takes one step forward from j = 0 to j = 1), the first two elements of this array (iceload files '0' and '1') overlap the last two elements of the binary template (see the top red box in 1b-2). Overlapping 170 elements are multiplied together to generate masked iceload history files for the sea-level model to read in at the time step. The sea-level model only reads in those ice files masked with binary value of 1.
However, to ensure that the solid Earth retains the memory of the initial loading, the sea-level model always reads the initial iceload file (see the dotted box and resulting masked iceload files in Fig. 1b-2).

175
At every simulation step j>1, the template marches forward by one element relative to the iceload file array, and the multiplication process repeats followed by the sea-level calculation. Our algorithm starts filling in the first internal time window to its prescribed length (L_ITW1), followed by the other internal windows in order (L_ITWk for k=2, 3, 4). By the end of the simulation, the time window grows to the full prescribed profile (see the j = 8 result in Fig. 1b-2).

180
Overall, this time window algorithm limits the increase in the amount of surface loading history with simulation length or temporal resolution, improving the computational efficiency of sea-level model calculations (compare Nj values in Fig. 1a and Fig. 1b-2). The time window algorithm also enables the sea-level model to capture both short-and long-term interactions between ice sheets, sea level and the 185 solid Earth in coupled ice-sheet -sea-level simulations.
In the next section, we perform a suite of sensitivity tests performing standalone sea-level simulations and coupled ice-sheet -sea-level simulations to test the sensitivity of model results to temporal resolution of the sea-level model. The sea-level model takes two inputs, an Earth model and an ice history model. 190 We adopt 1-D Earth model in all simulations; the elastic and density profile of the Earth structure are given by the seismic model PREM (Dziewonski & Anderson, 1981). For mantle viscosity, we adopt a lithospheric thickness of 120 km and upper mantle viscosity of 5×10 20 PaS and lower mantle viscosity of adopt the best-fitting radially varying Earth model from Barletta et al. (2018), characterized by a 195 lithospheric thickness of 60 km and upper mantle viscosities of ~10 18 -10 19 PaS. Ice history inputs are described in each section with reference to corresponding figures. In all simulations, we perform sea-level calculations up to spherical harmonics degree and order of 512.

Sensitivity of sea-level model outputs to temporal resolution
Before exploring the new time window algorithm in sections 3.2 and 3.3, we begin by performing a series of experiments with a standalone sea-level model adopting the classic algorithm to demonstrate how 205 predicted topography changes (that is, negative of sea-level changes) at near and far-field locations vary with the (uniform) temporal resolution of the ice history (as in Fig. 1a). We generate an idealized axissymmetric input ice history based on the equation for an equilibrium ice surface profile for viscous ice provided in Cuffey and Paterson (1969, The Physics of Glacier), and the ice sheet grows and retreats centered at the South pole (Fig. 2). The initial topography for the buildup-phase simulations is idealized 210 such that its elevation is 1000 m between latitude 60-90 S degrees and -1000 m (ocean with a depth of 1000 m) everywhere else. In Fig. 2, we consider predicted changes in topography at three locations: A location at the center of loading (latitude 90 S degrees), near the periphery of the ice sheet at its maximum extent (latitude 65S degrees, where the peripheral bulge formed around the ice sheet is largest), and in the far-field near the equator at latitude 0 degrees. We begin by discussing the general behaviour of 215 topography at these sites in benchmark simulations performed at a temporal resolution of 1 ky (shown by the black dots in Fig. 2). sites experience a decrease in sea level (increase in topography) as ice becomes locked up on land (middle panel of Fig. 2c).

225
Figures 2d-f show results over a 20-ky long retreat phase. The initial topography for these simulations is adopted from the final topography modelled in the benchmark buildup-phase simulation with dt = 1 ky (i. e., black dots in the middle panel of Fig. 2a). As the ice sheet retreats, the topography uplifts at the center of the ice load (middle panel of Fig. 2d), subsides at locations peripheral to the ice (middle panel loading. For example, the loading change for the first 10 ky and the 5 ky resolution simulation is applied in two steps at 5 and 10 ky, while in the 10 ky resolution simulation, the full load is applied once at 10 ky. The latter thus does not capture the viscous signal due to the loading that takes place before 10 ky. These maximum differences ("errors") and the spread of the errors decrease gradually towards zero with Note that at the equatorial site, the rate at which the error decreases towards zero is slower than the near-250 field sites (e.g., see the pink line after reaching its peak in Figs. 2c and f). This is because there is active water loading occurring at this site even after the ice has stopped evolving, which prolongs the differences in deformation of the lower resolution simulations compared to the benchmark simulation ("water-loading effects", e.g., Han et al., 2018).

255
The timing of the maximum errors at each site corresponds to the size of topography changes at the site, which in turn depends on the distance to, and size of an evolving ice sheet. For example, at the near-field sites during the buildup phase, the peak differences in simulations occur as soon as the ice starts loading at the centre-of-loading site (bottom panel of Fig. 2a) while at the peripheral-bulge site the peak differences occur at 20 ky when the ice sheet reaches its maximum volume and extent (bottom panel of 260 Fig. 2b). During the retreat phase, the peak differences at the peripheral-bulge site occur as soon as the ice starts retreating. The differences then get smaller at this site as the ice sheet and its peripheral bulge retreat further away from the site towards the pole (bottom panel of Fig. 2e). The ice sheet's centre, on the other hand, experiences the peak differences at 20 ky when ice loading directly above this site disappears (bottom panel of Fig. 2d).

265
While the timing of the maximum differences in the near-field sites is most sensitive to the size and proximity of ice to the sites, at the equatorial site it is the most sensitive to the ice volume change. During the buildup phase, the greatest sea-level changes at the equatorial site (and thus the maximum errors) occur at 20 ky (bottom panel of Fig. 2c). During the retreat phase, the errors peak as soon as the ice starts 270 retreating and the rate of change of ice volume is largest (bottom panel of Fig. 2f). For both phases, the timing of maximum error is related to the maximum ice volume change. We note that the ice thickness at the centre of loading (as shown in the top panels of Fig. 2a and d) changes linearly, but the actual volume change is nonlinear because of the changes in the ice sheet's extent; the volume change across one time step is greater when the ice sheet is more extensive.
Overall, the idealized-loading simulations show that sea-level model outputs are sensitive to the model temporal resolution. This is because the timing of ice loading is different with different temporal resolutions. The sensitivity at a location depends on its setting (above or below sea level), the size of ice loading changes and the distance (near-field or far-field) to the changing ice load. However, the sensitivity 280 at all sites decreases with time after the ice loading event. These results suggest (as expected from the literature on the viscoelastic response of the Earth to surface loading, e.g., Peltier, 1974) that higher resolution information about ice cover changes is required for the ice history immediately prior to the current time step in a simulation, and lower resolution will suffice for earlier ice cover changes. The specific temporal resolution required will depend on both the rates of change of the ice cover and the 285 Earth's viscosity structure, which we explore in two contrasting examples in Section 3.3.
In this section, we have highlighted the sensitivity of predicted sea-level changes to the temporal resolution of the inputted ice history in classic, standalone sea-level simulations. In the following section, we explore how the differences in sea level predicted with different temporal resolution influence the ice 290 sheet evolution in coupled ice sheet -sea level simulations and how the time window algorithm can reconcile the errors. 295 This section explores how the differences in predicted sea-level change due to temporal resolution of the input ice history discussed in the previous section impact ice dynamics in coupled ice-sheet-sea-level model simulations. To do this, we perform a suite of coupled simulations over the Northern Hemisphere through the last glacial cycle (125 ky) incorporating different sizes of uniform time steps with the standard algorithm ( Fig. 3a) and nonuniform time steps applying the time window algorithm (Fig. 3b). We employ 300 the PSU 3D dynamic ice-sheet model by Pollard and DeConto (2012) and adopt the same set of ice model parameters (e.g., climate forcing, basal slipperiness, spatial and temporal domain and resolutions) used in the simulations from the main text of Han et al. (2021). The ice-sheet model has a standalone time step of 0.5 yr. We note that the coupling interval over which the ice-sheet model and sea-level model exchange their outputs (i.e., ice thickness and topography, respectively) corresponds to the size of the most recent 305 time step within the sea-level model (i.e., dt1 in Fig. 1b).  Fig. 3a). Spatially, the differences occur mainly in the Laurentide Ice Sheet in North America (we don't show this in the figure). As illustrated in Fig. 2, a lower temporal resolution of the ice history during 315 the Laurentide Ice Sheet retreat before 80 ka leads to less uplift of the bedrock beneath the ice sheet, keeping the ice surface at a lower (and thus warmer) evelation in the atmosphere. This lower ice elevation causes more intense deglaciation of the Laurentide Ice Sheet. It also prohibits the ice sheet from growing large during buildup phases later on (the role of deformational effects on ice sheet dynamics is discussed in detail in Han et al. 2021). Furthermore, The NHIS volume fluctuation becomes less smooth when the 320 coupling time interval is increased to dt = 5 ky and dt = 10 ky (red and grey lines in Fig. 3a). We presume that this is because a large change in bedrock over a long coupling time causes the ice-sheet model to respond unstably. These results suggest that Northern Hemisphere coupled simulations over the last glacial cycle require a coupling time of hundreds of years or less to accurately capture the interactions between the ice sheets, bedrock elevation and sea level (compare the black-dotted, black and blue lines 325 in Fig. 3a).

Sensitivity of modelled ice sheet dynamics to temporal resolution with a coupled sea-level model
While we might expect that the compute time would always increase with higher temporal resolution in the case of uniform time stepping (Fig. 1a), it is interesting to note that the 10-ky time-step simulation took an hour longer than the 5-ky case. This is because in the former simulation, the ice model took longer 330 to converge to a solution because of infrequent and dramatic bedrock changes provided by the sea-level model (as hinted by the unstable fluctuation in the ice volume -the grey line in Fig. 3a). Finally, while there are very small differences in predicted ice volume between the dt = 0.2 ky and dt = 0.1 ky simulations (black and black-dotted lines in Fig. 3a), CPU time increases by greater than a factor of two from ~45 to ~98 hr, suggesting that dt = 0.2 ky is a suitable choice of coupling time for glacial cycle simulations.

335
In Fig. 3b, we apply the time window algorithm to the coupled glacial-cycle simulation rather than adopting uniform temporal resolution in the ice history. We perform three simulations with different time window profiles illustrated in the schematics shown in  In this section, taking the last glacial cycle as an example, we have shown that a coarse temporal resolution (e.g., dt=1 ky or longer) causes less precise coupled ice-sheet -sea-level simulation results. We have also demonstrated that how the time window algorithm can be used in coupled simulations to maintain the precision of the modelled topography changes and ice sheet dynamics while significantly reducing the computational cost compared to simulations with the standard algorithm. In the same way, the time 355 window algorithm can be applied to other coupled simulations that are otherwise infeasible. In the next section, we derive time window profiles that are suitable for two-glacial-cycle global ice-sheet simulations and the future rapidly retreating Antarctic Ice Sheet simulations.

Derivation of time window profiles for different applications
In this section, we apply the time window algorithm in the sea-level model to two contrasting examples.
First, we consider sea-level changes in response to the evolution of global ice sheets over the last two glacial cycles (240 ky) modified from Han et al. (2021). Then, we consider a simulation of the rapid future retreat of the West Antarctic Ice Sheet in the coming centuries taken from DeConto et al. (2021). West Antarctica is known to have an upper mantle viscosity up to several orders of magnitude lower than on average (Barletta et al., 2018;Lloyd et al., 2020;Nield et al., 2014). For each scenario, we perform a suite of simulations in which we vary the time window parameters (i.e., L_ITWk and dtk; Fig. 1) and compare them to a benchmark simulation with uniform high-resolution time stepping to arrive at an optimal choice 370 of a time window profile. Here we note that the experiments are done in standalone sea-level simulations (i.e., ice cover is prescribed rather than provided by a dynamic ice-sheet model) because the benchmark coupled simulations for these scenarios become less feasible without the time window algorithm.

380
(2021) that we adopt here as an input to the sea-level model. The original simulation covers the last glacial cycle (125 ka). It includes ice-cover changes predicted with the dynamic PSU ice-sheet model (Pollard & DeConto, 2012) in the Northern Hemisphere and Antarctic ice cover changes taken from the ICE-6G_C ice history model (Argus et al., 2014;Peltier et al., 2015). To extend our ice history input to cover two glacial cycles, we first take the ice history for 120 ka from the original simulation (Han et al. 2021), then 385 repeat this ice history to cover an additional glacial cycle going back to 240 ka. We replace the ice history between 125-120 ka with the ice history between 120-115 ka. This is to make the ice volume curve continuous at the last interglacial. We note that the goal of this experiment is not to produce an accurate glacial history but to produce a sample long timescale, global ice history that contains the spatiotemporal detail provided by a dynamic ice-sheet model.

390
To explore the choice of time window parameters for this global glacial-cycle scenario, we first perform a standard sea-level simulation in which we assign a uniform temporal resolution of 0.2 ky throughout the 240 ky simulation. We take this simulation as our benchmark, and then perform a suite of simulations in which we systematically vary the temporal resolution (dt) of internal time windows (ITW) that cover 395 periods 240-120 ka, 120-50 ka and 50-20 ka in the simulation (see the internal time windows marked by dashed lines in Fig. 4a). We choose these internal time windows based on the timing of ice volume variations and the results of our idealized tests in section 3.1. That is, our first internal time window covers the last deglaciation, the next covers the preceding growth phase in the simulation, then the rest of the glacial cycle back to the last interglacial, and finally the entire previous glacial cycle. We note that in the 400 absence of knowing the specific details of the ice cover changes a priori (as in coupled model simulations), the internal time windows may also be set based on the timing of the climate forcing that serves as input to the model. Sensitivity tests (not shown here) varying the internal time windows lengths to account for potential offsets between the timing of climate forcing and ice sheet response indicate that the timing of these internal time windows need not be set very precisely, with less sensitivity for earlier ice history.

405
When we explore each internal time window in turn by varying the internal time step (i.e., temporal resolution dt), starting from the earliest, and fixing the temporal resolution at 0.2 ky for all periods beyond the internal time window. Then, we compare the total CPU time (Fig. 5d) and the precision of our results by calculating the root mean squared errors (RMSE) in predicted topography from these simulations to 410 results from the benchmark simulation. The RMSE are calculated based on the following expression: Once we choose an optimal temporal resolution for the internal time window based on the calculated RMSE and CPU time, we move on to explore the next internal time window and we repeat the same 420 procedure.
We start by exploring the internal time window covering the earliest period between 240-120 ka (see the purple bar in Fig. 5a). Varying the internal time step between 5-40 ky for this period (Figs. 5a-d), the RMSE in predicted topography is zero for the first 120 ky (Fig. 5b) 5h). Hence, we adopt a temporal resolution of 1 ky for this period.
Finally, Figs. 5i-l explores the effects of varying the size of the internal time step for the period between 50-20 ka (purple bar in Fig. 5i). Based on the above discussion, the size of the internal time steps for the 460 periods 240-120 ka and 120-50 ka are fixed at dt = 10 ky and 1 ky, respectively. We arrive at the final optimal time window profile by identifying dt = 0.4 ky (red line in 5h-i) as an optimal internal time step for this internal time window. This profile keeps the RMSEs in output topography below 0.4 m throughout the entire simulations. The total CPU time is reduced to ~26.9 hr, a 54% reduction compared to the benchmark with uniform time stepping of 0.2 ky. We note that the CPU times shown in Fig. 5 are based 465 on standalone sea-level simulations only. This time window algorithm is designed for sea-level calculations performed within a coupled ice sheet -sea level simulation, and compute times will be similarly reduced in this context. Moreover, the reduction will grow for longer simulations as the CPU time in the standard simulation will increase quadratically whereas the time window simulation will suppress the rapid growth of rate of increase.

Application to future Antarctic Ice Sheet changes
In this section, we develop a time window profile for application to simulating future West Antarctic Ice  We start by finding an optimal internal time step (i.e., temporal resolution) for the period between 0-200 500 yr simulation time ( Fig. 6a; also marked as a purple bar in Fig. 7a) during which the WAIS starts retreating, and the rate of retreat starts accelerating just after 100 yr. Fig. 7b shows that RMSEs in predicted topography compared to the benchmark simulation start increasing after 350 yr (this is because the temporal resolution is set to be the same as that in the benchmark simulation during the first 350 yr).
The RMSEs remain smaller than 10 cm for simulations with an internal time step of 10 yr and 25 yr (red 505 and magenta lines in Fig. 7b) and below 0.85 cm with an internal time step of 5 yr (black line in Fig. 7b).
The number of ice history steps (Nj) that the sea-level model considers at time step tj starts diverging after 350 yr. At the last time step of the simulation, the benchmark simulation considers 550 ice history steps while Nj considered for the simulations with an internal time step dt = 5-50 yr are reduced by ~ 29-36 %, respectively (Fig. 7c). The total CPU times are reduced by ~ 4-8 %, from ~17 hr with the standard 510 simulation to between ~15.6-16.3 hr for the others (Fig. 7d). We choose dt = 5 yr as an appropriate internal coarser temporal resolution, keeping the RMSEs in predicted topography below 0.07 m throughout the simulation (black line in Fig. 7j), without a significant increase in compute time and ice history steps compared to the other coarser simulations (Fig. 7k). The total computing time for the 5 yr simulation is ~ 8.5 hr (Fig. 7l), which is a ~ 50 % reduction from the computing time of the benchmark simulation (~ 17

Discussion and Conclusions
We have developed a new time window algorithm that assigns nonuniform temporal resolution to the inputted ice cover changes in a forward sea-level model (Gomez et al., 2013, and restricts the linear 555 increase in the number of ice history steps that a sea-level model has to consider at each time step. Our algorithm allows coupled ice-sheet-sea level models to capture short-term O (≤10 2 yr) interactions between ice sheets, the solid Earth and sea level within simulations across a range of timescales. The algorithm improves computational feasibility while maintaining the precision of the sea-level (and thus coupled ice-sheet -sea-level) simulations.

560
In benchmarking the algorithm, we first tested the sensitivity of sea-level model outputs (i.e., predicted topography) to the temporal resolution adopted in idealized simulations (Fig. 2). Our results show that sea-level simulations with coarser temporal resolution do not accurately capture the timing and geometry of ice loading, and this leads to missing viscous signals and thus an underestimation of topography 565 changes. We also found (as suggested in earlier literature, e. g., Peltier, 1974) that there is a stronger sensitivity to more recent loading, indicating that higher temporal resolution is required close to the current time step in a simulation. We then performed coupled ice-sheet -sea-level simulations through the last glacial cycle over the Northern Hemisphere with varying temporal resolution. Our results demonstrated that the underestimated magnitudes in predicted topography and infrequent topography 570 updates in the coupled simulation with a lower-temporal resolution lead to smaller and sometimes unstable ice volume fluctuations (Fig. 3a). Our results also identify that 0.2 ky is the optimal coupling time interval for glacial-cycle simulations with broad spatial scale on a model of Earth structure representative of the global average Earth structure typically adopted in ice-age sea-level studies (e.g., Lambeck et al., 2014). When we utilize the time window algorithm and capture short-term, recent Overall, our results demonstrate that capturing short-term responses during a period including and temporally close to ongoing surface loading changes is important. At the same time, a coarser temporal resolution can be used for past loading changes. This is expected based on normal mode theory where the 590 solid Earth signals comprised of normal modes with shorter decay times associated with the loading changes would have already relaxed out after simulations have proceeded (Peltier, 1974). and this has the potential to have a significant impact on future evolution of marine ice in the region (Gomez et al., 2015). Furthermore, recent work by Larour et al. (2019) suggests that high spatial resolution and short time-stepping may be required to capture the elastic component of deformation in this region. This work together suggests that annual to decadal scale coupling time is likely needed to capture the short-term interactions in a coupled model that may play a significant role in the stability of 630 marine-based WAIS. In Section 3.3.2, we have performed a benchmark sea-level simulation with a future WAIS evolution at 1 yr temporal resolution. We introduced a set of time window parameters that allows us to keep a coupling interval of 1 yr while improving the total CPU time by 50 % (Fig. 7) and maintaining the RMSE of predicted topography below 0.24 % across the grounding line in West Antarctica (Fig. 8).
We have adopted the shortest temporal resolution suggested in the literature to date (Larour et al., 2019) for the benchmark sea-level simulation in our analysis, but given the complexity of Earth structure and ice dynamics in this region, further exploration with a coupled ice-sheet -sea-level model will be needed to rigorously assess the necessary coupling time interval needed to simulate ice-sheet evolution in marine sectors of the AIS.

640
In this study, we have presented a new time window algorithm in a global sea-level model and provided sample time window parameters for applications to, global glacial-cycle ice-sheet evolution and rapid marine ice sheet retreat in a region with weaker Earth structure. In addition to these applications, the time window algorithm has the potential to unlock opportunities to tackle a range of questions using coupled ice-sheet -sea-level modelling, such as evaluating shorelines during and since the warm mid-Pliocene (3

645
Ma; Raymo et al., 2011;Pollard et al., 2018), investigating the effects of short-term interactions between ice sheets, sea level and the solid Earth on the dynamics of the marine-based portion of Eurasian Ice Sheet during the last deglaciation phase (e. g., Petrini et al., 2020) and the associated impact on abrupt or episodic global sea-level events such as MWP-1A (e.g., Harrison et al., 2019) and understanding the dynamics of ice sheets during past warm interglacial periods (e. g., Clark et al., 2020). Finally, the 650 improved computational feasibility with the time window could allow for ensemble simulations of coupled ice-sheet -sea-level dynamics for the future under different warming scenarios, which will provide useful insight into projected future sea-level hazard.     30 ky and 90 ky with dt1 = 0.2 ky, dt2 = 1 ky and dt3 = 5ky, and TW profile 3 also applies three internal time steps, each of which covers 5 ky, 60 ky and 60 ky with dt1 = 0.2 ky, dt2 = 1 ky and dt3 = 5ky, respectively. Note that all three profiles all assign the first internal time window (L_ITW1) to a 5-ky length with the internal time step (dt1) Figure. 7h), standard simulation with dt = 5 yr (red lines), dt = 10 yr (magenta lines), and dt 760 = 50 yr (blue lines). Note the change in the y-axis in (e). The sea-level code and data can be found in https://osf.io/8ptfm/.

Author contribution
HKH has developed and implemented a novel sea-level model time window algorithm, performed numerical simulations and analysis and wrote the manuscript. NG provided motivation and vision for the 770 work and contributed to manuscript revisions, JXW contributed to discussions of results during the COVID-19 pandemic and isolation. All authors contributed to designing the numerical experiments.