Constraints on glacier flow from temperature-depth profiles in the ice . Application to EPICA Dome C .

A one-dimensional (1-D) ice flow and heat conduction model is used to calculate the temperature and heat flux profiles in the ice and to constrain the parameters characterizing the ice flow and the thermal boundary conditions at the Dome C drilling site in East Antarctica. We use the reconstructions of ice accumulation, glacier height and air surface temperature histories as boundary conditions to calculate the ice temperature profile. The temperature profile also depends on a set of poorly known 5 parameters, the ice velocity profile and magnitude, basal heat flux, and air-ice surfaces temperature coupling. We use Monte Carlo methods to search the parameters’ space of the model, compare the model output with the temperature data, and find probability distributions for the unknown parameters. We could not determine the sliding ratio because it has no effect on the thermal profile, but we could constrain the flux function parameter p that determines the velocity profile. We determined the basal heat flux qb = 49.0± 2.7(2 )mW m , almost equal to the apparent value. We found an ice surface velocity of 10 vsur = 2.6±1.9(2 )m y 1 and an air-ice temperature coupling of 0.8±1.0(2 )K. Our study confirms that the heat flux is low and does not destabilize the ice sheet in east Antarctica.


Introduction
Climate change, as a consequence of human activities, has become a focus of attention of the scientific community.Climate change can be seen as the consequence of an imbalance between the emission and absorption of energy by the planet.Such global energy imbalance generates exchanges of energy among climate subsystems as the Earth follows a path towards a new equilibrium.The cryosphere, and in particular the polar ice sheets, play a key role in the climate system, influencing ocean temperatures, sea level, thermohaline circulation and global albedo, responding to the energy imbalance on time scales ranging from decades to millennia (Clark et al., 1999).Therefore, understanding ice sheet dynamics, controlled by basal heat flux (Oppenheimer, 1998), is of major importance for modelling future climate change.
The analysis of air bubbles trapped in the ice at the time of its formation allows for the reconstruction of past atmospheric concentrations of CO 2 , CH 4 and N 2 O (e.g.Barnola et al., 1987;Spahni et al., 2005).Past temperatures can be estimated through the analysis of the stable isotope ratios of 18 O and deuterium in the ice molecules (Pol et al., 2010), which also yields past snow accumulation rates (Jouzel et al., 2007).
The interpretation of paleoclimatic records from ice cores, as well as the interaction of the ice sheet with the climate system and the solid earth are not straightforward.Ice sheets are complex systems subject to significant variability in their dynamics and structure, which depend on the thermal and mechanical basal boundary conditions (e.g. the basal heat flux and the sliding ratio of the ice sheet) (e.g.Marshall, 2005), and surface conditions (e.g.snow accumulation rates and atmospheric temperatures).
Analysis of age versus depth of the ice cores provides some insight on ice flow mechanics (Parrenin et al., 2004(Parrenin et al., , 2007b)), through the determination of the ice thinning function, the vertical velocity, as well as the basal melting rate and the sliding ratio.These age depth profiles also provide some constraints on the basal heat flux which is a essential parameter for the dynamics of ice sheet models (Fisher et al., 2015).
High resolution temperature profiles measured in drill holes also retain a record of the time dependent boundary conditions and ice sheet dynamics.These temperature-depth profiles provide another means to determine basal heat flux.In this paper, we use a numerical, one-dimensional (vertical) advective-conductive model which includes the key parameters driving icesheet dynamics.The model is forced by the local histories of atmospheric temperatures, snow precipitation and ice thickness, to determine temperature profiles for the ice-sheet.We then use a Monte Carlo inversion scheme to find set of values of the parameters that fit the observed temperature profile, determine their probability distributions and their most likely values.
We have applied this method to the ice core drilled at EPICA Dome C in Antarctica, one of the longest ice core records presently available.The temperature profile through the entire thickness of the ice-sheet was measured in the drill hole (Pol et al., 2010).The forcing is provided by the precipitation and atmospheric conditions for the past 800ka that have been reconstructed by analyzing the ice core extracted from the drill hole (Jouzel et al., 2007).The method used in the present study could be applied to any ice temperature profile extending through most of the ice-sheet thickness and provide independent constraints for ice dynamics parameters.

