Regional modeling of the Shirase drainage basin , East Antarctica : full Stokes vs . shallow ice dynamics

A hierarchy of approximations of the force balance for the flow of grounded ice exists, ranging from the most sophisticated full Stokes (FS) formulation to the most simplified shallow ice approximation (SIA). Both are implemented in the ice flow model Elmer/Ice, and we compare them by applying the model to the East Antarctic Shirase drainage basin. First, we apply the control inverse method to infer the distribution of basal friction with FS. We then compare FS and SIA by simulating the flow of the drainage basin under present-day conditions and for three scenarios 100 years into the future defined by the SeaRISE (Sea-level Response to Ice Sheet Evolution) project. FS reproduces the observed flow pattern of the drainage basin well, in particular the zone of fast flow near the grounding line, while SIA generally overpredicts the surface velocities. As for the transient scenarios, the ice volume change (relative to the constantclimate control run) of the surface climate experiment is nearly the same for FS and SIA, while for the basal sliding experiment (halved basal friction), the ice volume change is ∼ 30 % larger for SIA than for FS. This confirms findings of earlier studies that, in order to model ice sheet areas containing ice streams and outlet glaciers with high resolution and precision, careful consideration must be given to the choice of a suitable force balance.


Introduction
The Shirase drainage basin (SDB) in East Antarctica covers an area of ∼ 2 × 10 5 km 2 .The basin is characterized by slow, vertical-shear-dominated ice flow inland, which converges into the Shirase Glacier, one of the fastest-flowing glaciers in Antarctica.As the glacier calves nearly 90 % of total ice discharged from the basin, the converging flow regime and the fast-flow feature play a crucial role in the mass budget of this region.The dynamic conditions of the drainage basin were extensively measured, with flow speeds at the grounding line between 2300 and 2600 m a −1 (Rignot, 2002;Pattyn and Derauw, 2002;Nakamura et al., 2007;Rignot et al., 2011) and an ice mass discharge of about 10 to 14 Gt a −1 (Fujii, 1981;Nakamura et al., 2010Nakamura et al., , 2016)).Pattyn and Naruse (2003) studied the complex ice flow in the Shirase Glacier catchment.They showed that it originates from the Mizuho triangulation chain (Naruse, 1978), and between there and the grounding zone the ice flow regime is mainly dominated by longitudinal stresses and basal sliding.
Due to the nature of the ice flow, the SDB is an interesting location for comparing different formulations of ice dynamics.Such comparisons have already been the subject of several studies in the past.In an early work by Mangeney and Califano (1998), the shallow ice approximation (SIA; Hutter, 1983) was used to model ice flow, including induced anisotropy, in a vertical plane with either flat or uneven bedrock under steady state conditions.Further simplifications were Newtonian (linear-viscous) ice rheology and isothermal conditions.Results were compared to respective full Stokes (FS) solutions.For flat bedrock, the SIA gave excellent results for both isotropic and anisotropic conditions.For uneven bedrock, the SIA was found to be still acceptable for isotropic conditions but not appropriate for anisotropic conditions.Hindmarsh (2004) and Gudmundsson (2008) performed further investigations of different approximations of ice dynamics for idealized geometries based on perturbation theory.Community-wide intercomparisons of higher-order and FS ice sheet models for idealized 2-D and 3-D geometries were conducted by the ISMIP-HOM and MISMIP projects (Pattyn et al., 2008(Pattyn et al., , 2012)).A particularly interesting finding of ISMIP-HOM was that the FS models produced very consistent results, while the spread across other higherorder models was larger.SIA models were not part of this investigation.
Comparisons of the impact of different formulations of ice dynamics on simulation results have also been carried out for real-world problems.Morlighem et al. (2010) simulated the present-day velocity field and the spatial pattern of basal drag for the West Antarctic Pine Island Glacier, using a finite element framework, a control method and three different representations of ice dynamics, namely FS, the shelfy stream approximation (SStA) (MacAyeal, 1989) and the Blatter-Pattyn (a.k.a.first-order) approximation (Blatter, 1995;Pattyn, 2003).They found that SStA and Blatter-Pattyn produce a large basal drag near the grounding line, while the FS solution has almost none.The discrepancy is attributed to neglected bridging effects by the simpler dynamics (SStA, Blatter-Pattyn) in an ice stream region of rapidly varying thickness.Seddik et al. (2012) applied the models Elmer/Ice with FS and SICOPOLIS with SIA to the entire Greenland ice sheet, forced by scenarios from the Sea-level Response to Ice Sheet Evolution (SeaRISE) project (Bindschadler et al., 2013).They found that Elmer/Ice with FS was more sensitive than SICOPOLIS with SIA for the dynamical forcing (doubled basal sliding) experiment but less sensitive for the direct global warming scenario.However, comparability between FS and SIA was limited by the use of the two different models because differences can also arise from different numerical techniques.Fürst et al. (2013), who also studied the centennial evolution of the entire Greenland ice sheet, overcame this disadvantage by using only one model.Their five approximations to the force balance were between SIA and Blatter-Pattyn, thus excluding FS.
Here, we apply the model Elmer/Ice (Gagliardini et al., 2013;Elmer/Ice Project, 2017) to the SDB.Elmer/Ice can be run in either FS or SIA mode.This allows us to compare FS and SIA results computed by the same model with the same mesh.To our best knowledge, this is the first time that such a comparison is done for these two "end-member" force balances within one model applied to a large, ice-sheet-sized area.We consider the results obtained with the FS dynamics to be more physically adequate and use them as a reference to infer the basal friction coefficient with the InSARbased Antarctica ice velocity map provided by Rignot et al. (2011) and a control inverse method.The optimized distribution of the basal friction coefficient is then used to compute the present-day flow of the SDB with both FS and SIA, thus allowing a direct comparison between the two methods.We also carry out three transient experiments 100 years into the future, of which the forcings are taken from the suite of SeaRISE experiments.The surface velocities and volume changes produced by FS for each experiment are compared to those produced by SIA.In addition, we investigate the differences in the slip ratios (the ratio of basal to surface velocity) of the FS and SIA solutions.

