Analysis of ice-sheet temperature profiles from low-frequency airborne remote sensing

Abstract

Ice internal temperature and basal geothermal heat flux (GHF) are analyzed along a study line in northwestern Greenland. The temperatures were obtained from a previously reported inversion of airborne microwave brightness-temperature spectra. The temperatures vary slowly through the upper ice sheet and more rapidly near the base increasing from ~259 K near Camp Century to values near the melting point near NorthGRIP. The flow-law rate factor is computed from temperature data and analytic expressions. The rate factor increases from ~1 × 10⁻⁸ to 8 × 10⁻⁸ kPa⁻³ a⁻¹ along the line. A laminar flow model combined with the depth-dependent rate factor is used to estimate horizontal velocity. The modeled surface velocities are about a factor of 10 less than interferometric synthetic aperture radar (InSAR) surface velocities. The laminar velocities are fitted to the InSAR velocities through a factor of 8 enhancement of the rate factor for the lower 25% of the column. GHF values retrieved from the brightness temperature spectra increase from ~55 to 84 mW m⁻² from Camp Century to NorthGRIP. A strain heating correction improves agreement with other geophysical datasets near Camp Century and NEEM but differ by ~15 mW m⁻² in the central portion of the profile.

Introduction

Intra-glacial temperature remains one of the most elusive of the various parameters required to understand and predict the behavior of the polar ice sheets. Temperature is one factor that moderates ice rheology and hence the rate at which ice flows over the bed. Knowledge of the temperature thus influences mass continuity calculations of ice flux through a volume and the prediction of ice flow rates under changing climate conditions.

Temperature is measured in the few and widely scattered deep boreholes located in both Greenland and Antarctica, but such measurements provide only a limited validation for the broader temperature fields and basal conditions deduced from coupled dynamic and thermodynamic models. Consequently, approaches utilizing remote-sensing data are being explored to fill this knowledge gap. In contrast to earlier approaches based on the use of ice-penetrating radar (MacGregor and others, 2015, 2016; Schroeder and others, 2016; Jordan and others, 2018), recent studies have introduced the measurement of low-frequency microwave emissions to directly sense physical temperature at depth in the polar ice sheets. By observing brightness temperatures (Tb) over a wide band of microwave frequencies, the variations in penetration depth over the band enable inversion of the measured Tb spectra into the physical temperature depth profile. As discussed by several authors (Macelloni and others, 2018; Johnson and others, 2021), new technologies for low-frequency wide-band radiometric measurements are enabling continued and future airborne and spaceborne remote sensing of the physical temperature field of polar ice sheets.

Yardim and others (2021) report an inversion of brightness temperature spectra from a 2017 airborne campaign over northwestern Greenland into a set of geophysical parameters necessary for calculating ice physical temperatures at depth. Here we take a next step by illustrating the glaciological applications of the Yardim temperature profiles and associated parameters. We first discuss the obtained temperature profiles and examine variations of basal temperature along the flight line, which proceeded along an ice divide and connected multiple deep borehole sites. We go on to use the temperature data to estimate the ice flow-law rate factor. The temperature-dependent rate factor relates the strain rate to the applied stress through Glen's Law for polycrystalline ice, and so is an important control on the ice motion (Glen, 1958). By adopting a simple, laminar flow model of ice motion, we estimate an enhancement to the rate factor through comparison with interferometric synthetic aperture radar (InSAR) measurements of surface speed. The enhancement factor incorporates the additional influences of ice grain-size and crystalline fabric (the distribution of crystal orientations) on the rate at which the ice will deform. The modified rate factor and the laminar flow model are then used to compute the integrated strain heating over the ice column. The strain heating is subsequently used to correct the remotely sensed geothermal heat flux (GHF) for comparison with independent estimates of GHF obtained using a geothermal model of the lithosphere beneath Greenland (Martos and others, 2018).

Temperature data

The Ohio State University developed the Ultra-WideBand software-defined microwave RADiometer (UWBRAD) to investigate the proposal by Jezek and others (2015) that low-frequency, wide-band radiometers could sense physical temperatures deep within polar ice. UWBRAD was designed to observe circularly polarized 0.5–2 GHz thermal emissions in a nadir viewing geometry (Johnson and others, 2021). Given the heavy use of this portion of the radio spectrum, radio frequency interference (RFI) originating on board the aircraft or from external sources is removed using multiple processing algorithms (Andrews and others, 2018). The observed 0.5–2 GHz spectrum is subdivided into 12 individual frequency sub-bands, typically between 80 and 100 MHz, within which RFI contaminated portions of the data in time and frequency are removed.

