Comparison of ice dynamics using full-Stokes and Blatter-Pattyn approximation: application to the Northeast Greenland Ice Stream

. Full-Stokes (FS) ice sheet models provide the most sophisticated formulation of ice sheet ﬂow. However, their ap-plicability is often limited due to the high computational demand and numerical challenges. To balance computational demand and accuracy, the so-called Blatter-Pattyn (BP) stress regime is frequently used. Here, we explore the dynamic consequences :: of :::: using ::::::::: simpliﬁed ::::::::: approaches : by solving FS and the BP stress regime applied to the Northeast Greenland Ice Stream. To ensure a consistent comparison, we use one single ice sheet model to run the simulations under identical numerical conditions. A 5 sensitivity study to the horizontal grid resolution (from 12.8 down to :: to : a :::::::: resolution :: of : 0.1 km) reveals that velocity differences between the FS and BP solution emerge below ⇠ 1 km horizontal resolution and continuously increase with resolution. Over the majority of the modelling domain both models reveal similar surface velocity patterns. At the grounding line of 79  North Glacier the simulations show considerable differences whereby BP ::: the ::: BP ::::: model : overestimates ice discharge of up to 50% compared to FS. A sensitivity study to the friction type reveals that differences are stronger for a power-law friction than a 10 linear friction law. Model differences are attributed to topographic variability and the basal dragsince , :::::: where neglected stress terms in BP become important. compared to a linear friction law (up to 17% compared to 8%, respectively). Furthermore, the weak implementation produces larger differences compared to the strong scheme. The P2P1 discretization reveals differences that are in between P1P1GLS-strong and P1P1GLS-weak results.

