Grounding line transient response in marine ice sheet models

Marine ice sheet stability is mostly controlled by the dynamics of the grounding line, i.e., the junction between the grounded ice sheet and the ﬂoating ice shelf. Grounding line migration has been investigated in the framework of MISMIP (Marine Ice Sheet Model Intercomparison Project), which mainly aimed at investigating steady state solu- tions. Here we focus on transient behaviour, executing short-term simulations (200 yr) of a steady ice sheet perturbed by the release of the buttressing restraint exerted by the ice shelf on the grounded ice upstream. The transient grounding line behaviour of four di ﬀ erent ﬂowline ice sheet models has been compared. The models di ﬀ er in the physics implemented (full-Stokes and Shallow Shelf Approximation), the numerical approach, 10 as well as the grounding line treatment. Their overall response to the loss of buttressing is found to be consistent in terms of grounding line position, rate of surface elevation change and surface velocity. However, large discrepancies ( > 100%) are observed in terms of ice sheet contribution to sea level. Despite the recent important improvements of marine ice sheet models in their ability to compute steady-state conﬁgurations, our 15 results question models’ capacity to compute reliable sea-level rise projections. amplitudes of perturbation are investigated with corresponding modiﬁed of C F

. However, implementing GL migration in ice flow models still represents a challenge to be faced by the community of ice sheet modellers (Vieli and Payne, 2005;Pattyn et al., 2012).
As mentioned above, Schoof (2007) developed a boundary-layer theory establishing the relation between ice flux and ice thickness at the GL, which can be implemented as 15 a boundary condition in ice-flow models. The boundary layer is a zone of acceleration, generally 10-20 km in extent (Hindmarsh, 2006), where the stress regime adjusts from being shear-dominated to extension-dominated. This theoretical development demonstrated the uniqueness of steady solutions of marine ice sheets resting on a downward sloping bedrock and their unstable behaviour on an upward sloping region. Based on 20 the Schoof (2007) results, an intercomparison effort compared the behaviour of the GL evolution on a flowline of 26 different models, as part of the Marine Ice Sheet Model Intercomparison Project (MISMIP, Pattyn et al., 2012), which was essentially designed to compare models with the semi-analytical solution proposed by Schoof (2007). However, Schoof's flux formula is derived on the assumption of near-steady-state, and its 25 ability to represent transient behaviour has not been fully investigated. This issue was briefly touched upon during the MISMIP experiments, but it was not the primary focus of investigation.

3905
The MISMIP experiments showed a broad range of behaviour of numerical implementations in response to an instantaneous global change of the ice rheology, with some quantitative consistency between different numerical formulations. The MISMIP experiments highlighted, along with Schoof's studies, the importance of obtaining high accuracy in the numerical solution in the boundary layer near the GL, which in practice 5 means the use of high resolution or high accuracy methods, which has the consequence that the numerical approach used is a significant issue.
Short term predictions of rapid change in the Antarctic Ice Sheet necessarily involve transient processes, and the ability of marine ice sheet models to represent these requires quantification. Therefore, we conduct a model intercomparison dealing with 10 rapid change in order to evaluate the transient behaviour of all models. A particular aim is to investigate the divergence of ice sheet models from the Schoof (2007) solution during these very short time scale processes. Furthermore, owing to the use of different physical approximations and numerical approaches, we expect that the same experiment carried out with different ice sheet models may give different results. Therefore, 15 another aim of this study is to quantify these differences and understand their origin. We choose to investigate the physically more reasonable transient forcing of a decrease in ice-shelf buttressing. This is implemented by means of a plane-flow model with grounded part and a floating ice shelf. As is common with previous studies (Nick et al., 2009;Price et al., 2011;Williams et al., 2012), buttressing is implemented by The first one is the finite element full-Stokes Elmer/Ice model, denoted FS-AG for Full-Stokes -Adaptive Grid, developed at CSC/LGGE. In this application, the adaptive grid refinement proposed by Durand et al. (2009b) is used. This model is computationally two dimensional in this plane-flow representation. The three remaining models solve the SSA, and are therefore vertically integrated and thus computationally one-10 dimensional. SSA-FG (for SSA-Fixed Grid) and SSA-H-FG (for SSA-Heuristic-Fixed Grid) use a fixed grid with a resolution of 50 m and 10 km, respectively. The GL migration of SSA-H-FG is computed according to the Pollard and DeConto (2009) heuristic rule that implements the Schoof (2007) boundary condition. The last model solves the SSA equations using pseudo-spectral method (Fornberg, 1996;Hindmarsh, 2012) on 15 a moving grid, and will be denoted SSA-PSMG for SSA -Pseudo-Spectral Moving Grid. For this model, grounded ice and floating ice shelf are solved on two coupled domains, with continuity of stress and velocity across the grounding-line guaranteed. The first two models approach the problem of modelling the flow in the boundary layer by increased resolution, the third model uses a coarse resolution and a heuristic rule 20 at the GL, and the last model addresses this issue by using high-accuracy spectral methods.
Details and numerical characteristics of the four models are summarised in