UWBRAD was flown aboard a Bassler aircraft in 2016 and 2017 across northern Greenland (Fig. 1). Both flights originated from Thule Airbase. The 2016 flight ended near Camp Century after experiencing a technical problem with the radiometer. The 2017 campaign proceeded toward the Camp Century and the North Greenland Eemian Ice Drilling (NEEM) drilling sites. The flight approached the North Greenland Ice Core Project (NorthGRIP) site but was diverted due to weather before the camp was reached. A second flight in 2017 proceeded northwards, and the microwave spectra obtained were used in an additional study of sea-ice remote sensing over the Lincoln Sea (Jezek and others, 2019).

Fig. 1. 2016 (blue) and 2017 (heavy black) flight lines over the Northwestern Greenland Ice Sheet. Data from 2017 discussed in this paper are identified by the red line, which is the portion of the flight located near the ice divide. Slope-derived flow lines are thin black lines (Thomas and others, 2000). Locations of deep boreholes are shown by white dots. Base map is a 16 July 2016 MODIS Terra 3-band optical image.

Figure 1 also includes ice-sheet flow lines derived from surface elevation data (Thomas and others, 2000) highlighting the location of the ice divides. The 2017 flight line diverged from the ice divide along the path from Camp Century to NEEM with a maximum deviation of ~50 km. The flight then closely followed the ice divide from NEEM to ~50 km north of NorthGRIP where weather conditions forced a deviation from the flight plan. The data acquired during the flight were averaged over six successive measurements during which the aircraft traveled an average of 1.1 km (varying between 0.8 and 1.4 km depending on the speed of the aircraft). The aircraft typically flew at an altitude of ~500 m over the ice-sheet surface.

The 2016 and the portion of the 2017 flight which overlapped the 2016 route overflew Greenland melt facies. Inverting the spectra for physical temperature in these regions remains a challenge because of the strong scattering that can occur within wet snow and percolation facies. Nevertheless, information about the size and density of ice lenses, which represent the dominant scattering sources at UWBRAD's frequencies, has been retrieved (Jezek and others, 2018).

Spectral data collected in 2017 over the dry snow facies between Camp Century and ~50 km northwest of NorthGRIP were sufficiently free of scattering to enable Yardim and others (2021 and hereafter referred to as Yardim) to invert the spectra into the estimates of physical temperature with depth (red line in Fig. 1). Yardim implemented a Bayesian inversion scheme using an ice-sheet radiative transfer model and a simple parametric model of ice-sheet temperatures first proposed by Robin (1958).

Upwelling microwave radiation from the ice sheet can be described approximately using a ‘cloud’ model (van der Veen and Jezek, 1993) in which brightness temperature contributions from morphologically homogeneous ice layers at varying depths and physical temperatures within the ice sheet are combined and then multiplied by the transmissivity at the upper boundary of the ice sheet. The model must further account for the impact of fine, sub centimeter layering in the firn that impacts the frequency dependence and amplitude of the firn transmissivity (Tan and others, 2015; Tan and others, 2020). Because individual firn layers are typically much thinner than a wavelength, the transmission across a single, low-loss layer is near unity, but the combined impact of many layers in the upper 50–100 m of the ice sheet produces a significant effect. Consequently, the model used in Yardim accounts for stochastic density variations in the upper firn. These density fluctuations typically damp the brightness temperature variation with frequency and reduce mean brightness temperatures.

Combining these effects, the model represents ice brightness temperatures as

(1)$$T_{\rm B} = T_{\rm B}^{{\rm cloud}} ( {T( z ) } ) \times \Gamma ( {\rho ( z ) } ) \times {\rm exp}( {-\tau_{\rm s}} ) , \;$$

where T _B is the modeled brightness temperature, Γ is the bulk transmissivity at the surface representing the stochastically layered firn, and τ _s represents any scattering loss from ice lenses (which are neglected for this study area). T _B^cloud is the ‘cloud model’ prediction given by

(2)$$T_{\rm B}^{{\rm cloud}} ( {z = H} ) = \mathop \int \limits_0^H ( \kappa _{\rm a} + \kappa _{\rm s}) T( z ) {\rm e}^{-\mathop \int \limits_z^H ( \kappa _{\rm a} + \kappa _{\rm s}) {\rm d}{z}^{\prime}}{\rm d}z + T_{\rm B}( {z = 0} ) {\rm e}^{-\mathop \int \limits_0^H ( \kappa _{\rm a} + \kappa _{\rm s}) {\rm d}z}. $$

Here, T _B^cloud(z = H) is the brightness temperature at the ice-sheet surface just prior to transmission across the surface and which consists of an integration of the emission from layers at each depth z plus the attenuated subglacial emission (T _B(z = 0)). The latter is taken as the physical temperature of the ice-sheet base multiplied by the transmissivity across the ice-sheet base. We assume also that the scattering coefficient (κ _s) within the ice is small compared to the ice absorption coefficient (κ _a). All parameters in Eqns (1) and (2) are frequency-dependent, excepting the ice thickness H and the physical temperature profile T(z).