FS formulation
Since ice is an (almost) incompressible material, conservation of mass entails that the velocity field is solenoidal.Furthermore, the acceleration (inertial force) is negligible.The FS equations are thus where T is the Cauchy stress tensor, g = −ge z the gravitational acceleration vector pointing downward, ρ the density of ice and v the velocity vector (e.g., Gagliardini et al., 2013).
The Cauchy stress tensor can be expressed as T = −pI+T D , where T D is the traceless stress deviator, p = −tr T/3 is the pressure and I is the unit tensor.Ice rheology is represented by a nonlinear Norton-Hoff-type flow law (commonly referred to as Glen's flow law): where D is the strain rate (stretching) tensor.The effective viscosity η is defined as where d is the effective strain rate, A is the rate factor and E is the flow enhancement factor (e.g., Gagliardini et al., 2013).
The temperature equation follows from the general balance equation of internal energy and is (e.g., Greve and Blatter, 2009) where κ and c are the heat conductivity and specific heat of ice, respectively (Table 1).

Shallow ice approximation
In contrast to the FS formulation, which neglects none of the stress components, SIA assumes that grounded ice flow is governed only by ice pressure and the vertical shear stresses.This yields the following equations for the velocities v = A ) and the pressure p (Greve and Blatter, 2009): where, in the Cartesian coordinate system (x, y, z) with z positive upward, z = h(x, y, t) denotes the free surface.The viscosity term η is computed with Eq. ( 4) and the SIA effective strain rate Equations ( 6) and ( 7) are solved by formulating a degenerated Poisson equation that is discretized with the finite element method (Appendix A).
The SIA formulation of the temperature equation (Eq.5) follows from the negligible horizontal heat conduction assumption.Then, by expressing the effective strain rate d in terms of the effective stress σ e in the dissipation term (Greve and Blatter, 2009), the equation is 2.2 Boundary conditions

Ice surface
For both FS and SIA, we assumed a stress-free ice surface (atmospheric pressure and wind stress neglected).In the case of FS, the evolution of the upper surface h(x, y, t) is governed by the kinematic boundary condition: where a s (x, y, t) is the accumulation-ablation function or surface mass balance.The free surface evolution equation employed in SIA results from the ice thickness equation (Greve and Blatter, 2009): where a b is the basal melting rate, z = b(x, y, t) is the ice base, and D is defined by the function where H is the ice thickness and β is the basal friction coefficient (Sect.2.3).
The mean annual surface temperature T ma and the mean summer surface temperature T ms for Antarctica were parameterized as functions of surface elevation h, latitude φ and time t.Following Fortuin and Oerlemans (1990), the current mean annual surface temperature is  Huybrechts and de Wolde (1999), the current mean summer surface temperature is where the temperature constant is d ms = 16.81 • C, the mean lapse rate is γ ms = −0.00692• C m −1 and the latitude coefficient is c ms = −0.27937• C ( • S) −1 .For climatic conditions different from the current, mean annual surface temperature was modified by the additive anomaly T ma (φ, h, t) and mean summer surface temperature by T ms (φ, h, t).