Basic assumptions and equations
The drilling at the Dome C site (75°6 0 6.35 00 S, 123°23 0 42.76 00 E), with an ice thickness of 3273±5m (Tabacco et al., 1998), was stopped ⇠ 15m above bedrock, allowing the measurement of the temperature profile through the entire glacier.The temperature profile (Fig. 1) is determined by both the thermal boundary conditions and the ice dynamics which controls heat transport by advection and produces shear heating.
We follow the theory of ice sheet dynamics of Paterson (1994).The ice sheet grows by accumulation of snow at the top, increasing its gravitational potential energy that drives the flow of the ice sheet toward the ocean.Due to basal friction, the velocity of the flow is maximum at the top of the ice sheet and decreases with depth, consequently thinning the ice layers and reducing the height of the glacier.The ice surface is in contact with the air above, absorbing heat from the ice sheet.As heat flows into the glacier from the underlying bedrock, the temperature gradient in the ice is positive downward, which could bring the base to melt.Meltwater reduces basal friction and allows the ice sheet to slide over the bedrock, a movement defined through the sliding parameter s, the ratio between the ice horizontal velocity at the base and at the top of the ice (Fig. 2).
We have developed a forward model to simulate the thermal processes that take place in the ice sheet, including both heat diffusion and advection by ice movement, and therefore defined as an advective-conductive model.This model calculates a temperature profile that is determined by heat conduction, the dynamics of the ice sheet, and their boundary conditions.Ice sheet dynamics is described through the field equations, i.e. the conservation of mass and momentum.The mass conservation for a material of density ⇢ is: where !u is the velocity of the ice.
For an incompressible material ( D⇢ Dt = 0), a reasonable assumption for ice, the mass conservation equation reduces to: !r • !u = 0.
(2) For steady state flow (without acceleration), momentum conservation requires the balance between the body forces acting on the ice volume, i.e., weight ⇢g (where g is the acceleration of gravity), and the internal forces described by the stress tensor (3) For ice, the constitutive equation has been established by experimental work (Glen, 1955) and corresponds to non linear viscous rheology.The empirical law that gives the relationship between the strain rate ✏ and the shear stress ⌧ is known as Glen's law: where A is a temperature-dependent empirical quantity (Paterson, 1994).
At the local scale it is possible to model ice flow in two dimensions by choosing the horizontal axis in the direction of ice flow.In addition, because the size of the ice sheet is much greater than its thickness and variations in crustal heat flux occur on a scale of tens of kilometers (Jaupart and Mareschal, 2015), we can assume that horizontal temperature gradients are negligible and therefore conductive heat flow q (W m 2 ) is vertical and given by Fourier's law in one dimension: where (W m 1 K 1 ) is the thermal conductivity and z is the vertical coordinate, defined positive upwards with the origin at the base of the ice.
The one dimensional heat equation is (Carslaw and Jaeger, 1959): where  = ⇢cp (m 2 s 1 ) is the thermal diffusivity, ⇢ = 916.2kgm 1 is the ice density, c p (J kg 1 K 1 ) is the specific heat, and u z (m s 1 ) is the vertical component of the ice velocity.The function ⌦ (µW m 3 ) is the rate of heat production per unit volume, due to shear heating.Absorption of latent heat at the base is not included in this term but introduced in the basal heat flux boundary condition.
The conductive heat flux profile obtained for Dome C (Fig. 1) can be used to estimate the basal heat flux.From a linear regression of the lowermost 300m of the heat flux profile, we obtain an apparent basal heat flux of q b,app = 49.4±1.1 mW m 2 .
However, internal heating caused by shear deformation affects the profile and its effect must be taken into account.Therefore, we consider q b as a free parameter, whose probability distribution is to be evaluated.

