Response to the Interactive comment on “ Basal drag of Fleming Glacier , Antarctica , Part A : sensitivity of inversion to temperature and bedrock uncertainty ”

In response to the reviewer 2's question about our choice of enhancement factor, we implemented a new sensitivity test. This was more thorough than our original test, and with a more up-to-date setup. And in fact, it reveals that our original choice was not optimal. So we added the sensitivity tests for various E values (0.5, 1.0, 2.0, 4.0) in Sect. 3.6 and Sect. 4.2, and the optimal value of 1.0 was chosen as the E in the CONTROL experiment. We redid all the simulations and modified the text and figures accordingly. Our conclusions have not changed.

Both reviewers commented on the choice of a non-zero sea level.The authors need at least a sentence to explain their choice of WGS84, rather than a geoid, as their reference surface, and to evaluate the magnitude of the problems this might cause.This is a potential problem in calculating the driving stress (which should be determined by the surface gradient with respect to the geoid) and in calculating heights above flotation (as is relevant to the companion paper).Given that the authors are using WGS-84 as their reference surface, it is not obvious that the height of the sea surface in BEDMAP-2 is the appropriate reference surface, so the authors should refer to the BEDMAP-2 documentation to see how what the ocean surface there is intended to represent: is it a mean sea-surface height or the height at some tidal phase (i.e.mean low-low tide)?
The topographic datasets in this study originate in different (geoid) reference frames.We chose WGS84 ellipsoid heights (for which convenient conversion data were available) rather than a geoid to make sure all the elevation datasets used in this study are in the same reference system.The SPOT DEM and ASTER DEM used the EGM96 geoid (Lemoine et al. 1998) as the height reference, while BEDMAP-2 used the GL04C geoid (Förste et al. 2008) as the absolute height reference.To clarify this, we added a sentence at Line 88-91 "Both SPOT and ASTER DEM products used the EGM96 geoid (Lemoine et al. 1998) as the height reference.However, the bed elevation data from Bedmap2 dataset (Fretwell et al. 2013) adopted the EIGEN-GL04C geoid (Förste et al. 2008) as its height reference, and we chose to convert all the elevation datasets to the WGS84 ellipsoid." We also modified the sentence at Line 308-310, "In addition to the ice front position, there are other sources of uncertainty in the vicinity of the ice front: ice thickness, bedrock depth, height conversion from geoid to ellipsoid, and backstress due to the presence of ice mélange.",and the sentence at Line 477-481, "These possible causes include uncertainty in local bedrock elevation (or equivalently ice thickness), uncertainty in the geoid-ellipsoid height conversion, uncertainty in observed sea level, uncertainty in exact ice front position and grounding line position, uncertainty in surface velocity, and uncertainty in potential backstress due to ice mélange and/or grounded icebergs in contact with the ice front.".
We converted all the BEDMAP-2 datasets used in this study from the GL04C geoid to the WGS84 ellipsoid based on the conversion file (gl04c_geoid_to_WGS84.tif)provided by Fretwell et al. (2013).We have clarified this at Line 108-110, "The first is from the Bedmap2 dataset (Fretwell et al. 2013) with a resolution of 1 km (hereafter bed_bm; Fig. 2b), which we converted from the EIGEN-GL04C geoid (Förste et al. 2008) to WGS84 ellipsoid heights." We have checked that both geoid-ellipsoid separation fields vary very slowly spatially compared to the surface elevation of the ice sheet, so that we do not expect any significant change in the computed surface slope that enters the driving stress calculations from mapping the geoid-based elevations into the ellipsoidal frame.Ice thickness is preserved in converting the datasets to the ellipsoid reference frame.One additional issue in using ellipsoidal heights rather than elevation relative to the geoid is that sea surface height is required for the water pressure boundary condition at the ice front, rather than being implicitly zero.This is already mentioned in Sect 3.6 in the description of the CONTROL experiment where we mention using the observed sea level height (Line 310-320) while the specification of other sea level values forms part of the sensitivity test to ice front boundary conditions (see Sect. 4.4).The sea level is obviously relevant to the heights above flotation in the companion paper -and the relevant equations there explicitly include the sea level.To clarify these matters, we added one sentence at Line 96-102, "Both geoid-ellipsoid separation fields vary very slowly spatially compared to the surface elevation of the ice sheet, so that we do not expect any significant change in the computed surface slope that enters the driving stress calculations from mapping the geoid-based elevations into the ellipsoidal frame.Ice thickness is preserved in converting the datasets to the ellipsoid reference frame (see Sect. 2.2).Clearly, the sea level height in the ellipsoidal reference frame enters the calculation of ocean water pressure on the ice front explicitly, as we discuss under experimental design in Sect 3.6 and Sect.4.4.".Also, while we found no relevant information characterising the sea surface height in the BEDMAP2 documentation, we remind the Editor that the basis of the BEDMAP2 reference frame is the EIGEN-GL0C4 geoid (Förste et al. 2008).As one would expect from the standard definition of the geoid, construction of EIGEN-GL0C4 involves mean sea surface height -to quote their Conclusions: A new satellite-only gravity field model, EIGEN-GL04S1, complete up to degree and order 150 has been inferred from GRACE and LAGEOS data.Using this satelliteonly model as a starting point, a new combined global gravity fieldmodel EIGEN-GL04C, complete up to degree and order 360, has been developed, incorporating surface gravity data including newly available or improved data sets of the Arctic, Antarctica and North-America and improved mean sea surface heights from altimetry processing at GFZ.
At least one referee commented on the linear sliding law.The authors should make clear in the text that the sliding law is only used as a numerically convenient tool for calculating the basal shear stress.
Thanks for the editor's suggestion.Reviewer 1 said: Line 163: What is your justification for using a linear sliding law?
We replied: Different sliding laws in inverse modeling will not change the inversed basal shear stress distribution, and it will just lead to different basal friction coefficients based on different sliding law.In diagnostic studies that invert to find the basal shear stress which gives the best agreement with observed surface velocities, the choice of sliding "law" is not relevant provided that the required stress can be generated by adjustments of the parameters in the sliding law -in this case the coefficient C. The inversion procedure modifies C to modify stress -adjusting the momentum balance.That solution of the Stokes equation provides an updated estimate of basal velocity -which enters the next cycle of the inversion search.The question does remain whether this is physically suitable relationship to apply when the system is evolving, but this is not relevant here.So we adopted the simplest sliding law here following Gillet-Chaulet et al. (2012), Gagliardini et al. (2013).We clarified this in the text.
To clarify this, we altered one sentence at Line 202-208, "Considering that in this diagnostic study the sliding law is only used as a numerically convenient tool for calculating the basal shear stress, a simple linear sliding law following Gillet-Chaulet et al. (2012), Gagliardini et al. (2013) is applied on the bottom surface: where C, the basal friction coefficient, is used as the adjustable parameter in the inversion scheme described below." We would be happy to provide more detail in the main text if the editor feels this is needed.
I recommend that the authors convert figure 5 to include only two versions of the temperature initialization (the one selected, and a second version, of their choice, that demonstrates unacceptable behavior), and only the first and last steps of the temperature iteration.The intermediate steps and the full range of models tested can be described, very briefly, in the text.
We can certainly comply with this recommendation if the editor insists, but we feel it undermines the presentation of the multiple cycle spin-up procedure.We think demonstrating the extent to which multi-cycle spin-up reduces dependence on assumed initial temperatures will be of broad interest to modellers, beyond our other location-specific sensitivity studies.If we must modify Fig. 5 as the editor suggests, we would want to include the original figure and relevant discussion in the supplementary material but we would prefer to have it more directly citable in the main paper.We could show the most disparate cases -simulated basal friction coefficient distributions for CONTROL and TEMP1 (T = -20 C) -in the main text, and the results of other two experiments TEMP2 and TEMP3 in the supplementary material but this fragmentation of the examples seems clumsy and the loss of discussion of the progression through the cycles regrettable.
I recommend that the authors remove the right-hand column of figure 9 (the relaxation height difference).This column is given only cursory attention in the text.
Thanks for the editor's suggestion.We have removed the right column of Fig. 8 (the original Fig. 9) and modified the relevant text.
For all figures (5-9) with multiple panels, the rows and the columns of the figures should be labeled, so that the reader does not need to refer to the letter key in the captions to interpret each panel.
Thanks for the editor's suggestion.We have labeled the rows and columns of those figures (Figs. 5-8, and figures in the supplementary).
These changes are somewhat substantial, and I would like to look over the authors' response to them before I consider whether it needs a re-review or not.