H. Seddik et al.: Modeling of the Shirase drainage basin
The surface mass balance a s (snowfall minus melting) was computed with the SeaRISE specifications for Antarctica (Sato and Greve, 2012;Bindschadler et al., 2013;Nowicki et al., 2013).For the current mean annual precipitation rate, P ma,present (λ, φ) (function of longitude λ and latitude φ), we used the data from Le Brocq et al. (2010) (newly compiled from Arthern et al., 2006).For climatic conditions different from the current, they were modified by multiplicative anomalies (factors) f P ma (λ, φ, t).Surface melting was parameterized by a positive degree-day (PDD) method (Reeh, 1991) supplemented by the semi-analytical solution for the PDD integral (Calov and Greve, 2005).The PDD factors were β ice = 8 mm ice equiv.d −1 • C −1 for ice melt and β snow = 3 mm ice equiv.d −1 • C −1 for snow melt (Huybrechts and de Wolde, 1999).Furthermore, the standard deviation of short-term, statistical air temperature fluctuations was σ = 5 • C (Huybrechts and de Wolde, 1999), and the saturation factor for the formation of superimposed ice was chosen as P max = 0.6 (Reeh, 1991).Conversion from precipitation to snowfall (solid precipitation) was done on a monthly basis with the empirical relation by Marsiat (1994).

Ice bed
On the short timescales of our simulations, the bed topography b(x, y) can be assumed to be rigid (isostatic compensation neglected) and thus at all times equal to the prescribed initial condition (the last term on the right-hand side of Eq. 11 therefore vanishes).For FS, basal shear stresses and basal velocities are related according to where β is the basal friction coefficient (Sect.2.3), n is the normal unit vector and t is the tangential unit vector aligned with the flow direction.For the SIA, Eq. ( 15) simplifies to a statement about the horizontal components of the basal velocity: In both cases, the normal velocity at the ice bed interface is given by where the basal meting rate a b is determined by the energy jump condition (Greve and Blatter, 2009): where L is the latent heat of ice and q geo the spatially varying geothermal flux (Shapiro and Ritzwoller, 2004).
Finally, the temperature equations (Eqs.5, 9) were solved assuming a temperate base everywhere (see Sect. 2.4).This yields the Dirichlet-type condition where T m is the temperature at the pressure melting point.The pressure p for SIA was assumed to be hydrostatic (p = ρg(h − z)) in the calculation of T m .

Side boundaries
Our computational domain (see Sect. 2.4 for details) has two different types of side boundaries, namely two lateral flow lines and a downstream grounding line.For the FS case, we assumed vanishing normal velocities and zero stress conditions for the flow lines: The additional condition (zero tangential velocity) was only applied to nodes with an ice thickness of less than 10 m, otherwise a tangential free slip condition was assumed.At the grounding line, the normal component of the stress vector is equal to the hydrostatic pressure exerted by the ocean for nodes under the water line and zero elsewhere, and the tangential component is zero: where ρ w = 1025 kg m −3 is the sea water density and l w is the sea level.
For the SIA case, neither of the dynamic conditions of Eqs. ( 20)-( 22) are required.Instead, a boundary condition for the free surface evolution (Eq.11) is needed.For the lateral flow lines, we assume a zero surface gradient in normal direction, which, in the SIA, is equivalent to the zero-normal-velocity condition (Eq.20 1 ).At the grounding line, we keep the ice surface h fixed at its observed, present-day distribution.

Control inverse method
The control inverse method, introduced by MacAyeal (1993), was implemented for the FS mode of Elmer/Ice by Gillet-Chaulet et al. (2012).It relies on the computation of the FS adjoint (Morlighem et al., 2010) and on the definition of a cost function expressed as the difference between the norm of the modeled and observed horizontal velocities (v h and v obs h , respectively): where s is the footprint of the ice surface in the x-y plane.The regularization procedure is constructed by adding a smoothness constraint to the cost function in the form of a Tikhonov regularization term: where b is the footprint of the ice base in the x-y plane.
The parameter α is used to avoid negative values of the basal friction coefficient β and is defined as The total cost function is then given by where λ is a positive ad hoc parameter.Minimizing the cost function J tot with respect to α was implemented using the quasi-Newtonian routine M1QN3 (Gilbert and Lemaréchal, 1989;Gillet-Chaulet et al., 2012).