Governing equations
The problem consists of solving a gravity driven flow of incompressible and isothermal ice sliding over a rigid bedrock noted b(x). The ice is considered as a nonlinear viscous material, following the behaviour of the Glen's flow law (Glen, 1955): where τ is the deviatoric stress tensor, D is the strain rate tensor defined as D i j = (∂ j u i + ∂ i u j )/2 and u = (u, w) is the velocity vector. The effective viscosity η is defined as follows: 10 where A and n are the Glen's law parameter and flow law exponent respectively, and D e is the strain-rate invariant defined as D 2 e = 2D i j D i j . The ice flow is computed by solving the Stokes problem, expressed by the mass conservation equation in the case of incompressibility and the linear momentum balance equation where σ = τ −pI is the Cauchy stress tensor with p = −trσ /3 the isotropic pressure, ρ i the ice density and g the gravity vector. Both the upper ice/atmosphere interface z = z s (x, t) and the lower ice/bedrock or 20 ocean interface z = z b (x, t) are allowed to evolve following an advection equation: where (u i , w i ) is the surface velocity (i = s) or the basal velocity (i = b). For this application, the mass flux at the surface is constant and uniform (a s (x, t) = a s , see Table 2) and a b = 0.

Boundary conditions
The geometry is restricted to a two-dimensional plane flow along the x-direction and 5 the z-axis is the vertically upward direction. The upstream boundary of the domain x = 0 is taken to be a symmetry axis (ice divide), where we impose the horizontal velocity u(x = 0) = 0. The downstream boundary, x = x f corresponds to the calving front. The position of the calving front x f is fixed, and the GL position x g is delimited by 0 ≤ x g ≤ x f . The upper ice surface z = z s (x, t) is in contact with the atmosphere, where pressure 10 is negligible with respect to involved stresses inside the ice body. This is a stress free surface, implying the following condition: where n is the outward pointing unit normal vector. The lower surface z = z b (x, t) is either in contact with the bedrock or with the ocean, 15 and two different boundary conditions will be applied for the Stokes problem on these two different interfaces, defined as: In Eq. (7), the water pressure p w = p w (z, t) is defined as: where ρ w is the water density and w is sea level, assumed constant and equal to zero in what follows. Where the ice is in contact with the ocean (first condition in Eq. 7), the following Neumann boundary condition applies for the Stokes equations: Where the ice is in contact with the bedrock (second condition in Eq. 7), a nonpenetration condition is imposed as well as a friction law, such as where τ b is the basal shear stress, t is the tangent vector to the bedrock, u b is the sliding velocity, C is the friction parameter and m is the friction law exponent (see Table 2 for the adopted values).

Shallow shelf/shelfy stream approximation (SSA)
As mentioned previously, three of the four models use the Shallow Shelf Approximation

15
(SSA) which is a vertically integrated approximation of the Stokes Eqs.
(3) and (4). The horizontal velocity u(x) is obtained by solving the following equations (Morland, 1987;MacAyeal, 1992): where h = h(x) is the ice thickness, τ xx = 2η∂ x u is the longitudinal deviatoric stress and 20 u is the horizontal velocity in the flow direction. The effective viscosity, η, is computed as in Eq.
(2), where D e ≈ ∂ x u. The parameter γ is defined as: According to the SSA approximation, ice deformation is dominated by membrane stresses and vertical shear within the ice is neglected. For the SSA model, the only boundary condition is u(x = 0) = 0 at the ice divide, whereas other boundary conditions are already implicitly included in the set of Eqs. (11). The lower surface z b is determined from the non-penetration condition and the float-5 ing condition: The upper surface z s = z b + h is deduced from the vertically-integrated mass conservation equation giving h as 10