The transmissivity Γ(ρ(z)) depends on frequency and the firn layer structure, the latter of which also impacts the spectral shape and magnitude. Because the level of detail required to model these structures exceeds conventional in situ firn sampling techniques, the Yardim retrieval approach used the physical temperatures measured in each borehole as a constraint on local firn density parameters. The information obtained from this process was then interpolated to allow temperature profile retrievals at other points along the study line.

The Bayesian inversion of the brightness temperature uses a Robin steady-state temperature model along with prior estimates for the Robin model's parameters. The Robin temperature model assumes a planar stratified medium parameterized by its surface temperature, GHF, accumulation rate and ice thickness. Vertical velocities vary linearly with depth with the surface velocity equal to the surface mass-balance rate and zero at a frozen base. Horizontal advection of upstream ice is also ignored in the model; this is reasonable given the proximity of the study line to the ice divide. Using an analytical solution of the heat flow equation for these conditions, the Robin solution for the temperature profile (T(z)) is obtained:

(3)$$T( z ) = T_{\rm s}-\displaystyle{{G\sqrt \pi } \over {2k_{\rm c}\sqrt {\displaystyle{M \over {2k_{\rm d}H}}} }}( {\rm erf}\left({z\sqrt {\displaystyle{M \over {2k_{\rm d}H}}} } \right)-{\rm erf}\left({H\sqrt {\displaystyle{M \over {2k_{\rm d}H}}} } \right)).$$

Here, erf is the error function, G is the geothermal heat flux, T _s is the surface temperature, k _c is the ice thermal conductivity (2.7 W m⁻¹ K⁻¹), k _d is the ice thermal diffusivity (45 m² a⁻¹), H is the ice thickness, M is the accumulation rate in ma⁻¹, and z is a depth variable here defined so that z = H is the ice surface. Note that the model ignores any climate signals of long duration, but these are generally on the order of the uncertainty of our analysis. Any strain that induced heating near the ice-sheet base can also be considered in the model as having been added to the ‘true’ GHF. Although the Robin model requires significant assumptions to be applicable, it is sufficient for the purposes of the current initial study near the ice divide. The inversion scheme developed by Yardim yields estimates of the surface temperature, accumulation rate and GHF parameters that are then used in the Robin model to predict ice temperature profiles. Ice thickness in the retrieval is obtained from Operation IceBridge multichannel coherent radar depth sounder (MCoRDS) data (Paden and others, 2010).

The physical temperature field estimated by Yardim is shown in Figure 2; selected temperature profiles along the flight are highlighted in Figure 3. The Robin temperature model governs the functional form of the temperature field that exhibits the characteristically slow change in temperature through the upper half of the ice sheet and more rapid warming in the lower half of the ice column (Fig. 3). This behavior is largely driven by surface temperatures that cool by ~7 K along the transect and relatively high (10s of cm) accumulation rates that drive the downward descent of cold ice at the surface. In the lower portion of the ice sheet, geothermal and strain heating drive the increase of temperatures. The UWBRAD-retrieved surface temperatures compare favorably (within ~1 K) to Regional Atmospheric Climate Model (RACMO) model results (van Meijgaard and others, 2008) and to Moderate Resolution Imaging Spectroradiometer (MODIS) surface temperature data provided by Hall and others (2012) (see Fig. 8 of Yardim and others (2021) for comparisons). The UWBRAD-derived accumulation rates also compare favorably with ancillary information available for the borehole locations.

Fig. 2. Physical temperature (K) estimated along the flight line from Yardim and others (2021). The longitudes of Camp Century (C) and NEEM (NE) are indicated. The study line terminated west of NorthGRIP and that station is not shown (see Fig. 1 for relative locations). Temperatures are slowly varying near the surface. Temperatures are warm near the bed and are warmest in the vicinity of the NorthGRIP drill site located near the end of the profile. White bars at top of figure correspond to locations of temperature profiles shown in Figures 3 and 4.

Fig. 3. Temperature profiles selected every 50 points or ~56 km along the profile – longitudes are indicated in the plot titles. Horizontal axes are temperature in Kelvin; vertical axis represents depth (m) within the ice sheet. The profiles are characteristic of locations where horizontal advection of upstream ice is largely absent, accumulation rates are on the order of 30 cm a⁻¹, and surface temperatures are ~244 K.

