Geothermal heat flux is the dominant source of uncertainty in englacial-temperature-based dating of ice-rise formation

. Ice rises are areas of locally grounded, slow-moving ice adjacent to floating ice shelves. Temperature profiles measured through ice rises contain information regarding changes to their dynamic evolution and external forcings, such as 10 past surface temperatures, past accumulation rates and geothermal heat flux. While previous work has used borehole temperature-depth measurements to infer one or two such parameters, there has been no systematic investigation of parameter sensitivity to the interplay of multiple external forcings and dynamic changes. A one-dimensional vertical heat flow forward model developed here examines how changing forcings affect temperature profiles. Further, using both synthetic data and previous measurements from the Crary Ice Rise in Antarctica, we use our model in a Markov Chain Monte-Carlo inversion to 15 demonstrate that this method has potential as a useful dating technique that can be implemented at ice rises across Antarctica. However, we also highlight the non-uniqueness of previous ice rise formation dating based on temperature profiles, showing that using nominal values for forcing parameters, without taking into account their realistic uncertainties, can lead to underestimation of dating uncertainty. In particular, geothermal heat flux represents the dominant source of uncertainty in ice-rise age estimation. For instance, in Crary Ice Rise higher heat flux values (i.e., about 90 mW m -2 ) yield grounding timing of 20 1400±800 years, whereas lower heat flux of around 60 mW m -2 implies earlier ice rise formation and lower uncertainties in the ice rise age estimations (500±250 years). We discuss the utility of this method in choosing future ice drilling sites and conclude that integrating this technique with other indirect dating methods can provide useful constraints on past forcings and changing boundary conditions from in-situ temperature-depth measurements.


25
Present-day englacial temperatures are the product of the millennial-scale histories of ice flow and thermal boundary conditions experienced by an ice sheet (Robin, 1955). Temperature measurements from boreholes drilled through ice sheets have been widely used to extract important paleoclimatic archives, such as surface temperature and accumulation history, as well as information about the conditions at the base of an ice sheet (i.e., glacial thermal regime and geothermal heat flux), both in Antarctica, the ice sheet contains ice rises -regions of slow-flowing, locally elevated, grounded ice embedded within or adjacent to fast-flowing, floating ice shelves; one way they form is through an ice shelf grounding on marine bed (e.g., Martin and Sanderson, 1980;Matsuoka et al., 2015;Wearing and Kingslake, 2020). This shift in boundary condition at the base of an ice-sheet results in a transient evolution of the temperature-depth profile within an ice rise (MacAyeal and Thomas, 1980;Bindschadler et al., 1990). Therefore, due to their proximity to the marine ice-sheet periphery and negligible horizontal flow, 35 ice rises can retain an imprint of past grounding line migration on millennial timescales, a record that is otherwise largely inaccessible beneath the ice sheet or its fringing ice shelves (e.g., Conway et al., 1999;Matsouka et al., 2015;Kingslake et al., 2018;Neuhaus et al., 2020).
Past work used temperature-depth measurements within ice rises in Antarctica to estimate the timing of ice-shelf grounding.
For instance, Bindschadler et al (1990) developed an advection-diffusion thermal model of Crary Ice Rise, West Antarctica 40 ( Fig. 1). The model calculated the initial steady-state ice-shelf temperature profile from the specified set of parameters, including ice thickness, surface temperatures and accumulation rates. The steady-state ice-shelf profile was perturbed using thermal properties of the bed (geothermal heat flux, diffusivity and conductivity of the bedrock) and a specified vertical ice velocity function to calculate transient thermal evolution after ice-shelf grounding. The modelled temperature profiles were then compared to in-situ borehole measurements at Crary Ice Rise, and minimizing the mismatch between the measured and 45 synthetic profiles yielded the best age estimate of the ice rise in its thickest part to be 1100 years. Recent work by Neuhaus et al (2020) built upon the model of Bindschadler et al (1990) to evaluate the timing of grounding at three sites in the Ross Sea sector of Antarctica, where previous measurements showed anomalously high basal temperature gradients (Engelhardt, 2004a).
These results largely corroborate hypotheses of late-Holocene re-advance in the region (Kingslake et al., 2019) and associated grounding at these sites between 1100 and 500 years ago (Neuhaus et al., 2020).