Best
Ben Smith

Introduction
In response to rapid changes in both atmosphere and ocean, glaciers in West Antarctica (WA) and the Antarctic Peninsula (AP) have undergone rapid dynamic thinning and increased ice discharge over recent decades, which has led to a significant contribution to global sea level rise (Cook et al., 2016;Gardner et al., 2018;Wouters et al., 2015).Understanding the underlying processes, especially for fast-flowing outlet glaciers, is crucial to improve modeling of ice dynamics and enable reliable predictions of contributions to sea level change, especially for fast-flowing outlet glaciers.
The high velocities of the fast-flowing outlet glaciers are determined byarise from both internal ice deformation and or ice sliding at the bed or both.Internal dDeformation is strongly dependent on gravitational driving stress, englacial temperature, the development of anisotropic structure at the grain scale in polycrystalline ice (e.g.Gagliardini et al. (2009)) and larger scale weakening from fractures (Borstad et al., 2013).Basal sliding is strongly dependent on the gravitational driving stress, bedrock topography and the basal slipperiness, which in turn is affected by the roughness of the bed, the presence of deformable till, or the basal subglacial hydrology.Therefore, one of the keys to modeling fast-flowing glaciers is accurate knowledge of the basal conditions: the bedrock topography and the basal slipperiness (Gillet-Chaulet et al., 2016;Schäfer et al., 2012).Inverse methods are commonly used in ice sheet models to infer the basal shear stressfriction coefficient, and basal velocities, and ice rheology from the glacier topography geometry and observed surface velocities (Gillet-Chaulet et al., 2016;Gladstone et al., 2014;Morlighem et al., 2010).
However, pPoorly constrained quantities, like the basal topography, and the distribution of internal temperature, have provided major challenges for modeling the basal shear stress (Vaughan and Arthern, 2007), especially for small-scale glaciers.However, iIn a study studies carried out on a fast-flowing outlet glacier draining from the Vestfonna ice cap in the Arctic (Schäfer et al., 2014;Schäfer et al., 2012), it was found that the Robin inverse methods did not depend strongly on the mesh resolution or uncertainties in the topographic and velocity data.In their case, sliding dominated the flow regime, and tThe impact of ice temperatureinternal deformation on ice internal deformationice velocity was relatively small compared to the important role of friction heating at the bed on the basal sliding (Schäfer et al., 2014;Schäfer et al., 2012).However, It is unclear whether this property is specific to the Vestfonna situation or if it also applies to other fast flowing glaciersno generalization on these findings to Antarctic outlet glaciers has been investigated.The motivation of this paper is twofold: to test the sensitivity of a variational inversion inverse method (MacAyeal, 1993;Morlighem et al., 2010)s for basal friction to basal geometry and to an assumed initial englacial temperature distribution for a different outlet glacier system, and to determine a robust basal shear stressfriction coefficient pattern for the Fleming Glacier, located in the southern AP, in 2008.
The Wordie Ice Shelf (WIS) (Fig. 1b) in the southern AP has experienced ongoing retreat and collapse since 1966, with its almost-complete disappearance by 2008 (Cook and Vaughan, 2010;Zhao et al., 2017).The Fleming Glacier (FGLFG) (Fig. 1b), as the main tributary glacier that fed the WIS, has a current length of ~80 km and is ~10 km wide near the ice front (Friedl et al., 2018).This glacier has recently shown a rapid increase in surface-lowering rates (doubling near the ice front after 2008) (Zhao et al., 2017), and the largest velocity changes (> 500 m yr -1 near the ice front) across the whole Antarctic ice sheet over 2008-2015 (Walker and Gardner, 2017).
In this study, we employed the Elmer/Ice code (Gagliardini et al., 2013), a three-dimensional (3D), finite element, full-Stokes ice sheet model, to invert for the basal dragfriction coefficient distributions over the whole WIS-FGLFG system in using a parallel computing environment.
Here, wWe deduce the distribution of basal shear stress using a control inverse method to assess its sensitivity to bedrock topographies, assumptions about the initial temperature distribution, bedrock topographies, ocean boundary conditions and other constraint parameters in the model.We introduce the data in Sect.2, present the ice sheet model, spinup scheme and experiment design in Sect.3, and discuss the results in Sect. 4 before we give the conclusions in Sect. 5.

Surface elevation data in 2008
The surface topography in 2008 (Fig. 2a) is combined from two SPOT DEM products acquired on 21 st Feb, 2007 (resolution: 240 m) and 10 th Jan, 2008 (resolution: 40 m) (Korona et al., 2009) and an ASTER DEM product ranging from 2000 to 2009 (resolution: 100m) (Cook et al., 2012).The surface elevation data for the Fleming Glacier is mainly from the SPOT DEM product acquired on 10 th Jan, 2008 (see masks of different DEM products in Fig. S1 in the supplementary material).Here we apply the SPOT DEM precision quality masks on the raw data to extract the DEM data with correlation scores from 20% to 100%.Areas with low correlation scores were filled with the ASTER DEM data.To remove noise from the DEM data, the combined DEM (resolution: 40 m) is resampled to 400 m with a median filter and a window size of 10×10 pixels.
Both SPOT and ASTER DEM products used the EGM96 geoid (Lemoine et al., 1998) as the height reference.However, the bed elevation data from Bedmap2 dataset (Fretwell et al., 2013) adopted the EIGEN-GL04C geoid (Förste et al., 2008) as its height reference, and we chose to convert all the elevation datasets to the WGS84 ellipsoid.The EGM96 geoid (Lemoine et al., 1998) and EIGEN-GL04C geoid (Förste et al., 2008) is are used to convert from the EGM96 Geoid geoid and EIGEN-GL04C geoid values to WGS84 ellipsoidal heights, respectively.We extract a median value of 15 m for the DEM data over Marguerite Bay (Fig. 1a) as the mean local sea level in the ellipsoid frame.
Both geoid-ellipsoid separation fields vary very slowly spatially compared to the surface elevation of the ice sheet, so that we do not expect any significant change in the computed surface slope that enters the driving stress calculations from mapping the geoid-based elevations into the ellipsoidal frame.Ice thickness is preserved in converting the datasets to the ellipsoid reference frame (see Sect. 2.2).Clearly, the sea level height in the ellipsoidal reference frame enters the calculation of ocean water pressure on the ice front explicitly, as we discuss under experimental design in Sect 3.6 and Sect.4.4.