The std dev. of the physical temperature (Fig. 4) estimated by the retrieval algorithm at each point ranges from ~1 K in the upper part of the ice column to ~3 K at the base (see Appendix in Yardim and others, 2021). Additional accumulated errors are potentially introduced as the calculation proceeds sequentially along the line.

Fig. 4. Temperature uncertainty estimated at depth for selected locations – longitudes are indicated in titles. Horizontal axis is temperature uncertainty (Kelvin); vertical axis is depth (m) within the ice sheet.

Figure 5 plots retrieved temperatures at the base of the ice sheet along the study line. These temperatures are coldest (~256 ± 3 K) at the east of Camp Century and then steadily rise to ~271 ± 3 K at the line's end. Comparisons with borehole basal temperatures are also shown, and match at the Camp Century (Robin, 1983) and NorthGRIP locations to within the estimated uncertainties. The borehole temperatures at NEEM (G. Clow, personal comm., 2018) are ~269 or 7 K warmer than the point of closest approach on the flight line, which is offset from NEEM by ~17 km to the northeast. Estimated basal temperatures are near or at the melting point at the end of the flight line ~50 km north of NorthGRIP, where basal temperatures extrapolated from measurements in the deep borehole are ~270 K (Dahl Jensen and others, 2003), and where basal water has been detected based on radio-echo sounding data (Oswald and Gogineni, 2012). Note again that measured temperature profiles from Camp Century, NEEM and NorthGRIP locations were used in part to constrain the integrated firn transmissivity in the retrieval process, so that the basal temperatures retrieved are not completely independent from those of the borehole measurements. That said, the estimated basal temperatures tend to be cooler than the borehole basal temperature particularly near NEEM. The Robin model assumes a linearly varying vertical velocity with depth. Rezvanbehbahani and others (2019a) applies a non-linear vertical velocity profile to investigate the effect on modeled temperatures. The differences at any depth between the linear and curved velocity models are small (several hundredths of ma⁻¹) and the magnitude of the difference depends on the choice of model parameters including the exponent in Glen's flow law. Nevertheless, Rezvanbehbahani concludes that the linear model leads to an overestimation of the downward advection of cold ice and hence an underestimation of the basal temperature by several degrees Kelvin for an ice-sheet configuration similar to north central Greenland.

Fig. 5. Basal temperatures (black line) along the flight line shown in Figure 1. Estimated std dev. about the mean temperatures are shown in dashed red lines and are ~±3 K at basal depths. The black solid circles are the basal temperatures obtained from the Camp Century (left circle), NEEM (center circle) and NorthGRIP (right circle) boreholes. The abrupt increases in basal temperature at ~−58 and −51° longitude are associated with spikes in the lowest-frequency T _b channel and are likely associated with brief bursts of RFI that are not completely mitigated.

Ice rate factor

Temperature is a key factor controlling the deformation rate of the ice sheet. The relation between strain rate and deviatoric stress is usually taken as a power law first proposed by Glen (1958) and generalized by Nye (1953):

(4)$$\dot{\varepsilon }_{ij} = A\tau _{\rm e}^{n-1} \tau _{ij}, \;$$

where $ \dot{\varepsilon }_{ij} $ are the components of the strain rate tensor, A is the rate factor, τ _e is the effective deviatoric stress and τ_ij are the components of the deviatoric stress tensor. The flow-law exponent n is typically taken equal to 3. The effective stress is defined as

(5)$$2\tau _{\rm e}^2 = \tau _{xx}^2 + \tau _{yy}^2 + \tau _{zz}^2 + 2( {\tau_{xy}^2 + \tau_{xz}^2 + \tau_{yz}^2 } ).$$

Several investigators have reported on the functional form of the rate factor for polycrystalline ice; many of these are summarized in van der Veen (1999, Chapter 2.4) and Paterson (1994, p. 96). Here we adopt the equations proposed by Hooke (1981):

(6)$$A = A_0{\rm exp}\left({\displaystyle{{-Q} \over {{\rm RT}}} + {\displaystyle{{3C} \over {{( {T_{\rm r}-T} ) }^k}}}_\;} \right).$$

The constants are taken as A ₀ = 9.302 kPa⁻³ a⁻¹, Q = 78.8 kJ mol⁻¹, R is the gas constant (8.321 J (mol K)⁻¹), C = 0.16612 K^k, T _r = 273.39 K, k = 1.17 and T is the ice temperature in Kelvin. Figure 6 plots A versus ice temperature. Comparisons with Paterson's recommended values (Paterson, 1994, p. 97) are also included, and show reasonable agreement with those computed from Eqn (6) at moderate temperatures but diverge at lower temperatures. Additional estimates are summarized in van der Veen (1999, Fig. 2.6, p. 17) and vary within a factor of 2 at 258 K.