Grounding line treatment
The implementation of GL treatment differs from one model to the other. In this section we define for each model the specificities regarding the treatment of the GL.
The FS-AG model solves the contact problem between the ice and the bedrock. During a time step, the contact condition (Eq. 7) is tested at each node of the mesh 15 and the bottom boundary conditions (Eq. 9) or (Eq. 10) are imposed accordingly. More details about this method and its implementation can be found in Durand et al. (2009a). The consistency of this GL implementation strongly depends on the grid resolution, and a grid size lower than 100 m is needed to obtain reliable results (Durand et al., 2009b). In order to reach this resolution while considering a reasonable number of mesh nodes, 20 an adaptive mesh refinement around the GL is applied: the horizontal distribution of nodes is updated at every time step, such that finer elements are concentrated around the GL. For the SSA-FG model the grid points are kept fixed in time and the last grounded grid point is determined through the flotation criterion, i.e. by solving the following equation: The GL position x g is given with sub-grid precision between the last grounded grid point 5 and the first floating point following the method proposed by Pattyn et al. (2006). The GL position is also determined with sub-grid precision following Pattyn et al. (2006) for the SSA-H-FG, but while SSA-FG uses the flotation criterion as a boundary condition at the GL, the SSA-H-FG model makes use of an additional boundary condition based on the semi-analytical solution of Schoof (2007). The ice flux at the GL q g 10 is calculated as a function of ice thickness at the GL h g : and is used in a heuristic rule to enable GL migration (Pollard and DeConto, 2009). This parameterization allows relatively coarse resolutions to be used (10 km in this study) and gives steady-state results of GL position that are independent of the cho- 15 sen resolution and agree well with the semi-analytical solution given by Schoof (2007) (Docquier et al., 2011). In Eq. (16), the coefficient θ accounts for buttressing and is defined as The numerical approach used by the pseudo-spectral SSA-PSMG model consists in 20 explicitly calculating the rate of GL migration,ẋ g , according to the following explicit  (Hindmarsh and LeMeur, 2001) x where F is given by Eq. (15). At each time step, a new position is computed and the grid moves accordingly, so that the GL coincides exactly with a grid point (Hindmarsh, 1993). Moving grids have the ability to ensure that a grid-point always coincides with 5 the GL, allowing easy representation of gradients at this location but, are not always convenient to implement.

Calving front boundary condition and the specification of buttressing
The experiments we propose are driven by changes in the buttressing force. One approach could have consisted of applying lateral friction on the ice-shelf following the 10 method of Gagliardini et al. (2010), but the total buttressing force would then have been function of the ice-shelf area and ice-shelf velocities, and therefore different for all models. In order to ensure the same buttressing force for all models, we follow the method proposed by Price et al. (2011), in which the inward force at the calving front is modified by a factor, noted C F in our study. 15 For vertically integrated models, the horizontal force acting on the calving front is entirely due to the hydrostatic water pressure and the longitudinal deviatoric stress at the front is given by (MacAyeal et al., 1996): where h f is the ice thickness at the calving front. In the case of the vertically integrated 20 models SSA-FG, SSA-H-FG and SSA-PSMG, a factor C F is then used to modify longitudinal deviatoric stress (Eq. 19), which becomes:

3913
A value of C F = 1 means that the longitudinal deviatoric stress at the calving front implies that ice extension is opposed solely by water pressure, corresponding to no buttressing. Values less than one induce a lower extensional longitudinal deviatoric stress at the front, simulating the effect of buttressing. Note that this procedure implies an additional force applied at the calving front; this results in a varying stress upstream 5 as the ice thickens. Moreover, for SSA-H-FG, the buttressing parameter C F is by construction incorporated in the boundary condition at the GL. This boundary condition relates the ice flux q g to the ice thickness h g at the GL and includes the buttressing factor θ as defined by Eq. (17). From the SSA equations in the ice shelf, we derive (see Appendix A) the 10 relation that links θ and C F through both the ice thickness at the GL h g and the ice thickness at the calving front h f : The other two SSA models solve for the longitudinal variation of τ xx in the shelf to compute the value at the grounding-line. 15 For the FS-AG model, the hydrostatic pressure p w (z) is imposed along the ice column in contact with the sea, so that the longitudinal Cauchy stress is not uniform on this boundary. This non-uniform stress induces a bending of the ice-shelf near the front. To avoid an increase of this bending when adding the buttressing, the stress condition at the front is modified by adding a uniform buttressing stress p b , such that Using Eqs. (22) and (20), and assuming the equality of the mean longitudinal Cauchy stress for both parameterisations, the buttressing stress to be apply at the front of the full-Stokes model is obtained as a function of C F (see Appendix B), such as Note that p b has to be computed at each time step since it depends on the ice thickness at the front, which is not constant.