Computational domain and meshing
The Elmer/Ice model was applied to the SDB.The presentday surface and bed topographies were extracted from the Bedmap2 data set by Fretwell et al. (2013).Our domain originates at Dome Fuji, is delimited by two approximate flow lines and converges into Shirase Glacier (Fig. 1a).The downstream end of the domain is the grounding line of Shirase Glacier, which means that we have removed the small ice shelf of the glacier flowing into Havsbotn (the bay that forms the head of Lützow-Holm Bay).The bed topography (Fig. 1b) shows by trend a decrease in elevation from the interior of the ice sheet towards the coast, and the bed of the most downstream regions including Shirase Glacier is below sea level.
For this domain, we solved Eqs.(1)-( 5) for the FS case and Eqs. ( 6)-( 9) for the SIA case with the finite element method.Since radar data indicate that the basal temperature in the area is at the melting point everywhere except for the highest peaks of the bed topography (Fujita et al., 2012), for simplicity we assumed a temperate base for the entire domain (Eq.19).The finite element mesh was constructed by using an adaptive meshing method with horizontal resolutions up to 200 m near the grounding line, 300 m near Dome Fuji (Fig. 1c) and coarser resolutions down to 15 km elsewhere.
We first constructed a 2-D footprint mesh and then extended it to give a 3-D mesh with 56 394 elements and 11 equidistant, terrain-following layers.
The nonlinearity of the model equations was dealt with by Picard iteration as in Seddik et al. (2012).Stabilization methods (Franca et al., 1992;Franca and Frey, 1992) were applied to the finite element discretization.The resulting system of linear equations was solved with a direct method using the parallel sparse direct solver MUMPS (Amestoy et al., 2001(Amestoy et al., , 2006)).The position of the ice front was assumed to be fixed.A minimum ice thickness of 10 m was applied everywhere and for all times.This means that locations that are initially glaciated were not allowed to become completely ice free but rather are always covered with at least 10 m of ice.The minimum thickness is required in order to avoid having 2-D (rather than 3-D) elements, which would be treated as degenerated elements during the finite element assembly so that the assembly would fail.The reason why we chose the value of 10 m is to avoid a too-low aspect ratio (thickness-to-width ratio) of the 3-D elements, which can cause the numerical solution to become unstable.

Present-day state
In order to infer an initial spatial distribution of the temperature field that contains a historical footprint of ice sheet evolution during glacial cycles, a spin-up of the whole ice sheet is generally needed.However, this procedure does not always produce a distribution in satisfactory agreement with present conditions, particularly at the ice base (Sato and Greve, 2012), and the simulation periods of a glacial cycle needed for a proper spin-up state imply the need to use an SIA model to reduce the computational cost.Instead, here we used a simple 1-D steady state approach to initialize the temperature distribution.Equation (5) simplifies to Following Greve et al. (2002), the vertical velocity v z in Eq. ( 28) was computed assuming a Dansgaard-Johnsen distribution (Dansgaard and Johnsen, 1969).This consists of a constant vertical strain rate ∂v z /∂z from the free surface down to z = 0.25H and a linearly decreasing vertical rate below.
Using the temperature field computed by solving Eq. ( 28) for each ice column of the domain, we applied the control inverse method (Sect.2.3) to infer in FS the spatial distribution of the basal friction coefficient β (starting with a spatially uniform value of β = 10 −4 MPa m −1 a as an initial guess) and the present-day 3-D velocity field.The observed velocities used for the inversion were given by the InSARbased Antarctica ice velocity map (Rignot et al., 2011, and  Fig. 1c), and the basal melting rate a b in the boundary condition (Eq.17) was neglected (thus the normal velocity at the ice base was set to zero).For comparison, we also computed the SIA velocity field based on the same distribution of β.