Fig. 6. Functional form of A using Hooke's relationship (black) and Paterson's recommended values (dashed).

Given these relationships, the UWBRAD pressure-corrected temperature field (a 1–2 K correction depending on ice thickness (Cuffey and Paterson, 2010, section 3.4.5.2)) and excluding corrected temperatures near the melting point, it is straightforward to calculate the rate factor A as a function of position and depth (Fig. 7) along the flight line and at selected locations (identical to those of Fig. 3) in Figure 8. At the surface, the rate factor decreases from ~2.4 × 10⁻⁹ kPa⁻³ a⁻¹ near Camp Century to ~9 × 10⁻¹⁰ kPa⁻³ a⁻¹ near NorthGRIP and corresponding to the cooling surface temperatures as the elevation along the transect line increased. Examining the vertical temperature profiles, A increases slowly with depth before more rapid increases begin in about the lower third of the ice column. The A parameter at the ice-sheet base is plotted along the study line in Figure 9, which also indicates the one std dev. uncertainty about the mean. The basal rate factor increases by almost an order of magnitude from ~1 × 10⁻⁸ kPa⁻³ a⁻¹ near Camp Century to ~8 × 10⁻⁸ kPa⁻³ a⁻¹ near NorthGRIP. The measurement gap at ~−43° longitude corresponds to locations where the pressure-corrected basal temperatures approach the melting point and so are ignored.

Fig. 7. Rate factor A (kPa⁻³ a⁻¹) inferred along the study line using the pressure-corrected physical temperatures and Hooke's functional relationship. The longitudes of Camp Century (C) and NEEM (NE) are indicated. The study line terminated west of NorthGRIP and that station is not shown (see Fig. 1). The color scale is logarithmic. White bars at top of figure show location of profiles presented in Figure 8.

Fig. 8. Rate factor A (kPa⁻³ a⁻¹) scaled by 10⁷ (horizontal axis) versus depth in m (vertical axis) for selected points along the flight line corresponding to the temperature profiles of Figure 3. Following the temperature gradients, the rate factors at the base of the ice sheet increase with distance along the flowline by about an order of magnitude.

Fig. 9. Rate factor A scaled by 10⁷ near the base of the ice sheet along the UWBRAD flight line. Rate factors calculated from physical temperatures are shown by the black line. Red dashed lines correspond to calculating A with one std dev. of uncertainty in the physical temperature. The gap at the easterly end of the profile corresponds to locations where the pressure-corrected temperature approaches the melting point and is considered unreliable for estimating A.

Temperature controls on ice flow

Ice flow within ~4 ice thicknesses away from an ice divide can be approximated using laminar flow theory (Raymond, 1983) for which the only important stress is the shearing across planes parallel to the bed. This condition is roughly satisfied along the portion of the flight line from Camp Century to NEEM along which the thickness varies from ~1300 to 2500 m. Under this model:

(7)$$\displaystyle{1 \over 2}\dot{\varepsilon }_{xz} = A\tau _{xz}^n , \;$$

which assumes that the driving stress is solely resisted by basal drag

(8)$$\tau _{xz} = \rho g{\rm sin}( \theta ) ( {H-z} ), \;$$

where θ is the surface slope, ρ is the density and g is the acceleration due to gravity. Because the rate factor varies with depth and position along the flight line, we calculate the horizontal velocity in the direction of the surface slope by integrating the strain rate (Eqn (7)) with depth. This speed (U(z)) is then given by

(9)$$U( z ) -U{\rm b} = 2 \times ( {\rho g{\rm sin}( \theta ) } ) ^n\mathop \int \limits_0^Z E( z ) A( z ) ( {H-z} ) ^n{\rm d}z, \;$$

where U _b is the basal sliding speed, ρ is the density, g is the acceleration due to gravity, H is the Operation IceBridge-derived ice thickness and n is a constant taken equal to 3. The slope θ is computed using the Greenland Ice Mapping Project (GIMP) digital elevation model (Howat and others, 2014, 2015) over 5 km intervals and corrected to elevations relative to the geoid (Forsberg, 2016). Finally, E is an enhancement factor that attempts to account for the effect of impurities, grain size and crystalline fabrics on ice hardness. Taking the sliding velocity to be zero and the enhancement factor equal to unity, selected velocity profiles along the study line are shown in Figure 10.

Fig. 10. Horizontal velocity profiles for selected sites along the dataset labeled with their longitudes. No enhancement factor is applied to the flow-law rate factor.