Experimental setup
We consider an ice sheet resting on a downward sloping bedrock, with the calving front fixed at 1000 km. The GL never advances as far as this in the experiments. The flow pa-5 rameters summarised in Table 2 are used by each model in order to calculate a steady state geometry. The steady state is obtained with a buttressed ice-shelf (C F = 0.4).
Computed steady surfaces are in good agreement between models, exhibiting only a slight difference in GL position of less than 20 km (see Fig. 1). We chose the simpler, stable case of a forward slope for the simple reason that computing comparable initial 10 starting conditions on the unstable reverse slope is a practical impossibility. Groundingline retreat rates are governed by the water depth and the buttressing, and we chose values that were physically acceptable and also produced physically reasonable retreat rates.
Ice-sheet geometry is subsequently perturbed by a release of the initial buttressing 15 force. This process, arising from increased melt of the ice shelf, appears to be responsible for the observed acceleration of Antarctic outlet glacier (Pritchard et al., 2012). Starting from the steady geometries obtained with initial factor C F = 0.4, the buttressing force is decreased at t = 0 ( i.e. C F increases) and kept constant during the simulation.

Transient behaviour of direct observable variables on actual ice sheets
We first evaluate the response of the various models regarding the variables that are currently observed over actual ice sheets, namely GL position (Fig. 2), surface elevation change (Fig. 3) and surface velocities (Fig. 4).

5
As expected, release of buttressing induces a GL retreat, and the greater the release, the higher the amount and rate of retreat (Gagliardini et al., 2010). Retreat can reach up to almost 100 km in 200 yr following a complete loss of buttressing restraint (C F = 1, see Fig. 2 and Table 3). The different models show a similar trend regarding the temporal evolution of GL position (left panels in Fig. 2). However, owing to the var-10 ious initial steady state profiles, the GL position differs between models. For the three perturbations, SSA-H-FG shows the highest GL retreat compared to the initial position, followed by SSA-FG, then SSA-PSMG, and finally FS-AG (Table 3). The evolution of the GL position of SSA-H-FG has a step-like behaviour due to the model grid size (10 km). 15 Rates of GL migration (right panels in Fig. 2) for SSA-PSMG and SSA-FG exhibit a very similar pattern, i.e. a high retreat rate value in the beginning of the perturbation and then a convergence towards a zero-value. Moreover, the greater the perturbation (higher value of C F ), the higher the retreat rates in the beginning of the perturbation. The smooth decrease of the migration rate computed by SSA-PSMG is due to the ex-20 plicit way the GL migration is computed (see model description above). Because the SSA-FG interpolates the GL position between the last grounded point and the first floating point (Pattyn et al., 2006), it also ensures a smooth description of GL migration rate. However, FS-AG and SSA-H-FG show discontinuous GL migration rate induced by numerical artefacts: both models give results that are affected by their grid size. The 25 stepped patterns obtained with FS-AG are due to high frequency oscillation between two successive nodes during GL migration: the GL retreats, then stays at the same position during one time step, then retreats, etc. so that the GL migration rate oscillates 3916 Introduction with an amplitude of 500 m yr −1 (i.e. grid size divided by time step). The numerical noise found in SSA-H-FG is due to a combination of both the grid size effect and single-cell dithering (Pollard and DeConto, 2012). As a general trend, the GL retreats by 10 km steps as a consequence of the model resolution (grid size effect). At some discrete GL positions (every 10 km), the rate of GL migration varies significantly due to the heuris-5 tic rule used in the model (flux imposed either upstream or downstream the GL), so that the GL slightly advances and retreats within the same grid cell (single-cell dithering). In summary, the GL retreats by 10 km (corresponding to the model resolution) and reaches a discrete position where it oscillates within the same grid cell, and then retreats before reaching another discrete position again, etc.