Boundaries
The heat equation requires two boundary conditions.At the top of the ice column, the temperature is the Glacier Surface Temperature (GST).It is obtained directly by assuming a constant offset (see subsection 3.3) from the Surface Air Temperature (SAT) history, which is known from the analysis of the deuterium content of the ice core (Pol et al., 2010), and is in excellent agreement with the reconstruction from 18 O record (Jouzel et al., 2007).At the bottom, the contact between the ice column and the bedrock, the boundary condition is a fixed and constant basal heat flux q b .If the basal temperature reaches the melting point, and basal heat flux is enough to maintain it, the excess heat flux is absorbed as latent heat (see subsection 4.1), and the melting temperature becomes the effective bottom boundary condition.

Ice thickness history
Ice thickness H varied significantly (3100 ⇠ 3300m) during the past 800ky that we model.We used a prescribed the history of ice thickness, obtained from the linear perturbation model developed in Parrenin et al. (2007b), which is consistent with a large and complex 3-D thermo-mechanical model (Ritz et al., 2001).The ice thickness history obtained from this model is consistently lower (⇠ 350m) than in that obtained by Parrenin et al. (2007b), but the relative variations of ice thickness are consistent.Therefore, we have used the relative variations of ice thickness defined by this model, using as reference the current ice thickness at Dome C H 0 = 3266m (Fig. 3).This history of ice thickness was used to calculate @H @t dynamically, so that H follows the prescribed history.

Surface accumulation and SAT histories
The accumulation rate a (m y 1 ), measured in ice-equivalent units, and the SAT are determined from deuterium content D (Jouzel, 2007;Pol et al., 2010):  where a 0 and T 0 are the accumulation rate and the surface temperature for a reference deuterium content of 396.5‰.
D cor is the corrected deviation (Lliboutry, 1979) from the current deuterium content of the ice.D smo is a 50-year smoothed version of D cor because the accumulation rate is only supposed to be related to the deuterium content over a certain time interval (Parrenin et al., 2007b). We

Temperature offset and snow-firn cover
Accumulation at the surface of the ice column is measured in ice equivalent units, but precipitations are first deposited in the form of snow.Before being transformed in ice, snow passes through an intermediate phase called firn (Arnaud et al., 2000;Goujon et al., 2003).This process is driven mainly by the atmospheric temperatures and the accumulation rates.The density of snow/firn increases gradually with depth, eventually becoming ice.We estimated the average firn density from Dome C data (Arnaud et al., 2000) and modeled the upper 80m as firn, with a uniform density 75% that of ice.
The temperature signal that propagates into the ice column is the GST, a filtered version of the SAT.Measurements at Vostok station over the last 50 years (Vostok, 1958(Vostok, -2016) ) show that the GST is ⇠ 5K warmer than the SAT in summer and ⇠ 2.5K colder in winter.We assume the mean temperature offset T o↵set (K) to be time-independent.
Because imprecisions on the estimate of the snow cover of the ice column make T o↵set difficult to estimate, it is set as a free parameter restricted within a range [ 2.5K, 5K].

Basal melting
The forward model calculates dynamically the basal melting rate.When the basal ice layer reaches melting temperature T m (K), the extra incoming energy is used to calculate melting rate m (m s 1 ).The melting temperature depends on the pressure as: where CC = 7.42 ⇥ 10 8 K Pa 1 is the Clausius-Clapeyron slope and P (Pa) is the pressure at the glacier base.
When the basal temperature reaches T m , the excess heat flux (incoming minus outgoing) at the base q b (mW m 2 ) is absorbed as latent heat of fusion by the ice, and the melting rate m is calculated as: where L f = 3.337 ⇥ 10 6 J kg 1 is the latent heat of ice.If q b < 0, the basal temperature decreases below the melting point, and melting stops.
The melted ice is subtracted from bottom of the ice column and reduces its thickness.The effect of melted ice on sliding is neglected.