Future climate experiments
The obtained present-day state of the SDB served as initial condition for runs into the future.We used a subset of the SeaRISE experiments defined in Bindschadler et al. ( 2013 -Experiment CTL ("constant-climate control run") starts at the present (more precisely, the epoch 1 January 2004 0:00 corresponding to t = 0) and runs for 100 years, with the climate held at the current state.
-Experiment S1 ("basal sliding experiment") is constantclimate forcing with increased basal lubrication.This was implemented in Elmer/Ice (for both FS and SIA) by halving the basal friction coefficient (approximately doubling the basal sliding) everywhere in the domain.
-Experiment C2 ("surface climate experiment") is 1.5 × AR4 surface climate forcing, starting with the current state as in CTL and S1, but with climatic forcing (changes in mean annual temperature, mean summer temperature and precipitation; see Sect.2.2.1).The changes were derived from an ensemble average of 18 AR4 models, run for the period 2004-2098 under the A1B emission scenario.For runs beyond 2098, conditions were held at those of 2098.
Note that, as defined by SeaRISE, the combination experiment R8 designed for the IPCC AR5 does not include increased basal sliding during the first 100 years and the increase in basal melting is applied only to the Amundsen Sea and the Amery Ice Shelf.Therefore, experiment C2 is equivalent to R8 for the SDB and the period considered here.All runs were carried out in both FS and SIA mode.For the FS runs, the numerical time step was t = 0.02 a for the first 5 years of model time and t = 0.05 a thereafter, while for the SIA runs it was t = 0.05 a for the entire 100 years.4 Results

Present-day state
As described previously (Sect.3.1), we used the control inverse method to compute the spatial distribution of the basal friction coefficient and the 3-D velocity field in FS.The parameter λ in Eq. ( 27) was chosen using the L-curve method (Hansen, 2001;Gillet-Chaulet et al., 2012), which plots the regularization, i.e., the term J reg as a function of the term J o in Eq. ( 27).The L curve obtained with the control inverse method is shown in Fig. 2. When we increased λ from 0 to 10 10 , the roughness of the basal friction coefficient distribution, represented by J reg , decreased by several orders of magnitude, while the mismatch between the observed and the computed velocities (J o ) remained essentially constant.For higher values of λ, the roughness J reg continued to decrease but at the cost of an increasing J o as the basal friction coefficient distribution became too smooth.We therefore chose λ = 10 10 as the optimal value.Figure 3a-c show the comparison between the observed velocities and the velocities computed using FS and the control inverse method.The overall agreement is good; in particular, the model reproduced well the fast-flowing ice towards the grounding line and the slower ice motion elsewhere in the drainage basin.The most significant mismatch occurs at the upstream end of the domain in the vicinity of Dome Fuji, where ice flow is slowest.This is due to the formulation of the control inverse method.The cost function works better in regions with high velocities than in slow-moving areas be-cause of the quadratic dependence of Eq. ( 24) on the velocity mismatch.This leads to a larger contribution of mismatches that occur at high velocities to the cost function, while the contribution of mismatches that occur at low velocities is limited.Another notable mismatch occurs in the right (eastern) part of the domain about one-half to two-thirds in downstream direction.
The spatial pattern of the basal friction coefficient obtained by the optimization (Fig. 3d) is characterized by a region with low friction at Shirase Glacier and the adjacent, downstream area of the drainage basin.Further upstream, the basal friction generally increases, with the exception of the vicinity of Dome Fuji, where the friction is again low.As stated previously, the latter result is not very relevant because the cost function is rather insensitive to velocity mismatches in this slowest-flowing area, which limits the quality of the obtained friction coefficients.Along the lateral boundaries of the domain, the basal friction is also low.This is likely related to the choice of the boundary conditions (Eq.20) and thus an artefact rather than a real phenomenon.
An initial velocity distribution was also computed with SIA, using the basal friction coefficient inferred from FS and the control inverse method (Fig. 3d). Figure 4a shows the resulting surface velocities.They were generally higher than those for FS (see also Fig. 4b).In particular, SIA did not reproduce so well the sharp contrast between the narrow fast-flowing region near the grounding line and the slowerflowing ice further upstream.Figure 4c, d show the surface velocity ratio FS / SIA vs. the FS slip ratio (panel c) and the SIA slip ratio (panel d).In both cases, the scatter of the surface velocity ratio increases as the slip ratio increases from near zero to unity.This is because the SIA is best suited for areas where the slip ratio is low and the flow regime is dominated by vertical shear, while in areas with a high slip ratio near-plug-flow conditions prevail for which the SIA is less appropriate (e.g., Blatter et al., 2011).

Future climate experiments
These experiments with evolving ice surface were carried out with the previously computed present-day state as the initial condition.For both FS and SIA, the distribution of the basal friction coefficient obtained by FS and the control inverse method was used (Fig. 3d), and it was kept constant with time.
Figure 5 depicts the simulated evolution of the ice volume (V ) of the Shirase drainage basin (panel a) and the volume relative to CTL (panel b) for the three scenarios.For CTL-FS (Fig. 5a, solid blue line), the average ice volume change (net mass balance) in the first year of the simulation is +2.0 × 10 −3 mm SLE a −1 (where SLE is sea level equivalent), and in the first 5 years it is −14.6 × 10 −3 mm SLE a −1 .So the mass balance is very slightly positive initially, but then turns to a negative value.Both values are consistent with the observed mass balance reported by Nakamura et al. (2016)  −1.9±3.3Gt a −1 , which is equivalent to −5.3×10 −3 ±9.1× 10 −3 mm SLE a −1 .By contrast, the initial ice volume change produced by CTL-SIA (Fig. 5a, dashed blue line) is much more negative than that of CTL-FS and by an order of magnitude larger than the most negative value in the range of uncertainty given by Nakamura et al. (2016).
For both FS and SIA, all scenarios produce a volume loss over the 100 years of model time.However, it is much smaller for FS than for SIA, so that the results are more dependent on the used model dynamics than on the scenario.The more rapid volume loss in the SIA experiments is certainly related to the larger flow velocities produced by the SIA for the initial state (see above).
The difference between FS and SIA becomes much smaller if we consider volume changes relative to the respective control run ( V ) rather than absolute volumes (Fig. 5b).It is evident that the basal sliding experiment S1 produces a much larger response than the surface climate experiment C2 over the entire model time.This finding agrees qualitatively with the simulated response of the entire Antarctic ice sheet (AIS) to these forcings (Bindschadler et al., 2013).For both FS and SIA, C2 produces first a slightly lower volume than CTL, but V goes through a minimum after ∼ 50 a and catches up thereafter.Apparently, after ∼ 50 a, the increased precipitation of the surface climate scenario C2 outweighs the initially dominant impact of warming and surface melting.After 100 a, the volume changes are -S1 (basal sliding exp.) -CTL (control): V FS ∼ −3.2 mm SLE, V SIA ∼ −4.4 mm SLE.
Figure 6 shows the surface velocities and slip ratios (ratio of basal to surface velocity) computed with FS dynamics for each scenario at t = 100 years.The distribution of the surface velocities obtained with the control experiment (CTL; Fig. 6a) is similar to the initial velocities computed with the control inverse method (Fig. 3b).The main difference occurs near Dome Fuji, where the ice flow has slowed down, which agrees better with the expected slow flow in the vicinity of topographic domes.The slip ratio (Fig. 6b) is characterized by a complex fine structure but increases generally from the interior ice sheet towards the grounding line.In the narrow channel of Shirase Glacier, the slip ratio is close to unity, which means that the fast flow of the glacier is controlled mainly by basal sliding (as expected).Experiment S1 (halved basal friction) produces a marked acceleration of the ice flow (Fig. 6c), particularly near the grounding line.The slip ratio distribution (Fig. 6d) indicates that basal sliding contributes significantly to the ice flow for practically the entire drainage basin.By contrast, experiment C2 does not differ greatly from CTL (the surface velocities and slip ratios are similarly distributed; Fig. 6e and f).The slip ratios of all three experiments exhibit boundary effects along the lateral margins of the domain.These are artifacts resulting from the low basal friction near the margins (see Sect. 4.1).
The surface velocities and slip ratios obtained with SIA dynamics are shown in Fig. 7.The differences compared to the respective FS solutions are notable.For all three scenarios, SIA produces higher surface velocities than FS.Further, the sharp contrast between the narrow fast-flowing region of Shirase Glacier and the slower-flowing ice further upstream is not so well reproduced; the zone of fast flow is smeared out over a larger area.This is the same behavior we found for the control inverse method (Fig. 4a).The slip ratio distributions of the SIA solutions are also very different from those of FS.Not considering the vicinity of Dome Fuji (where the distribution of the basal friction coefficient β obtained by the control inverse method is not well constrained; see Sect.4.1), the slip ratios show generally lower values.The only exception is the channel of Shirase Glacier where the slip ratios approach unity, similar to the FS solutions.However, as we have seen previously, the SIA solution is more sensitive to the S1 perturbation (halved basal friction) than the FS solution.Therefore, the smaller SIA slip ratios are due not to less basal sliding but rather to more internal deformation.The latter results from the fact that, in the SIA, the flow is controlled by the local surface topography only and does not experience any resistance due to longitudinal stresses and horizontal shear, which is accounted for by the FS dynamics.

Discussion
Our finding that the volume change of the Shirase drainage basin relative to CTL is small for the surface climate experiment C2 is consistent with the SeaRISE findings.For both the entire AIS (Bindschadler et al., 2013) and the extended Queen Maud Land basin that includes the SDB (Nowicki et al., 2013), the SeaRISE results for six different models showed a generally small sensitivity of the simulated ice volume to this forcing with an even unclear sign (some models predicted a volume increase, others a decrease).By contrast, all SeaRISE models produced a significant volume decrease for the S1 basal sliding experiment, which also agrees with our results.For the entire AIS, the range of SeaRISE S1 results after 100 model years was ∼ −300 to −75 mm SLE (Bindschadler et al., 2013).Scaling this down to the SDB by simply using the volume ratio V SDB /V AIS ≈ 1.148 m SLE/58.3m SLE ≈ 0.02 (V SDB : our value, see Fig. 5a; V AIS by Vaughan et al., 2013) yields a range of ∼ −6 to −1.5 mm SLE.The volume changes of S1 relative to CTL after 100 years, as we show in Fig. 5b, are well in the center of this range for both FS and SIA.Especially in the fast-flowing region near the grounding line, the Shirase Drainage Basin is characterized by a complex stress regime (Naruse, 1978;Pattyn and Naruse, 2003) that the SIA force balance cannot capture fully due to its neglect of longitudinal stresses and horizontal shear.We investigate the difference between FS and SIA further by considering the stress ratio computed at the ice base, where the effective stress σ e is and the τ ij and τ D ij are the components of the stress tensor T and deviator T D , respectively (see Sect. 2.1.1).Since the longitudinal stresses τ D xx and τ D yy as well as the horizontal shear stress τ xy are neglected in the SIA, the SIA stress ratio is equal to unity everywhere.Thus, the deviation of the FS stress ratio from unity reveals where the stress components neglected in the SIA method play an important role for the dynamics of the drainage basin.
Figure 8 depicts the stress ratio at the base of the SDB that results from the FS solution for the control run CTL (stress ratios for S1 and C2 are not shown because they are very similar).The stress ratio is close to unity everywhere except for two regions, namely in the vicinity (several tens of kilometers) of Dome Fuji and in the narrow trough near the grounding line where Shirase Glacier is located.This numerical finding agrees perfectly with the theoretical understanding that the ice flow regime deviates most from simple, bed-parallel shear (and thus the SIA is violated) in the vicinity of ice domes and for fast-flowing ice streams with high slip ratios (e.g., Blatter et al., 2011).
When starting transient numerical simulations of ice sheets, depending on the employed initial conditions, spurious noise in the computed velocity field (e.g., Herterich, 1988) and/or undesirable transients in ice geometry and other state variables (e.g., Perego et al., 2014) may occur.Such initialization shocks can be due to an unsmoothed, slightly noisy initial geometry or an imbalance of the initial thermodynamic condition with the applied surface mass balance.Strategies to avoid them are either "relaxation" or "flux correction" approaches (e.g., Edwards et al., 2014).However, in our case, the topography constructed as described in Sect.2.4 was sufficiently smooth to render a relaxation or flux correction unnecessary.This becomes evident in Fig. 5: all ice volumes computed by both FS and SIA evolve smoothly out of the initial topography without any sign of an initialization shock, and, as discussed in Sect.4.2, the initial ice volume change for FS is even consistent with the observed mass balance of the drainage basin.
As described previously (Sect.4.1 and 4.2), for both the FS and SIA experiments, we used the distribution of the basal friction coefficient obtained by inverting the FS problem (Fig. 3d).The idea behind this was to keep the set-up of the FS and SIA experiments as similar as possible, which facilitates comparison between the respective results.However, the price to pay for this approach is that FS is clearly favored because the obtained basal friction is custom-tailored for FS but not for SIA.It is therefore not surprising that the CTL run produces a much smaller, and more realistic, volume change in FS than in SIA (Fig. 5a).Investigating the results of the C2 and S1 runs in an experiment-minus-control way, as we did previously, eliminates this bias towards FS to a large extent.
The Cryosphere, 11, 2213-2229, 2017 The alternative would have been to invert the SIA problem for basal friction separately and run the SIA experiments with the result of this inversion.We also attempted to do so.Due to the prescribed ice geometry, temperature field and the local nature of the SIA flow field (Eq.6), this is a straightforward exercise that does not require the control inverse method or anything similarly sophisticated.However, the inversion failed to produce meaningful results.For most parts of the domain, even for no-slip conditions at the base the SIA produces surface velocities that exceed the observed ones (Fig. 3a), so that a minimization of the misfit between modeled and observed surface velocities by choosing an optimal distribution of the basal friction coefficient could not be achieved.
We have seen in Sect.4.2 that, when looking at the ice volume changes relative to the control run CTL, the results of the scenarios C2 and S1 are not overly different when FS and SIA are compared.For the surface climate scenario C2, they are almost identical; for the basal sliding scenario S1, SIA shows a ∼ 30 % larger volume change than FS.However, it must be noted that the comparability of the results is somewhat limited by the fact that we applied different boundary conditions at the grounding line (Sect.2.2.3):For FS, the hydrostatic pressure exerted by the ocean on the ice front (grounding line) provides some buttressing, and the thickness at the front can evolve, whereas for SIA the thickness is kept fixed at its present-day distribution.This different treatment results from the mathematical nature of the respective field equations and is thus unavoidable.The FS equations are a set of 3-D partial differential equations that allow prescribing stresses and/or velocities at all boundaries.By contrast, in SIA, the pressure, vertical shear stresses (all other stress components are neglected), horizontal and vertical velocities can be computed for each ice column individually, which only allows boundary conditions at the top and bottom.At the sides, the stresses and velocities are part of the SIA solution so that a horizontally applied buttressing cannot be incorporated in the formulation.Due to the absence of longitudinal stress transfer in the SIA, it would have no effect anyway.Further, we compared only the flow dynamics for grounded ice and assumed a fixed grounding line for Shi-rase Glacier.Therefore, we did not account for the potentially important impacts of grounding line migration and ice shelf buttressing on the dynamics of the system.Such effects are beyond the scope of the SIA and require at least some flavor of higher-order dynamics.

Conclusion
We compared two approaches to represent ice flow dynamics for the Shirase Drainage Basin, namely the FS formulation and the SIA, implemented within the same dynamic/thermodynamic ice flow model (Elmer/Ice).The complex nature of the stress regime in the drainage basin allows a good characterization of the differences in the evolution and The Cryosphere, 11, 2213-2229, 2017 dynamics of the area resulting from the two approaches.In the first step, we applied an inverse method to infer the distribution of the basal friction coefficient with FS.We then compared FS and SIA by assessing the respective response of the drainage basin different climatic and dynamic and forcings.
There were evident differences in the computed surface velocities between the two approaches.The surface velocities computed with FS showed a distinct, well-defined fastflowing area near the grounding line.A similar flow feature is observed in current surface velocities, essentially coincid-ing with Shirase Glacier.In contrast, the SIA produced a less well-defined contrast between the narrow fast-flowing region of Shirase Glacier and the surrounding, slower-flowing ice; the zone of fast flow is distributed over a larger area.In general, the SIA overpredicted the surface velocities everywhere in the domain, which is a consequence of the neglected longitudinal stresses and horizontal shear stress that can generate an efficient resistance to ice flow.Consequently, in transient scenarios, SIA runs consistently produced smaller ice volumes than FS runs.However, when considering ice volume evolution relative to a control run, the difference between the FS and SIA results was not overly large.Nevertheless, our findings show clearly that FS is superior to the SIA in modeling the ice flow in the area, in particular in fast-flowing regions with high slip ratios.
In this study, we considered grounded ice only and kept the grounding line fixed.A desirable extension would be to include floating ice (the small ice shelf attached to Shirase Glacier) and compare FS to coupled SIA and shallow shelf approximation dynamics within the same model.This would reveal whether the complex interactions between grounded and floating ice, including grounding line dynamics, lead to further differences in the response of the system to external forcings.

Figure 1 .
Figure 1.Shirase drainage basin: (a) location in East Antarctica.The two red lines delimit the drainage basin.They originate at Dome Fuji and follow the steepest descent of the free surface, hence approximately representing flow lines.Underlying map generated with data by Bamber et al. (2009).(b) Bed topography by Fretwell et al. (2013).(c) Surface velocities by Rignot et al. (2011), finite element mesh of the computational domain superimposed.Note the mesh refinement towards the grounding line.

Figure 2 .
Figure 2. L curve obtained with the control inverse method: cost function J o and Tikhonov regularization term J reg , parameterized by the regularization parameter λ (which determines the contribution of J reg to the total cost function, J tot = J o + λJ reg ).
Figure 3. (a) Observed InSAR-based ice surface velocities (Rignot et al., 2011), (b) surface velocities computed using FS and the control inverse method, (c) absolute error of the surface velocities |v h − v obs h | at the end of the inversion procedure and (d) computed basal friction coefficient β.

Figure 4 .
Figure 4. (a) Surface velocities computed using SIA, applying the basal friction coefficient β obtained with the control inverse method (Fig. 3d), (b) surface velocity ratio (FS / SIA), (c) scatter plot of the surface velocity ratio vs. the FS slip ratio and (d) scatter plot of the surface velocity ratio vs. the SIA slip ratio.The slip ratio is the ratio of the basal to the surface velocity for any given position.

Figure 5 .
Figure 5. (a) Ice volume (V ) simulated with FS and SIA for experiments CTL (constant-climate control run), S1 (basal sliding experiment) and C2 (surface climate experiment) and (b) ice volume relative to CTL ( V ) for experiments S1 and C2.The volumes are expressed in sea level equivalents (SLE).t = 0 corresponds to the year 2004.

Figure 8 .
Figure 8. Stress ratio ς = τ 2 xz + τ 2 yz /σ e computed at the ice base for the constant-climate control run CTL with FS.