The surface velocities are all <2 ma⁻¹, or about a factor of 5 smaller than those measured with InSAR systems (see below). Consequently, a multiplicative enhancement factor is incorporated into Glen's Law to empirically account for the additional effects of grain size, crystalline fabric and impurities (Cuffey and Paterson, 2010, section 3.4.5.8). Layers composed of finer ice grains, as found in dust-particle-rich ice-age ice, tend to be softer (Cuffey and others, 2000). Fabrics (the relative orientation of crystalline c-axes for an ensemble of grains) develop by rotation and recrystallization that evolve as the ice is strained to favor basal glide (van der Veen and Whillans, 1994) that again leads to softer ice. Both grain size and fabric change down the ice column. The average grain size at Camp Century decreases by almost an order of magnitude at the transition from Holocene to Wisconsin age ice at ~1140 m depth or roughly the lower 18% of the ice column (Herron and Langway, 1982). There is a related change in fabric at the ice-age transition from a broad, single maximum to a strong, single maximum. Montagnat and others (2014) compare crystal fabrics at NEEM and NorthGRIP. At NorthGRIP, they find a progressive strengthening of the fabric below 2300 m depth or about the lower 25% of the ice column and well below the ice-age transition (~1500 m depth). In contrast at NEEM, there is an abrupt change in fabric between 1400 and 1500 m depth or about the lower 40% of the ice column and corresponding to the ice-age transition. Montagnat and others (2014) go on to report a strong c-axis clustering in the range of 1800–2000 m depth or about the lower 25% of the ice column.

Based on the ice core properties described above and minimizing the root-mean-square deviation between the measured surface and InSAR velocities, the enhancement factor (E(z)) in Eqn (9) is taken to be:

$$E( z ) = 1\;{\rm for}\;z < 0.75\;H, \;$$

$$E( z ) = 8\;{\rm for}\;z > 0.75\;H.$$

Horizontal velocity profiles using this ‘softer’ basal ice are shown at selected sites in Figure 11. Also shown is U(z) at Camp Century as reported in Figure 7 of Gundestrup and others (1993) and located ~1 km downstream of the nearest airborne data point. Gundestrup and others (1993) computed the velocities using borehole deformation data and integrating the horizontal motion gradients up the ice column. These borehole data show that the largest velocity gradients corresponding to 85% of the deformation occur in the lower 250 m of the ice column corresponding to the glacial stage ice. The borehole-deformation surface velocities are ~0.3 ma⁻¹ faster than those predicted by the laminar flow model, but the latter are highly dependent on the choice of enhancement factor which here is selected to provide a reasonable fit to all of the InSAR velocities.

Fig. 11. Laminar-flow horizontal velocities with enhancement factor applied to the lower part of the ice sheet. Dashed line is from Figure 7 of Gundestrup and others (1993) for Camp Century for comparison with the laminar flow velocities located at the point of the closest approach of the flight line.

For a comparison of enhancement factors, Shoji and Langway (1985) found that the enhancement factor E at DYE-3 increased irregularly from ~1 to 5 between 500 and 1500 m depth, with values between 5 and 17 within the lower 250 m glacial ice. Cuffey and Paterson (Table 3.5, p. 77) propose values ranging from 1 to 12 for cases ranging from polar ice shelves to debris laden basal ice. Moreover, the presence of folded, Eemian-age ice in the lower 200 m of ice at NEEM (Dahl Jensen and others, 2013) introduces additional complications including debates about the proper form of the ice flow law. Kuiper and others (2020) argue that a modified, composite flow law better describes the deformation of stratigraphic features within the lower, Eemian-age layers, and contends that strain rates between folded fine and coarse grain layers can vary by about an order of magnitude. This contrasts with the Glen model which is independent of grain size when the enhancement factor is neglected.

Figure 12 compares surface velocities obtained using the proposed enhancement factor to those from two-dimensional InSAR velocity data (Joughin and others, 2010, 2015). The overall trends are in general agreement. Some of the rapid changes in UWBRAD speeds are partly due to errors in consecutive UWBRAD data. Because the estimation is done sequentially and uses results at the previous point as an approximation to the prior in a Bayesian sense, this affects the results at the current step and beyond. Other along track speed variations in both the InSAR and UWBRAD velocities are driven by variations in surface slope, ice thickness and the hardness. On the one hand, the latter two parameters fall within the velocity integral (Eqn (9)) and so do not provide a convenient normalization for removing known behaviors. On the other hand, normalizing the speed by the slope is practical and the result shown in Figure 13 where the normalization factor is (sin(θ))³. Figure 13 illustrates good agreement between the two speed estimates up till ~−53° longitude. There as in Figure 12, the UWBRAD speeds and normalized speeds exceed the InSAR estimates. Larger changes in basal topography and basal slope occur at this location (Fig. 14), which may indicate that neglected longitudinal and/or lateral stresses are important. That said, it is not possible to address reasons for the difference as other than speculations absent additional independent information on the ice physical properties.