Bed elevation data
The bed topography plays an very important role in the basal sliding and distribution of fastflowing ice (De Rydt et al., 2013).However, high-resolution observations of bedrock elevation for this regionthe WIS-FG system are still not available.To explore the sensitivity of the basal friction coefficientshear stress distribution to the uncertainty in the bedrock topography, we adopt three basal topographies.The first is from the Bedmap2 dataset (Fretwell et al., 2013) with a resolution of 1 km (hereafter bed_bm; Fig. 2b), which we converted from the EIGEN-GL04C geoid (Förste et al., 2008) to WGS84 ellipsoid heights.
The other two are derived using the equations below: where S2008 is the 2008 surface DEM in 2008 mentioneddescribed in Sec.2.1, and Sbm is the surface elevation data from Bedmap2 (Fretwell et al., 2013), again relative to the WGS84 ellipsoid.Sbm is downscaled to 500 m with a bilinear interpolation method., and Hmc (where "mc" refers to "mass conservation") is the ice thickness data with a resolution of 450 m combined fromcovering three sources regions shown in Figs.2e: .Hmcdata for the yellow area is computed from the Center for Remote Sensing of Ice Sheets (CReSIS) ice thickness measurements using the Ice Sheet System Model's a mass conservation method (Morlighem et al., 2011;Morlighem et al., 2013), based on ice thickness measurements from ; the Center for Remote Sensing of Ice Sheets (CReSIS), using ice surface velocities in 2008 from Rignot et al. (2011b), surface accumulation from RACMO 2.3 ( van Wessem et al., 2016) and 2002-2008 ice thinning rates from Zhao et al. (2017).The thickness data for the grey area is interpolated from Bedmap2 (Fretwell et al., 2013),; while the data in the red area thickness ensures a smooth transition between the two regionsis interpolated from CReSIS and Bedmap2 data.The yellow area is indicates the Fleming Glacier system with ice velocity >100 m yr -1 .The uncertainty of Hmc (Fig. 2f) ranges from 10 m to 108 m.For the calculation of Hmc, we assume the that the ice elevation changes over 2002 to 2008 (Zhao et al., 2017) were small compared to the uncertainty uncertainties in ice thickness (Fig. 2f) and could be ignored in the ice thickness measurements which span a wider time frame.Both bed_mc (Fig. 2c) and bed_zc (Fig. 2d) have a higher resolution of 450 m while bed_bm (Fig. 2b) has a resolution of 1 km.The uncertainty of bed_bm for the fast-flowing regions of the Fleming Glacier (yellow and red area in Fig. 2e) ranges from 151 m to 322 m (Fretwell et al., 2013), while the uncertainty of bed_mc and bed_zc ranges from 10 m to 108 m (from uncertainties in Hmc).
The bed topography data (Fig. 2b) indicates the essentially marine character of the Fleming Glacier, showing two basins featuring retrograde slopes, both located underneath the main trunk of the Fleming Glacier's fast flow region.The region basin further upstream (hereafter "FG upstream basin") has a steeper retrograde slope than the one closer to the grounding line of those basins (hereafter "FG downstream basin").For the FG downstream basin, elevation differences between bed_bm and the other two datasets (Figs 2c,2d) show that bed_bm has a generally steeper retrograde slope.The sensitivity of basal shear stressfriction coefficient distributions to the three bed datasets is discussed in Sect.4.2.

Surface velocity data in 2008
The surface velocity data used for 2008 (Fig. 1b) were obtained from MEaSUREs InSARbased Antarctic ice velocity (from the fall 2007 and/or 2008) produced by Rignot et al. (2011b) (version 1.0) with a resolution of 900 m and with uncertainties ranging from 4 m yr -1 to 8 m yr -1 over the study area.For the regions without data (grey area in Fig. 1b), we prescribe the surface speed to be 0. We do not use the finer (450 m) resolution MEaSUREs velocity here since the coarser (900 m) resolution data have been subjected to some postprocessing, including smoothing and error corrections.

Method
All the simulations are carried out using the Elmer/Ice model (Gagliardini et al., 2013).These simulations are used to solve the ice momentum balance equations with an control inverse method to determine the basal dragfriction coefficients, and the steady state heat equation for to model the ice temperature distribution.The ice rheology is given by Glen's flow relation (Glen, 1955): with where  is the deviatoric stress and ̇ is the strain rate tensor.The viscosity  is computed as: using where E is an overall flow enhancement factor, A is a temperature-dependent rate factor calculated using an E, and a function of the ice temperature relative to the pressure melting point according to the Arrhenius equationLaw (Gagliardini et al., 2013), .  ̇= √(̇2) 2 ⁄ is the effective strain rate, and n is the exponent in Glen's flow law.Table 1 lists the parameters used in this study.

Mesh generation and refinement
We used GMSH (Geuzaine and Remacle, 2009) to generate an initialthe 2-D horizontal footprint mesh with the boundary defined from the grounding line data in 1996 (Rignot et al., 2011a) and the catchment boundary of the feeding glacier system (Cook et al., 2014), with the assumption that the ice front position in 2008 was coincidednt with the grounding line position in 1996 (Rignot et al., 2011a).This assumption is tested as part of the sensitivity tests to various ice front positions.
To reduce the computational cost without reducing the accuracy, we refined the mesh withusing the anisotropic mesh adaptation software YAMS (Frey and Alauzet, 2005) using the local Hessian matrix (second order derivatives) of the surface velocity data in 2008 from Rignot et al. (2011c) as a metric for the mesh density.The resulting mesh is shown in Fig. 3 with and has the minimum and maximum element sizes of approximately 250 m and 4 km, respectively.The 2-D mesh is then vertically extruded using 10 equally spaced, terrain following layers.Sensitivity tests have been done on the Vestfonna ice cap (Schäfer et al., 2014;Schäfer et al., 2012) to prove demonstrate the robustness of inverse simulations to the vertical mesh resolution.In the current study an experiment with 20 extruded layers (not shown) gives very similar results as with 10 layers, confirming those findings also apply to It would be useful to know whether the WIS-FGLFG system shows same robustness to the vertical resolution, but this is beyond the scope of current study.Experiments with various horizontal resolutions (1 km, 500 m, 250 m, and 125 m) show that 250 m are sufficient for the simulations on the WIS-FG system.