Glacier movement
We calculated the vertical movement of the glacier with the 1-D ice flow model used in Parrenin et al. (2007b), which gives vertical velocity u z (m s 1 ) as: where a (m s 1 ) and m (m s 1 ) are respectively the accumulation (at the surface) and the melting (at the base) rates, and is the height of the glacier, expressed as ice equivalent.![⇣] is the flux shape function (Parrenin et al., 2006), which depends on the non-dimensional vertical coordinate ⇣ = z/H and includes contributions of the basal sliding term and the shear deformation term: where s is the sliding ratio (varying between 0 and 1) and !D [⇣] is the vertical profile of deformation (Lliboutry, 1979), given by: where p is the parameter determining the shear deformation component of the flux shape function.Lliboutry (1979) suggested that it is approximately given by: where n is the exponent of Glen's law from Eq. ( 4), Q = 6 ⇥ 10 4 J mol 1 is the activation energy, R = 8.3145J mol 1 K 1 is the gas constant, T b is the basal temperature, and G 0 (K m 1 ) is the vertical temperature gradient at the bottom of the ice column.
In uncertainty on the parameter p, we decided to set p as a free parameter within a range [0, 9] for which we will try to obtain a probability distribution.
The value of the sliding parameter s at Dome C is not well constrained, though Parrenin et al. (2007b) determined that it is less than 0.3, and less than 0.1 with 50% confidence.However we do not include sliding as a free parameter, because sliding has a little effect on the thermal model and the calculated thermal profiles are much less sensitive to sliding than to the other parameters.This was verified in preliminary tests of the model taking sliding as a free parameter.For simplicity, we set the sliding parameter as a fixed value s = 0.

Shear heating
Preliminary attempts to model the temperature profile at Dome C showed that the profile could not be fitted without some internal heat sources, that we assumed to be shear heating.The heat production by shear deformation ⌦ (µW m 3 ), or simply 'shear heating' is given as (see appendix C): where ✏zz is given by Eq. ( A1) and ✏xz is given by: where Ūx is the average horizontal velocity.We calculate it as Ūx = v sur /! 0 (1), where v sur is the ice velocity at the surface, relative to the bedrock.Estimates of ice velocity from geodetic surveys (Vittuari et al., 2004) at Dome C give a horizontal ice velocity of a few mm per year at the topographical dome, and 15 ± 10mm y 1 at the drilling site by tying it to a 25kmaway point.However, satellite measurements of ice velocity (Mouginot et al., 2012;Rignot et al., 2011) with a resolution of 450m gave velocities of v 13m y 1 at the four closest points to the drilling site.Tests with surface velocity in the order of 15mm y 1 yield a negligible shear heating, while the heating rate for a surface velocity on the order of 13m y 1 seems excessive and results in higher than observed curvature of the profile.In absence of reliable estimates of the ice surface velocity, we set it as a free parameter.
The influence of shear heating can be clearly seen in the heat flux profile (Fig. 5b).The heating rate is maximum at the bottom of the ice column and decreases upwards.For a given basal heat flux, shear heating gives a characteristic shape to the heat flux profile.The vertical distribution of shear heating is ultimately dependent on the flux shape function in Eq. ( 13), therefore dependent on the value of p (Fig. 5b).This characteristic shape can be observed in the Dome C heat flux profile (Fig. 1), which shows a small but non negligible shear heating.

Free parameters
We have introduced four parameters that we consider free: the basal heat flux q b , the parameter p, the ice surface velocity v sur , and the ice-air temperature offset T o↵set .For each set of parameter values, we calculate a temperature-depth profile with our advective-conductive model.
We select randomly the values of the parameters within a range of plausible values and calculate a temperature profile that we compare with the one measured at Dome C. For each combination of parameters, we calculate the misfits to temperature and heat flux profiles as the root mean square (RMS) difference between calculated and measured profiles.We retain a set of parameters if both misfits are below two respective maximum values or cutoffs for temperature and heat flux.
We set these cutoffs in such a way that we eliminate as many trials as possible, while keeping enough to obtain a meaningful statistic.We set the temperature cutoff to 0.5K and the heat flux cutoff to 2mW m 2 .The procedure to obtain these cutoffs is discussed in the appendix B. With both cutoffs as acceptance condition, we retained ⇠0.18% of all trials, obtaining histograms of accepted values for each free parameter.The number of accepted experiments within a range of values of the free parameters defines a probability density, from which we derive the mean and standard deviation to obtain the most likely values and a 2 confidence interval.