50
Thus, the methods used in these studies have potential if future boreholes are drilled at Antarctic ice rises in locations suspected of undergoing significant dynamics change. Yet, the uncertainties inherent in these approaches must be carefully assessed to target drilling, maximize the utility of borehole drilling, and increase the accuracy of ice-dynamic and paleoclimatic inferences.
Previous work has included sensitivity tests where some predefined variables, such as accumulation, ice thickening, and melt rate, were assigned several different values to examine how they affected the final temperature profile, and their relation to 55 inferred timing of grounding (Bindschadler et al., 1990;Neuhaus et al., 2020). Yet, there is a lack of systematic investigations of temperature profile sensitivity to the cumulative effects of multiple external forcings and dynamic changes, particularly given that some parameters (e.g., geothermal heat flux) have considerable uncertainties (e.g., Fudge et al., 2019). In addition, time-variable parameters, such as ice thickness, accumulation and surface temperature may significantly increase the dimensionality of the problem, solutions to which need optimized inversion methods, as opposed to exhaustive global search algorithms where highly dimensional inversion tasks become computationally unfeasible (Mosegaard and Tarantola, 1995). accumulation rate, heat flux and thickness history) affect the englacial temperatures. Using previous measurements from the Crary Ice Rise, West Antarctica ( Fig. 1), we also implement a Markov Chain Monte-Carlo inversion to explore the contributions of multiple uncertainties to the inferred timing of ice-rise formation. Arrow on (a) shows the location of Crary Ice Rise, where previously timing of grounding was inferred using borehole temperature measurements (Bindschadler, 1990

80
In this section we outline the numerical forward model and inversion method. The forward model builds upon and extends the models used by Alley and Koci (1990), Orsi et al (2012) and Neuhaus et al (2020). The inversion method is similar to that previously used to reconstruct past surface temperatures in Greenland and Antarctica (e.g., Dahl-Jensen et al., 1998;. Extending previous work, we include the effects of different vertical velocity functions, we use an optimised set of some parameters (e.g., pressure melting/freezing point of ice/seawater), and we introduce temporal variability to ice thickness,

85
surface temperature and accumulation rates, as well as multiple phases of grounding and ungrounding and corresponding changes to the boundary conditions.

Model equations
We use the following form of the vertical diffusion-advection equation to simulate the time-and depth-evolution of temperature 90 T: where t is time, z is height above the bed, w is vertical ice velocity (positive upwards) and α is thermal diffusivity: Where ρ is density, k is the thermal conductivity, c is the specific heat capacity. We assume that α is uniform and constant, and 95 the internal heat production can be neglected due to insignificant horizontal ice flow (Dahl-Jensen et al., 1999), typically on the scale of metres per year (e.g., Matsouka et al., 2015;Kingslake et al., 2016).
For sensitivity experiments, we implement three analytical approximations for vertical velocity w within grounded ice. The Dansgaard and Johnsen (1969) vertical flow approximation can be formulated as: where ws is time-varying vertical velocity at the surface, z = 0 represents the ice-bedrock interface, H is the ice thickness and zk is a free parameter (Martín & Gudmundsson, 2012). Lliboutry (1979) provides another commonly used approximation (Wearing & Kingslake, 2019, Fudge et al., 2019: where n is a rheological parameter (Glen's flow law exponent, n=3).
Finally, a simplified approximation where vertical velocity varies linearly with depth from zero at the base to a maximum value at the surface, following Bindschadler et al (1990), was implemented for sensitivity experiments and estimation of associated uncertainties: 110 When time-varying ice thickness is introduced, the difference between accumulation and thickening rates determines the vertical velocity at the surface: For temperature calculations within the underlying bedrock, no vertical advection is assumed (w=0), reducing Eq. 1 to only account for heat diffusion.

Spatiotemporal domain and boundary conditions
We define a one-dimensional spatial domain that extends vertically from the bedrock base to the ice surface. The vertical coordinate is z, which increases upwards, and zb and zs are the elevations of the ice base and ice surface, respectively. We implement a simple finite-difference scheme to solve Eq. 1. We discretize the spatial domain with a minimum of 50 nodes within the ice column and 50 nodes with fixed 10 m vertical spacing in the subjacent bedrock (Neuhaus et al., 2020). We used In the ice-shelf scenario (Fig. 2a) boundary conditions at the ice surface are set by time-variable surface temperatures T(t) and accumulation rates a(t). The temperature at the top of bedrock/base of ice column T(z=zb) is forced to equal the freezing point 125 of seawater (calculated as a function of ice thickness, to account for its pressure dependence), and the temperature gradient in the bedrock is calculated based on the geothermal heat flux G (Millero, 1978;Determann and Gerdes, 1994). We assume that the rates of accumulation and basal melting balance the depth-uniform vertical velocity (Holland and Jenkins, 1999).
When modelling temperatures within a grounded ice sheet (i.e., an ice rise; Fig. 2b), the boundary condition at the base of the ice column is set such that the vertical gradient in T corresponds to a diffusive heat flux (there is no advection because w(zb) = 130 0) that balances the geothermal heat flux G through the underlying bedrock. Thickness history H(t) and the vertical velocity function (Eqs. 3-6) within the ice rise are also used. The boundary conditions at the top of the ice column are prescribed based on the temperature and accumulation rate histories; prior values of parameters that dictate these boundary conditions are discussed further in Section 2.1.3.

Input parameters
Parameters characterising the properties of ice and the bedrock are assumed to be constant ( Table 1). Values of geothermal heat flux at specific locations are sampled from the Antarctic-wide heat flux data compilation derived from spectral analysis of airborne magnetic data (Martos et al., 2017). Surface temperature evolution, accumulation rate histories and ice thickness histories are sampled from distributions provided by an ensemble of simulations using the Parallel Ice Sheet Model (PISM) Kingslake et al., 2018;Albrecht et al., 2020aAlbrecht et al., , 2020b. PISM is a three-dimensional, thermomechanical ice-sheet model that solves a hybrid shallow approximation of Stokes flow. PISM produces continental-scale (16 km grid cell size), long-term (multi-millennia length) simulations of ice-sheet thickness. Associated with these fields are surface-temperatures reconstructed from the West Antarctic Ice Sheet (WAIS) divide ice-core , scaled by modelled ice-surface elevation and accumulation patterns simulated by the regional climate model RACMOv2.1 (Ligtenberg et al., 2013), scaled 150 by 2% per degree of temperature change from present (Kingslake et al., 2018). Other specifications and parameters used for the chosen PISM model output, including mantle viscosity and flexural rigidity are described in detail in Kingslake et al. (2018). Despite inherent uncertainties associated with large-scale numerical models, which have relatively coarse resolution compared to the length-scales of ice rises, PISM outputs provide first-order estimates of prior information about the timevariable parameters (i.e., temperature, accumulation and ice thickness) at any specified location on the Antarctic Ice Sheet

Experimental design and computational environment
Prior to imposing time-variable forcings, we allow the simulation to reach a steady state using surface temperatures, thicknesses, and accumulation rates equal to their values at the beginning of the chosen simulation period ( then compared to examine how altering our prior parameters affected the final temperature profiles. In order to evaluate the fit between two temperature profiles, we used a root-mean squared mismatch, where and are the measured and predicted temperatures at grid point i, and n is the number of grid points. To increase computation efficiency in our forward simulations where Monte-Carlo approach is used (i.e., where parameters 180 are repeatedly randomly sampled from a range of their prior distributions), we use Dask, a flexible library for parallel computing in Python, which is implemented using cloud-based computing clusters provided by the Pangeo project (Odaka et al., 2019). Pangeo is a developing ecosystem of scalable, open-source tools for cloud-based parallel computation and interactive large-scale computation and data analysis. Automatic parallelization on tens to hundreds of workers significantly increases computation performance, as compared to a standard approach using a desktop computer, and is particularly 185 applicable to tasks that are easily parallelized, such as forward Monte-Carlo simulations.

Inverse model
Where borehole temperature-depth measurements are available, they can be used to infer the history of dynamic change and evolution of boundary conditions (e.g., Dahl-Jensen et al., 1999) via inversion (e.g., Mosegaard and Tarantola, 1995). In this section, we outline what observational data and inversion methods we use to infer past ice-rise evolution.

Temperature-depth data
Among numerous ice rises mapped across the AIS only a few have been sampled by borehole temperature measurements (Matsouka et al., 2015). These sites include Siple Dome and Crary Ice Rise of the Ross Ice Shelf, as well as Mill Island of the Shackleton Ice Shelf in East Antarctica (e.g., Koci and Bindschadler, 1989;Engelhardt et al., 2004b;MacGregor et al., 2007;Roberts et al., 2013). We digitized Bindschadler et al.'s (1990) temperature observations from Crary Ice Rise and used them 195 as input data for our Bayesian inversion method (i.e., probabilistic data analysis that involves using the prior information and computing the posterior probability distribution for the parameters of the model; Section 2.2.2). Due to the quality of available data, digitisation implies inherent uncertainties (up to 0.1ºC RMS difference between independently digitised temperature profiles). Furthermore, as uncertainties associated with measurement calibrations can reach up to 0.1°C (e.g., Orsi et al., 2012), we consider 0.1ºC as an upper threshold value for estimating the degree of fit between measured and predicted temperature 200 profiles. In addition, we used the forward model outlined in Section 2.1 with predetermined parameters to produce synthetic temperature-depth data. These synthetic profiles were then used as an input for validation of our inverse method.

Markov Chain Monte Carlo inverse method
The Markov Chain Monte Carlo (MCMC) method tests randomly selected combinations of prior variables using a random walk through a high-dimensional parameter space. The variables are assumed to be 'unknown' parameters and are prescribed prior distributions (or simply a range of realistic prior values). The forward model uses these parameters from their prior distributions as inputs, and its output is compared to the measured (or synthetic) temperature profiles. In each step of random walk with predefined length, a perturbed model of the current model is proposed. The MCMC then uses a likelihood function to estimate the agreement between the modelled and measured profiles (Mosegaard and Tarantola, 1995;Dahl-Jensen et al., 1998): where and are the measured and predicted temperatures for each point of the grid, and Terror is the uncertainty of measurements.
The perturbed model is either 'rejected', in which case a new random perturbation is applied from the same starting location in parameter space, or 'accepted', in which case this location becomes the new starting location for the next random 215 perturbation. Even when a new perturbation yields better fit to the observations, the model may be rejected, or if the new perturbation produces larger misfit, the model may be accepted. Whether a model is accepted is based on an acceptance probability (see Dahl-Jensen et al., 1999), which introduces a degree of stochasticity, ensuring that the random walk avoids entrapment in local likelihood maxima (Mosegaard and Tarantola, 1995). Eventually, the paths converge towards the regions corresponding to model parameters that yield the lowest misfit between modelled and observation temperature profiles.

220
In practice, MCMC can be used to infer any inputs to the forward model. For instance, in the Dye 3 borehole drilled through the 3-km thick Greenland Ice Sheet, the unknown temperature history has been divided into 125 intervals, which, together with heat flux (also assumed to be unknown) yielded a 126-dimensional parameter space (Dahl-Jensen et al., 1999). Here, we focus on inferring past dynamic changes (i.e., timing of grounding), but also assume heat flux, temperature, accumulation, and thickening rates to be unknown in order to explore the possible combinations of realistic parameter values that could provide 225 a close fit between the field observations and modelled temperature-depth profiles. To perform the MCMC we used Emcee, a Python-based software package which implements the affine-invariant ensemble sampler (Goodman & Weare, 2010;Foreman-Mackey et al., 2014). Prior to using actual borehole measurements, we tested the MCMC inversion method on a synthetic temperature profile calculated using a simplified forward model that uses the Lliboutry vertical velocity function and assumes a hypothetical 1000m thick, 1 k.y. old ice rise (i.e., formed when an ice shelf grounded 1000 years ago), forced by a heat flux of 50 mW m -2 , a constant surface temperature of -30º C and an accumulation of 0.1 m yr -1 . This experiment also allows us to evaluate some 235 aspects of temperature-depth variability that may result from changing forcing parameters (Fig. 3). The variables are prescribed the widest range of possible values as their prior distributions. From sampling of 100,000 different combinations of these variables from their prior range, over 9,000 yield an RMS mismatch with the synthetic temperature profile of less than 0.1ºC, with best fit samples (RMS ~ 0.01ºC) yielding parameters that very closely match the parameters used to produce the synthetic profile (Fig. 3). Nevertheless, even in this idealised setup, a wide range of ice-rise age estimations (mean±standard deviation 240 of ~1200±580 years,) closely match the synthetic profile (i.e., RMS <0.1ºC). Moreover, while derived accumulation rates and surface temperatures are tightly clustered (i.e., within a few percent) of their prescribed values, heat flux values are more widely distributed 50±20 mW m -2 (Fig. 3d). In addition, we ran a synthetic inversion while making the assumption that thermal diffusivity of ice and bedrock are unknown (James, 1968;Fuchs et al., 2021). Our results show that closely matching temperature profiles can be reproduced using a wide range of thermal diffusivity values, although their effect on inferred 245 grounding time is negligible (Fig. 1 in the Supplement). In summary, these synthetic tests show that our inverse approach effectively recovers surface parameters (accumulation and temperature), but uncertainty inherent in the system may introduce unavoidable uncertainty in estimated geothermal heat fluxes and ice-rise grounding times.

Crary Ice Rise borehole measurements
To implement our MCMC inversion with englacial temperature data from Crary Ice Rise, we prescribe a wide range of prior 260 values for model parameters. This allows us to probabilistically evaluate the possible combinations of realistic parameter values that could closely match the measured temperature-depth profiles (Fig. 4). For example, instead of assigning a heat flux of 77 mW m -2 , as used by Bindschadler et al (1990), The distributions of the colored points in the scatter plots in Fig. 4 provide insight into the dependence of each parameter and its uncertainty on the other parameters. For example, as with synthetic data inversion (Section 3.2.1), the inferred age of ice rise formation is strongly dependent on heat flux (Fig. 4c). Higher heat flux (i.e., about 90 mW m -2 ) yields grounding timing of 1400±800 years, whereas lower heat flux of around 60 mW m -2 implies earlier ice rise formation and lower uncertainties in 275 the ice rise age estimations (500±250 years). For comparison, using their deterministic least square approach with fixed parameter values, Bindschadler et al (1990) estimated an age of ice rise formation of 1100 years with a standard error between measured and simulated temperature profiles of 0.08ºC. Our best fit (i.e., RMS of 0.1ºC) temperature profile is inferred when the timing of grounding is 1300 years (with heat flux of 86 mW m -2 , accumulation of 0.14 m yr -1 and initial thickness of 460 m being the estimated best fit values). Here, we focus on the interplay of all parameters and do not assess the inferred initial 280 thickness required to explain grounding at a particular Crary Ice Rise site. However, more extra information onpast sea level, as well as the history of bedrock isostatic fluctuations could be easily incorporated into our model to improve initial thickness estimation. Our MCMC inversion experiments highlight that the inferring timing of ice-rise formation from englacial temperature measurements requires accurate knowledge of heat flux. However, Figures 3 and 4 show that the inverse approach does not in full illustrate the sensitivity of the temperature profile to model parameters, and we discuss this in Section 3.2.

295
To explore temperature profile sensitivity to model parameters and to identify which parameters are particularly important for ice-rise dating and borehole thermometry, a series of forward simulations (Figs. 5a-c, Fig. 6 increases the sensitivity of basal temperature to surface temperatures (Fig. 5a). Similarly, Figure 5E shows that the temperature profiles are much more sensitive to grounding time for younger ice rises. Finally, changing temperature-dependant thermal conductivity value from 2 W m -1 K -1 (for ice close to the pressure melting point) to 2.6 W m -1 K -1 (for ice at -30ºC) yields up to 0.3ºC RMS difference between the resultant temperature profiles for ice that is 500 m thick. We also investigate how different velocity functions and temporal variability of thickness, accumulation and surface temperatures may affect the distribution of temperature within grounded ice of various thicknesses (Fig. 6). For these 320 simulations, variable temperature and accumulation for the last 40 k.y. were obtained from PISM simulation outputs (Fig. 1g-i). Advection effects of vertical velocity approximations on englacial temperature profiles are illustrated in Figure 6a, d.
Temperatures are affected most in the lower part of the ice column (Fig. 6a). Deviations between profiles produced with Dansgaard and Johnsen (1969) and Lliboutry (1979) velocity approximations also increase with ice thickness. For example, in 500-m-thick ice the RMS difference between two profiles is 0.18ºC, while in 1800-m-thick this difference is 2.55 ºC (Fig.  6d).
Inversion experiments described in Section 3.1 assumed constant surface temperatures, Ts, and accumulation rates, a, over the simulation period. To examine how temporal variability of these parameters impacts the temperature-depth data, we ran forward simulations and compared resultant profiles for four scenarios (Ts and a constant, Ts constant and a variable in time, a constant and Ts variable in time, both variable in time; Fig. 6b). These experiments reveal considerable deviations between 330 profiles forced only by one time-variable parameter, and while these effects typically do not exceed 1ºC RMS for ice thinner than 1 km, they result in as much as few degrees RMS difference between two profiles within ice about 3 km thick (Figs. 6b,e). The effect of the ice thickening rate is most significant in relatively thin ice (up to 2ºC RMS change for ice 500 m thick, Figs. 6c,f), but decrease rapidly as ice becomes thicker.
Overall, forward sensitivity experiments outlined in this section provide insight into what model parameters have the strongest 335 effect on the ice column, how these effects manifest in the shape of temperature profile, and how this is affected by ice thickness. For thick ice, more likely to be found in the ice-sheet interior, parameters such as accumulation rates and vertical velocity profile stronger affect temperature profile, compared to thinner ice, more characteristic for ice-sheet periphery, including ice rises (Fig. 6). These results highlight the importance of careful vertical velocity parametrization for borehole thermometry experiments in thick ice-sheet settings. In contrast, knowledge of ice-sheet thickening history is more important for dating shallow ice rises (Fig. 6c).

Inversion of ice-rise age
Bayesian inversion of englacial temperature-depth profiles indicates that inferences of ice rise age (i.e., timing of grounding) may significantly vary depending on the values of other forcing parameters. Among these, geothermal heat flux may have a 355 significant effect on the inferred timing of grounding, with lower heat flux yielding earlier age estimates, along with much smaller corresponding uncertainties (that is, narrower range of solutions acceptable within a prescribed degree of misfit, Figs. 3, 4). In case of the Crary Ice Rise inversion experiment, we infer a range of possible ice rise age that encompasses the value of 1100 years previously reported by Bindschadler et al (1990), but increases from 500±250 to 1400±800 years ago, as the corresponding values of heat flux decrease (assumed here to vary between 60 mW m -2 and 90 mW m -2 ).
The effect of heat flux may be significant when prior knowledge about its values is poor or unconstrained. For instance, in the synthetic data inversion experiment with prescribed grounding timing of 1000 years, heat flux of 30 mW m -2 yields ice rise age estimation of 520±80 years, whereas a value of 80 mW m -2 corresponds to ice grounding 4030±610 years ago, with both combinations falling within misfit of 0.1ºC RMS (Figs. 3, 7). This result is important because previous deterministic approaches, for example the one implemented by Bindscahdler et al (1990), assumed a fixed value of heat flux. However, 365 recent studies have shown that despite its importance to understanding ice sheet evolution (Pollard et al., 2005;Seroussi et al., 2017), heat flux beneath ice sheets remains relatively poorly constrained, with discrepancies between continental-scale reconstructions and targeted heat flux measurements, and associated uncertainties of up to ±15 mW m -2 (e.g., An et al., 2015;Martos et al., 2017;Fudge et al., 2019). Furthermore, heat flux at the base of an ice sheet may show substantial lateral variations (i.e., up to 200 mW m -2 ) over relatively short distances of less than 100 km, as for example observed below the Whillans Ice

370
Stream (Begeman et al., 2017), as well as in other locations beneath the Antarctic and Greenland ice sheets (e.g., Cuffey et al., 1995;Dahl-Jensen et al., 2003;Schroeder et al., 2014;White-Gaynor et al., 2019). Careful estimation of heat flux uncertainties and their incorporation in temperature-depth inverse models are therefore essential for accurate estimation of past dynamic changes to ice rises.
The uncertainties in ice rise age estimates are also determined by the degree of accepted misfit between predicted and measured 375 temperature profiles, which in turn relies on the accuracy of englacial borehole thermometry. Depending on the instrumentation used, previous temperature-depth data have been collected with typical uncertainties of around 0.1ºC (Orsi et al., 2012;Yang et al., 2017, Fudge et al 2019, and occasionally up to 0.05ºC and even 0.02ºC, as demonstrated by work on the Antarctic Ice Sheet and Himalayan glaciers (Van Ommen et al 1999;Miles et al., 2017;Talalay et al 2020). Improvements in accuracy of measurements can put tighter constraints on inverted parameters: for example, in our synthetic data inversion experiments, Higher accuracy of englacial thermometry implies that the grounding-induced evolution of the temperature profile can be 390 detected with two measurements separated by a few decades, which could potentially be utilised in previously drilled boreholes ( Fig. 7b). In the case of Crary Ice Rise, a series of forward models shows that 30-year increase in the timing of grounding (e.g., difference between temperature profiles within 1030 and 1000 years old ice rise) may yield a difference in englacial temperatures that could be detected in the lower part of the ice column as well as the upper section of the underlying sediment/bedrock (Fig. 7b). Therefore, if a sufficiently deep borehole is drilled, measurements through the underlying bedrock 395 may provide useful additional constraints on both the heat flux and the timing of ice-rise formation. However, rapid attenuation of temperature differences with time since grounding implies that this method could only be applicable for ice rises that experienced recent grounding (i.e., less than around 1 ka) assuming measurement uncertainties of 0.02ºC, as reported from a borehole drilled through the Dome Summit South in East Antarctica (Van Ommen et al., 1999). This loss of temporal resolution also implies increasing uncertainties associated with inferring the timing of grounding for ice rises that are older than 400 approximately 4 k.y. (Fig. 5e). Additionally, we do not take into account variations of thermal conductivity through the ice column. Temperature-dependant conductivity may yield a slight but measurable effect on the shape of englacial temperature profile, thus contributing to cumulative uncertainty in ice-rise dating. Furthermore, our model does not account for the heating effects that may occur due to basal friction upon regrounding, as well as the freezing of groundwater and its corresponding effect on cooling. A more sophisticated forward model that takes this into account would be required to investigate this effect in water-saturated bedrock or sediments with high porosity.

Inversion of time-variable forcings
The Bayesian inversion method presented here has potential use in providing information about other model parameters, including time-variable thickness, temperature and accumulation rates, as previously done in multiple locations in Greenland and Antarctica (e.g., Dahl-Jensen et al. 1998;Waddington, 2005;Orsi et al., 2012;Cuffey et al., 2016). Since the addition 410 of temporal variability to external forcings significantly increases the dimensionality of the inverse problem, with associated exponential growth of the computation cost, here we focus largely on thermal effects of dynamic ice-rise evolution. However, our forward simulations show the variable impact of these forcings on englacial temperature distribution and their variability with ice thickness and depth within the ice column (Figs 5,6). For example, the effects of relatively small perturbations to accumulation rates are generally greater for thicker ice, whereas the opposite is true for timing of grounding and thickening 415 rates, which play a more important role when ice is thin. Anomalies associated with changes to most parameters under consideration typically increase with depth (insets of Figs.5a-c and Figs. 6a-d), with the exception of surface temperatures, which exert a strong influence on the temperature-depth distribution in the upper part of the ice column (e.g., Dahl-Jensen et al. 1998). This implies that the upper part of the ice column is likely to contain a more recent and less diffuse record of past surface temperatures (e.g., Orsi et al., 2012).
420 Matsouka et al.'s (2015) Antarctic ice-rise inventory shows that ice rises rarely exceed 500 m in thickness (Fig. 1d), suggesting that these features may store a stronger signal of past dynamic changes and thickening history, whereas thick, permanently grounded ice domes may retain more information about accumulation and temperature histories and are thus more appropriate locations for paleoclimatic inferences from englacial temperature measurements (e.g., Dahl-Johnsen et al., 1999;Engelhardt, 2004;Orsi et al., 2012;Cuffey et al., 2016). Forward model experiments with different vertical velocity approximations show 425 that the impact of vertical ice flow parametrization becomes more significant for thicker ice (Fig. 6a,d). Yet, previous studies from Greenland and Antarctica have also shown that analytical ice flow approximations from Dansgaard and Johnsen (1969) and Lliboutry (1979) cannot fully capture the nonlinearity of vertical velocity profiles in ice divide/ice rise settings, and phasesensitive radar measurements can provide useful additional constraints on vertical ice-sheet velocities (e.g., Gillet-Chaulet et al., 2011;Kingslake et al., 2014;Buizert et al., 2021). Therefore, integration of these techniques would help improve inferences 430 of external forcings from borehole temperatures, in particular surface temperature and accumulation histories from deep boreholes. where direct geological sampling methods are inaccessible (e.g., Bentley et al., 2010;Spector et al., 2018). Integrating this technique with other methods, such as: (1) indirect estimates of timing of grounding from radar observations and modelling (e.g., Schroeder et al., 2014;Kingslake et al., 2016;Wearing and Kingslake, 2019), (2) parametrization of vertical velocities (Kingslake et al 2014), and (3) adoption of more tightly restricted, informative prior constraints from geochemical ice-core data (e.g., for past temperature proxies; Cuffey et al., 2016), will allow for more accurate inferences of dynamic ice-rise/ice-440 sheet evolution and grounding line migration (e.g., Orsi et al 2012). Moreover, we have demonstrated an approach for better quantifying uncertainties in these inferences. Borehole measurements through the upper tens of meters of underlying sediment/bedrock could place additional important constraints, both on the geothermal heat flux and ice-rise evolution. This technique could even provide insights into dynamic ice-sheet evolution if future boreholes are drilled through floating ice and sediment in the vicinity of the grounding line, in places where recent ungrounding has left a pronounced vertical temperature 445 anomaly within both ice and sediment/bedrock columns. Our results prompt the question of what characteristics make a location favorable for borehole drilling and measuring temperature-depth data within an ice rise. Due to the diffusive nature of the englacial thermal signal, and as synthetic data 455 experiments have shown, the temporal resolution decreases with time since grounding (Figs. 3, 5b,e, 7a). Therefore, ages of ice rises that are over 4 k.y. old may be difficult to determine accurately, subject to other parameters like heat flux and thickening rates. Areas that are located close to the present-day grounding line (and thus more likely to have been formed relatively recently) with well-constrained, low values of heat flux and low thickening rates could represent optimal locations for implementation of this method and could yield accurate (i.e., on the order of 10%) ice-rise age estimations. Kingslake et al 460 (2018), Venturelli et al. (2020) and Neuhaus et al. (2020) have shown evidence of potential regrounding across large areas of West Antarctica. Preliminary investigations, juxtaposing these maps of potential grounding line migrations (Kingslake et al., 2018;Albrecht et al., 2020b) with Matsouka et al.'s (2015) ice-rise inventory (Figs.1a-d) and Martos et al. (2017) heat flux model (Fig. 8) shows several ice rises where ice is relatively thick and accumulation rates are relatively low for example Korff Ice Rise in the Ronne Ice Shelf area (Kingslake et al., 2016) and small ice rises to the east of the Weddell Sea. These ice rises 465 could prove to be optimal locations for application of this technique. Future work could systematically quantify the suitability of these locations following this approach.

Conclusions
In this paper, we combined Bayesian inversion and forward modelling to make an evaluation of uncertainties inherent in inferences of ice-rise dynamic evolution from temperature-depth profiles. Tested with both synthetic datasets and borehole 470 temperature measurements from Crary Ice Rise, Ross Sea Embayment, our method explores the interplay of surface temperature, rates of accumulation and thickening, geothermal heat flux and parameterized vertical velocities. We show that depending on the accuracy of borehole thermometry, the same temperature profile (within the accuracy of measurements) may result from a range of forcing parameters, of which geothermal heat flux through underlying bedrock plays a particularly important role. The key implication is that careful model parametrization and evaluation of uncertainties are essential to infer 475 dynamic ice-rise evolution from borehole thermometry. We highlight that uncertainties in inferred ice-formation time may increase significantly with ice-rise age. Accuracy of inversion relies on the low measurement uncertainties (i.e., <0.05ºC) and can be high (i.e., uncertainties <10%) for relatively young ice rises (i.e., formed <4 ka) that are grounded in areas where heat flux is low, and its value is well-constrained.

480
Code availability. The code related to this article is available on-line at: https://github.com/sashamontelli/borehole_temperature_models/blob/a85d7b81cea89c413884ce3af6f167ff828e4ed6/Annotated%20temper ature%20profile%20MCMC%20inversion.ipynb Author contributions. A.M. co-designed this research, performed the analysis, and wrote the manuscript. J.K. co-designed this research and 485 contributed to the writing and editing of the manuscript.
Competing interests. The authors declare that they have no conflict of interest.

Acknowledgements.
A.M. is grateful to the Schmidt Science Fellowship for funding of this research. We thank Thorsten Albrecht for granting access to the PISM model outputs and Nicholas Holschuh for helpful discussions of temperature measurements from Crary Ice Rise.
Financial support. This research has been supported by the Schmidt Science Fellowship to A.M.