Boundary Conditions
For transient simulations (surface relaxation, section 3.3), the stress-free upper surface is allowed to evolve freely, with a minimum imposed ice thickness of 10 m over otherwise icefree terrain.For inverse and temperature simulations, the upper surface height and temperature are fixed.
The surface temperature is defined by the yearly averaged surface temperature over 1979-2014 computed from the regional atmospheric climate model RACMO2.3/ANT27(van Wessem et al., 2014).The geothermal heat flux (GHF) at the bed is obtained from Fox Maule et al. ( 2005) using input data from the SeaRISE project, and the GHF is interpolated with bilinear interpolation method from the standard 5 km grid onto the anisotropic mesh.A basal heat flux boundary condition combining GHF and basal friction heating is imposed for temperature simulations.
At the ice front, the normal component of the stress where the ice is below sea level is equal to the hydrostatic water pressure exerted by the ocean.The uncertainties of ice thickness and bedrock topography, the low accuracy of ice front and grounding line locations, and the possible buttressing on the ice front by partly detached icebergs and ice mélange (see Fig. 1c) would affect the calculation of ocean forcing there.Accordingly, wWe will discuss the sensitivity to the ice front boundary condition in Sect.4.34.On the lateral boundary, which falls within glaciated regions, the normal component of the stress vector is set equal to the ice pressure exerted by the neighboring glacier ice while the tangential velocity is assumed to be zero.
The bedrock is regarded as rigid, impenetrable, and temporarily temporally fixed in all simulations.The present-day solid Earth deformation rate in the Fleming glacier region (Zhao et al., 2017) is negligible compared to the uncertainty of the bedrock data.Assuming that basal melt is negligible under grounded iceSo, the normal basal velocity is assumed setto be zero at the ice/bed interfacehere.The sliding relation relates the basal sliding velocity ub to basal shear stress   .Considering that in this diagnostic study the sliding law is only used as a numerically convenient tool for calculating the basal shear stressHere, a simple linear sliding law following Gagliardini et al. (2013); Gillet-Chaulet et al. ( 2012) is applied on the bottom surface: where C, is athe basal sliding friction coefficient, is used as the adjustable parameter in the inversion scheme described below.During the initial surface relaxation, and at the start of the inversion, C is initialized to a constant value of 10 -4 MPa m -1 yr (following Gillet-Chaulet et al. ( 2012)), which is replaced with the inverted C in subsequent following steps.The surface temperature is defined by the yearly averaged surface temperature over 1979-2014 computed from the regional atmospheric climate model RACMO2.3/ANT27( van Wessem et al., 2014).The geothermal heat flux (GHF) at the bed is obtained from Fox Maule et al. ( 2005) using input data from SeaRISE project, and the GHF is interpolated with bilinear interpolation method from the standard 5 km grid onto the anisotropic mesh.A basal heat flux boundary condition combining GHF and basal friction heating is imposed for temperature simulations.

Surface relaxation
There may be non-physical spikes in the initial surface geometry, caused for example by observational uncertainties of the surface or bedrock data and/or by the resolution discrepancy between mesh and geometry data.To reduce these features, we relaxed the free surface of this domain during a short transient simulation of 0.2 yr length with a timestep of 0.01 yr.This is long enough to remove the non-physical spikes, but too short to significantly modify the geometry of the fast flowing regions of the Fleming Glacier.

Inversion for basal shear stressfriction coefficient
After the surface relaxation, we used a variationalthe control inverse method (MacAyeal, 1993;Morlighem et al., 2010) implemented in Elmer/Ice (Gagliardini et al., 2013;Gillet-Chaulet et al., 2012) to constrain the basal sliding friction coefficient C in Eq. ( 35).To avoid non-physical negative values, we used a logarithmic representation of the basal dragfriction coefficient, C  = 10  , where  can take any real value.
The inverse method is based on adjusting the spatial distribution of the basal dragfriction coefficient to minimize a cost function that represents the mismatch between the magnitudes of the simulated and observed surface velocities.: where Γ  is the upper surface of the domain, u and   are the simulated and observed surface velocities, respectively.We do not try to fit velocity directions.
To avoid over-fitting of the inversion solution to non-physical noise in the observations, a regularization term is added to the cost function as: where J0 represents the square of the magnitude of the mismatch between the simulated and observed surface velocities, Jreg is the a regularization term imposing a cost on spatial variations in the control parameter ,  is a positive regularization weighting parameter, and Jtot is the total cost (following for example Gillet-Chaulet et al. ( 2012)).Thus, the minimum of this cost function is no longer the best fit to observation but a compromise between fit to observations and smoothness in .An L-curve analysis (Hansen, 2001) has been carried out for inversions in the current study to find the optimal  by plotting the term Jreg as the function of J0 (Fig. S2 in the supplementary material).The optimal value of 10 8 is chosen for  to minimize J0.
With  = 10 8 , we compute the total cost Jtot with different values of flow enhancement factor E (0.7, 1, 2.5, 5, 10).It is found that inversions with smaller E gave a better-simulated surface velocity for slow ice-flow regions while greater E gave a better velocity for fast ice-flow regions.The optimal value of E = 2.5 is chosen for the current study.

Steady-state temperature simulations
In the absence of a known englacial temperature distribution for the Fleming Glacier system, the steady state ice temperature is solved heat transfer equation is solved using an iterative method as described in Gagliardini et al. (2013) to provide temperatures for use in the inversion process.The ice velocity and geometry are held constant for this part of the simulation.Steady-state temperature simulations for a non-steady-state glacier system will result in the estimations of the temperatures that deviate from actualityreality.However, sSimilar experiments on the Greenland Ice Sheet indicated that the simulated steady-state temperature field could present provide a reasonable thermal regime for calculation of basal conditions (Seroussi et al., 2013).Heat sources and internal energy transfer determine the temperature distribution within the ice.The heat transfer equation is solved using an iterative method as described in Gagliardini et al. (2013).