10
Rates of surface elevation change through time and distance from the ice divide are presented in Fig. 3 for the various models and perturbations. The horizontal velocity on the surface velocity is similarly plotted (see Fig. 4). The largest perturbation (C F = 1) exhibits rates of surface elevation change of a few meters per year in the beginning, with horizontal velocities above one kilometer per year. Together with GL migration 15 rates of the order of a kilometer per year (Fig. 2), those are in general agreement with the obervation for currently recessing glaciers of West Antarctica, and Pine Island Glacier in particular (Rignot, 1998;Rignot et al., 2011). That confirms the relevance of the amplitude of the perturbations applied. Rates of surface elevation change are quite similar between the four models (Fig. 3). The highest thinning rates appear in the 20 vicinity of the GL at the beginning of the perturbation. Similarly, the surface velocities steadily decrease during the simulation (Fig. 4). High frequency and small amplitude numerical noise in FS-AG appear not to significantly affect the surface response. However, with SSA-H-FG the high frequency and amplitude variabilities drastically affect the surface thinning rate and velocities over short time scales (i.e. about a decade). 25 We deliberately chose a low spatial resolution (uniform 10 km along the flowline) for the SSA-H-FG model compared with other models. As shown by Docquier et al. (2011), such models perform well at low spatial resolution, which is the main motivation for applying such parameterizations in large-scale ice sheet models. Increasing the TCD Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | resolution (down to 250 m) allows removal of high frequency numerical artefacts and a better match to the three other models (data not shown). However, this significantly increases its numerical cost, so that the major advantage of this model is lost, as well as its applicability to large-scale ice-sheet models. 5 Despite the numerical noise exhibited by SSA-H-FG and FS-AG models, the evolution of the geometry during the simulations appears very similar for all four models. However, the boundary layer theory implemented in the SSA-H-FG model hypothesizes near-steady conditions and its ability to represent transients requires evaluation. In Fig. 5, the flux at the GL is plotted as a function of the instantaneous ice thickness at the GL for all models and simulations. By construction, SSA-H-FG essentially follows the boundary layer prescription. This can most clearly be seen for the case C F = 1 (see the bottom of Fig. 5) where the close correspondence of the curves of Schoof (2007) and SSA-H-FG is evident. This correspondence is not evident for the other perturbations, since the SSA-H-FG boundary condition for the flux now relies on a parameterization of θ, which in turn depends on the quantity h f /h g (see Eq. 21). Since this ratio varies in time, the steady-state condition of the Schoof condition is not fulfilled.

