Formulation, calibration and validation of the DAIS model (version 1), a simple Antarctic Ice Sheet model sensitive to variations of sea level and ocean subsurface temperature

15 The Dcess Antarctic Ice Sheet (DAIS) model is presented. Model hindcasts of Antarctic Ice Sheet (AIS) sea level equivalent are forced by reconstructed Antarctic temperatures, global mean sea level and high-latitude, ocean subsurface temperatures, the latter calculated using the Danish Center for Earth System Science (DCESS) model forced by reconstructed global mean atmospheric temperatures. The model is calibrated by comparing such hindcasts for 20 different model configurations with paleoreconstructions of AIS sea level equivalent from the last interglacial, the last glacial maximum and the mid-Holocene. The calibrated model is then validated against present estimates of the rate of AIS ice loss. It is found that a highorder dependency of ice flow at the grounding line on water depth there is needed to capture the observed response of the AIS at ice age terminations. Furthermore it is found that a 25 dependency of this ice flow on ocean subsurface temperature by way of ice shelf demise and a resulting buttressing decrease is needed to explain the contribution of the AIS to global mean sea level rise at the last interglacial. When forced and calibrated in this way, model hindcasts of the rate of present day AIS ice loss agree with recent, data-based estimates of this ice loss rate. 30


Introduction
The Antarctic Ice Sheet is a major player in the Earth's climate system and is by far the largest depository of fresh water on the planet. Ice stored in the AIS contains enough water to 35 raise sea level by about 58 m and ice loss from Antarctica contributed significantly to sea level high stands during past interglacial periods (Vaughan et al., 2013;Kopp et al., 2009;Naish et al., 2009). There is considerable uncertainty as to the amount the AIS will contribute to future sea level change in response to ongoing global warming (Church et al., 2013). 40 A broad hierarchy of AIS models have been developed and applied to try to understand the workings of the AIS and to form a robust basis for future projections of the AIS contribution to sea level change (e.g. Huybrechts, 1990;Huybrechts and de Wolde, 1999;Oerlemans, 2003;De Conto, 2009: Whitehouse et al., 2012). In some cases, AIS models have been coupled to global climate models (e.g. Pollard and DeConto, 2005;Vizcaino et al, 45 2010). A common feature of many of these models is an increase in AIS ice mass in response to warming. This is a consequence of increased snow fall as warming leads to more precipitable water in the atmosphere that falls as snow for cold Antarctic temperatures. However, observations show that the AIS is losing mass at present (Vaughan et al., 2013).
Two recent advances in our understanding of ice flow at the grounding line have been paving 5 the way for AIS model improvements. First, a detailed study comparing results from boundary layer theory with high resolution numerical modeling showed that ice flux at the grounding line increases sharply with ice thickness there (Schoof, 2007). Second, observations show that ice stream flow increases substantially when an adjoining ice shelf disintegrates, removing associated buttressing of the ice stream . Ice 10 shelf disintegration appears to be mainly associated with increased basal melting from increasing subsurface temperatures of the adjoining ocean (Shepherd et al., 2004). This chain of processes has been proposed as a trigger for Heinrich events, explaining why they occur during cold phases of millennial-scale climate variations when the North Atlantic warms at intermediate depths due to shutdowns of the Atlantic Meridional Overturning Circulation 15 (Shaffer et al., 2004). There is considerable support for this interpretation in the ocean sediment record (Marcott et al., 2011).
Here I take a simple modeling approach with recent work by Johannes Oerlemans as my point of departure (Oerlemans 2003(Oerlemans , 2004(Oerlemans and 2005. First I consider the mass balance 20 formulations in these publications whereby I correct some errors and make several parameter value adjustments that follow from the consequences of the corrections. Then I consider formulations within this modeling context of the two recent advances discussed above. The DAIS model is then forced by reconstructed time series of Antarctic temperature, global sea level and ocean subsurface temperature over the last two glacial cycles. Values for the 25 parameters used in the model formulations of the effect on ice flux of grounding line ice thickness and basal melting are then calibrated by comparing model hindcasts with reconstruction targets from the last interglacial period, the last glacial maximum and the mid-Holocene. Finally the calibrated DAIS model is validated against observations of recent AIS ice loss and the future applicability of the calibrated and validated model is discussed. 30

Model formulation and characteristics
The DAIS model builds upon the simple Oerlemans AIS model (Oerlemans 2003(Oerlemans , 2004(Oerlemans and 2005, referred henceforth as O3, O4 and O5). The point of departure here is O5; a more 35 detailed description is found in O3.

Mass balance
The Oerlemans model considers the mass budget of an axi-symmetrical ice sheet with ice 40 sheet radius R resting on a bed with a constant slope, s, before ice loading. The undisturbed bed profile, b, is where r is the radial coordinate and b 0 is the (undisturbed) height at the center of the continent (all model parameters and their standard values are listed in Table 1). Ice loading depresses the bed to immediate isostatic equilibrium. Constant stress is assumed at the ice sheet base leading to a parabolic profile for the ice sheet surface height h in this perfectplasticity limit: 5 where μ is a profile parameter related to ice stress (O3). Ice sheet evolution follows from conservation of mass: 10 where V is ice volume, B tot is the total mass accumulation rate on the ice sheet and F is the total ice flux across the grounding line. The mean annual air temperature reduced to sea level 15 and averaged over Antarctica, T a , and sea level, SL, (relative to its 1961-1990 mean) are the only forcings in the original model. As in O5 I take present day T a = -18 °C (in the following present day refers to a 1961-1990 mean). The second term on the right hand side of Eq. (3) is only considered for a marine ice sheet, that is for R > r c where r c is the distance from the continent center to where the ice sheet enters the sea. From Eq.
(1) , r c = (b 0 -SL)/s, I adopt 20 the O5 approximation of equating the distance from the continent center to the grounding line with the ice sheet radius. Fig. 1 shows a cross-section of a steady state solution of the model with standard parameter values and present day forcing.
The mass balance B at any height on the ice sheet surface is specified as 25 where h R (T a ) is the height of the runoff line above which precipitation, P(T a ), is assumed to accumulate as snow and β(P) is a rate of mass balance increase with height. Sublimation has 30 been disregarded here. The value of β is depends on a combination of (1) height variations of precipitation and (2) height variations in summertime melting coupled to the atmospheric lapse rate (Oerlemans, 2008). The problem can be easily reformulated in terms of an equilibrium height, h e , where the yearly-mean mass balance is zero: from Eq. (4) Likewise total accumulation and total ablation can be calculated by 35 integrating B over the ice sheet surface above and below h e , respectively. Results from such a calculation will be presented below. However, integration of B over the whole surface to obtain B tot is more straightforward for obtaining problem solutions (O3).
The mass balance part of the problem is closed by choosing appropriate expressions for h R , P 40 and β. For the runoff line height, where h 0 is the runoff line height for T a = 0. The values in Table 1 for h 0 and the coefficient c were taken from mass balance studies (O4). For precipitation, where P 0 is the (ice equivalent) precipitation at 0 °C and the value of the coefficient κ is chosen assuming a relationship between columnar water content (that increases exponentially with temperature) and precipitation (O5). No attempt is made here (nor in O5) to deal with decreasing precipitation toward the center of the continent. This feature was included in O4 but at the cost of one extra free parameter. For the mass balance gradient β, 10 2 1 P ν β = where this relationship and the value of the parameter ν in O5 were also based on mass balance observations that reflect in part larger vertical precipitation gradients for greater 15 precipitation (Oerlemans, 2008). Integration of B over the ice sheet surface yields 2 3 2 This correction to the original work leads to much reduced runoff and, for the O5 parameter 30 values, an AIS model much less sensitive to climate at warm temperatures. In the corrected model, complete Antarctic deglaciation occurs at a T a of about 8 °C higher than in the original model (O4) and about 4 °C higher than in a 3-D thermomechanical model ( Fig. 2; Huybrechts, 1993). One possible explanation for this behavior is the β value of ~ 0.001 m ice m -1 yr -1 for O5 parameter values (calculated at 0 °C). This is a low value for the balance 35 gradient compared to observations, even for dry polar climates (Oerlemans, 2008). To address this I doubled the O5 value of ν from 0.006 to 0.012 m -0.5 yr -0.5 yielding a climate sensitivity at warm temperatures similar to the 3-D model (Fig. 2). Detailed agreement of the idealized DAIS model with the 3-D model cannot be expected since the latter includes a representation of real Antarctic topography. Further improvements/adjustments of the model 40 runoff formulation/calibration will be undertaken in future work. For the present I will concentrate on revised formulations of the other main component of the model: ice flux across the grounding line. This has been the dominant process for ice loss from Antarctica over ice age cycles and will continue to be so in the near future (Pollard and DeConto, 2009). 5

Ice flux at the grounding line
In the following I retain the above O5 mass balance treatment (but with the sign corrections) as well as all O5 parameter values except for the revised value of ν and a slight increase of b 0 from 760 m to 775 m. As in O5, the parameter values are chosen to reproduce present day 10 AIS volume, area, mean surface elevation and mass throughput ( Table 1). The ice flux at the grounding line, F, is, where ρ w and ρ i are water and ice densities, H is the water depth, and S is the ice speed. H and S should be thought of as some average around the ice sheet periphery over all ice streams. To the grounding line/ice sheet radius approximation mentioned above, The ice flux problem was closed in O5 by assuming that the ice speed is related linearly to the water depth: where the value of constant of proportionality f 0 was chosen to reproduce present day AIS throughput (Table 1). 25 The two recent advances in our understanding of ice flow at the grounding line discussed in the Introduction call for a revision of the model ice speed formulation. First, I take ice flux at the grounding line to depend on water depth there raised to a power (Schoof, 2007). Second, I take melting at the base of a marine ice shelf to depend on the temperature difference between subsurface temperature of the adjacent ocean and temperature at the ice shelf base, 30 the latter anchored to the freezing temperature of sea water. Modeling studies have yielded a range of such dependencies from linear to quadratic (e.g. Williams et al., 2002;Holland et al., 2008). Below I am guided mainly by the comprehensive, 3-D ocean general circulation model study of Holland et al. (2008) who find a quadratic dependency on this temperature difference for a wide range of shelf-slope topographies: They find the melt rate is found to be 35 proportional to the product of ocean flow speed and ocean temperature beneath the ice shelf, both of which increase linearly with ocean warming. However I also touch upon results for a model version with a linear dependency on this temperature difference.
With the above motivations I choose to formulate the ice speed as follows : 40 [ ] where T o is the ocean subsurface temperature adjacent to the AIS, α is a partition parameter and γ is the power for the relation of ice flux to water depth. Furthermore, T f is the freezing temperature of sea water, T o,0 is a present day reference T o and R 0 is a reference R. In the following, the parameter α will be taken to vary from 0 to 1 (but values become unphysical toward the upper limit of this range as there will be flow at the grounding line even when T o 5 approaches T f ). The parameter γ will be taken to vary from 1/2 to 19/4. This can be compared to an O5 value of 1 and a value of 15/4 from the results of Schoof (2007; his preferred value for p is 19/4 where F H p , therefore from Eq(s).. (9) and (11), γ = p-1). The value for T o,0 was obtained as described in Appendix A and the value for R 0 was taken from the steady state solution of the O5 model (α = 0; γ = 1) with present day forcing (T a = -18 °C; SL=0; see 10 below and Table 1). The scaling in Eq. (11) by present day ocean-ice shelf temperature difference and water depth at the grounding line anchors the present day ice speed values to the reference solution ice speed for all α and γ values. In the spirit of the original work (O3; O4; O5), the ice speed formulation in Eq. (11) is meant to capture the bulk effect of all individual ice streams and associated embayed ice shelves around the AIS periphery. 15 The conservation of mass is now The ice sheet is now also forced by ocean subsurface temperature T o that is determined by a combination of local and remote conditions. Model formulation in terms of the independent variable R is completed by relating ice volume to R, taking into account immediate isostatic adjustment and the effect on this adjustment of the displacement of sea water by ice (O3): 25 where ε 1 is ρ i (ρ m -ρ i ) -1 , ρ m is rock density and ε 2 is ρ w (ρ m -ρ i ) -1 . The term multiplied by ε 2 in Eq. (13) is only considered for a marine ice sheet (R > r c ). Finally from Eq. (13), 30 where the terms multiplied by ε 2 are only considered for a marine ice sheet. Note the sign correction in the last term in curly brackets as well as the extra term at the end of the equation when compared with the original derivation in Eq. (13) of O3. 35 . Present day ice volumes for transient model solutions slightly exceed this steady state value since these solutions are still responding to past temperature and sea level rises. This will be discussed in more detail in Section 3. Fig. 3a shows the steady state distribution of AIS ice volume for the original (but corrected) 5

Steady state solution properties
O5 model with Table 1 parameter values. Each of the solutions used in the figure was obtained by a 100 kyr integration of the above time dependent model equations. Each such integration of this semi-analytical model with a one year time step takes a fraction of a second on a personal computer. Steady state ice volumes increase for warming up to 5 -7 °C above present day, due to increased snow fall and accumulation. Volumes decrease rapidly 10 for still warmer temperatures as melting becomes important. For very warm temperatures, ice volume isolines become horizontal as the ice sheet recedes out of the ocean, eliminating the dependency on sea level. Ice volume is greater during maximum glacial conditions with T a about 10 °C colder and SL about 130 m lower than present. Due to the reduced snow fall, increased ice sheet area is needed to balance the ice flux at the grounding line. 15 Fig. 3b shows the steady state distribution of AIS ice volume for a preferred model configuration from the calibration in Section 4 (case 4 of Table 2; γ = 2, α = 0.35 in Eq. (11)). This configuration exhibits enhanced ice flux at the grounding line from (1) a higher-order ice flow dependency on water depth there and (2) ice flow increase from ice shelf demise by 20 way of basal melting. In this case high latitude, ocean subsurface temperature, T o , also comes into play. For the calculations upon which this figure is based I related T o to T a by applying a second-order polynomial fit to values of these temperatures from reconstructions over the past 240,000 years (Appendix A). A best fit with a RMSE of 0.16 °C was found for: For example, this yields T o = -0.50, 0.72 and 3.32 °C for T a = -28, -18 and -8 °C, respectively. During colder periods, greater model dependency on sea level leads to less relative dependency of ice volume on temperature (isolines more vertical). Furthermore, ice 30 volume now decreases with warming above present day values. The increase in ice flux at the grounding line from warmer ocean subsurface temperatures, increased basal melting and increased ice shelf demise balances increased snow fall and accumulation from warmer atmospheric temperatures without the need for ice sheet growth. Rather, ice sheet contraction is now needed for equilibrium as temperatures warm. 35 Fig. 4 provides a closer look at ice conservation terms for steady state solutions of this preferred model configuration. The sum of total accumulation (a), total ablation (b) and total ice flux at the grounding line (c) add up to zero for each T a , SL combination. As mentioned above, there is a balance between total accumulation and total ice flux for T a less than -15. 5 40 °C. Likewise, there is a balance between total accumulation and ablation for T a greater than about -1 °C, at which point the ice sheet has receded out of the ocean. Total accumulation is greatest for T a of about -5.5 °C, the temperature around which ablation and ice flux are of comparable importance. Ablation and ice flux are greatest for T a of about -2 °C and -13 °C, respectively. All these terms increase as sea level decreases illustrating the effect of increased 45 ice sheet area and circumference for a lower sea level.

Model calibration and validation
I calibrate the DAIS model by comparing sea level equivalent (SLE) hindcasts from growing and shrinking of the modeled AIS over the last two glacial cycles to paleoreconstructions that provide SLE constraints. In particular I consider the period from 240 kyr BP up to the year 5 2010 AD (I take 0 BP to be 2000 AD). The first step is to construct credible time series over this time period for the three model forcings, T a , SL and T o . The forcing time series used here are shown in Fig. 5. A detailed description of how these series were constructed is given in Appendix A. The second step is to choose time slices for which suitable paloereconstructions are available. For this first calibration I choose to work with the last interglacial (LIG), the 10 last glacial maximum (LGM) and the mid-Holocene (HOL) at about 6000 BP.
Maximum LIG sea level was considerably higher than present. Recent estimates converge on a range of about 6 -9 m above present for this maximum (Kopp et al., 2009;Dutton and Lambeck, 2012 (Vaughan et al., 2013); their contribution to LIG sea level rise was probably less than this given that, during the LIG, global mean temperatures were no more than 2 °C above present 20 day (Fig. A1). Taken together, these two sources can explain a LIG sea level rise of at most 1 m. The Greenland contribution to LIG sea level rise may have been in the range 2.0 -2.5 m (Kopp et al, 2009;NEEM, 2013). The sum of these contributions leaves a remaining 2. 5 -5.5 m to be explained by ice loss from Antarctica. I adopt this range as my LIG constraint. 25 The AIS was larger during the LGM. Ice sheet and glacial isostatic adjustment (GIA) modeling a decade ago indicated a range of 14 -21 m SLE for this size increase (Clark and Mix, 2002;Peltier, 2004;Huybrechts, 2002). However, more recent modeling and GIA studies using local GPS observations show lower values of 8 -10 m (Ivins and James, 2005;Whitehouse et al. 2012). This issue needs to be resolved but for present purposes I assume 30 the range of 8 -17 m SLE for my LGM constraint. By the mid-Holocene, the Northern Hemisphere ice sheets had melted and mean sea level had risen to about 2 -3 m below present (Lambeck et al., 2010). Since temperatures had been slightly warmer than present for thousands of years by the mid-Holocene (Marcott et al., 2013), the ocean may have been warmer and less ice than present may have been found in mountain glaciers and ice caps and 35 perhaps also on Greenland (Vinter et al., 2009). These effects would have raised sea level to or slightly above present levels. Taken together with the mean sea level estimate, this implies an AIS size then about 2 -4 SLE above present. I adopt this range as my HOL constraint.

Hindcasts over the last two glacial cycles 40
Here I integrate the time dependent equations for ice sheet radius, R, from 240 kyr BP to 2010 AD for the forcing in Fig. 5 and for five different model configurations (cases 1 -5 in Table 2). The integrations start from the initial condition of R = R 0 (1.8636 x 10 6 m), the radius of the steady state model solution for present day T a , SL and T o . This is a reasonable 45 initial condition for 240 kyr BP during an interglacial period not unlike the present one. Tests with other reasonable initial conditions varying by up to 10% from the above value produced identical hindcasts from about 220 kyr BP onward. With the use of Eq. (13), ice sheet volume was then calculated and converted to SLE: an SLE of 57 m was taken to correspond to the ice volume of the steady state model solution, 24.78 x 10 15 m 3 (assuming ice and sea water densities from Table 1 and seawater replacing ice below sea level). 5 The model hindcast with the original (but corrected) Oerlemans model (case 1) exhibits a slow, low amplitude response to the forcing whereby maximum SLE occurs about 30 kyr after the LIG (red lines in Fig. 6). In this configuration the AIS continues to respond significantly at present to sea level rise over the last deglaciation and would continue to lose mass equivalent to more than 1 m SLE if present day temperature and sea level were 10 maintained (Table 2). A model hindcast with increased sensitivity of ice flow to sea level rise (case 2) shows a more rapid, higher amplitude response, more in accord with the data constraints (maroon lines in Fig. 6). In this case, the timing but not the amplitude of the SLE target at the LIG is achieved. As the model now responds more rapidly to sea level change, much less sea level rise remains from continuing model response to the last deglaciation 15 ( Table 2). A model hindcast with increased sensitivity of ice flow to ocean subsurface temperature (case 3) shows a still higher amplitude, a good agreement with the SLE target at the LIG but a slow response to the last deglaciation (green lines in Fig. 6). Here almost all the AIS ice loss occurs after 10 kyr BP, the AIS ice volume is too large during the mid-Holocene and there is a remaining sea level rise of nearly 2 m for present day forcing ( Fig. 6; Table 2). 20 A model hindcast with both increased sensitivity of ice flow to sea level rise and to ocean subsurface temperature (case 4) meets all three reconstruction targets (blue lines in Fig. 6). The responses at both the LIG and the last deglaciation are now sufficiently fast and large. Maximum SLE during the LIG was 3.07 m. When rerun using a double-peak structure for 25 LIG sea level rise (Kopp et al., 2009), case 4 yields a very similar maximum with a 500-1000 year extension of high SLE. When rerun using the original Waelbroeck et al. (2002) sea level curve across the LIG (Appendix A; Fig. 5), case 4 yields a slightly lower maximum of 2.31 m. The results of cases 3 and 4 show that the key model feature for simulating an LIG-like ice loss is ocean subsurface warming leading to ice flow acceleration. AIS ice loss during the 30 last deglaciation now occurs mainly after 15 kyr BP and includes a SLE ice loss of less than one meter at melt water pulse 1A, forced by ocean subsurface warming and sea level rise across the pulse. This result is consistent with a mainly Northern Hemisphere source of this melt water event (e.g. Gregoire et al., 2012). The case 4 hindcast is also consistent with a reconstruction of the 220 kyr BP sea level high stand several meters below present day sea 35 level (Waelbrocke et al., 2002). A model hindcast with higher sensitivity of ice flow to sea level rise and to ocean subsurface temperature (case 5) barely meets all three reconstruction targets (cyan lines in Fig. 6). This rapid response of this configuration to forcing leads to a more rapid AIS shrinking at the last deglaciation. Furthermore, by the Medieval Warm Period around 1000 years ago, the model AIS has receded to a smaller size than present day. 40 Fig. 7 shows the results of a more systematic search for the degree to which increased sensitivity of ice flow to sea level rise (as characterized by the parameter γ) and to ocean 45 subsurface temperature (as characterized by the parameter α) lead to model hindcasts that satisfy all three reconstruction targets. From top to bottom panes, the figure shows isolines of Antarctic Ice Sheet SLE for the last interglacial, the last glacial maximum and the mid-Holocene, respectively, for model hindcasts spanning the ranges of γ and α. The five cases considered in section 3.1 are plotted as colored dots with the coloring scheme as in Fig. 6. The isoline ranges of SLE defined by the reconstruction targets are shaded. Of the 336 hindcasts tested, 29 meet all three reconstruction targets (black dots in the figure). 5

Comparison with past and present day estimates
Acceptable values for α are 0.25 -0.45 with the lower values constrained mainly by the LIG target and the higher values mainly constrained by the LGM target. Acceptable values for γ are 1 -3.75 and are mainly constrained by the HOL target. However most of the acceptable solutions are grouped in the γ range of 1.75 -3. These values are somewhat lower than expected from boundary layer theory (3.75; see above) but it should be remembered that they 10 represent some mean of all ice streams around Antarctica with varying degrees of ice shelf buttressing. The low value outliers for γ are coupled to low value outliers for α.
The above results show some promise for the use of such well-calibrated versions of DAIS in other future applications. After all, the model captured well the timing and amount of ice loss 15 during the warmer-than-present LIG and evidence has been growing for the importance of the process by which this loss occurred in the model: ice flow increase driven ultimately by ocean subsurface warming (Pritchard et al., 2012). However, the well-calibrated model versions should also pass verification tests whereby their hindcasts are compared with reconstructions outside of the calibration interval. Much recent effort has gone into 20 quantifying the rate of ongoing sea level rise due to ice loss from Antarctica. The most recent IPCC estimate of this for the period 1993 -2010 AD is 0.27 ± 0.11 mm yr -1 (Vaughan et al., 2013). Fig. 8 shows isolines of this rate for this period as calculated from the same model hindcasts as in Fig. 7. The isoline range defined by the recent IPCC estimate is shaded and the well-calibrated model configurations from Fig. 7 are plotted in Fig. 8. There is a large 25 overlap between the rate of ongoing sea level rise from AIS ice loss for these model combinations and the IPCC estimate. For example, the rise for case 4 is 0.24 mm yr -1 (Table  2). Only a few well-calibrated configurations with low values of α and high values of γ fall outside the IPCC range. 30 I also tested a model version identical to the present one except for a linear, rather than quadratic, dependence of basal melting on the temperature difference between the subsurface temperature of the adjacent ocean and the temperature at the ice shelf base. With this version it was also possible to satisfy all three reconstruction targets but in a much narrower parameter space with acceptable values for α of 0.65 -0.7 and for γ of 1.5 -2. These few 35 acceptable hindcasts produced rates of ongoing sea level rise from AIS ice loss in the range of 0.24 -0.31 mm yr -1 , also in good agreement with the recent IPCC estimate.

Discussion 40
Here I formulated a simple Antarctic Ice Sheet model by first adopting and making corrections to a published mass balance approach (Oerlemans 2003(Oerlemans , 2004(Oerlemans and 2005 and then by advancing a new treatment of ice flow at the grounding line based on recent advances in our understanding of the controls on this flow. It was then possible to calibrate the resulting DAIS model, as forced by reconstructed Antarctic temperature, sea level and ocean 45 subsurface temperature, to simultaneously satisfy reconstructions of AIS contributions to sea level during the last interglacial, the last glacial maximum and the mid-Holocene. Finally it was shown that the well-calibrated model hindcasts also reproduced best estimates of the present rate of ongoing sea level rise from AIS ice loss. These results lend some support for the use of the present DAIS model for hindcasts farther back into the past as well as for projections into the future. However the question arises of how much confidence can be placed in the results of such a simple model. On one hand, well-calibrated simple models 5 have proven useful in many contexts in the past (Shaffer et al. 2009;Meinhausen et al. 2009). On the other hand there are key questions regarding the AIS that cannot be addressed in detail with the simple DAIS model. Perhaps the most emblematic of such questions regards the collapse of the West Antarctic Ice Sheet. 10 Analyses of sea level rise during the last interglacial period indicate that during this period AIS ice loss compared to present day was 2.5 -5.5 m in sea level equivalent (Kopp et al., 2009;Dutton and Lambeck, 2012). The lower end of this estimate coincides with a recent estimate of how much sea level would rise from the collapse of the present day WAIS (3.3 m; Bamber et al., 2009.). The upper end of this estimate would require some additional 15 contribution from East Antarctica. Basal melting in combination with a high order dependency of ice flow on grounding thickness facilitate marine ice-sheet instability, such as would likely be active in a WAIS collapse, but would also increase flow in East Antarctic ice streams (Schoof, 2007;Joughin et al., 2012). Both basal melting and the high-order ice flow dependency were included, albeit schematically, in the DAIS model and made it possible for 20 this simple model to simulate observed AIS ice loss during the LIG when forced by reconstructed ocean subsurface temperatures. Marine ice-sheet instability also requires complex bed geometry with the bed sloping downwards in the inland direction. As shown in Fig. 1, the depressed bed of the DAIS model slopes in this manner at the periphery of the ice sheet. However, assumed immediate isostatic adjustment in this highly-idealized model 25 immediately restores the initial upward slope. While the model has been calibrated to produce observed LIG ice loss of the size of a WAIS collapse when forced with some realism, it can therefore not be expected to capture detailed timing of a marine ice sheet instability nor high, short-term ice loss rates that would be associated with one. Still, when calibrated with paleoconstraints the model was successful in reproducing present, short term 30 rate of ice loss from the AIS.
Recently a much more complex and complete model of the AIS was used to address past WAIS collapses (Pollard and DeConto, 2009). A long hindcast of this model showed considerable skill in capturing the timing of such collapses when compared to the ocean 35 sediment record. However, this hindcast also showed a large collapse in the next-to-last interglacial and a smaller one in the last interglacial, in conflict with sea level reconstructions (Waelbroeck et al, 2002: Kopp et al, 2009). The well-calibrated hindcasts of the DAIS model presented above (Fig. 6)  applications. This very fast, semi-analytical model is well suited for very long hindcasts that would not be feasible with 3-D models. One such possibility would be a study of AIS emergence and subsequent evolution over Cenozoic time. The DAIS model is also well positioned to be a component in integrated assessment modeling by allowing extensive sensitivity studies, using Monte Carlo-type analyses for example (Hargreaves and Annan, 2002;Applegate et al., 2012). After all, more than 10 million, 1000 year DAIS model simulations can be run in a day on a single computer processor. Such a study would 5 determine the relative importance of and provide Probability Density Functions for the parameters in the model (Table 1). It should be remembered that the model has been calibrated and validated here for conditions no warmer than the last interglacial for which ice flux at the grounding line is essentially the only ice loss term. If future warming becomes sufficiently strong, ablation will also become an important factor in AIS mass balance. Work 10 is underway to apply the model to the Greenland Ice Sheet for which ablation has been an important ice loss term during interglacial periods (Applegate et al., 2012). Such an application will allow improved calibration of the model parameters controlling ablation.

Code availability 15
Matlab codes for the DAIS model, the 240 kyr forcing data and instructions for using the model can be downloaded from www.dcess.dk under "DCESS model".
Acknowledgments. I thank the following people for discussions and/or for supplying data: For annual and area mean Antarctic temperature reduced to sea level, T a , I use temperature anomaly reconstructions adjusted to a present day mean temperature of -18 °C (present day refers to a mean of the period 1961-1990 AD). For the period 1500 -240000 BP, I adopt 30 temperature anomaly estimates from the Dome C ice core (Jouzel et al., 2007) as referenced to 1.2 times the mean global temperature anomaly in the period 500 -1500 BP from Mann et al. (2008). The factor 1.2 is an estimated interglacial polar amplification factor for Antarctica. This referencing led to a small adjustment of -0.1 °C relative to the published anomalies. For the periods 500-1850 AD and 1851 -2010 AD I used the global mean 35 temperature anomaly from Mann et al. (2008) and Morice et al. (2012), respectively, multiplied by 1.2. These two time series have already been referenced to present day global means. The composite time series for T a was then interpolated to one year time steps (Fig. 5).
There are no detailed and reliable time series for sea level around Antarctica so I fall back 40 upon global mean estimates. Reconstruction of a specific Antarctic sea level time series would require knowledge of the sources of sea level change over the two glacial cycles and application of non-eustatic corrections for changes in ice mass (e.g. Mitrovica et al, 2001). This is beyond the scope of the present paper. For the period 21 -240 kyr BP, I adopt SL values from Waelbroeck et al. (2002) but with an adjustment for the last interglacial period 45 (LIG). In this work, the transition of LIG sea level to levels above present day took place at about 124 kyr BP. However, subsequent work put that transition earlier, at 126 ± 1.7 kyr BP (Waelbroeck et al., 2008) and around 130 kyr BP (Dutton and Lambeck, 2012). I adjust the LIG transition to occur at about 127.5 kyr, based on the more recent work, and also adjust the Waelbroeck et al. (2002) curve by a corresponding amount back to the preceding glacial maximum (Fig. 5). For the period 7000-21000 BP, I adopt SL values from Clark et al. (2012) that include representations of melt water events during the last deglaciation. For SL in the 5 period 6000 BP to 1869 AD, I use the curve of Lambeck et al. (2010). For 6000 -7000 BP, I interpolated linearly between values from the above two sources. Finally, for 1870 -2010 AD, I adopt SL values of Church and White (2011) adjusted to a present day value of zero. The composite time series for SL was then interpolated to one year time steps (Fig. 5). 10 High latitude, ocean subsurface temperature, T o , enters the model as a forcing that ultimately modulates ice flow at the grounding line (see Section 2). It is a major challenge to construct a realistic time series for T o since it will depend on some combination of local and remote processes like wind-driven upwelling/downwelling and Atlantic Meridional Overturning Circulation strength. Furthermore, T o will vary around Antarctica. For this task I take a very 15 simple approach in tune with the scope of the present paper and model. First I construct a global mean atmospheric temperature anomaly (GMATA) time series from 240000 BP until 2010 AD. This time series is relative to the mean temperature of the reference period 1961-1990 AD and 0 BP is taken to be 2000 AD. The GMATA series is based mainly on Antarctic and Greenland ice core data before 1500 BP (see details below), adopts a global 20 reconstruction afterward (Mann et al., 2008) and finally adopts a reconstruction from direct temperature observations after 1850 AD (Morice et al., 2012;Fig. A1a). Then I force the low-mid (0 -52°) and high (52 -70°) latitude DCESS model ocean with corresponding atmospheric temperature constructed from the global series using amplification factors (0.928 and 1.266, respectively) and mean reference period temperatures (20.55 and -4.39 °C, 25 respectively), all taken from DCESS model simulations (Shaffer et al., 2008;Fig. A1b). In the model, sea ice extent is diagnosed from zonal profiles of atmospheric temperature and only the sea ice-free part of the high latitude ocean exchanges heat with the atmosphere. The DCESS model ocean with 100 m vertical resolution includes parameterized overturning, horizontal mixing, small-scale vertical mixing in the low-mid latitude zone and enhanced 30 vertical mixing in the high latitude zone. Parameter values were calibrated by fitting to ocean temperature and carbon 14 observations (Shaffer et al., 2008). The model has only one high latitude zone but the poleward extent of that ocean zone matches the equatorward extent of Antarctica. The T o used for the forcing of the DAIS model is then taken to be the average temperature in the depth range 200 -800 m of the high latitude ocean zone ( Fig. A1b and  35 Fig . 5). The value for T o,0 in Table 1 is the mean T o for the period 1961-1990.
For the period 1500 -122400 BP temperature estimates are available from Antarctic and Greenland ice cores (Jouzel et al., 2007;Masson-Delmotte et al., 2006). The GMATA time series for this period is calculated as a mean of the two estimates after each is (1) referenced 40 to the 1961-1990 period by referencing to the mean temperature anomaly in the period 500 -1500 BP from Mann et al. (2008) and (2) divided by appropriate polar amplification factors. These factors were found to be 3 and 6 for Antarctica and Greenland, respectively, by fitting to the Shakun et al. (2012) temperature reconstruction for the period 11000 -22000 BP after referencing that reconstruction to the 1961-1990 period by referencing to a common period 45 with the results of Marcott et al. (2013) (Fig. A1a). Lower amplification factors of 1.2 and 2.4, respectively, were adopted for the interglacial part of the 1500 -122400 BP period.
The temperature estimates for the Greenland ice core used above are not available before 122400 BP. GMATA time series were constructed as above for earlier times but using other estimates for Greenland temperature. For the period 128700 -240000 BP, this temperature was estimated from a proxy based on methane data from an Antarctic ice core (Spahni et al., 5 2005). Global atmospheric methane increases during glacial times are due largely to enhanced emissions from boreal wetlands in response to boreal warming (Fischer et al., 2008). A linear regression of Greenland temperature on Antarctic methane for the period 150 -122400 BP results in the relation T = -51.5 +0.0802 [CH 4 (ppb)]. I also tried several more sophisticated approaches included non-linear regression and applying the simple DCESS 10 model module for atmospheric methane (Shaffer et al., 2008) but found no significant improvement over the simple linear regression in terms of modeled versus observed Greenland temperature for this calibration period.
A different approach was required for the remaining period of 122400 -128700 BP 15 encompassing Termination II and much of the last interglacial as little guidance can be found in complex methane -Greenland temperature relationships at Termination I into the present interglacial (Spahni et al., 2005;Masson-Delmotte et al., 2006). A detailed study of Termination II from a variety of climate archives indicates that Greenland warming lags that of Antarctica with rapid warming commencing around 128500 BP in the northern North 20 Atlantic and reaching full interglacial levels by about 127000 BP (Masson-Demotte et al., 2010). Guided by these results I extended the interglacial, Greenland ice core temperature at 122400 BP back to 127000 BP and applied a linear interpolation between that value at that time and the (methane-based) temperature at 128700 BP. 25  (Huybrechts, 1993).  (Table 1; case 1 in Table 2). b. Solutions for a 5 preferred model configuration (case 4 of Table 2) that also includes sensitivity to high latitude, ocean subsurface temperature (T o ) as calculated from Eq. 15. The black crosses mark solutions for present day T a , SL and T o (Table 1). Figure 4. Sensitivity of total equilibrium volume fluxes (10 12 m 3 ice yr -1 ) to Antarctic temperature, sea level and high latitude, ocean subsurface temperature for the preferred model configuration (case 4 in Table 2). The volume fluxes are a. Accumulation, b. Ablation 5 and c. Ice flux at the grounding line. The black crosses mark solutions for present day T a , SL and T o (Table 1). 15 Figure 6. Hindcasts of sea level equivalent (SLE) from changes in AIS ice volume for the five model configurations in Table 2 and for the period 240 kyr BP to 2010 AD. The hindcasts are plotted relative to present day values (means for 1961-1990 AD). Shown are 5 results for the original (but corrected) Oerlemans (2005) model (case 1; red lines), for increased sensitivity of ice flow to sea level (case 2, maroon lines), for increased sensitivity of ice flow to ocean subsurface temperature (case 3, green lines), for increased sensitivity of ice flow to sea level and ocean subsurface temperature (case 4, blue lines) and for high sensitivity of ice flow to sea level and ocean subsurface temperature (case 5, cyan lines). 10 Frames bd are enlarged sections of the full hindcasts shown in a. Paleoreconstruction targets for the last interglacial, the last glacial maximum and the mid-Holocene are shown as vertical bars in ac, respectively (see text in section 3 for details). Also shown is reconstructed global mean sea level from 6000 BP to the present (black dashed lines in c and d; see Appendix A for details). 15 Figure 7. Model hindcasts in specific time slices of sea level equivalent (m) from AIS ice volume changes, as functions of model parameters γ and α. These parameters enter in the dependency of ice speed on water depth at the grounding line and on ocean subsurface 5 temperature, respectively (see Eq. 11). The time slices are a. the last interglacial period, b. the last glacial maximum and c. the mid-Holocene. Reconstruction targets for each time slice are shaded. The black dots mark hindcasts that meet all three paleoreconstruction targets. The colored dots mark the five hindcasts of Table 2 and Fig. 6 with the Fig. 6 coloring scheme. 5 10 Figure 8. Model hindcasts of the mean rate of sea level rise (mm yr -1 ) from ice loss from the AIS for the period 1993 -2010 AD as functions of model parameters γ and α. The present IPCC best estimate for the range of this mean rate is shaded (Vaughan et al., 2013). The black dots mark hindcasts that meet all three paleoreconstruction targets (from Fig. 7). The colored dots mark the five hindcasts of Table 2