Initial conditions
The initial temperature profile of Dome C at 800ka is impossible to determine from available data.However, because of the character of heat diffusion, the sensitivity of the temperature profile to the initial condition decreases with time.We performed several tests to confirm that the outcome of the simulations is independent of the initial temperature condition.We have thus used the same initial temperature profile for all simulations.This initial condition is obtained by applying the forward model to the present profile of Dome C, using as boundary conditions the time reversed history of ice thickness, surface accumulation and SAT.The values of the free parameters are unknown for this backward in time simulation, but are not of critical importance 10 because of the insensitivity to initial conditions.We chose basal heat flux q b = 49.4mW m2 (the apparent heat flux), while for the parameter p, the surface velocity and the temperature offset we used p = 9, v = 0m y 1 and T o↵set = 0K.

Results
The histograms of retained values for each parameter are shown in Fig. 6.The peaks corresponding to the most likely values of the parameters are well marked suggesting that the parameters are well constrained, except for the parameter p.We can see how high values of q b are accepted, when in combination with low values of surface velocity (and therefore shear heating), but even very high (⇠ 60mW m 2 ) are accepted.This is happens due to melting, as high values of basal heat flux are able to keep the base of the glacier at fusion point through the whole simulation.Further increase in basal heat flux increase melting while the base of the glacier stays at melting temperature, having little effect on the resulting profiles.
The most probable value for the basal heat flux is q b = 49.0 ± 2.7mW m 2 , slightly lower than the apparent heat flux of q b,app = 49.4±1.1mW m2 obtained from thermal conductivity and the thermal gradient at the bottom 300m of the measured profile, stopped 60m above the bedrock.We expected to obtain a slightly lower value than that of the apparent heat flux which includes a contribution from shear heating.The parameter p is not tightly constrained but low values of p are more likely, with a ⇠ 90% probability that p is < 2.5, and ⇠ 50% probability that p is < 1.
The most likely surface velocity is v sur = 2.6 ± 1.9m y 1 .This range is compatible with the surface velocity v 13m y 1 inferred from satellite measurements (Mouginot et al., 2012;Rignot et al., 2011), and results in substantial heat generation by shear deformation.The ice velocity measurements of 15 ± 10mm y 1 from geodetic surveys (Vittuari et al., 2004) is too small to generate enough shear heating.As we assume a sliding ratio s = 0 in our simulations, the surface velocity is the component corresponding to the shear deformation movement.The real surface velocity, sum of deformation and sliding components, is higher than the one obtained.The most likely temperature offset between the GST and the SAT is T o↵set = 0.8 ± 1.0K.This result could be affected by errors on the thickness of the firn layer, as well as on the SAT reconstruction itself.
We have recalculated the temperature profile with the most likely values, and a value p = 1 for the parameter p (Fig. 7).The good match of the measured and calculated temperature profiles supports that the basal heat flux and the shear heating are in the correct range of values.