Experiment design
Gong et al (2017) adopted a four-step spin-up scheme (Gladstone et al., 2014) in inverse modelling using Elmer/Ice (Gagliardini et al., 2013), without testing the effect of assumptions about the initial englacial temperature distribution on the inversion resultsThe four-step spinup scheme (Gladstone et al., 2014) has been adopted in inverse modeling using Elmer/Ice (Gong et al., 2017), without testing the effect of initial temperature assumption on the inversion results.To explore the sensitivity of inverse modeling to initial temperature assumptions, we proposed a spin-up scheme with more cycles (three cycles in this study as presented in Fig. 4).For each cycle, we followed the spin-up scheme of from Gladstone et al. (2014): 1. surface relaxation; 2. inversion foof ther basal friction coefficient using the relaxed surface geometry; 3. a steady state temperature simulation using the simulated velocity velocities from that inversion; 4. another inversion with the previously obtainedsimulated steady-state temperature.
The surface relaxation for each cycle starts from the same initial geometry described in Sect.3.3.For cycle 1, the surface relaxation and first inversion are implemented with an initial temperature assumption (described below) and uniform basal dragfriction coefficient of 10 −4 MPa m −1 a (following Gillet-Chaulet et al. ( 2012)).For cycles 2 and 3, the surface relaxation and inversion are implemented initiated with the simulated steady-state temperature and an initial distribution of basal dragfriction coefficient C from the final state of the previous cycle.
To explore the sensitivity of our inverse methods to assumed initial englacial temperature distribution, enhancement factor (E), basal topography, ice front positions, and the ice front boundary conditions, we carried carry out the experiments summarized in Table 2.
An assumed iInitial englacial temperature distribution is used in the first cycle of the scheme above and would affect the surface relaxation, would affect the modelled ice deformation and the ice velocity field, especially for fast-flowing regions, and consequently lead to a difference in the relaxed upper surfaceaffect the steady-state temperature calculation, which might affect the subsequent inversion process.To explore the impact of initial temperatures on inversion results with the three-cycle spin-up scheme, we proposed experiments with different initial temperature assumptions for the surface relaxation and initial inversion in Cycle 1. TEMP1: a uniform temperature of -20 ℃; TEMP2: a uniform temperature of -5 ℃; CONTROL: a linearly increasing temperature from the upper surface values (see also Sect. 3.2) to the pressure melting temperature at the bed.To test the sensitivity of basal dragfriction to the relaxed geometry, we also added another experiment experiment -"TEMP3": surface relaxation in the first cycle using the linear temperature, followed by inversion with a uniform temperature of -20 ℃.Experiments TEMP1, TEMP2 and TEMP3 differ from CONTROL only in the temperature fields imposed before the first temperature simulation.Ma et al. (2010) tested the influence of ice anisotropy on the ice flow through various enhancement factors, and found that ideal appropriate E-values for the grounded ice are usually >1.0.To find out the most appropriate value of E (in Eq. 4) in this study, we evaluate inversion carried out with different values of E (EF1: E = 0.5, CONTROL: E = 1.0,EF2: E = 2.0, EF3: E = 4.0; Table 2).Experiments EF1, EF2 and EF3 differ from CONTROL only in terms of the value used for E.
As described in Sec.2.2, we generated three different bed topography datasets to explore the sensitivity of the inverse modelling.The three-cycle spin-up scheme is carried out for each bed dataset using the linear (described above) initial temperature distribution described above.These experiments are referred to as CONTROL, BEDZC, and BEDMC (Table 2).
Experiments BEDZC and BEDMC differ from CONTROL only in terms of the bedrock data set used.
In our standard model domain we assume the 2008 ice front is coincident with the 1996 grounding line, which has an error of several km on fast-moving ice (Rignot et al., 2011a) and might have changed since 1996.The frontal surface elevation is from the SPOT DEM data in Jan 2008, which shows the ice front position is ~1.5 km downstream of the 1996 grounding line position.Since such a narrow residual ice shelf is considered unlikely to have a major influence, we construct the model geometry to have the ice front coincide with the 1996 grounding line for simplicity, i.e. all ice is considered grounded.To test the sensitivity of inverse modelling to the ice front positions, we implement two further scenarios with different ice front positions: downstream (experiment IFP1) and upstream (experiment IFP2) of the 1996 grounding line position (CONTROL).In IFP1, we assume the ice front position is coincident with the frontal boundary of SPOT DEM data (~1.5 km downstream).In IFP2, we artificially put the ice front position ~1.5 km upstream of the 1996 grounding line position for ~1.5 km.IFP1 and IFP2 differ from CONTROL only in their ice front position.
In addition to the ice front position, there are other sources of uncertainty in the vicinity of the ice front: ice thickness, bedrock depth, height conversion from geoid to ellipsoid, and backstress due to the presence of ice mélange.These uncertaintiesUncertainties from the ice thickness and bedrock datasets may have an significant effect on the pressure boundary condition applied to the ice front, which conventionally balances the normal stress in the ice against the ocean water pressure.In view of the ice thickness uncertainty (ranging from 10 m to 100 m) and hence bedrock depth around the grounding line, and given the possibility of increased additional pressure buttressing force due to floating icebergs and ice mélange as indicated in many previous studies (e.g.Amundson et al. (2010); Krug et al. (2015); Robel (2017); Todd and Christoffersen (2014); Walter et al. (2017)) and clearly seen in Fig. 1c, we vary the ocean pressure this boundary condition by varying the sea level used to calculate ocean water pressure.This approach directly represents some small uncertainty in the actual exact sea level itself, but is also a proxy for pressure errors variations due to bedrock elevation/ice thickness uncertainty and mélange back stress.Firstly, in the CONTROL experiment, we assume an ocean pressure at the ice front computed using the observed sea level of 15 m, as mentioned in Sec.2.1.We adjust the sea level by 10 m from hydraulic equilibrium to test the sensitivity of the inverse modeling to the ice front boundary condition.Firstly, we assume an ocean pressure at the ice front computed using the sea level mentioned in Sec.2.1.We further simulate two alternative scenarios for the sea level used in the simulations to calculate ocean pressure: IFBC1 with a sea level of 5 m and IFBC2 with a sea level of 25 m.Another extreme scenario (IFBC3, Table 2) is adopted here by setting the ice front pressure to the ice overburden: where   () is the pressure at the ice front as a function of height ,   is ice density (Table 1),  is the gravitational constant (Table 1), and   is the height of ice upper surface at the ice front.This is the pressure that would be imposed by a hypothetical undeforming continuation of the advancingneighboring glacier, and imposes zero normal strain rate at the ice front.The ice surface elevation   at the front is ~115 m, approximately 100 m above actual sea level.The total vertically integrated pressure imposed by this condition is equivalent to a sea level of ~60 m, although the vertical distribution of pressure is different todiffers from an ocean pressure condition.Experiments IFBC1, IFBC2 and IFBC3 differ from CONTROL only in their ice front boundary condition.
In our model domain we assume the 2008 grounding line is consistent with the 1996 grounding line, which has an error of several km on fast-moving ice (Rignot et al., 2011a) and might have changed since 1996.The frontal surface elevation is from the SPOT DEM data in Jan 2008, which shows the ice front position is ~1.5 km downstream of the 1996 grounding line position.Since such a narrow residual ice shelf was considered unlikely to have a major influence we constructed the model geometry to have the ice front coincide with the 1996 grounding line for simplicity, i.e. all ice is considered grounded.

Results and discussions
The main focus of the current study is the sensitivity of the inversion to the variations of three five factors: temperature initialization, enhancement factor, bed topography, ice front positions, and ice front stress balanceoceanic pressure boundary condition.The evaluation criteria are the robustness of simulated basal dragfriction coefficient distribution to experiment design and the mismatch between the simulated and observed surface velocities.