Fig. 12. Laminar surface velocities (blue line) using the depth-dependent enhancement factor compared to InSAR surface velocities (red dashed line) (Joughin and others, 2015). Estimated InSAR errors vary between 0.5 and 1 ma⁻¹ along the profile. Black error bars on laminar model velocities span the associated uncertainty range of the physical temperatures.

Fig. 13. Slope normalized speed along the profile line. The normalization factor is (sin(θ))³. Blue line is the normalized UWBRAD speed; red line is the normalized InSAR speed.

Fig. 14. Retrieved geothermal heat flux (blue line) with error bars computed as part of the parameter retrieval. Black line is the basal topography. Red line is GHF from Martos and others (2018). Black circles are average GHF values at the Camp Century, NEEM and NorthGRIP borehole sites taken from Martos and others (2018).

A more complete description of the force balance equations is required beyond NEEM where the flight line is much nearer to the ice divide (van der Veen, 1999, p. 45). Numerical solutions to the necessary equations are beyond the scope of this analysis. However, the stronger variations in A(z) along the path from NEEM to NorthGRIP (Fig. 9) and amplified by enhancements associated with grain size, shape and orientation suggest increasingly important rheological controls on the strain rate field and hence the velocity field in this region.

Geothermal heat flux estimates

An estimate of GHF is obtained in the Bayesian inversion of the brightness temperature spectral data as discussed by Yardim. This inversion is independent of geological constraints saved for the influence of basal topography on ice thickness. Because the inversion is based on the use of the Robin model, sources of heat within the ice sheet are neglected. Hence it is assumed that strain heating is included as part of the estimated GHF, given the concentration of strain near the ice-sheet base.

GHF values retrieved along the flight line are shown in Figure 14. These values are ~55 mW m⁻² from Camp Century to ~−55^° longitude. GHF increases on average along the remainder of the transect that more closely follows the ice divide and reaching ~85 mW m⁻² at the end of the dataset. Bedrock topography from BedMachine is also shown (Morlighem and others, 2017, 2021), which on average negatively correlates with GHF.

Figure 14 also includes GHF estimates by Martos and others (2018 and 2018 supplement and hereafter referred to as Martos) for comparison. Martos applied spectral analysis to the World Digital Magnetic Anomaly Map 2 to estimate Curie depths on a 15 km grid. They then used solutions to the heat conduction equation to compute GHF. They assume that the temperature at the Curie depth is 580°C, that the temperature at the ice-rock interface is 0°C, and that the thermal conductivity and radioactive heat production at the surface are consistent with the Precambrian geology of Greenland, and further refine their model based on comparisons with Greenland borehole temperature data. Compared to earlier heat flux maps, Martos report a reduced amplitude and range of heat flux variations but consider their range of GHF values to be consistent with the cratonic nature of the continent. Rezvanbehbahani and others (2019b) compared GHF estimates based on radar-identified regions of basal meltwater with those of Martos and suggest that additional, more localized mechanisms not included in standard GHF models are important. These include modified heat fluxes near topographic edges (van der Veen and others, 2007; Colgan and others, 2021), interplay between glacial loading and the underlying crust, and discharges from groundwater aquifers.

The eastward generally increasing trend in GHF is similar for both the Martos and Yardim data although the Yardim data increase at a faster rate. The microwave-derived estimates are ~6 mW m⁻² higher than those of Martos near Camp Century, similar at NEEM, and show an increasing divergence farther along the flight line reaching a maximum offset ~24 mW m⁻² at the end of the profile. Yardim GHF estimates show much greater variability. Some of the GHF fluctuations, as with the physical temperatures, are due to the retrieval algorithm recovering from poor brightness temperature data that can be caused by brief periods of RFI associated with other instruments and radios on the platform. There is a possible exception located at −44.6^o longitude where the transect intercepts a portion of the 100s of km long valley that extends northwards toward Petermann Glacier (80.2^oN, 58^oW). At that point, the relative topography drops ~250 m and the GHF decreases by ~10 mW m⁻². Topography controls local GHF though modification of isotherms by sloping walls (van der Veen and others, 2007) resulting in local maxima and minima in GHF (Fig. 1 of Colgan and others, 2021). That said, it is difficult to attribute the anomaly to either physical mechanisms or to data uncertainty without additional modeling and geologic control.

Martos also summarized GHF at the Camp Century, NEEM and NorthGRIP deep drilling sites (see their Table S1 of supporting information) and concluded that local GHF values are between 41 and 50 mW m⁻² at Camp Century, near 58 mW m⁻² at NEEM, and between 55 and 160 mW m⁻² at NorthGRIP. Average values for these sites are shown as circles in Figure 14. The borehole site estimates at Camp Century are lower than both the Yardim and Martos data, similar to both at NEEM, and significantly higher at NorthGRIP where the Yardim data fall about halfway between the Martos and the averaged borehole site values.