Discussion and conclusions
In this study, we have modelled the ice flow and heat transport at Dome C with a 1-D thermal and mechanical model, using the histories of air temperature, ice thickness and snow accumulation.We have tuned key parameters of ice flow and thermal boundary conditions so that the calculated temperature-depth profile fits the measured thermal profile at Dome C, in order to obtain the most likely values of these parameters.
We have found that shear heating has a strong effect on our thermal model.Attempts to fit the temperature and heat flux profiles without taking shear heating into consideration yielded profiles determined by the thermal boundary conditions at the surface and the base of the glacier, i.e.SAT and basal heat flux.These factors alone failed to explain the shape of the temperature-depth profile.For the highest estimates of basal heat flux (obtained from the temperature gradient near the base) the temperature profiles were systematically ⇠ 5K colder than the one measured.Further tests showed that both temperature and heat flux profiles could be successfully fitted with a depth-uniform heating rate of ⇠ 3.5µW m 3 , but a uniform heating rate has no physical justification.Clim. Past Discuss., doi:10.5194/cp-2016-116, 2016 Manuscript under review for journal Clim.Past Published: 30 November 2016 c Author(s) 2016.CC-BY 3.0 License.Shear heating is the only possible source of internal heating, but it varies with depth.The variation of shear heating with depth depends on the flux shape through Eq. ( 12).Therefore it is possible to extract information about the flux shape function because the goodness of the fit of the model to the profiles depends on it.The flux shape function is controlled by the parameter p, for which the theoretical formula (Eq.15) predicts values on the order of p ⇠ 9. Initial tests with this theoretical value failed to fit the Dome C profile satisfactorily, because shear heating was excessively concentrated in the lower part of the glacier.Parrenin et al. (2007b) argued that this formula is not be applicable to Dome C, and found the most likely value for this parameter p = 1.97 ± 0.93, based on a best fit of ice chronology to age markers.Although we have not been able to precisely determine the parameter p, our reconstruction suggests that the value of p is likely to be closer to the value obtained by Parrenin et al. (2007b) than to the theoretical formula's range of values.
The deformation component of ice surface velocity, the parameter that defines the magnitude of shear heating, is reasonably well constrained.Accounting for the sliding ratio s < 0.3 at Dome C (Parrenin et al., 2007b), the value of the ice surface velocity is of the same order as but lower than the values obtained by satellite measurements of ice velocity (Mouginot et al., 2012;Rignot et al., 2011) in the vicinity of Dome C. The estimates of ice velocity from geodetic surveys at Dome C (Vittuari et al., 2004) are too small to produce relevant shear heating.
The temperature offset is also well constrained, but it is an operational parameter that cannot be measured directly.
The basal heat flux value 49mW m 2 is almost equal to the apparent heat flux near the base of the glacier, as expected because, without melting, the heat flux must be continuous.Because the apparent heat flux is estimated over 300m, it may be slightly affected by shear heating near the base of the glacier.There is a very good fit between the Dome C profile and the thermal profile determined by the most likely values we obtained for our parameters, suggesting that the values we obtained are correct.The crust below East-Antarctica is believed to be made up of Archean and Proterozoic terranes that were welded together in the mid-Proterozoic (Dalziel, 1992;Harley, 2003).The value of 49mW m 2 obtained for the heat flux below Dome C is well within the range of values characteristic of Archean or Proterozoic terranes (Jaupart and Mareschal, 2015).In stable continents, local variations of surface heat flux do occur on different scales due to variations in crustal heat production (Jaupart et al., 2016).The heat flux value obtained is plausible and consistent with an estimate of heat flux beneath east Antarctica based on shear wave velocity profiles in the upper mantle (Shapiro and Ritzwoller, 2004).The latter estimate of average regional heat flux is insensitive to local heat flux variations.Although the upper part of the glacier may move over the bedrock by maybe as much as 100km per 100,000 years, the lowermost oldest part of the glacier does not move if there is no slip and it has experienced a constant heat flux for a time much longer than the characteristic heat conduction time.
The thermal model allows us to calculate basal melting dynamically.While basal melting does not have a direct feedback on our model, it is an important basal parameter that controls sliding.Solutions were obtained with high values of basal heat flux (more than 65mW m 2 ) that lead to constant melting temperature at the base.In this situation, increased basal heat flux can be absorbed by melting without affecting significantly the temperature profile.This has not been an issue in Dome C, but the method will not work for determining basal heat flux if it is high enough to keep the base of the glacier at melting temperature, unless we introduce another constraint on the melting rate.For Dome C we can determine an upper limit to melting, as both glacier height and accumulation follow prescribed histories.This upper limit is never reached through the simulations.through the Glen's law: where u i (m s 1 ) are the components of the ice velocity. 0x is defined as: ✏xx + ✏yy + ✏zz = @u x @x + @u y @y + @u z @z = 0. (C6) In Eq. (C2) and Eq.(C4), n = 3 is the Glen's exponent and A (Pa n s 1 ) is a quantity whose units depend on the Glen's exponent and that depends only on temperature: where A 0 = 3.985 ⇥ 10 13 Pa n s 1 , Q 0 = 6 ⇥ 10 4 J mole 1 is an activation energy, R = 8.3145J mole 1 K 1 is the ideal gas constant, and T (K) is the absolute temperature.
The quantity ⌧ is the effective shear stress, defined as: For the sake of simplicity, we assume a state of plane strain, with all components of velocity, strain and stress are independent of y, and their y-components u y = 0, ✏yy = 0 and y = 0 are zero.Consequently, we have from Eq. (C4) and Eq.(C5): while Eq.(C6) and Eq.(C8) are reduced respectively to to: ✏xx + ✏zz = @u x @x + @u z @z = 0, (C10) With these simplifications, the heat production Eq.(C1) becomes: We can solve ⌧ , ⌧ xz and 0 x from the Eqs.(C2), (C3) and (C11).It gives: The value of ✏zz is obtained from Eq. (A1).Assuming the ice flow to be laminar ( @uz @x =0), ✏xz , Eq. ( C2) is reduced to: The horizontal velocity u x is unknown, but it varies with depth (Parrenin et al., 2006) as: where Ūx is the average horizontal velocity at the surface of the glacier.Therefore, Eq. (C14) can be written as:

Figure 1 .
Figure 1.(a): Temperature-depth profile measured in Dome C. (b): Conductive heat flux profile, calculated from the temperature profile by Fourier law with a temperature-dependent thermal conductivity.Heat flux has been truncated at a depth of 225m because it is very noisy near the surface of the ice sheet.

Figure 2 .
Figure 2. Sketch of flow in an ice sheet.The velocity of ice particles, both vertical and horizontal components, is highest near the surface and decreases at depth.The ice motion thins the ice layers and reduces the height of the glacier, while accumulation of snow at the surface increases its height.Temperature in the ice sheet increases downward, as heat flows out of the bedrock, and melting of ice could happen at the base.The ratio of the horizontal components of velocities at the bottom and at the surface of the ice defines the sliding ratio s.

Figure 3 .
Figure 3. Variations of thickness in Dome C for the last 800ky, simulated with the linear perturbation model developed in Parrenin et al. (2007b).
Dome C, using the values n = 3, G 0 = 21.62Kkm 1 , H = 3266m and the current basal temperature T b = 270.2K,we obtain a value p ⇡ 9.All the parameter values except T b and H are assumed constant, and p remains almost constant during the simulation.However, initial tests of the model using the theoretical formula from Eq. (15) were not able to produce thermal profiles close enough to the measured Dome C profile.Parrenin et al. (2007b) pointed out several reasons why this theoretical value of p could be invalid at Dome C, and estimated its value by inverse method obtaining p = 1.97 ± 0.93.Due to this Clim.Past Discuss., doi:10.5194/cp-2016-116,2016   Manuscript under review for journal Clim.Past Published: 30 November 2016 c Author(s) 2016.CC-BY 3.0 License.

Figure 5 .
Figure 5. (a) Temperature-depth and (b) heat flux-depth profiles calculated with different vertical distributions of shear heating, determined by the value or the parameter p, for which we use the theoretical value (red) obtained from Eq. (15) and the value p = 2 (blue) obtained by Parrenin et al. (2007b).Note how shear heating influences the shape of the heat flux profile.The other parameters for these calculations are q b = 49.4mW m2 , vsur = 8m y 1 , and T o↵set = 0K.The no-heating experiment (black) uses the theoretical value p ⇡ 9 with shear heating set to 0.

Figure 6 .
Figure 6.Histograms of retained values for the free parameters: (a) basal flux, (b) parameter p, (c) surface velocity, and (d) GST-SAT temperature offset.

Figure 7 .
Figure 7. (a): Temperature-depth profile of Dome C (black) versus a simulation with parameter p = 1 and the most likely values of the other parameters q b = 49.03mW m2 , vsur = 2.59m y 1 and T o↵set = 0.815K (blue).(b): Heat flux profile determined from the temperature profile, measured at Dome C (black) and calculated for the most likely values (blue).

Figure A1 .
Figure A1.Change in the position of the q b peak, as function of the temperature cutoff T cuto↵ .Heat flux cutoff was set to q cuto↵ = 2mW m 2 to obtain this graphic.

Figure A2 .
Figure A2.Change in the position of the q b peak, as function of the heat flux cutoff q cuto↵ .Temperature cutoff was set to T cuto↵ = 0.5K to obtain this graphic.