Sensitivity to initial temperature
We present the results for the inferred basal friction coefficients of from the CONTROL and four three TEMP experiments (Sect.3.76, Table 2) for the WIS-FGLFG system in Fig. 5.The 2008 ice velocity contours are added as visual references for comparing the basal dragfriction coefficient patterns in the regions of fast flow, since the largest observed ice velocity changes occurred in fast flowing outlet flow regions (Mouginot et al., 2014;Walker and Gardner, 2017).
All those experiments showed thatIn each cycle, the absolute differenceroot-mean-square deviation (RMSD, sometimes also called root-mean-square error) between the relaxed and the observed surface was < 30 25 m (see Table S1 in the supplementary material), smaller than the ice thickness uncertainty (> 50 m) used in this study.However, the systematic changes generated at the ice front during the surface relaxation may have an effect on the inversion, and this is further discussed in Sect.4.4.However, we think some of the systematic changes generated by surface relaxation are not correcting real errors in the surface topography data, as discussed later in Sect.4.3.After the first cycle (left column, Fig. 5), results showed different patterns of basal dragfriction coefficient for each experiment, especially in the fast-flowing regions with surface velocity higher thanexceeding 1000 m /yr -1 (yellow contour in Fig. 5).The basal dragfriction coefficients from TEMP2 (Fig. 5d5g) and CONTROL (Fig. 5g5a) share a very similar rib-like patternsticky spots around the ice front, and some isolated sticky spots ~3-5 km upstream of the ice frontand another rib close to the yellow contour in the fast flow regions (> 1000 m yr - 1 ), but the TEMP1 (Fig. 5a5d) and TEMP3 (Fig. 5j) display different patterns, indicating dependence on the initial temperature assumption.The RMSDs of key properties are computed to evaluate the consistency of these experiments (Table S2-S5).This is in contrast to a similar inverse study on the Vestfonna ice cap (Schäfer et al., 2012), which showed little impact of temperature distribution on the basal sliding coefficient.That was due to a low contribution of ice deformation to ice motion compared to the basal sliding (Schäfer et al., 2012).We return to this contrast after considering the effect of the second and third cycles of our spin-up.
To remove reduce the dependence on initial temperature and achieve a consistent equilibrium thermal regime with respect to the given slip friction coefficient distribution for surface relaxation, we carried out the second cycle shown in Fig. 4. The basal dragfriction coefficients from the final step of Cycle 2 (the middle column in Fig. 5) shows greater similarity across all the temperature experiments.However, for experiments CONTROL and TEMP2, the isolated sticky points ~3-5 km upstream of the grounding lineice front (with horizontal scale around ~1 km and peak basal friction coefficient of around 6×10 -5 MPa m -1 yr) in the downstream basin show a trend ofmostly decreasing decrease and or disappearing from the first cycle (left column of Fig 5Figs.5a, 5g) to the second cycle (middle column of Figs.5b, 5h) for experiments CONTROL" and TEMP2.Therefore, a third cycle was implemented for all temperature assumptionsto test whether a two-cycle spin-up scheme was enough to reduce the dependence on the initial temperature assumptions.After the third cycle, all the scenarios depicted a similar basal dragfriction coefficient pattern (right column in Fig. 5).These differences in basal friction coefficients between the TEMP simulations can also be analyzed through Table S2 and Fig. S4.While these statistics and visualizations confirm the similarity between CONTROL, TEMP2 and TEMP3, it is evident that TEMP1 still shows notable differences to these simulations, even after three cycles (see also Table S3 for basal velocity RMSD).The CONTROL simulation, starting with a linear interpolation of temperature from upper to lower surfaces, seems to be the best option for several reasons: the choice of temperature value for upper and lower surfaces is physically motivated, which is not true for the other assumptions; it shows the lowest RMSD between simulated and observed upper surface velocity of the temperature sensitivity simulations (Table S5); and it shows the least change in the temperature distribution over the three cycles (Table S4).Given this choice of preferred temperature initialization (CONTROL), and the significant difference between this and the cold initialization (TEMP1), we argue that TEMP1 likely deviates furthest from an ideal temperature initialization, and that such a large initial deviation would require more than three cycles to converge on a basal friction coefficient distribution.In other words, we postulate that the three cycles are likely sufficient to provide a robust inversion only for initial temperatures moderately close to reality, with the linear interpolation in the vertical providing the most appropriate initial guess amongst our tests.HenceThe differences between the simulated and observed surface speed for the above experiments (Fig. 6) also prove that the three-cycle scheme could provide relatively robust inversion results with little sensitivity to the initial temperature.Considering the linear temperature is likely closer to a realistic temperature distribution, we adopted the scenario with initial linear temperature for the experiments described hereafter.
The present study is focused on exploring the effects of uncertainties and their control, and while the dynamics of the FGLFG system will be discussed in more detail in a companion paper (Zhao et al., companion paper).However, a few comments are in order regarding the contrast with the previousan earlier study on the Vestfonna ice cap.The low impact of temperature distribution profile on the basal sliding friction coefficient distribution in that study was due to a lower contribution of ice deformational to ice motion compared to the basal sliding (Schäfer et al., 2012).Internal ice deformation, and hence temperature, may be especially important for the WIS-FGLFG system due to steep surface slopes and corresponding high driving stresses in the region between the downstream and upstream basins (~8-12 km upstream of the ice front in Fig. 7aS5a).The patterns of basal dragfriction coefficient (right column of Fig. 5) all indicate substantial differences spatial variation in basal dragfriction over the fast flowing part of the FGLFG.For example, in the region flowing at overfaster than 1000 m yr -1 (inside the yellow contour), we see very low dragfriction over the downstream basin, but higher dragfriction coefficients over the upstream bedrock high, and in a narrow band along the ice front.The nature of the basal shear stress is further complicated by substantial variations in the contribution of basal sliding velocity and of vertical shear deformation to the flow.A comparison between the simulated basal and surface velocities (Fig. 7bS5b) shows that internal vertical deformation shear dominates the ice dynamics in the some of the fast-flowing regionsregion of high slope between the downstream and upstream basins, where the driving stress is relatively high.This alone would suggest a high sensitivity of modelled sliding velocity and basal dragfriction to the englacial temperature.
The threemulti-cycle iterative spin-up scheme is suggested as an effective set-up for inverse modelling of fast-flowing glaciers that have high surface slopes and vertical shear strain rates and therefore are sensitive to the internal vertical ice temperature fielddistribution.In the present application to the Fleming system, three cycles were sufficient, except in the case of an unphysically cold initialization.In other cases, the inversion process is not so heavily dependent on the temperature field, for example for reproducing the shear margins of the outlet glacier of Basin 3 on Austfonna ice cap, Svalbard (Gladstone et al., 2014).

Sensitivity to enhancement factor
Sensitivity of inverse modelling to the flow enhancement factor has been explored by experiments EF1-3 and the results (after three-cycle procedure) are shown in Fig. 6.The simulated basal friction coefficients (left column in Fig. 6) show different patterns with different E values.Recall that from Eq. ( 4), smaller E means higher ice viscosity.The local high friction coefficient sticky spots near the ice front expanded both upstream and along the ice front with increased E values, forming a band across the ice front for E = 4.0 (EF3).Conversely, inversions with smaller E give a better-simulated surface velocity at the ice front (middle column in Fig. 6), and also lead to smaller differences between the observed and relaxed surface elevation after the surface relaxation (right column in Fig. 6), whereas for EF3 the surface relaxation generates a considerable steepening of the surface slope towards the ice front (Fig. 6l).However, the computed RMSD of the surface velocity mismatch for the fast flowing regions (> 1500 m yr -1 , middle column in Fig. 6 and Table S5) indicates that the experiment EF1 (E = 0.5) (Fig. 6e) shows greater underestimation of surface velocity than CONTROL (Fig. 6f).Therefore, the optimal value of E = 1.0 is chosen as the most suitable enhancement factor for the Fleming system.

3 Sensitivity to bedrock topography
Figure 8 S6 summarizes results from the three experiments using different bed topographies (Sect.3.7, Table 2).The 2008 ice velocity contours are added as visual references for comparing the basal drag coefficient patterns in the regions of fast flow, since the largest observed ice velocity changes occurred in fast outlet flow regions (Mouginot et al., 2014;Walker and Gardner, 2017).As shown in Fig. S68, the simulated basal friction friction coefficient C varies slightly with bedrock geometry and its distribution shows greater similarity between BEDZC and BEDMC, compared with CONTROL.CONTROL (using "bed_bm"Bedmap2 bedrock data; Fig. 8aS6a) shows slightly smaller basal dragfriction coefficients than BEDMC (Fig. 8bS6b) and BEDZC (Fig. S8c6c) in the fast-flowing region (>1500 m yr -1 , cyan contour in Fig. 8S6).and tThe pattern in the region between the yellow 1000 and cyan 1500 m yr -1 contours also differs compared toin the CONTROL case, which might be caused by the deeper bedrock of bed_bmBedmap2 in the FG downstream basinthis region (Fig. 8gS6g), compared to the other two datasets (Figs. 8hS6h,8iS6i).However, all three cases feature aindicate similar regions with low basal dragfriction coefficient in fast flow regions (>1500 m yr -1 , cyan contour in Fig. 8S6), which is approximately consistent coincident with the boundary of the FG downstream basin.
The simulated surface velocities from BEDZC (Fig. 8eS6e) and BEDMC (Fig. 8fS6f) match the observed surface velocities better than those from CONTROL (Fig. 8dS6d) in the regions around the ice front/grounding line and more broadly for velocity exceeding 1000 m yr -1 .This point is supported by the computed RMSD of surface velocity mismatches (Table S5).One possible cause of the different basal friction coefficient distributions in these inversions might be the changed surface topography during the surface relaxation, especially near the ice front (Figs.S7).
The deeper retrograde bed in the CONTROL simulation may indicate increased vulnerability to marine ice sheet instability, and more overestimation of surface velocity is found around the grounding line (Fig. 8d).One possible cause of the different basal shear stress in these inversions might be the increased slope caused by the surface relaxation.However, we find the inversion process is not sensitive to the surface relaxation, and this is further discussed in Sect.4.3.It means a high-accuracy bedrock topography data is very important for inverse modeling owing to the fact that the bedrock resolution around the grounding line determines the ice dynamics .Comparisons of the distributions of velocity mismatch and of C between BEDZC and BEDMC does not provide a clear direct insight into which is the best more accurate basal geometry for modelling this the Fleming system.We compute tThe computed root mean square errors (RMSED) of the velocity mismatch for the regions with velocity >1500 m yr -1 (Table S5) is only slightly higher for , and find the RMSE of BEDMC (62.60 m yr -1 ) is than for marginally larger than the RMSE of BEDZC (61.78 m yr -1 ), and both are much lower than CONTROL.Both BEDMC and BEDZC use the 2008 surface DEM and this improvement over the Bedmap2 surface DEM in CONTROL appears significant, even before turning to the matter of ice thickness.While bBoth cases use the ice thickness extracted using the mass conservation mechanismapproach (which is independent of surface geometry) and the bed geometries are accordingly more similar to each other than they are to CONTROL (see Fig. 2b-d).However, BEDZC maintains better internal consistency with the 2008 surface elevation, since it results in the mass conserving ice thickness Hmc being employed, whereas, by the construction of bed_mc (Eq.2), the ice thickness in BEDMC is not entirely consistent with mass conservation, although still a more physically motivated interpolation than bed_bm in CONTROL.The BEDMC and BEDZC ice thicknesses clearly differ by the difference between the Bedmap2 and 2008 DEMs, which should be greatest in areas of greatest lowering, and as we see BEDMC provides a useful sensitivity test case.
Since bed_zc is extracted from the accurate and contemporary DEM2008, it should also incorporate into the bed geometry (via Hmc) more detail from the then current surface, compared to bed_mc, extracted from Bedmap2's surface DEM, which was generated over a longer time range data used in current study than BEDMC.. Therefore, bed_zc is suggested as the best current bedrock elevation data for further ice sheet modelling of the WIS-FGLFG system.