Martos suggest that the wide variation of GHF estimates near NorthGRIP could be due to the presence of water at the bed which is an unaccounted-for latent heat contribution. Microwave-derived GHF may also be overestimated because strain heating at the base of the ice sheet is incorporated into the GHF parameter in Robin's formulation. Given the laminar flow model assumption that τ _xz is the only operative stress, strain heating over the ice column is given by

(10)$${\rm SH} = \mathop \int \limits_0^H \dot{\varepsilon }_{xz}\tau _{xz}{\rm d}z.$$

Assigning most of the strain heating to the bottom of the ice, a correction to the GHF can be estimated by subtracting the integrated strain heating calculated using Eqn (10) from the Yardim GHF as shown in Figure 15, for which uncertainties are also plotted based on the range of rate factor values. Between Camp Century and NEEM, there is now a pronounced local minimum in the Yardim GHF. The minimum also is present in the Martos data but relative to the westerly end point the minimum is only ~1 mW m⁻² for the Martos data but almost 10 mW m⁻² for the Yardim data. The difference may in part be attributed to the average 1.2 km sampling of the airborne data as compared to the 15 km Curie depth grid used by Martos that could in turn result in attenuated GHF variations.

Fig. 15. Blue line is the GHF estimated by Yardim; red line is from Martos. The black line is the strain heating corrected GHF and the gray circles are average GHF values at borehole sites. The error bars reflect the maximum and minimum values of the GHF using the bounds of the rate factor at selected sites.

Summary

Airborne measurements of Greenland ice sheet internal temperatures were used to calculate the rate factor in Glen's flow law. The rate factors at the base increase by about an order of magnitude along the flight line. The largest values are found northwest of the NorthGRIP deep drilling site where basal temperatures approach the melting point.

Horizontal velocities for measurements ~4 ice thicknesses from the ice divide were computed as a function of depth using the estimated rate factor and a simple laminar flow model. Modelled surface velocities were found to be smaller by about a factor of 8 from surface velocities derived from satellite InSAR data. Consequently, an enhancement factor was introduced that improved agreement between predicted and InSAR observed surface velocities between Camp Century and NEEM. The estimated enhancement factor, which empirically accounts for the additional influence of grain size and crystalline fabric on ice deformation, was applied to the lower 25% of the ice sheet. The depth corresponds roughly to the transition where grain size reduces and fabrics cluster to a single maximum of crystalline c-axis orientations. Because the flux of ice across a vertical plane within the ice is proportional to the depth-integrated velocity, mass continuity calculations are improved with better knowledge of the rate and enhancement factors.

GHF is one of the controlling parameters in glacier thermodynamics as expressed in Eqn (3). Yardim GHF increases along the transect and the trend compares favorably with other geophysical estimates reported in the literature. However, the magnitude of the Yardim data is offset compared to magnetically constrained estimates by Martos. To account for the difference, a strain heating correction of ~15 mW m⁻² was subtracted from the Yardim GHF between Camp Century and NEEM. The corrected GHF emphasizes a local minimum in the Yardim data and that is more pronounced than in the Martos data possibly due to differences in sampling intervals.

The results presented here illustrate the potential of remotely sensed ice-sheet temperature profiles for glaciological studies. Broader applicability awaits several developments. First and foremost, the analysis must go beyond the parameters estimated from the restrictive Robin ice-sheet temperature model. To that end, a more complete thermodynamic model (e.g. Rezvanbehbahani and others, 2019a) that captures the different glaciological regimes away from the ice divide is needed to initialize the retrieval process. Alternatively, Eqns (1) and (2) explicitly relate the physical temperature to the brightness temperature. In principle, direct inversion of this equation bypasses entirely the need for a thermal model. However, experiments to date indicate that methods must be found to reduce errors on the computed temperature.

Additional improvements in temperature sensing may be possible through combined active and passive microwave observations as well. As shown by MacGregor and others (2015), radar sounding data can be used to estimate average temperature with depth. Oswald and Gogineni (2012) and Oswald and others (2018) show how radar reflectivity at the base of the ice can be used to infer locations where the base is at the pressure melting point. These observations can constitute another important, independent constraint on estimated temperature profiles. Equally important, Xu and others (2020) show how the integrated effect of density fluctuations on the transmissivity through the upper layers of the ice sheet can be measured with an ultra-wideband radar operating in the same frequency range as the radiometer. The result can be used to compensate the brightness temperature for firn transmissivity without resorting to in situ observations. The importance of temperature information in predicting future ice-sheet evolution will continue to motivate research in these areas.