Brief communication: Estimating the ice thickness of the Müller Ice Cap using an inversion of the shallow ice approximation

The Müller Ice Cap will soon set the scene for a new drilling project. Therefore, ice thickness estimates are necessary for planning since thickness measurements of the ice cap are sparse. Here, two models are presented and compared, i) a simple inversion of the shallow ice approximation (SIA inversion) by the use of a single radar line in combination with the glacier outline, surface slope, and elevation, and ii) an iterative inverse method using the Parallel Ice Sheet Model (PISM). The two methods mostly agree about a good drill site candidate. However, the new semi-empirical SIA inversion is insensitive to mass 5 balance, computationally fast, and provides better fits.

. The black and red dots are long-term HIRHAM5 SMB averages from 1980-2016, with the red dots corresponding to the SMB within the White Glacier polygon. qualify as good drill site locations due to limited horizontal ice flow. Furthermore, from ITMIX1 it is not possible to determine 25 what model approaches are to be preferred on ice caps. Thus, to narrow down areas of conducting ground based radar measurements, the ice thickness is modelled using two different techniques that differs greatly in computational demands, input data and model setup. One is a simple inversion of the shallow ice approximation (SIA inversion), only relying on the surface elevation, the ice cap outline and a single Operation IceBridge (OIB) flight line with thickness measurement. This is a fast new method that differs from existing SIA approaches due to its semi-empirical nature which makes it less sensitive to mass 30 balance, steady state assumptions and ice flow physics. The second is an iterative inverse method using the Parallel Ice Sheet Model (PISM) as a forward model, which has to be forced with both climate data and an initial guess of ice cap geometry.