Divergence from the boundary-layer solution
Interestingly, and despite their very different physical and numerical approaches, all the other models show very similar behaviour, with the boundary layer theory result 20 attained after some time. This is most obvious for the largest perturbation (C F = 1) but also clearly visible for the weaker perturbations (C F = 0.8 and 0.5). However, during the highly transient phase, for a given ice thickness at the GL, the ice flux is substantially overestimated by the boundary layer theory, consequently overestimating the outflow during the whole period of 200 yr. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Changes in Volume Above Flotation (∆VAF)
From the perspective of projecting the future contribution of Antarctica to sea level rise, the change in Volume Above Flotation (∆VAF = VAF(t) − VAF(t = 0)) is certainly a pertinent variable to investigate. Indeed, ∆VAF has the advantage of integrating through time both the contribution coming from outflow at the GL and the consequence of 5 grounding-line retreat in terms of ice release. In our case, this further allows investigating the spread in the transient behaviour of the various models in response to similar perturbations. We arbitrary choose to express the evolution of ∆VAF for each model relative to ∆VAF computed by FS-AG. Corresponding evolutions are presented in Fig. 6. It is first of all striking that response in terms of relative ∆VAF from one model 10 to the others are extremely similar, irrespective of the amplitude of the perturbation. As anticipated, SSA-H-FG shows the greatest change in VAF compared with other models. Relative to FS-AG, SSA-H-FG overestimates the contribution to SLR by more than 100 % during the first 50 yr of the simulation, which decreases to a 40 % overestimation after 200 yr. SSA-FG shows a similar pattern with a smaller overestimation 15 (about 15 % after 200 yr). On the other hand, SSA-PSMG briefly underestimates the change in VAF relative to FS-AG at the beginning of the perturbation, but after 20 yr the contribution of the models to SLR is remarkably similar to the one computed by FS-AG, with relative difference below 5 %. While this intercomparison has not given definitive values for the GL migration rates, it strongly suggests that models prescrib-20 ing flux at the GL according to the boundary layer theory most probably overestimate ice discharge. It also clearly shows that the rate of contribution to SLR significantly differs from one model to the other, even for a relatively simple and constrained experiment. When extrapolated to the current imbalance of the Antarctic ice sheet, this would have important consequences. According to Rignot et al. (2011), the Antarctic Assuming that ice-sheet models were 3919 capable of describing exactly the ice dynamical conditions in 2000, and also assuming the parameters forcing enhanced ice discharge to be properly known, the above uncertainties in modelled ice discharge would lead to a erroeneous contribution between 3.1 and 18 mm in 2010, due to the over-estimation of SSA-H-FG of more than 300 % in the first 10 yr after a given perturbation. Furthermore, as ice sheets are still in a transient 5 phase (i.e., perturbations are sustained through time) the discrepancy of the models would eventually increase with time integration. Of course, these assertions have to be moderated by the fact that the complexity of actual 3-D geometries could mitigate the discrepancy between model results, which is the focus of future research. 10 We have computed the transient response of four flow-line ice-sheet models to a reduction in the buttressing force exerted by an ice shelf onto the upstream grounded ice sheet. The intensity of buttressing perturbations was chosen in order to reproduce changes in geometry that are comparable to those observed on current ice sheets. Compared with MISMIP, we investigated the transient response in more detail and ap- 15 plied a perturbation that reflects direct mechanical forcing. The physics are implemented in a different way in the different models (from SSA to the solution of the full-Stokes equations), while the models differ in their numerical treatment as well (finite difference and finite element). One of the models includes the heuristic rule of Pollard and DeConto (2009), i.e. the flux-thickness relation proposed 20 by Schoof (2007) is imposed at the GL. All models have successfully participated in the MISMIP benchmark (Pattyn et al., 2012), exhibiting unique stable positions on downward sloping beds, unstable GL positions on retrograde slopes and related hysteresis behaviour over an undulated bedrock.

Conclusions
Surprisingly, and despite the different physics and numerics implemented, all models 25 give consistent results in terms of change in surface geometry and migration of the GL. A major divergence is found with the SSA-H-FG model which directly implements the 3920 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | boundary layer theory proposed by Schoof (2007). Here, the prescription of flux at the GL introduces high frequency, large amplitude numerical noise deteriorating the surface change signal over decadal time scales. Moreover, it seems that at least in these experiments that the boundary layer theory overestimates the discharge during the transient evolution. As a consequence, models that prescribe the flux at the GL 5 should be used with particular caution when dealing with small spatial and temporal scales. Estimation of the contribution to SLR through numerical modelling still exhibits large uncertainties, with results from different models showing > 100 % spread on a decadal time-scale and still around 40 % two hundred years after the initial change in buttress- 10 ing. There may be a large uncertainty in models that are seeking to establish reliable projection of coming contribution of Antarctic ice sheet to SLR. Further model intercomparison must be pursued to better constrain the rate of discharge, and intercomparisons on specific Antarctic outlet glaciers should be encouraged in the near future.

15
In this Appendix the relation between the buttressing factors θ in Eq. (17) and C F in Eq. (20) is derived. The ice-shelf equation is where h is the ice thickness along the ice-shelf. The longitudinal deviatoric stress within the ice shelf is then obtained as

Appendix B
In this appendix, we demonstrate how is obtained the buttressing pressure p b (t) in Eq. (22) giving the front-stress for the FS-AG model. We need to find p b (t) such that the 3922 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | mean longitudinal Cauchy stress be the same for all models. This equality is expressed as follows: whereσ SSA xx andσ F S xx are the longitudinal Cauchy stress of SSA models and FS-AG model,= respectively.

5
The mean longitudinal Cauchy stress for SSA models reads: whereσ zz = − ρ i gh f 2 andτ xx is given by Eq. (20). The longitudinal Cauchy stress for FS-AG model, given by Eq. (22), and once integrated over the ice column gives: Using Eq. (B2) for SSA models and Eq. (B3) for FS-AG, Eq. (B1) leads to Using the flotation condition ρ i h f = ρ w z b , and after simplifications, p b can be isolated and deduced as