4 Sensitivity to ice front position and boundary condition
All the inversions presented so far feature both a bandsome sticky spots of with high basal dragfriction coefficient near the ice front of the Fleming Glacier (right column of Fig. 5 and left column of Fig. 8S6) and a similar localized overestimate of upper surface velocities at the ice front (right column of Fig. 6 and middle column of Fig. 8).We now consider causes for possible uncertainties about in the force applied to the ice front, and whether the high basal friction near the ice front is likely to be a feature of the real system or emerges from the inversion process as a compensating response to incorrect boundary forcing by the inversion process.These possible causes include uncertainty in local bedrock elevation (or equivalently ice thickness), uncertainty in the geoid-ellipsoid height conversion, uncertainty in observed sea level, uncertainty in exact ice front position and grounding line position, uncertainty in surface velocity, and uncertainty in potential backstress due to ice mélange and/ or grounded icebergs in contact with the ice front.The sensitivity to various bedrock uncertainty datasets has been discussed in Sec.4.23.In our model domain we assume the 2008 grounding line is consistent with the 1996 grounding line, which has an error of several km on fast-moving ice (Rignot et al., 2011a) and might have changed since 1996.The frontal surface elevation is from the SPOT DEM data in Jan 2008, which shows the ice front position is ~1.5 km downstream of the 1996 grounding line position.Since such a narrow residual ice shelf was considered unlikely to have a major influence we constructed the model geometry to have the ice front coincide with the 1996 grounding line for simplicity, i.e. all ice is considered grounded.In this frameworkBy assuming the ice front position to coincide with the 1996 grounding line, uncertainty about the bedrock depth at the ice front feeds in to significant uncertainty in the total restraining force from ocean pressure.Regarding velocities, Friedl et al. (2018) presented evidence that an acceleration phase occurred on the Fleming Glacier around between MarchJan-April 2008, but we are not sure the specific month of the surface velocity data used in this the current study was extracted from measurements in Fall 2007 and 2008 (Rignot et al., 2011b).ThisIt means the surface velocity data, which is provide the target to be matched by the control inversione process, might not be consistent with the DEM data used here (acquired in Jan 2008).
To explore the influence of these different sources of uncertainty, we adopt different ice front positions and effective sea level heights within our vertical reference frame to apply a range of ocean pressures to the ice front as described in Sect.3.6 (IFBC1-3 and IFP1-2, Table 2).
Experiments with different ice front positions (IFP1-2 in Table 2) directly affect the ice thickness and bed elevation at the ice front, which affects the ice front pressure condition.The simulated basal friction coefficients (left column in Fig. 7) show that the high sticky spots near the ice front migrate with the ice front position but with different patterns.The experiment IFP1 with a seaward shifted ice front position shows a decrease in magnitude of the high friction spots (Fig. 7b) and a better match with the observed velocity (Fig. 7e), while the IFP2 with a retreated ice front shows an increased C (Fig. 7c) and worse surface velocity match (Fig. 7f) compared with CONTROL experiment (Figs. 7a,7d).In experiment IFP1, thinner ice at the ice front leads to a relatively smaller ice velocity compared with CONTROL, so the model does not need to increase C to match the observed surface velocity.This does not mean that ice front position in IFP1 is more accurate than CONTROL, since the time inconsistency of surface DEM data, ice front and grounding line position, and surface velocity data is the obstacle to obtaining a reliable basal friction pattern.Therefore, we speculate that some of the high basal friction spots near the ice front are artefacts.However, we do not exclude the possibility of high basal friction spots caused by the pinning points located at the 1996 grounding line, which is also proposed by Friedl et al. (2018).An accurate location of the ice front and grounding line is clearly important for inverse modelling of fast flowing glaciers like the Fleming Glacier.
A higher sea level in the ice front boundary condition imposes a higher pressure at the ice front, i.e. a higher total retarding force, and we impose these different boundary conditions as a proxy for the sources of uncertainty discussed above.
Basal dragfriction coefficients C simulated from the IFBC1-2 and CONTROL experiments (Figs.8a-c) present different similar patterns but differ systematically around the ice front regions of the FGL (within ~1 km of the grounding line).Experiments with higher sea levels display smaller C there (Fig. 98, left column) and provide a better match between modeled and simulated surface velocities (Fig. 98, rightmiddle column), which is consistent with the computed RMSD of the surface velocity mismatch (Table S5).If the applied ice front boundary condition underestimates the real world forcing, the inversion process will compensate by increasing the basal dragfriction in this region.However, the large vertical shear strain rate imposes a limit to how much increasing basal drag can reduce the surface velocity, which could explain why the mismatch between the modeled and observed velocity is still large in the narrow band near the ice front (Fig. 9, middle column).For the fastflowing region (velocity > 1500 m a -1 ), the decreased basal shear force from IFBC2 to IFBC3 (~1.1×10 11 N) roughly matches the increased the ice front pressure over a 6 km length of ice front (~2.8×10 11 N).
Experiment "IFBC3", with an extreme assumption of applying ice pressure corresponding to a neighbouring column of ice matching the ice front, shows very small basal dragfriction for the ice front area around the grounding line, and also resulted in lower drag over the downstream basin (Fig. 9d8d).However, "IFBC3" introduces a much greater mismatch to the observed surface velocities (Fig. 9h), with lower simulatedunderestimated velocities over a substantial region extending upstream from the ice front and greater overestimate of velocities further upstream.This is only a sensitivity test but implies a potentially suitable ice front pressure may lie between "IFBC2" and "IFBC3".This set of experiments also suggests that moderate changes influence only a limited area.It is hard to decide the best ice front boundary condition here owing to the lack of precise bedrock data (as seen above) and difficulty of estimating the additional pressure from the partly detaching icebergs and ice mélange.As an indicator, the simulated ice mélange depth-integrated back stress (~1.1×10 7 N m -1 ) required to prevent the iceberg rotation at a calving front (Krug et al., 2015) would be comparable to an additional ~2.3 m in sea level in terms of ice front boundary condition for the Fleming Glacier.The thickness and density of mélange may affect this estimation.But it is certainly clear that the ice front boundary conditions can have a significant effect on the inversion results near the grounding line.
The different ice front boundary conditions also lead to minor significant differences in the surface relaxation at the ice front, with lower sea levels leading to slightly greater lowering and corresponding steepening of the surface adjacent to the ice front (Fig. 9, right columnfor example, ~8 m lowering from IFBC2 to CONTROL and from CONTROL to IFBC1 at the ice front).The differences in surface elevation are localized to the ice front zone, with the relaxation over the rest of the domain essentially unaffected, even except for the most extreme forcing.This The lowered surface at the ice front in experiments IFBC1 and CONTROL is apparently the consequence of rapid deformation due to its own weight spreading(longitudinal extension with locally high vertical shear) of an ice cliff, which is over 100 m higher than the control sea level at the ice front due to its own weight.Thus, the inversions are potentially influenced both by the ice front condition directly in the overall momentum balance and also by the increased local driving stress due to this artificially steeper surface slope.The band of higher basal drag near the ice front may be partly a response to these issues.However, an additional simulation (not shown), in which a high sea level was used for the surface relaxation step and a low sea level for the inversion, gave a relaxed surface very close to observations (no steepening) and stillthe sticky spot located ~1 km upstream the ice front is a persistent feature except for the experiment IFBC3 showed the high basal friction band near the ice front.This implies that the high friction near the ice front is directly sensitive to the boundary condition at the ice front but not to associated artifacts in the surface relaxation.Friedl et al. (2018) Based on the experiments IFP 1-2 and IFBC1-3, we cannot be sure whethersuspect the high friction near the ice front is likely a real feature, an arteifact due to errors in the ice front boundary condition but we cannot rule out the possibility that this may be a real feature, or a combination of both.However, the impact diminishes rapidly with distance inland for moderate sea level shifts, which do not affect the general pattern of basal friction coefficients or the quality of the velocity matching more than ~2 km upstream of the grounding line.