Data
Ice thickness measurements taken from the Multichannel Coherent Radar Depth Sounder from Operation IceBridge (OIB) (Paden et al., 2010(Paden et al., , updated 2019 are used to validate all models and to perform the SIA inversion. Only a few flight lines The surface elevation is obtained from the Arctic Digital Elevation Model (ArticDEM, Porter et al. (2018)) and is shown as contours on Fig. 1. To ensure continuity over the entire ice cap, the OIB surface elevations are not used when present, thereby only using the ice thicknesses from OIB.

40
As climate inputs to initialise PISM, SMB and 2-m temperatures are obtained from the Regional Climate Model HIRHAM5 (HIRHAM5) Mottram et al., 2017). Since steady state is assumed in the models presented here, the longterm average of these fields are used. They are calculated from the yearly model results of HIRHAM5 from 1980HIRHAM5 from -2014HIRHAM5 from and 1980HIRHAM5 from -2016 for the 2-m temperature and the SMB, respectively. The long-term average SMB is shown in Fig. 1, and compared to in-situ measurements from White Glacier (Cogley et al., 1996). Furthermore, a uniform geothermal heat flux of 0.055 Wm −2 45 is used based on Minnick et al. (2018).
The initial guess of ice cap geometry for PISM is obtained from ArcticDEM and the bedrock topography from  All data are interpolated onto a 100 m grid and 900 m grid for the SIA inversion and PISM, respectively. Thus, gaps in the before mentioned bedrock topography are filled in the interpolation process. The initial PISM ice thickness is calculated from the interpolated versions of surface and bedrock elevations.

55
The SIA is one of the oldest ice flow models, which has proven to show good results on ice sheets and glaciers in areas with no to little basal sliding. Combining this with the simplicity of it and how computationally cheap it is, still makes it a good and often used zero order approximation. Thus, a simple inversion of the SIA to estimate ice thicknesses from sparse data is tested.
From theory the SIA without sliding reads

60
where Q is the horizontal ice flux, A the rate factor, n the creep exponent, H the ice thickness and τ b is the basal shear stress given by Here ρ is the density of ice and α is the surface slope. Substituting eq. (2) into eq. (1) and introducing the constant c, implies that eq. (1) can be reduced to Isolating the ice thickness thus results in Linearity can be obtained by introducing the three tuning parameters; k = log c 1 n+2 , a = 1 n+2 , and b = n n+2 , resulting in Hence, the assumption is that if the ice thickness is known on parts of the ice cap and the surface slope and ice flux are known everywhere on the ice cap, one can perform a least squares regression to obtain the tuning parameters a, b, and k in the areas with known ice thicknesses and apply those in the areas with unknown ice thickness. Thereby, assuming that the areas with known ice thicknesses are representative of the entire ice cap. Note that the empirical regression is insensitive to global 75 multiplicative errors in Q as that will be accounted for by adjustments to k. The approach is semi-empirical because the form of eq. (6) is justified by the SIA theory, but the parameters are tuned empirically. The least squares estimates of a and b are not necessarily consistent with the exponents in eq. (5), which was derived for isothermal ice with no sliding. However, the empirical calibration has freedom to adjust a and b so that it better matches the physical reality.

Surface slope 80
The surface slope is obtained from ArcticDEM by smoothing the surface elevation using a Gaussian filter with a standard deviation of 250 m, from which the surface slope is calculated. The smoothing is done to prevent surface depressions on the ice cap, thus ensuring that there is a downward slope from all top points of the ice cap all the way to the ice margin. To prevent H → ∞ in low sloping areas a minimum slope of 0.01 is introduced on the ice cap. This corresponds roughly to 0.6 • , and is relatively lower than the 2 • -5 • minimum slope used in other similar models (Farinotti et al., 2009(Farinotti et al., , 2017. Furthermore, the 85 surface slope is set to be zero outside of the present day ice cap margin.

Ice flux
The ice flux of the MIC, and most ice caps and glaciers in general, is not known, why the balance flux is calculated. This is done based on the top model approach described in Quinn et al. (1991).
The ice then flows from cell i,j to cell i,j+1 by where nbh is the neighbourhood of cell i,j . At the very top of an ice cap the ice flux is equal to the SMB, and all areas further down the ice both have a flux input from the SMB at that given point, and from the ice flowing in from higher elevations. In this model we use a uniform SMB of ones in every ice covered cell, as it has been shown that the SMB only has little influence on the modelled ice thickness (Zinck, 2020). This is due to the fact that any multiplicative errors in SMB can be captured in k in eq. (6). Do notice, that this might not hold in areas with negative SMB. According to in-situ measurements from White 100 Glacier and HIRHAM5 SMB (Fig. 1), the areas with negative SMB are restricted to the outlet glaciers and the glacier tongues at the margin, which are anyway outside the area of interest as surface melt at the ice core site should be avoided.
Once again, to prevent H → ∞ in eq. (6), a minimum ice flux of 10 −3 m s −1 per grid cell area is introduced.

PISM
To estimate the ice thickness using PISM (stable version 1.1.4, Bueler and Brown (2008) shelf approximation (SSA) (Bueler and Brown, 2008). In the first of a total 10 iterations, PISM is forced with an initial guess of geometry, with surface elevation from ArcticDEM (S ref ) and bedrock elevation from Farinotti et al. (2019). The surface elevation is kept constant when initialising every iteration assuming that ArcticDEM represents the ground truth. However, 110 the bedrock topography is adjusted after each iteration (n) by adding bedrock in areas with too low surface elevation (S n ) as compared to S ref , and vice versa removing bedrock in areas with too high surface elevation. Thus, the bedrock topography used in iteration n + 1 (B n+1 ), is given by where B n is the bedrock used to force PISM in iteration n. K is a relaxation parameter which ensures that the bedrock 115 topography is not overcompensated (van Pelt et al., 2013;Koldtoft et al., 2021). In this study we use a varying relaxation parameter with elevation, to prevent overcompensation in bedrock and ice build-up on the outlet glaciers west of the ice cap (see Fig. 1). K is set to 0 below 500 m and increases linearly to 0.5 at 1000 m above sea level, from where it is kept constant.
Thus, under the assumption of steady state, the idea is that the bedrock should come closer to the true bedrock after each iteration. It is important to note that this is not a guarantee.

120
Since the properties of the bedrock below the ice cap are unknown, a bigger parameter space of till friction angles (φ) and enhancement factors (E) are tested using the following possible combinations φ ∈ {10, 20, 30, 40} To test convergence and find the most suitable parameters for φ and E, the root mean squared error (RMSE) between the 125 following misfit metrics are evaluated: -Modelled surface elevation and ArcticDEM. -Modelled ice thickness and OIB ice thickness (both the entire flight line and the part in between the triangle and the square on Fig. 1).
-Modelled surface velocity and median surface velocities obtained through feature tracking of Landsat 8 images (Zinck,130 2020).
For a further description and analysis the reader is referred to Zinck (2020).

Results
After solving the least squares regression of the SIA inversion using the OIB ice thicknesses, the parameters a, b, and k are applied in combination with the surface slope and ice flux over the entire domain to move from one to two dimensions. The 135 modelled ice thickness is also transformed into bedrock elevation by using ArcticDEM. The OIB cross section with corresponding SIA inversion bedrock is shown in Fig. 2, where the surface elevation is from ArcticDEM (blue-white intersection), the ice thickness in white is from OIB, resulting in the bedrock given by subtracting OIB from ArcticDEM in brown. The SIA inversion is only tuned to the non-shaded part of the ice cap, from the red triangle to the red square in Fig. 1 and 2, since it is the main ice cap which is of interest in this study and because the uniform SMB used in the SIA flux calculation is not valid at 140 low elevations. Fig. 3 (left) shows the modelled bedrock topography (bottom) and ice thickness (top) of the entire ice cap. The RMSEs of the ice thickness compared to OIB are 145 m and 126 m by including and excluding the outlet glaciers, respectively.
Further, the mean absolute deviation (MAD) of the entire OIB flight line is 109 m.
In the PISM method, 12 different combinations of till friction angles and enhancement factors are tested. Based on the before mentioned convergence metrics described in Sect. 3.2 a till friction angle of 10 • and an enhancement factor of 6 is the 145 parameter combination with the best results. The ice thickness after the tenth iteration is used as the main PISM result, from which the bedrock topography is calculated using ArcticDEM. The PISM bedrock of the OIB cross section is shown in Fig.