Many thanks for providing this new version of your manuscript. I went through the manuscript and have provided final feedback that I would like you to consider when re-submitting your files. The list may seem somewhat long at first, but most comments should be very easy to incorporate. Once these final suggestions are addressed, we should normally be able to proceed to the acceptance of the manuscript. l.3: we explore the dynamics consequences -> we explore the dynamic consequences of using simplified approaches by : i.e. suggest being a bit more specific here. Done l.9: whereby the BP model overestimates Done l.11: and the basal drag, where neglected stress terms in BP become important Done l.16: computationally intensive -> computationally expensive Done l.17: because FS effects (vs. sim models) only occur at higher resolutions Done. We deleted compared to simpler models . l.18: in the last decades allowed for ice flow models to directly rely on FS equations or approaches with only few simplifications Done l.20: require an accurate representation of stresses, e.g . Done l.21-22: necessary for long time (i.e. remove e.g. here) Done l.25: neglects several components of the full system : would be good to mention which ones here Done. We added (e.g., bridging stresses) . l.32: IPCC, 2013: true, but bit outdated. I guess a similar statement was made in SROCC or in AR6? Refer to one of these reports instead? Well, the second reference Oppenheimer et al. 2019 is the reference to SROCC (chapter 4). However, we added also the ref to SROCC chap 3 (Meredith et al. 2019). l.34: was explored: e.g. in ISMIP-HOM (and remove second bracket on l.37) Done l.37: focused on Done l.38: show a smaller spread (i.e. remove, much ) Done l.38: and with available analytical solutions (remove one ) Done l.39: was not possible Done l.43: cost when running FS models Done l. 45-46: simpler models was tested on Done l.59: on FS and BP Done l.59: more frequently used : compared to what? FS/SIA? Would be good to be more specific here to avoid possible confusion Done. We have rewritten the sentence to: Our comparison focuses on FS and BP, as the BP is used as a compromise for FS to balance computational costs and accuracy. l.62: FS and BP intercomparison, stepping away from synthetic scenarios, and performing a highresolution Done l.66: Our study does not treat the Done l.68: 12.8 km to a resolution of (i.e. remove down ) Done. We also replaced down in the abstract (line 6). l.71: given the importance of the full-Stokes (FS) equations in your storyline, would it make sense to also explain shortly how this relates to navier-Stokes, which may be more familiar to those outside the glaciological community? Done. We added: In the case of ice sheets, the Reynolds number is very low (e.g., Fowler and Larson, 1978) and therefore the inertial term in the Navier Stokes equations can be neglected. This approximation is often referred to as Stokes flow. l.104: not sure if can have colours in equations in final document. Needs to be seen with copyediting. Suggest adding a note to flag this. Thanks, we will take care of this. If colours are possible, we will think about an alternative way to highlight the differences (e.g., underline, subscript etc.). l. 117-118: van der Veen and Whillans Done l.122: reduces the computational demand : could you give an indication by how much this is the case? A value was already given in line 43. l.126: alternative way to align with the BP stress regime, thereby directly allowing to compare FS and HO model simulations : the second part of the sentence that I suggest to add is not a must, but think it would be nice to explicitly mention this here to conclude this section. Done l.130: remove bracket after 2 Done l.132: in order to align with the BP stress regime Done. l.133: That means that -> Through this approach, Done l.130-140: nicely formulated. Very clear, also for non-expert! Thanks! l.159: increases -> increase Done l.163: which requires more computation time and working (or less computation, in case I did not understand this correctly. There s some room for misinterpretation in actual formulation, therefore good to change also) Done. The sentence is changed to which requires less computation time and working memory compared to MUMPS l.183: Exp. C is a parallel-sided slab Done l.187: In the parallel direction, 15 layers : in general, many instances throughout the text where a , could/should be added. Could check for this, although will normally be done through copy-editing phase. : reveals distinct differences between FS and lower order approximations : not entirely clear here: is this FS vs. HO? Or vs. SIA? Would be good to be more specific here. Done. Sentence is rewritten to: reveals distinct differences between FS and all other approximations on short length scales. l.195: On those : shorter time scales? Would be good to specify what those refers to here Done. We deleted this sentence as it is basically a repletion of the previous sentence. l.197: we found a good agreement with the original ISMIP-HOM contributions and within our own model versions : not sure if with the latter within the sentence expresses what you wanted to say. But thought this d be clearer than original formulation. If disagree, feel free to leave as is of course Done. l.203: focuses on Done l.204: The second region also includes the Done l.207-208: bed and surface topography : what is source for surface topography in BedMachine? Done. Sentence is rewritten to: Bed and surface topography are taken from BM and the GIMP data set (Howat et al., 2014), respectively. caption figure 1: overlay on each other -> overlap ? Done l.214-15: when referring to the shock : suggest to also directly refer to the ISMIP initialization experiments (Goelzer et al., 2018, TC) here, as they clearly make this point Done l.218: allow for a better comparison Done l.220: is assumed to be: Done l.229: density of the ocean water. Possibly add a reference for taking this value? We think it is not necessary to add a reference. It is just a model choice and will most likely not influence the results (i.e., the differences between BP-like and FS). l.236: slip-free condition ? No. The BC is known as free-slip condition. However, we added a hyphen. l.238: For those outflow boundaries ? Not sure, but would be good to specify what the those refers to here Done l.242: ZI; where Done l.253: from a few seconds to several hours : also the upper limit does not seem to be a lot.
How does this align with the statement made in the previous line which refers to this being expensive? The longest simulation in Tab. 1 is about 11 hours on 1440 computational cores and huge RAM requirement. We consider this as an expensive model run. Note, that this just P1P1-GLS. l.260: is computed by Cuffey and Paterson -> is computed following Cuffey and Paterson + also further in sentence when referring to Lliboutry and Duval (1985). Otherwise sounds like Cuffey and Paterson calculated this themselves... Done For the reference to the value by Cuffey and Paterson, suggest to also include the Table where this comes from in the book (if I m not mistaken there s several values suggested there, no?) Done. We added Tab. 3.4 therein . Table 1: o MPI)) are kept constant between o hybrid mode, meaning the process o Each production node has 362 GB RAM and Done l.268: ice-stream region for l=6400 m and 100 m is provided in Fig. 3 Done l.269-270: (Fig. 3b,e) at higher resolution reveals a rapidly varying bed (remove in the finer resolution at end of the century) Done l.278: could certainly be Done l.285: The largest impact of using a BP-like solution results from the surface flow velocities, which are up to 43 ma-1 faster ( Done l.290: spread in differences is more pronounced for Done. for E = 0.1, : remove or Done l.311-312: split in two sentences: (Fig. 6c,f). However, Done l.316: difference is or differences are Done Fig. 6 caption: not easy to follow. Suggest putting the references to the panels before the explanation: i.e. (a) slope of bed topography, (b) basal friction coefficient for E=1 and (c) simulated surface Done Fig. 8  we strive to give a flavour of the origins of the model differences by . Figure 10: Lower shows -> Panels (c) and (d) show the Done l.365: Figure 11 and 12 shows -> Figure 11 and 12 show Done l.370: lower in some larger areas Done l.374: in the majority larger -> is generally larger for Done l.379: weak and strong setup: is not entirely clear what this is here. The weak and strong setups are introduced in Line 152. However, we have slightly rewritten the sentences here to clarify the weak and strong setup: In the outlet domain, the basal drag difference is more pronounced for the weak implementation of the basal boundary condition with m=3 compared to the strong implementation with m=1. At the 79NG and ZI area, the maximum difference in basal drag is 316 kPa (99%) and 121kPa (73%) for the weak implementation with m=3; for the strong setup with m=1 it is 227kPa (72%) and 96kPa (55%), respectively. l.386: , in the BP-like -> , the BP-like Done l.387-388: Similarly, the difference between FS and BP is most pronounced in region with a high slip ratio , while with a higher resolution : not sure about the second part ( while with ), as the sentence was not entirely clear to me. If not correctly reformulated, try to opt for slightly different formulation to better reflect what you mean here. Done. Sentence is rewritten as suggested. l.391: 0 to 0.05, the higher bed undulations lead to larger errors at higher resolutions; however it : not entirely clear to what the last it refers here: (the neglection of) the bed slope? Done. Sentence is rewritten to: By neglecting the bed slope bin from 0 to 0.05, the higher bed undulations lead to larger errors at higher resolutions; however, the bed slope seems to play a secondary role, compared to e.g. the slip ratio, as fewer grid nodes are involved. l.393: would be good to also briefly refer to work on (mountain) glaciers that has shown this also. Particularly thinking of Le Meur et al. (2004Meur et al. ( , https://doi.org/10.1016Meur et al. ( /j.crhy.2004.001) for this. Done. We added a sentence: A similar conclusion for mountain glaciers was drawn in Le Meur et al. (2004) by comparing FS and SIA. l.397: are to date rarely used : would be good to specify this, as will change in the (near) future Done l.398: we found distinct differences + l.399: these differences : not entirely clear, would be good if you could specify differences in what? Done. In both cases, we added velocity . l.400: that at a higher resolution the disagreements would be more substantial, but Done. l.400: bedrock topography database: how is this linked to the disagreements you mention before? Not entirely clear. Suggest reformulating. Done. We deleted the but a bed topography database to investigate this is not available to date. l.401: computationally expensive, it Done l.402: that are currently being run within the ice sheet modelling community : see comment above also Done Figure 11 caption: o Pas -> Pa s ? Well, in the tex file a half space is already typed (i.e. Pa\,s). We will check in the proofs. o that range from E=1 to E=6 Done. Figure 12 caption: equivalent fields but for m=3 , correct? Yes, thanks. l.409-410: Figure 15 (panels a and c) shows the Done l.414: that is for most areas higher Done l.416: made to BP : not entirely clear here. Maybe simply omit? Done l.432: using the geometry from Done l.434: respond somewhat similarly (see Done l.435: coverage or the fact that ZI has Done l .438: and analysing how the Done l.467: must be investigated in future studies Done l .468: within a thermo-mechanical model? Done l.469: ice flow. Additionally, there : as is not really in contrast with what you have said before it seems. Seems more like an extra suggestion. Done l.480-481: suggest omitting An intriguing effect was identified at start of the sentence and also reformulate in general: In the ice stream region, when considering a rheology with much softer ice (E=6) than typically considered, and enhanced spread and BP-like was modelled : than typically considered -> why is this? Done. However, we deleted typically considered . When looking into the literature, there seems a large range of values for E considered. The choice is dependent on whether E is used as tunning parameter to account for missing physics or the choice of E is supported by e.g. anisotropic ice flow modelling. However, we think typically is subjective and therefore we deleted it. l.481: In contrast, stiffer ice (E=0.1) leads to : was not clear what the reference is for 10 time stiffer. Therefore suggested removing this. Done. Ten times stiffer was meant compared to the classical case (E=1). l.484: as important surface velocity differences occur (up to of 79NG). Numerical Done l. 485-486: to the choice of friction type: model differences are particularly important when a power-law friction is applied (as opposed to a linear approach) Done l. 491: results from non-FS models : otherwise a huge range of models is possible (e.g. SIA, SIA-SSA,..etc) Done l.492: software that is not freely Done l.494: remove ( before Larour Done Acknowledgments: please acknowledge the input from the reviewers Done Best regards, Harry

Introduction
The most comprehensive description of ice sheet flow is given by the full-Stokes (FS) equations. Such a model is considered 15 the most accurate available, capable of describing highly dynamic ice sheets, including ice streams, ice shelves, and grounding line migration. But this formulation is also the most computationally intensive :::::::: expensive, both from the numerical perspective and because FS effects can only occur at higher resolutionscompared to simpler models.

35
In earlier ISMIPs the performance in simulating idealized ice stream and marine ice sheet dynamics by ISMs was explored( : : was : not possible due to the limited number of FS models 1 . The recent MISMIP+ effort reflects that the choice of the friction law plays a larger role compared to, e.g. ice flow approximations (two FS models contribute to this exercise). From these synthetic scenarios, it yet needs to be assessed if discrepancies between FS and simpler models are minor compared to other uncertainties in ice sheet modelling given the huge additional computational cost by :::: when ::::::: running FS models (about 10 times 45 slower than BP (Larour et al., 2012, Table 3 therein)).
In contrast to the idealized ISMIP experiments, the influence of FS on ice dynamical behaviour compared to simpler models were :::: was tested on realistic problems (e.g. Leysinger Vieli and Gudmundsson, 2004;Le Meur et al., 2004;Morlighem et al., 2010;Seddik et al., 2012;Favier et al., 2014;Seddik et al., 2017). These exercises provide important insights into the transient response and indicate that FS produces different results compared to simpler models. For example, in most of the conducted 50 projection experiments, FS tends to contribute less to sea-level rise than simpler models under identical projection scenarios (Seddik et al., 2012;Favier et al., 2014;Seddik et al., 2017). Morlighem et al. (2010) concluded that treating ice flow of the fast-flowing Pine Island Glacier (PIG, Antarctica) with a steeply rising bed near the grounding line with FS is essential.
However, assessing whether FS is urgently needed compared to simpler models is complicated because of many interacting processes: (i) comparability is often limited by the use of different ISMs as it is not entirely clear how much of the differences 55 is due to numerical treatments as, e.g. different horizontal grids resolve the bed differently, or discretization schemes. (ii) 1 Compared to the ISMIP-HOM experiments, the MISMIP experiments require a long transient integration, which is costly to perform for an FS model.
Differences are also subject to e.g. assumptions in basal flow conditions or the initialization approach. (iii) Moreover, when transient simulations are forced with climate scenarios, the response involves numerous processes and interactions that make it difficult to separate causes.

Boundary conditions
The upper surface is assumed to be traction free. The ice sheet base is subject to basal sliding according to a friction law for the basal shear stress t D b in the tangential plane with the unit normal vector n pointing out of the ice, and v b is the velocity in the tangential plane at the base. Basal refreezing or melting is neglected. The basal drag parameter 2 often includes the effective pressure to account for the presence of water lubricating the ice-bed interface, as well as bed properties such as roughness. For floating ice we use the same boundary condition but prescribe free slip ( 2 = 0). The parameterization for 2 and the remaining lateral boundary condition are given 100 below for each experiment setup.

Model classification
The mass and momentum balance (Eq. 1 and 2) can be written in terms of velocity components and pressure p @v x @x + @v y @y 105 as well as the effective strain rate occurring in Eq. 5 The red terms are omitted in the BP approximation as explained below. The four equations (Eqs. 1 and 8) with the four un- knowns v x , v y , v z and p form the FS equation system. Various approximations to this set of equation exist ranging from shallow ice approximation (SIA) to different types of one layer and multilayer approximations for non-horizontal shear components of the stress tensor (Hutter, 1983;Hindmarsh, 2004;Bueler and Brown, 2009;Ahlkrona et al., 2013). The FS problem has saddle point character which is difficult to solve (Benzi et al., 2005). In order to reduce the computational demand and derive a better posed problem, Blatter (1995) and Pattyn (2003) developed the so-called BP approximation scheme which is valid for most of the ice sheet domain and widely used in ice flow models. The main assumption of their model is that the vertical component of the momentum balance is approximated as hydrostatic (i.e. @t zz /@z = % i g). This assumption is used to eliminate 115 the pressure variable p from the FS equation system and reduces the FS problem with the two horizontal velocity components as unknown field variables. The further assumption that horizontal gradients of the vertical velocity are small compared to the vertical gradient of the horizontal velocity (i.e. @v z /@x ⌧ @v x /@z and @v z /@y ⌧ @v y /@z) leads to a closed problem for the two horizontal velocity components. This approximation is equivalent to dropping the red coloured terms in Eq. 8 and 9; it is termed LMLa in Hindmarsh (2004).
The simplifications made by BP on the FS system reduce :::::: reduces : the computational demand by deriving a better posed problem. By re-arranging the BP equations this system leads to a closed problem for the two horizontal velocity components v x and v y . In doing so, the BP approximations decouple the vertical velocity v z and the pressure p from the full system. Both unknowns are recovered diagnostically by integrating the incompressibility equation (Eq. 1) and the vertical component of the 130 momentum balance equation (Eq. 8), respectively. However, in our approach we follow an alternative way to realize :::: align :::: with the BP stress regime : , :::::: thereby ::::::: directly ::::::: allowing :: to :::::::: compare ::: FS ::: and ::: BP ::::: model :::::::::: simulations : (see Sect. 3.1 and 3.2).

Ice flow model
To solve the set of equations presented above (Eq. 1 and 2 ) and boundary conditions), we use the COMmercial finite element 135 SOLver (COMSOL) Multiphysics © (version 5.5, www.comsol.com, last access May 2021). COMSOL provides a FS interface which is additionally manipulated in order to realise :::: align :::: with : the BP stress regime. This is accomplished by discarding the corresponding terms (red terms in Eqs. 8 and 9) in the weak formulation but maintaining identical FS numerical details.
That means that ::::::: Through :::: this :::::::: approach : the BP stress regime does not further modify the default FS module, i.e the BP version implemented in this study is not forming a multilayer model with the two horizontal velocity components as field 140 variables. Therefore, our implemented BP approximation is hereafter termed BP-like. In doing so, both flow modes (FS and BP-like) are running through the same (FS) solving algorithm. Hereby we ensure that differences in the computed results due to numerical issues are largely eliminated. However, the disadvantage is that we do not benefit by the reduced complexity of the BP approximation. As a consequence of the avoided decoupling @v z /@z in Eq. 9 is not replaced by (@v x /@x + @v y /@y) as commonly done in HO approximations (e.g., Pattyn, 2003).

Numerics
The mass balance and momentum equations constitute a nonlinear equation system. Independent of the employed flow regime, a damped Newton method (which is COMSOL's default solver) is applied in a fully coupled manner to solve for the velocity vector v and the pressure p. The unstabilized Stokes equation (Eq. 8) is subject to the Babuska-Brezzi condition, which states that the basis functions for p have to be of lower order than for v. This issue could be circumvented with suitable stabilization 150 techniques. Here, we test the numerical robustness of our results by employing two different discretization schemes: (1) We use P2+P1 Lagrange elements which consists of quadratic basis functions for the velocity and linear for the pressure.
(2) The second scheme uses linear elements for both the velocity components and the pressure field (P1+P1 Lagrange elements) and numerical stabilization is achieved with streamline diffusion (Galerkin least-squares (GLS), Hauke and Hughes, 1994). Due to the equal-order interpolation, the latter is computationally more efficient but less accurate and sensitive to the GLS stabilization 155 parameters (Helanow and Ahlkrona, 2018). In the following, the schemes will be referred as P2P1 and P1P1GLS, respectively.
To further test the numerical robustness we compare two different implementations (called 'strong' and 'weak') for the Dirichlet condition (Eq. 7) at the ice base. The strong imposition applies the constraints pointwise at each node and, therefore, modifies the underlying variational formulation. The weak method is rather a component of the variational formulation by adding exterior facet integrals. For the latter, we use the so-called Lagrange multiplier method (Babuska, 1973;Verfürth, 1986; increase rapidly with increasing spatial resolution, and a direct solver would be uncommon for such large numbers of DOFs. For 165 the larger problems we rely on the iterative GMRES solver (Saad, 2003) which is accelerated with appropriate preconditioners.
Simulations that make use of the strong imposition rely on a Domain Decomposition solver with an overlapping additive Schwarz method (ASM, Widlund and Toselli, 2005) which is much superior in terms of the required ::::::: requires ::: less : computation time and working memory compared to MUMPS. Unfortunately, the ASM preconditioner shows a very high computational demand for simulations that employ the weak imposition. Since the involved Lagrange multiplier induces a zero on the diagonal 170 of the system matrix, several common preconditioners are not applicable. Therefore, we employ the Vanka algorithm (Vanka, 1986;John and Matthies, 2001) which is specifically designed for large indefinite problems with saddle point character. Based on our simulations, we found Vanka to be very memory efficient and computationally fast for large problems although it requires more Newton and GMRES iterations compared to ASM (see Sect. 4 and Tab. 1 and 2).
Regarding the desired consistency between the simulation results, we must clarify the different usage of the linear solvers: 175 the coarser models show a poor performance with the ASM and Vanka preconditioner. At the same time, the higher resolutions run out of memory with the MUMPS solver. Depending on modelling domain and employed discretization, MUMPS runs out of memory around 5 million DOF's. For the largest DOF model that could be solved with the direct solver MUMPS on our cluster system, we also performed a simulation with the iterative linear solvers ASM and Vanka. A comparison of both simulations reveals differences well below 10 7 m a 1 in the surface velocity, v s = |v| s ; differences between MUMPS and 180 Vanka are below 10 3 m a 1 . Consequently, we assume that all applied linear solvers provide comparable results.

Verification by benchmarks
The ISMIP-HOM benchmark targets comparison for three dimensional ice flow models arising from basal topography undulations, from basal sliding and by changing the length scale. From the ISMIP-HOM setup (Pattyn et al., 2008) (Pattyn et al., 2008). Please note that the BP-like P1P1GLS and 'rhi2' BP solutions in Exp. A overlay ::::: overlap : on each other; in Exp. C the 'rhi2' BP solution is overlaid by BP-like P2P1 strong and and BP-like P2P1 weak.
(it is termed 'ice-stream', see red line in Fig. 2). The second region comprises also ::: also :::::::: includes the western branch (WB) of NEGIS and the outlet regions of 79NG and ZI (it is termed 'outlet', see cyan line in Fig. 2). For the ice-stream experiment a 210 region upstream of the grounding line and downstream of the ice divide is selected to avoid grounding line treatment and to focus on an ice-stream like flow regime. The outlet region captures the hinge-zone :::::::: grounding :::: zone : of 79NG and the marine terminus of ZI. As we rely on the BedMachine v3 data-set (BM, Morlighem et al., 2017), a small area at 79NG becomes afloat.
In both setups, the geometry of the NEGIS regions is laterally defined along flow lines that were deduced from the MEaSUREs surface velocity data-set (Joughin et al., 2016(Joughin et al., , 2018.

8
we do not perform an initial relaxation of the free surfaces. A relaxation run based on e.g., the BP-like scheme and subsequent simulation runs with relaxed geometry with BP-like and FS would likely experience a 'shock' ::::::::::::::::: (Goelzer et al., 2018) for the FS model because the geometry is consistent to the BP-like physics (and in the same way also for a FS relaxation run).

220
Performing relaxation simulations for BP-like and FS individually would result in different initial geometries and also different initial conditions (e.g., velocity). This would be an undesired effect for our consistency approach which would make a clear comparison difficult. Here, we followed the strategy to have physical parameters and the geometry input equal in the FS and BP-like to allow :: for a better comparison. White lines in (a) indicate flux gate locations at 79NG and ZI (see Fig. 9 and related text below). WB and SB indicate the western and southern branch of the NEGIS. Symmetry, free-slip, and stress in (b) indicate the boundary condition types used a lateral margins. Background image is a RADARSAT Mosaic (Joughin, 2015;Joughin et al., 2016) In the NEGIS experiments we consider a Budd-like friction law (Budd et al., 1979) with a linear (m = 1) and non-linear 225 relationship (m = 3). The basal drag parameter 2 in the basal boundary condition (Eq. 6) is assumed as :: to ::: be: where k is the basal friction coefficient found by an inversion technique (see below and Appendix A). The effective normal pressure is the difference between the ice pressure P i and the basal water pressure P w N = P i P w .
The basal water pressure is computed in marine parts, i.e. where the ice base is below the sea-level (h b < z sl ), according to Huybrechts (1992) where % w = 1023 kg m 3 is the density of the ocean water and h b the bed elevation. Here, P i is the normal stress instead of 235 the approximated ice overburden pressure. The assumptions made about the water pressure imply that the base is perfectly connected to the ocean at any location in the domain that is below sea-level. Overall, this might be incorrect but it forms an appropriate assumption in the absence of an additional hydrological model.
Solving a subset of an ice stream poses unknown boundary conditions in the interior of the ice sheet. A simple approach would be prescribing the measured surface velocities as a depth-averaged velocity profile. However, we choose boundary 240 conditions that are free to adjust during the solution process (Fig. 2b). Inflow boundaries are chosen to have a symmetry boundary condition, which implies v · n = 0 and vanishing shear stresses. A free slip ::::::: free-slip condition is chosen at lateral along-flow boundaries and a normal stress condition at the outflow boundaries. Dependent on the setup, the outflow boundaries are land-or marine terminating fronts of the glacier or located within the ice. For those :::::: outflow :::::::::: boundaries we chose stress boundary condition t·n = F 0 n. For outflow boundaries located within the ice (e.g., at 79NG) we choose F 0 = % i g(z s z)+ ;

245
for land-or marine terminating fronts (e.g., ZI) we set F 0 = min(% i g(z s z), 0) + , where z s is the surface elevation and z the vertical coordinate. Here, is an additional stress that is tuned to match the observed velocities at those boundaries. The additional stress is set to zero except at 79NG and ZI; ::::: where : is 0.5 MPa and 0.3 MPa at 79NG and ZI, respectively. Overall, computational times varied between a few seconds and several hours. A broad overview of computational costs is given in Tab. 1 and 2. If not mentioned otherwise all setups employ 11 vertical layers (refined to the base) which is sufficient 260 compared to higher resolutions (Fig. S1).
For the parameters A(T, W ) (Eq. 5) and k 2 (Eq. 10) we utilize external products and capacities: 1. To rely on a realistic thermal state but avoiding an intensive thermal spin-up, we make use of an enthalpy field from a paleo-climatic spin-up simulated by the ice sheet model SICOPOLIS with SIA (Rückamp et al., 2019). The field is do not vary much among other E and stress regimes settings. Solver and cluster settings (OpenMP (OMP) shared-memory threads and distributed-memory tasks (MPI)) are maintained identical ::: kept ::::::: constant : between each resolution and stress regime. Computational time includes the solver procedure (e.g. assembling the stiffness matrix, load vector) and outside routines (e.g. construction of mesh, saving results). Generally, COMSOL is run in hybrid mode, that means ::::::: meaning the process is distributed to a user-defined number of OMP threads and MPI tasks. Without testing the performance thoroughly we achieve a good performance by OMP parallelism per physical node and MPI communication between physical nodes (the code is simply run in OMP mode if MPI is equal to one). The listed RAM usage refers to one MPI task. Please note, that the cluster settings (MPI and OMP) have been chosen in a way that the job returns a solution within the cluster 'wall-time' but do not indicate any performance scaling. bilinearly interpolated to the COMSOL meshes. The temperature-dependent part of the rate factor in Eq. 5 for cold ice is 265 computed by Cuffey and Paterson (2010) :::::::: following :::::::::::::::::::::::::::::::::::: Cuffey and Paterson (2010, Tab. 3.4 therein), and the water-contentdependent part for temperate ice is by :::::::: computed :::::::: following Lliboutry and Duval (1985).
2. We make use of the inversion capability in the Ice-sheet and Sea-level System Model (ISSM, Larour et al., 2012) to obtain a spatially varying basal friction coefficient k 2 ( Fig. A1 and A2, see Appendix A). The ISSM inversion is performed with l = 0.8 km, with BP and for the individual E values. Similarly, the SICOPOLIS temperature field is interpolated to 270 the 3D ISSM mesh. Subsequently, the inferred fields for k 2 are bi-linear interpolated to the COMSOL meshes.

Ice-stream region
An overview of input parameters and FS simulation results of the ice-stream region are compared for l = 6400 m and 100 m :: is ::::::: provided : in Fig. 3. The slope of the bed topography (Fig. 3a, d) and the friction coefficient (Fig. 3b, e) of the finer resolution 275 reveals more small scale features and stronger amplitudes; particularly a rapidly varying bedin the finer resolution. Despite both compared resolutions relying on the same inferred basal friction coefficient and input geometry, the interpolation to the computational grids produces striking differences. As a consequence the simulated surface velocity fields differ quantitatively  ice (E = 0.1) the BP-like solutions exhibit higher flow speeds which is almost confined to the fast flow region (Fig. 4a); this 290 area partly coincides with a smooth bed and surface topography (see Fig. 3). Overall, discrepancies between FS and BP-like increase with higher velocities. The largest impact of :::: using :: a BP-like on the solution is by enhancing ::::::: solution ::::: results ::::: from the surface flow velocities : , ::::: which ::: are : up to 43 m a 1 ::::: faster (relative difference is about 19%) compared to FS. For the 'classical case' (E = 1) the pattern is similar but with reduced magnitudes. In this setup, the BP-like velocities are up to 19 m a 1 higher (relative difference is about 6%). With softer ice (E > 1), the pattern in surface velocity differences changes and, eventually, 295 for very soft ice (E = 6) differences are more pronounced in regions with higher topographic variability (compare Fig. 3) and exhibit a rippled pattern. Velocity differences are scattered around the FS solution by ±80 m a 1 ; the spread in differences is much stronger :::: more :::::::::: pronounced : for softer ice (see scatter plots in Fig. 4).
Based on the spatial differences, we compute a spatially-averaged surface velocity for all parameter settings to obtain an integrative overview (Fig. 5). The analysis generally reveals that both the absolute spatially-averaged surface velocity (Fig. 5a- resolutions, solution differences increase between the stress regimes. For the finer resolutions BP-like produces higher mean surface velocities compared to FS. At the finest resolution (below the nominal resolution of the BM data-set), the mean velocity seem to reach a converged state for the BP-like stress regime (less pronounced for E = 0.1), but FS reveals a small dip towards 305 lower velocities. The tendency of the diverging relative differences at the finest resolution is therefore likely a result from the still changing FS solution under grid refinements. The obtained tendency below 150 m is intriguing as the differences solely arise from grid-refinements and not due to the input geometry.
Relative maximum differences between FS and BP-like are obtained for the stiff-ice case (E = 0.1). A discussion on why stiff ice is more sensitive will come in section 5.5. The spatially-averaged ice flow velocity difference is up to 7 m a 1 for 310 E = 0.1; the relative error between FS and BP-like is up to ⇠ 5.8%. The resulting difference for E = 1 shows the same trend but is of lower magnitude (up to ⇠ 1.5%). With softer ice (E > 1) the differences decrease. Interestingly, both stress regimes reveal an intermediate peak in mean surface velocity at l = 1600 m over all E values.

Outlet region
As for the ice stream region, an overview of input parameters and FS simulation results of the outlet NEGIS region are 315 compared for l = 6400 m and 150 m with P1P1GLS-weak and m = 3 (Fig. 6). Similarly, the slope of the bed topography ( Fig. 6a, d) and the friction coefficient (Fig. 6b, e) of the finer resolution reveals more small scale features and stronger amplitudes. In both resolutions, the simulated surface velocity fields agree fairly well and reproduce the NEGIS (Fig. 6c, f), however : . ::::::: However : in shear margins and gradients in along flow direction are more pronounced in higher resolution.
Differences of surface velocities between FS and BP-like are displayed for the 150 m grid resolution, E = 1 and P1P1GLS in 320 Fig. 7. Overall, spatial differences reveal a similar pattern among the employed friction exponents and friction implementations.
Largest differences occur in the vicinity of the 79NG grounding line (up to 570 m a 1 , see Tab. 3). For a power friction law (m = 3) the difference :::::::: differences : are larger compared to the linear friction law (m = 1). The P1P1GLS-weak setup produces slightly larger differences than P1P1-strong.  We compute relative spatially-averaged surface velocity differences for all parameter settings to obtain an integrative overview 325 (Fig. 8). The analysis discloses the clear trend that BP-like simulations lead to higher velocities than FS on average. Overall, the differences increase with higher resolution independent of the employed numerical characteristics. Interestingly, differences between FS and BP-like are more pronounced when using a power friction law compared to a linear friction law (up to 17% compared to 8%, respectively). Furthermore, the weak implementation produces larger differences compared to the strong scheme. The P2P1 discretization reveals differences that are in between P1P1GLS-strong and P1P1GLS-weak results. Differences still tend to increase further below the highest employed resolution. Notably, the area-averaged differences are higher in the outlet region than in the ice stream region (see Fig. 5f).

Impact on ice discharge
To estimate the impact on ice discharge, we calculated the ice mass flux along two flux gates of 79NG and ZI (see white lines in Fig. 2a). The flux gate of 79NG corresponds to the grounding line of the BM data set while the flux gate of ZI was shifted 335 approx. 3 km upstream from the calving front. The ice mass flux is calculated as Q = R vH% i ds, wherev indicates the depth- . For those :::: these ::: two ::::: stress :::::: regimes, relative differences are calculated between :::::: relative : to : FSand LTSML-like or FS wo bridging stresses.

Free surface evolution
In this section, we investigate how the velocity differences are represented in the solution of the ice surface position. The evolution of the free surfaces is given by the kinematic boundary condition 345 @z s @t whereȧ s (x, y) is the accumulation-ablation function at the surface. Here, we assume no melting or accumulation, i.e.ȧ s (x, y) = 0, and the so-called emergence velocity (v · n) determines the surface elevation change, @z s /@t. Since we do not run a forward model the computed free surface evolution is representative for the first time step. Figure 10 shows the emergence velocity for FS and the difference between FS and BP-like. The FS emergence velocity 350 reveals a somewhat noisy pattern. Due to :: the : hyperbolic nature of Eq. 14 numerical stabilization is necessary (Donea and Huerta, 2003;Riviere, 2008) for solving the free surface equation which certainly will smooth ::::::: smooths out sharp gradients The emergence velocity differences between FS and BP-like reveals a higher emergence velocity in FS, particularly at the grounding line at 79NG. According to Eq. 14, BP-like would therefore compute a surface elevation that is lower than the FS surface within the first time step. Although we are not running a time dependent simulation and, therefore, not capturing grounding line migration we evaluate the buoyancy balance from the stationary setup. The contact condition, i.e. P i > P w , at 365 the ice base (e.g. Durand et al., 2009). BP-like reveals a grounding line that lies further upstream (⇡ 0.5 km) than the one calculated from FS. However both are still located within the hinge zone of 79NG. Background image is a RADARSAT Mosaic (Joughin, 2015;Joughin et al., 2016).

Relevant physical processes
The reason for the difference in the simulated stress regime is a complex composition :: of : ice dynamic processes and how they are resolved in the ice flow model, such as the basal boundary condition and the internal deformation. In addition, basal 370 drag is coupled to internal deformation via the effective viscosity, making it difficult to separate the two entangled processes.
Nevertheless, we strive to give a flavour of the origins : of ::: the :::::: model ::::::::: differences : by analyzing key processes of ice dynamics below. Figure 11 and 12 shows :::: show : the difference of characteristic quantities between FS and BP-like for the ice-stream and outlet domain, respectively. First of all we would expect, that FS would lead to larger ice velocities than BP due to the effective 375 20 deformation rate (Eq. 9). FS captures the entire stress tensor, which is equivalent to making ice softer on average for the same flow parameters (assuming that strain-rate components shared by FS and BP-like have the same magnitude). Indeed, the depthaveraged viscosity,⌘, is basically lower for FS compared to BP-like (Fig. 11g-i and 12c,f). For E = 0.1 the depth-averaged viscosity is lower at :: in : some larger areas, while for E = 6 the pattern changes towards several smaller patches; the magnitude of⌘ among the E values changes by one order.

380
However, since FS generally reveals lower velocities than BP-like, the softer ice in FS appears to be compensated by another process. A fundamental control on ice flow is the basal drag indicating the glacier's behaviour for fast flow. The spatial differences of basal drag inferred from the two models reveal it is in the majority : is :::::::: generally : larger for FS compared to BP ( Fig. 11d-f and 12b, e). The spatial patterns generally coincide with the differences in surface velocity ( Fig. 4a and 7). In the ice stream region, the difference is largest for E = 0.1 (up to 140 kPa between FS and BP-like). In this setup, it seems that the 385 basal drag is an essential control of the ice flow. With softer ice, the difference in the basal drag decreases. In the outlet domain, the basal drag difference is more pronounced for the weak setup ::::::::::::: implementation :: of ::: the ::::: basal :::::::: boundary :::::::: condition with m = 3 compared to the strong setup :::::::::::: implementation : with m = 1. At the 79NG ::: and :: ZI : area, the maximum difference in basal drag is 316 kPa (99%) and 121 kPa (73%) for weak setup :: the ::::: weak ::::::::::::: implementation with m = 3; for the strong setup with m = 1 it is 227 kPa (72%) and 96 kPa (55%), respectively.

390
Inspired by the findings of Gudmundsson (2003) and Hindmarsh (2004), deduced from idealized setups, we explore whether the criteria when a BP solution is invalid (BP invalid in regions of high slip ratio, high aspect ratio and high topographic variations) hold in our realistic setup. In Fig. 13 and 14 we present two-dimensional histograms of the relative velocity difference between FS and BP-like for discrete bins of FS slip ratio v b /v s , the local aspect ratio, ✏ = H/l, and the bed slope, |rh b | (see detailed line plots Figs. S4-S7)). Simulated differences emerge with increasing aspect ratio, ✏. With increasing aspect ratio 395 neglecting terms of O(✏ 2 ), i.e. @v z /@x and @v z /@y, in the BP-like (and also BP) model becomes problematic and the solution inaccurate. Similarly, the error :::::::: difference between FS and BP occurs in regions with : is ::::: most :::::::::: pronounced :: in ::::: region :::: with :: a high slip ratio (i.e., near-plug-flow conditions prevail)and with : , ::::: while :::: with : a : higher resolution a larger number of grid nodes are involved in producing this error. At a first glance, the distribution of the relative error with respect to the bed slope unveils larger errors at smaller slopes. However, the distribution is certainly biased by large errors occurring in the fast flowing and 400 smooth trunk of the NEGIS region (see Fig. 2). By neglecting the bed slope bin from 0 to 0.05it turns out that : , ::: the higher bed undulations lead to larger errors in the finer :: at ::::: higher : resolutions; however, it :: the :::: bed :::: slope : seems to play a secondary role, compared to e.g. the slip ratio, as fewer grid nodes are involved. The analysis underlines that the source of stress regime deviations occurs in regions of high slip ratio, high aspect ratio and high topographic variability. :: A :::::: similar ::::::::: conclusion ::: for :::::::: mountain :::::: glaciers :::: was ::::: drawn :: in :::::::::::::::::: Le Meur et al. (2004) :: by :::::::::: comparing :: FS :::: and :::: SIA.
However, once the outlet regions of NEGIS are included (i.e., 79NG and ZI), the analysis shows that differences between FS and BP-like are large in the 79NG grounding zone. This finding is consistent with the work by Morlighem et al. (2010) at PIG.
The overestimation of BP-like surface velocities to FS in our study agrees with their necessity of a lower inferred basal drag for 420 FS to reproduce observed velocities compared to BP. They attribute this behavior to the developing bridging stresses in FS (t zz about 2% larger in FS compared to BP) due to the rising bed towards the grounding line, which in turn causes the reduction in basal drag. However, in our simulations a clear connection to bridging stresses could not be identified. Figure 15 ::::: (panels : a and c : ) shows the relative differences of t zz between FS and BP-like. Although the differences are in a similar range (±2%) as in Morlighem et al. (2010) there is no clear trend that could explain the lower velocities in FS. However, at PIG the bed rises 425 towards the grounding line, while at 79NG the bed mainly slopes down towards the grounding line. This geometric setting might control how the stresses develop. In our case, the differences mainly stem from the basal drag ( Fig. 11d-f and 12b, e) that is in the majority ::: for :::: most ::::: areas higher in FS compared to BP-like. The regions where BP-like simulates a lower basal drag than FS overlap with regions where the assumption that horizontal gradients of the vertical velocity are small compared to the vertical gradient of the horizontal velocity made to BP is invalid (Fig. 15b, d). These terms are of similar order and therefore, 430 dropping stress terms of O(✏ 2 ), i.e. @v z /@x and @v z /@y, lead to a lower basal drag in the BP-like case. Consequently, we assume that basal sliding is overestimated in BP-like simulations.
In order to highlight that bridging stresses play a minor role in our setup, we have run two other stress regime versions (executed with P1P1GLS-weak and up to resolution of l = 400 m). The first version keeps the bridging stresses in the third Figure 15. Comparison of (a, c) Relative difference of tzz at the ice base between FS and BP-like and (b, d) ratio of vertical gradient of the horizontal velocity to horizontal gradients of the vertical velocity (|(@vx/@z)/(@vz/@x)| + |(@vy/@z)/(@vz/@y)|) in absolute terms from the FS solution. Upper row (a, b) shows results for l = 150m, E = 1, m = 1 and P1P1GLS-strong. Lower row (c, d) shows results for l = 150m, E = 1, m = 3 and P1P1GLS-weak. Background image is a RADARSAT Mosaic (Joughin, 2015;Joughin et al., 2016 (2004)); the second version is a FS model without ('wo') bridging stresses. The area-averaged results reveal that LTSML-like is close to BP-like while 'FS wo bridging stresses' is rather similar to FS (Fig 8b).
In our diagnostic simulations, we compute lower ice mass fluxes with FS compared to BP-like. This supports the tendency that the sea level rise estimate of FS models is lower than the one of simpler models. (Seddik et al., 2012;Favier et al., 2014;Seddik et al., 2017). However, we observe a different response at 79NG and ZI. At 79NG, ice discharge differences are up 440 to 50%, while at ZI, ice discharge differences are about 5%. Also, at 79NG, differences seem to increase below the highest resolution employed here, while at ZI, the anticipated discharge difference increase appears to be moderate. A thorough analysis of why these two regions behave differently unveiled no significant differences between ice dynamic processes (e.g., basal drag shows a similar decrease). However, a potential source might be the geometric setting, e.g., the data coverage of available in-situ retrieved data leads to different topographic variability in both regions. To test whether the response behaviour is an 445 effect of topographic variability, we re-run a simulation with P1P1GLS-weak, m = 3 and l = 200 m but using ::: the geometry from a re-gridded l = 1600 m simulation. Re-gridding from a coarser geometry acts as smoothing and leads to a less variable geometry. The simulated discharge reveals that 79NG and ZI respond somewhat similar ::::::: similarly : (see triangles in Fig. 9b).
However, whether the different response at 79NG and ZI with the high-resolution bedrock stems from in-situ data coverage or :: the :::: fact ::: that : ZI has generally a smoother geometry than 79NG must be examined in further studies. 450 We tested the numerical robustness of the underlying finite element method by employing different discretization schemes and :::::::: analysing : how the non-penetration condition at the ice base is enforced. In general, the different schemes show a similar velocity difference pattern between FS and BP-like. For m = 1 and m = 3 the area averaged velocity differences for l = 150 m are within 6.1% to 7.6% and 13.1% to 16.2%, respectively. Therefore, the choice of the employed friction law plays a larger role than the choice of numerical characteristics tested in this study.

455
However, some tests regarding the numerical robustness are missing. Because we keep the number of layers constant over the entire modelling domains, large variations in the aspect ratio occur when refining the mesh. The element aspect ratio is critical for the numerical error and the classical BP system would suffer less under bad element aspect ratios (Larour et al., 2012, Sect 5.1.3 therein). In order to assess the impact of element aspect ratios, we carried out a study where we varied the number of vertical layers from 3 to 11. For this test we re-run the simulations up to a resolution of l = 400 m with the P1P1GLS-strong 460 discretization scheme and m = 3 applied to the outlet domain. The area-averaged surface velocities difference between FS and BP-like are very small (Fig. S8). This demonstrates that the comparison is not suffering under bad element aspect ratios and the numerical error behave similar in FS and BP-like.
The basal boundary conditions are retrieved by inverting for the basal friction coefficient with BP employing a resolution of l = 800 m. That means, that the boundary conditions are not consistent with FS or with other resolutions. In order to study 465 whether this inconsistency influences the results, we repeat a set of experiments with two different fields for the friction coefficient. To this end, we run two ISSM friction inversions with (1) E = 1, l = 200 m, m = 3 and BP and (2) E = 1, l = 800 m, m = 3 and FS and feed the inferred friction fields to COMSOL. The area-averaged differences between FS and BP-like of the three setups show very minor differences (Fig. S9).
By implementing FS and BP-like solvers using the COMSOL finite element package, we have developed a useful tool for 470 exploring differences between FS and BP-like, without complications due to differing numerics. The developed model is also very flexible in order to test various ice flow components (e.g., by simply including bridging stresses to convert the BP-like model to LTSML-like). In our model setup, we solve for the same physical problem as in the 'classical' BP scheme but without taking advantage of the reduced number of equations to be solved. Regarding the model physics, our results compare well with other BP implementations, but our BP-like model is more expensive. Therefore, the developed BP-like model is not intended 475 for long time integrations or large ensembles.
One limitation of our study is that the simulations are not prognostic, i.e. we have not investigated how the solution differences propagate to or interact with other components of an ice sheet model, e.g. by coupling to the ice thickness evolution. The calculation of ice discharge and the emergence velocity demonstrates that the stress regimes show distinct discrepancies in a key process, but whether these differences are accumulated in the glaciers mass balance over transient runs, eventually leading 480 to substantially different estimates of ice mass changes, or fading out must be postponed to :::::::::: investigated :: in future studies.
A process not considered here is the thermo-mechanical coupling. For instance, within a thermo-mechanical coupling ::::: model, the feedback of temperature on the ice viscosity, in turn, changes the ice flow. However :::::::::: Additionally, there are several feasible processes how the stress regime differences can modify the thermal regime: (1) the different representation of the basal drag presumably generates a different amount of frictional heating, (2) the different vertical deformation produces a different internal 485 heat source rate, and (3) the internal advection scheme leads to different heat transport. Though the role of these feedbacks is not explored here, there are indications that small initial differences become much larger over long time integrations (Zhang et al., 2015).

Conclusions
We compared two approaches to represent ice flow dynamics, namely the FS and the BP-like formulation, for two different 490 regions of the NEGIS. Both stress regimes are implemented within the same dynamic ice flow model to enable a consistent comparison where numerical differences are largely eliminated. Our comparison experiments between FS and BP-like stress regimes unveil that BP-like tends to overestimate surface velocities and the resulting ice discharge. The models disagreements start at a horizontal resolution of around 1000 m and continuously increase with finer horizontal resolution (down to 100 m).
In contrast, ten times stiffer (E = 0.1) ice leads to enhanced ice flow velocities in BP-like compared to FS. The classical rheology case (E = 1) reveals very moderate differences between FS and BP-like. However, once the outlet regions of NEGIS are included, i.e. 79NG and ZI, ice flow treatment with FS is essential as ::::::: important : surface velocity differences of :::: occur :: (up to 570 m a 1 at the grounding line of 79NGoccur). Numerical characteristics play a minor role compared to the choice of the 500 friction type, i.e. linear versus power-law friction; : : model differences are greater :::::::::: particularly :::::::: important : when a power-law friction is applied :: (as ::::::: opposed :: to :: a ::::: linear ::::::::: approach). The basal drag and topographic undulations are identified as the general causes for the model disagreements.
Whether FS will reduce uncertainties of future sea-level projections must be examined further in subsequent transient studies. Although computer power is increasing, transient simulations in FS with a horizontal resolution of about 100 m in dynamic, 505 relevant regions are still challenging. Our diagnostic simulations and previous studies indicate that FS and BP differ significantly in regions with a grounding zone, : and results from simpler :::::: non-FS models should therefore be viewed with caution.
. Code and data availability. COMSOL Multiphysics© is a commercial software that is not freely available. The models used here are accessible for COMSOL users upon request to the authors. The ice flow model ISSM is open source and freely available at https://issm.jpl.nasa.gov/ 510 (last access: April 6, 2022; Larour et al. (2012). Simulations results are available on Zenodo with digital object identifier: https://doi.org/10.5281/zenodo.640 (Rückamp et al., 2022). : .

Appendix A: Inversion with ISSM
The simulations presented make use of a basal friction coefficient, k 2 , that is inferred by an inversion method. For the inversion 515 of the basal friction coefficient, we operate the Ice-Sheet and Sea-level System Model (ISSM, Morlighem et al., 2010;Larour et al., 2012), an open source finite element flow model appropriate for continental-scale and outlet glacier applications. The setup is as described above, i.e. (1) domain outline and geometry data set is similar as in the COMSOL simulations (see Fig. 2), (2) ice rheology is taken from SICOPOLIS spinup (Rückamp et al., 2019) and interpolated onto the computational mesh and (3) the friction law follows the form in Eq. 6. But for the ease of computational time model calculations are performed on 520 a structured finite element grid with a horizontal resolution of 1 km and the BP-stress regime. The inversion is conducted separately for each enhancement factor E = 0.1, 1, 3 and 6 and the sliding exponents m = 1 and m = 3.
Within the inverse problem a cost function (J), that measures the misfit between observed, v obs = (v obs x , v obs y ), and simulated velocities, v = (v y , v y ), is minimized. We use the observed velocities from the MEaSUREs project (Joughin et al., 2016(Joughin et al., , 2018 as target within the inversion. The cost function is composed of two terms which fit the velocities in fast-and slow-moving 525 areas. A third term is a Tikhonov regularization to avoid oscillations due to over fitting. The cost function is defined as follows: where " is a minimum velocity used to avoid singularities and s and b are the ice surface and ice base, respectively. An Lcurve analysis was performed to pick the Tikhonov parameter t (not shown). We obtained a good agreement to the observed 530 velocities by choosing with 1 = 1 ⇥ 10 3 , 2 = 1 ⇥ 10 4 and t = 1 ⇥ 10 4 (see Fig. A1 and A2).
Author contributions. MR setup the COMSOL FS and BP model and the ISSM model, conducted the simulation and wrote large parts of the manuscript. AH designed the study. MR, TK and AH analyzed the results. All author contributed in writing the manuscript.