Conclusions
We have obtained a basal dragfriction coefficient distribution for the Wordie Ice Shelf-Fleming Glacier system in 2008, using an iterative spin-up scheme of simulations, observed surface velocities and a detailed surface DEM.We explored the sensitivity of the inversion for basal dragfriction to three four inputs to the modelling process.Within the approximation of using simulated steady-state ice temperatures, we showed that three multiple temperatureinversion cycles are necessary of iterationto removed the influence of initial englacial temperature assumptions, at least for plausible initial temperature assumptions, and that a poor initial assumption will lead to a requirement for a greater number of cycles.In contrast to the observed low sensitivity to the englacial temperature of outlet glaciers from the Vestfonna Ice Cap (Schäfer et al., 2014;Schäfer et al., 2012), the first cycle of our iterative process showed that the inferred basal stress of the Fleming Glacier system is highly sensitive to the englacial temperature distribution.This conclusion is expected to also apply to other fast-flowing glacier systems with a significant dependence onthat feature high rates of the internal deformation.For such glacier systems, a multiple-cycle spin-up scheme is likely to be necessary.
For oOur inversionmodel of the Wordie Ice Shelf-Fleming Glacier system, our is highly sensitiveity tests to different the choice of ice flow enhancement factors and basal elevation datasets indicate a high dependence of basal inversion on the accuracy of bed topography.
The "bed_zc" bed topography, which used ice thickness determined using the mass conservation method for the fast-flowing regions, using contemporary velocities and ice thinning rates, and applied to the then current DEM, is suggested as the best current bed topography for further simulations in this region.
For the Wordie Ice Shelf-Fleming Glacier system, which we treated as grounded adjacent to the ice front, the inferred basal dragfriction coefficient near that grounding lineice front is sensitive to the ice front position and ocean pressure boundary condition, emphasizing the importance of the normal force on the ice front and the accuracy of ice front positions.These factors have a very low impact on basal friction coefficients more than a few kilometers upstream of the grounding line, but may still be important when using inversion to initialize transient simulations, due to the high sensitivity of transient ice dynamic behavior to grounding line dynamicsThis finding, combined with the sensitivity of surface relaxation to ice front boundary condition, implies that an accurate representation of the ice front boundary will be important for inverse modeling and transient simulations of the Wordie Ice Shelf-Fleming Glacier system.Table 2 Experiment lists.n/a is short for "not applicable".EF and SL are short for 980 "enhancement factor" and "sea level", respectively.IF1 and IF2 represent the ice front positions located downstream and upstream of the 1996 grounding line position (Rignot et al., 2011a), respectively.

Figure 1 .
Figure 1.(a) The location of the Wordie Ice Shelf-Fleming Glacier system in the Antarctica Peninsula (pink polygon).(b) Surface speed in 2008 with a spatial resolution of 900 m obtained from InSAR data (Rignot et al., 2011c) for the study regions.Colored lines represent the ice front position in 1947 (red), 1966 (brown), 1989 (green), Apr 2008 (blue), and Jan 2016 (magenta) obtained from Cook and Vaughan (2010), Wendt et al. (2010), and Zhao et

Figure 3 .
Figure 3. (a) Mesh structure of the domain in the current study with surface velocity in 2008 (Rignot et al., 2011c) and the zoomed-in map for (b) the Fleming Glacier and (c) the Prospect Glacier.

Figure 4 .Figure 5 .
Figure 4. Flow chart of simulation spin-up with three cycles.

Figure 6 .
Figure 6.Distribution of basal friction coefficient C (MPa m -1 yr) (left column), mismatch between the observed and modeled surface velocity (observed minus simulated; middle column), and the difference between the observed initial surface and relaxed surface elevation (observed minus relaxed; right column) from experiments: (a, e, i) CONTROL (first row), (b, f, j) EF1 (second row), (c, g, k) EF2 (third row), and (d, h, l) EF3 (forth row)Mismatch between the observed and simulated surface speed in 2008 (observed minus simulated) from experiments: (a-c) TEMP1, (d-f) TEMP2, (g-i) CONTROL, and (j-l) TEMP3.The left (a, d, g, j), middle (b, e, h, k) and right columns (c, f, i, l) are the inferred basal drag coefficients from Cycle 1, Cycle 2 and Cycle 3, respectively.The black, and yellow and cyan solid lines represent observed surface speed contours of 100 m yr -1 , and 1000 m yr -1 , and 1500 m yr -1 , respectively.

Figure 7 .Figure 98 .
Figure 7. Distribution of basal friction coefficient C (MPa m -1 yr) (left column), the mismatch between the observed and modeled surface velocity (observed minus simulated; middle column), and the difference between the observed initial surface and relaxed surface elevation (observed minus relaxed; right column) from experiments: (a, d, g) CONTROL (first row), (b, e, h) IFP1 (second row), and (c, f, i) IFP2 (third row).The black, yellow, and cyan solid lines

Table 1 .
List of parameter values used in this study.