Borehole temperatures reveal a changed energy budget at Mill Island , East Antarctica , over recent decades

J. L. Roberts1,2, A. D. Moy1,2, T. D. van Ommen1,2, M. A. J. Curran 1,2, A. P. Worby3, I. D. Goodwin4, and M. Inoue2,5 1Department of Sustainability, Environment, Water, Population and Communities, Australian Antarctic Division, Channel Highway, Kingston, Tasmania 7050, Australia 2Antarctic Climate and Ecosystems Cooperative Research Centre, University of Tasmania, Private Bag 80, Hobart, Tasmania 7001, Australia 3CSIRO Marine and Atmospheric Research, Castray Esplanade, Hobart, Tasmania 7000, Australia 4Marine Climate Risk Group, Department of Environment and Geography, Macquarie University, Eastern Road, Macquarie University, New South Wales 2109, Australia 5Institute for Marine and Antarctic Studies, University of Tasmania, Private Bag 129, Hobart, Tasmania 7001, Australia


Introduction
Palaeoclimate records provide an essential context for present-day climate, and they help us to understand the drivers of climate change.Ice core records from Antarctica have provided insights into the cycle of glacial and in-terglacial periods as far back as 800 000 yr (Jouzel et al., 2007) and show how the ice sheets have responded to past changes in temperature.However a lack of high-resolution climate data, especially for the Southern Hemisphere (Mann and Jones, 2003;Neukom and Gergis, 2012), still limits our understanding of climate processes over the past several millennia.In particular, recent Northern Hemisphere temperatures have, on average, been warmer than anytime in at least the last 1300 yr (Mann and Jones, 2003;Mann et al., 2008).In contrast, the relatively large uncertainties (likely due to data sparsity) in Southern Hemisphere palaeoreconstructions (Mann and Jones, 2003;Mann et al., 2008) are of comparable size to the recent warming, so that the Southern Hemisphere situation is more ambiguous.Records from sites such as Mill Island help address this data sparsity and provide insights into the rate of change over the past century.
While spatially isolated temperature reconstructions are useful for investigating local changes in climatic conditions, consistent features of regional-scale climatic conditions can be reconstructed with a spatially distributed network of observations.Previous regional studies on late Holocene paleoclimate have demonstrated a relationship between Law Dome summer temperature in Wilkes Land and evaporation in Ace Lake in the Vestfold Hills, Princess Elizabeth Land (Roberts et al., 2001).This indicates that the East Antarctic region spanning Wilkes Land, Queen Mary Land and Princess Elizabeth Land may experience a contiguous air temperature history on climatological time scales.In order to explore the regional sensitivity to temperature variability further, we focused our research on the most northerly limit of the East Antarctic Ice Sheet, in Queen Mary Land.There are very limited instrumental temperature data records available for this broad region of East Antarctica.The closest long-term record is from Mirny Station with observations since the 1957-1958 International Geophysical Year.Hence, we employed palaeo-climate reconstruction methods using ice cores and borehole measurements to investigate temperature variability.The reconstructions were obtained from Mill Island (65 • 30 S, 100 • 40 E, see Fig. 1), a small ice cap (∼ 40 km radius) which is detached from the East Antarctic Ice Sheet and lies at the northern edge of the Shackleton Ice Shelf.The Mill Island summit has an elevation of ∼ 500 m above mean sea level.Furthermore, this is the most northerly region of Antarctica outside the tip of the Antarctic Peninsula, and therefore a Mill Island temperature reconstruction represents the most northerly temperature record for East Antarctica.
Mill Island experiences a polar maritime climate with four primary influences: (i) the passage and decay of circumpo-lar low pressure systems (cyclonic eddies and polar frontal depressions) that transport moist and relatively warm air masses from the Southern Ocean, and result in orographic precipitation over the ice cap; (ii) high pressure ridging associated with the circumpolar longwave connecting the Antarctic and the Indian Ocean subtropical high; (iii) strong katabatic wind drainage from the Denman Glacier valley, drains over the Shackleton Ice Shelf to Mill Island transporting relatively, cold, dry air mass; and (iv) localised summer season sea-breezes associated with sea-ice breakout that result in low level cloud, fog and rime formation over the ice cap summit.In the absence of long-term meteorological observations, the general climate overview for Mill Island can be approximated from Mirny Station at 66 • 33 S, 93 • 01 E that indicates prevailing east-southeasterly winds (Turner and Pendlebury, 2004).

Temperature observations
In the summer of 2009/10 the Australian Antarctic program drilled a 120-m-deep ice core as part of an ongoing program to obtain palaeoclimate records from a network of East Antarctic sites for the purpose of reconstructing regionalscale climate, including temperature, precipitation and circulation indices (Jones et al., 2009).The high temporal resolution, due to the high snow fall rates of 1312 kg m −2 yr −1 at this site, and other high accumulation coastal sites in East Antarctica, including Law Dome, makes these sites particularly valuable for obtaining detailed climate records.
In summer of 2010/11 a field party returned to Mill Island to measure the temperatures in the dry borehole (65 • 33 25.84 S,100 • 47 11.44 E and ice surface elevation of 503 m above mean sea level), which was capped in the summer 2009/10.The dry-hole temperature measurements were made using a Leeds and Northrup resistance bridge (model number 8078) to measure the resistance of a four-wire platinum probe.Measured resistances were converted to temperature using the standard relation for industrial platinum temperature sensors.The platinum probe was coupled to a rope marked with 1-5 m graduations, and weighted to ensure the sensor tip was in contact with the ice wall of the borehole.Temperature readings are shown in Fig. 2 and Table 1.The borehole temperature measurements have an instrumental accuracy of approximately 0.02 • C. All measured depths are relative to the 2010/11 surface, and ages are calculated from this datum.
Surface temperature variations are known to propagate into ice sheets, with the attenuation of this temperature signal with depth being strongly influenced by the frequency of the surface variations (Patterson, 1994).Seasonal-scale variations attenuate to 5 % of the surface amplitude at around 10 m below the surface and are "undetectable below a depth of 20 m" (Patterson, 1994).Due to the lack of discernible seasonal temperature variations for most of the depth of the  fore we discarded the upper three measurements when assessing the trend in Mill Island surface energy budget.These upper three temperature measurements correspond to an advection time of 6 yr, while the shallowest temperature measurement included in this analysis (at a depth of 19.07 m below the surface) corresponds to an advection time of 8.5 yr (see Fig. 3).The reconstructions of surface temperature history extrapolate any recent temperature trend over these final 8.5 yr.
Analysis of the temperature data below 15 m suggests a break in the gradient of temperature as a function of depth occurs at 49 m below the surface.Specifically, the BREAK-FIT algorithm (Mudelsee, 2009) was used to fit a two linear segment model to the temperature data using a discrete brute force search approach that constrains the depth of the gradient break point to correspond to the depth of one of the temperature measurements.The temperature gradient is 0.0165 ± 0.0007 K m −1 above the break point and 0.0034 ± 0.0003 K m −1 below it.
The stratigraphy of the Mill Island ice core indicates past surface melting events as shown by the presence of melt layers in the 120 m ice core record.The current stratigraphy suggests there is no discernible change in the thickness and concentration of melt layers throughout the entire length of the ice core.

Numerical model
To reconstruct a surface temperature history that is consistent with the observed temperature profile down the borehole, a  Two inverse models that optimise surface temperature histories to minimise (in a least squares sense) the mismatch between the simulated and observed down borehole temperature profiles are described.Specifically, the two inverse models minimise the unweighted RMS error: One of these inverse models, the least squares QR (LSQR) model (see Sect. 3.2.1), is computationally very efficient and produces a single optimal solution with minimal variance, while the second method, the particle swarm optimisation (PSO) method (see Sect. 3.2.2),requires significantly more computational resources, but produces a distribution of likely temperature reconstructions.
Such an approach of using both a forward model and a method of optimising surface temperature histories has been used previously (see, for example Johnsen, 1977;MacAyeal et al., 1991;Johnsen et al., 1995;Cuffey et al., 1995;Cuffey and Clow, 1997;Dahl-Jensen et al., 1998, 1999;Barrett et al., 2009;Muto et al., 2011;Orsi et al., 2012).The forward models used in these studies all included the heating associated with firn densification with the exceptions of MacAyeal et al. (1991); Johnsen et al. (1995); Barrett et al. (2009) and Muto et al. (2011).These studies all select which possible surface temperature history reconstructions to accept based on minimising the mismatch between the simulated and observed borehole temperature profiles.However, the methods of generating possible surface temperature histories (to select from) differ and fall into three broad categories (Orsi et al., 2012): (i) optimisations based on relationships between surface temperature and stable water isotope ratios (e.g.Johnsen, 1977;Cuffey et al., 1995;Cuffey and Clow, 1997;Johnsen et al., 1995); (ii) Monte Carlo based approaches (e.g.Dahl-Jensen et al., 1998, 1999;Barrett et al., 2009); and (iii) generalised least-squares solution of a linearised version of the problem (e.g.MacAyeal et al., 1991;Muto et al., 2011;Orsi et al., 2012).The methods used herein most closely resemble the latter two methods for the PSO and LSQR methods respectively.

Forward model
The equation for the evolution of temperature (T ) with time (t) in the absence of horizontal advection is (see Table 2 for nomenclature) where the terms on the right-hand side are the rate of temperature change due to conduction, advection, firn densification (Patterson, 1994, Chapter 10, Eq. 32) and the correction to the conduction term due to spatially varying thermal conductivity respectively.We use a linear relationship for the specific heat capacity, namely C = 152.5 + 7.122T (Patterson, 1994, Chapter 10, Eq. 1).The thermal conductivity in ice is taken as κ = 9.828 exp −5.7 × 10 −3 T , and as per van Ommen et al. (1999) the thermal conductivity of the firn is the arithmetic average of the empirical fits of Van Dusen (Patterson, 1994, Chapter 10, Eq. 3) and Schwerdtfeger (Patterson, 1994, Chapter 10, Eq. 4): The firn densification term in Eq. ( 2) accounts for the work required to strain the firn and is given by (Patterson, 1994, Chapter 10, Eq. 38).
We applied Sorges law (time invariant firn density profile) to piecewise exponential plus linear or dual exponential (cross-over at 57 m) empirical fit to the measured density profile (see Fig. 4a), extrapolated beyond the depth of the Mill Island ice core using densities from the deep ice core at Law Dome summit (van Ommen et al., 1999).This particular functional form (piecewise exponential plus linear or dual exponential) was chosen based on previous analysis on the Law Dome ice core density profile (van Ommen et al., 1999) and the similarity between the density profiles at Mill Island and Law Dome (Fig. 4a).
A finite difference implementation of Eq. ( 2) was used to simulate the temperature profile with depth under the influence of a time-varying surface temperature.The finite dif-ference scheme was second order in space and first order in time with a forward time discretisation used to simulate 130 yr, with the initial 10 yr discarded to minimise the impact of transients associated with inconsistent initial conditions.The influence of the size of time step used for the forward time discretisation on the calculated temperature profile with depth was minimised by using an adaptive time step.In particular, the time step was successively halved until the absolute difference in simulated temperatures at all simulated depths calculated using the current time step and a time step twice as large was less than 5 × 10 −3 K.

Initial and boundary conditions
The thermal boundary conditions applied to the model were a time-varying, prescribed surface temperature (the optimisation of which is the purpose of the inverse models), and a zero steady-state heat flux boundary condition at depth (the depth at which the upward conduction of the geothermal heat flux is balanced by the downward advective heat flux).The optimal surface temperature history was found to be essentially independent of the location of this bottom boundary condition for depths in excess of 180 m below the surface (see Sect. 3.3).
The vertical velocity distribution with depth is unknown, but because the borehole is relatively shallow compared to the likely icecap thickness, the variation of age with depth for velocity profiles spanning the likely range of icecap thicknesses is small (see Fig. 3 and Sect.3.3).Therefore the resulting reconstruction of surface temperature history is essentially independent of reasonable assumed velocity profiles.We assume that the surface vertical velocity is equal to the recent snow accumulation rate of 2.93 m yr −1 at a surface snow density of 448 kg m −3 .We assume a Nye style velocity profile (constant vertical strain rate) and correct for firn density.Assuming that the base of the ice sheet is approximately at sea level gives a maximum depth of around 400 m, and finally the depth that the Nye style vertical velocity profile acts over (herein referred to as the "effective Nye depth") is adjusted by a factor of 5/6 to ensure consistency between the vertical strain rate and the vertical velocity at the surface, yielding a maximum effective Nye depth of 333 m.This was the default velocity profile used; a second profile with an effective Nye depth of 200 m was used to test the model sensitivity to this parameter (see Sect. 3.3).Figure 4b shows these two velocity profiles.This assumed velocity profile is consistent with initial analysis of annual layer thickness derived from chemical concentration measurements of the ice core, and in particular the thinning of the annual layers near the bottom of the ice core is not suggestive of a significantly thinner icecap.
The initial temperature profile down the borehole is chosen to be isothermal, with a value equal to the prescribed surface temperature for the start of the simulation period.This initial temperature profile was chosen to simulate steady conditions prior to the simulation period.Furthermore, due www.the-cryosphere.net/7/263/2013/The Cryosphere, 7, 263-273, 2013 to the relatively (by East Antarctic standards) high accumulation rate and warm temperatures, the problem is strongly dominated by the advective component of the heat flux.The advection parameter is of order 20 (the absolute value of which is the Péclet Number, Patterson, 1994), indicating that the system is relatively insensitive to downstream conditions and therefore also insensitive to the initial temperature profile down the borehole.Furthermore, this initial temperature distribution down the borehole varies for each simulation as the optimal solutions for both the LSQR and PSO methods include optimisation of the initial temperature distribution down the borehole.The initial choice of the surface temperature history differs in both form and values for the LSQR and PSO methods, which are detailed in Sects.3.2.1 and 3.2.2respectively.

LSQR method
The LSQR method uses an iterative greedy algorithm which reduces the largest temperature residual in each iteration by adding a time-varying update T j to each point in the surface temperature time history T j .This residual minimisation is based on a local linearisation of the problem.Specifically, the sensitivity matrix S ij is given by where the response of the temperature is Tmodel (i) at the depth i of the largest residual to the surface temperature history T j at time j is calculated using the complex-step derivative approximation (Martins et al., 2003).
Due to the diffusive nature of the thermal regime, the minimisation problem is underdetermined, with multiple time points in the surface temperature history (i.e.multiple T j 's) contributing to the simulated temperature at any depth ( Ti ).Therefore multiple updates to the surface temperature history that minimise E rms are possible, and the least squares QR method (Paige and Saunders, 1982) is used to select the particular update with minimum variance of the update T j (minimises the Euclidean norm j T 2 j 0.5 ).As the calculation of S ij is time intensive, this matrix is kept constant for multiple iterations.This results in some updates actually increasing E rms .Rather than discarding such updates, the surface temperature history corresponding to the lowest E rms is recorded.The initial choice of the surface temperature history was a linear interpolation of the measured temperature profile with the depth mapped to time using the assumed vertical velocity profile (see Fig. 3).

PSO method
Particle swarm optimisation is an optimisation method well suited for non-linear problems with multiple minima and a large dimensional subspace to search (i.e. a large number of variables to optimise).It works in an iterative manner by evaluating a number (or swarm) of candidate solutions (or particles) against a selection criterion and then updating each candidate's location in the search subspace based on a separate selection criterion for the entire swarm.Each particle is initially assigned a random location in the search subspace and also a randomised velocity (change in subspace location for each iteration) with an associated inertia which acts to keep the particle moving along its current velocity trajectory.After the evaluation of the selection criterion for every particle at the current iteration, the velocity of each particle is nudged to move the particle toward the location where the best selection criterion has been found for the entire swarm.This combination of particle inertia and a force attracting each particle toward the location corresponding to the best swarm selection criterion results in a broad investigation of the search subspace with a focus on areas corresponding to good overall selection criteria.
Specifically, this study used a swarm of 15 particles that traverse the parameter space being drawn towards the best solution (lowest E rms ) that any of the particles has visited (Pedersen and Chipperfield, 2010).In order to suppress unresolved high frequency variations, a piecewise linear surface temperature history was used, with four linear segments resulting in an eight-dimensional search subspace (an initial and final surface temperature and three intermediate time points where both temperature and time are free to evolve).Once any of the solutions has achieved an RMS error below a specified threshold, the best ground surface temperature history corresponding to the best RMS error is saved, and a new swarm generated with a new randomised initial state.Five hundred such swarms were simulated for varying specified RMS thresholds, with a RMS threshold of 0.05 • C being an optimum (see Sect. 3.3), and the median and distributions of these swarms analysed.

Numerical convergence
Numerical convergence was tested using an inverse model to calculate the required surface temperature history corresponding to the observed temperature profile, specifically using the the LSQR solution method (see Sect. 3.2.1).In particular the four key assumptions of the numerical scheme (numerical spatial resolution, effective Nye depth, the depth of the bottom zero heat flux boundary condition and the parametrisation of the upper firn thermal conductivity) were tested for the influence on the reconstructed surface temperature history.For these tests we will exclude the more recent decade of the reconstruction, where we have a lack of data to constrain the solution for this period (see Sect.Grid independence was tested using a 180 m computational domain with either 2 m or 5 m discretisation.The RMS difference between the reconstructed temperature histories was 0.057 • C. Most of the variation was in the upper part of the profile, where the approximately 2 m spacing of the measurements meant that, for the 5 m discretisation case, a sensitivity matrix could not be generated for each of the observational depths independently. The influence of the location of the bottom zero flux boundary condition was investigated using a 2 m spatial grid with domains of 180 m, 250 m and 350 m depth and the zero heat flux boundary condition applied at the bottom of the computational domain (see Fig. 5).The RMS temperature difference between the 180 m and 350 m domains is 0.040 • C, while between the 250 m and 350 m domains this is reduced to 0.015 • C.So we conclude that there is a slight dependence on the location of the zero heat flux boundary condition for shallow computational domains, but for domains of at least 250 m extent the reconstructed ground surface temperature becomes independent of the computational domain size.
The influence of the effective Nye depth was studied using a 2 m spatial grid with a computational domain of 180 m and varying the velocity profile to have an effective Nye depth of either 200 or 333 m.The shallow computational domain was used to avoid unrealistic vertical velocities at depth as a result of the 200 m Nye depth.The respective surface temperature histories differ by 0.039 • C, with the deeper effective Nye depth being slightly cooler throughout most of the history.
The sensitivity of the PSO solution to the specified RMS threshold is shown in Table 3.As expected, the RMS error of the median solution decreases with a decreasing RMS threshold, but for low RMS thresholds the resulting distribution from the swarms is noticeably skewed and the maximum error in the median solution increases for these skewed distributions.Therefore the optimal value for the RMS threshold was chosen as 0.05 • C. The thermal conductivity of snow is quite variable and dependent on the snow micro-structure (Strum et al., 1997) and evolves over time (Zhang, 2005;Domine et al., 2007).Such variability in thermal conductivity in the upper firn column should have negligible influence on the reconstructed surface temperature history, as such variations only alter the spatial distribution of heat rather than the total heat content of the forward model (Sect.3.1).Such spatial variations will have diffused away along with the seasonal temperature cycle by the time the snow has reached 19.07 m, the depth of the first temperature measurement used herein to constrain the surface temperature history reconstructions.This insensitivity to the details of the upper firn thermal conductivity was confirmed by LSQR simulations where the thermal conductivity (κ modified ) in the upper 10 m is given in Eq. (3) (κ orig ) by where 583 is the density of the firn at 10 m.The RMS difference in the reconstructed surface temperature history between the control simulation using κ orig and the simulation using the modified thermal conductivity (κ modified ) was only 4.9 × 10 −3 K.

Surface temperature reconstruction
Two inversion methods (see Sects.3.2.1 and 3.2.2 for details) were used to select ground surface temperature histories that minimise the RMS error in the simulated temperature profile.Due to highly diffusive nature of the problem, no high frequency information is retained, and therefore it is not possible to reconstruct a unique solution.
The heating associated with firn densification was found to be a small but not insignificant term at this site.In particular, this term contributed a warming of up to 0.06 K, which is similar to the approximate 0.1 K in the top 50 m for the Crete ice core in Greenland (Johnsen, 1977).The least squares QR (LSQR) method produces a single solution that is in some sense a minimum variance solution (actually the sum of minimum variance additions to an initial solution), whereas the probabilistic particle swarm optimisation (PSO) method produces a distribution of possible solutions.The ground surface temperature history reconstruction using each of these methods is shown in Fig. 6.
As shown in Fig. 6, the reconstructed surface temperature history is consistent with a change in surface energy balance at the site, and there is around a 0.75 K temperature increase between 119 and 19 m (see Fig. 2).Key details of this site, such as the vertical velocity profile, have not been measured to date, but reasonable estimates of these values yield a timing of this surface energy budget change of circa 1980/81 AD ±5 yr.This result is robust to changes in the vertical velocity distribution to values which are probably at the bounds of what is likely.Specifically, there is little difference in the time required for ice to move from the surface to around 60 m depth for either of the velocity profiles considered; the resulting age profiles down the borehole only start differing significantly at greater depths (see Fig. 3).
In general the inverse surface ground temperature reconstructions from the two methods (LSQR and PSO) agree with the observations (see Fig. 7) and with each other (see Fig. 6), although the LSQR method tends to have smoothed out the transition around 1980 AD and overestimated the very recent (last 5 yr) temperature increase.This is, at least in part, due to the least variance (in a 2-norm sense) nature of a LSQR solution for under-determined problems (Golub and Loan, 1990).
The PSO solution (Fig. 6) shows a fairly constant median surface ground temperature history up until circa 1980 AD and then shows an approximately linear increase post 1980 AD at a rate of 0.37 K decade −1 .The PSO method   used piecewise linear reconstructions of the surface temperature history (with 4 linear segments) and was free to evolve to include large magnitude but short duration variations.Such variations would be fairly quickly diffused away.The spread of the PSO surface temperature reconstruction (as shown by the percentile bands in Fig. 6) is a reflection of this potential for short-term variability, the additional time for thermal diffusion to act deeper into the borehole, and the relative lack of constraining observations deeper into the borehole.This last point is due to a combination of both the physical spacing of the temperature observations and the non-linear mapping between depth and age, so that the same physical distance represents longer epochs deeper into the borehole.

Discussion
Several processes both local or regional are likely to have contributed to the changed surface energy budget at the borehole site.In particular, possible mechanisms contributing to the recently increased surface energy budget include an increase in air temperature due to local-or regional-scale processes, a change in the radiative balance, a change in the annual duration of fresh snow cover and latent heat processes.Roberts et al. (2001) demonstrated regional-scale coherence between coastal Antarctic sites that span Mill Island, in particular between multi-decadal averaged Law Dome summer temperatures and a similarly averaged moisture balance for Ace Lake in the Vestfold Hills.Furthermore, van Ommen and Morgan (2010) show a similar spatial-scale climate teleconnection, in this case between Law Dome annual accumulation and southwest Western Australia rainfall.In addition, the multi-decadal high accumulation at Law Dome starting circa 1970 AD is unique in the last 750 yr (van Ommen and Morgan, 2010).Together, this shows that suitable regional-scale spatial climate coherence exists over this region, and that current climate conditions are, in at least some respects, unusual over the last several hundred years.Therefore, while local processes (see below) might be directly responsible for the changed surface energy budget, ultimately regional-scale changes are likely to be contributing to these changed local processes.
More local-scale processes that will directly influence the air temperature are changes in air/sea exchanges, which are strongly moderated by the presence of either sea ice or a more substantial ice shelf.Mill Island is on the periphery of the Shackleton Ice Shelf, so any changes in its extent have the potential to significantly alter the local atmospheric conditions.Specifically, there has been an approximate 17 % decrease in the area of the Shackleton Ice Shelf between 1956 and 2006, with most of the changes occurring in the region bounded by Mill Island and Bowman Island (Young and Gibson, 2013), where the ice shelf has retreated to the Antarctic coast between these two islands.This local driver is not, in itself, a complete explanation as it still requires a cause for the ice shelf retreat.
Changes to the amount of cloud cover at the site or the timing of cloud coverage on both diurnal and annual time scales have the potential to contribute to the heat budget.In particular, an increase in either nighttime cloud cover or atmospheric water content would decrease the radiative losses at this site and would contribute to annual average warming.Due to the limited period of availability of satellite data, which does not extend to before the observed change in the surface energy budget, this factor can not be directly assessed, although analysis of atmospheric reanalysis products would offer some insight.
Observations made by the field party drilling the ice core over the 2009/2010 austral summer indicate that summer daytime heating of the open ocean surface produced evaporation and a sea breeze that resulted in a marine fog layer at near surface level over the summit that persisted during the afternoon into the evening.A rime coating covered the surface, hence raising the liquid water deposition at the site.However, the marine fog reduced daytime incident radiation and surface melting and early evening radiative loss.In contrast, on days when the sea breeze not occur, a combined geostrophic wind from the southerly quadrant prevailed and kept the surface cloud-free.Under these conditions, surface melting occurs in summer.
Key physical properties of snow alter as the snow ages.Of direct relevance to the surface energy budget are decreases in both surface albedo and thermal conductivity as snow ages (Zhang, 2005).A decrease in albedo results in more direct absorption of incoming radiation, while changes to thermal conductivity influence the insulating effect between the atmosphere and the deeper firn column.Changes in the frequency of snow fall events (and hence the annual duration of fresh snow cover), and other processes influencing the surface characteristics of the snow may therefore change the surface energy balance.
Finally, latent heat processes related to the freezing of liquid water deposited at the site have the potential to significantly alter the energy budget.In particular due to the relative magnitudes of the latent heat of freezing and the specific heat content of firn, freezing of liquid water has the capacity to raise 160 times its own mass of firn by 1 K (Patterson, 1994).Therefore the observed 0.75 K temperature rise over the bottom 100 m of the borehole could be accounted for by an additional 6 kg m −2 yr −1 of liquid water deposition at the site.Potential sources for the liquid deposition include, but are not limited to, liquid deposition associated with the marine fog (indicated by corresponding rime coatings) or alternatively from wind blown sea-spray.Salinity levels at the site are relatively high compared to other Antarctic sites, with mean salinity levels around 5 ppm.Therefore small changes in the mode of delivery of this salt to the site (i.e.direct sea spray compared with salt dissolved in snowfall) have the potential to contribute significantly to the altered surface energy budget.
It is interesting to consider these local changes in the context of wider regional changes.The borehole temperature trend of 0.37 K decade −1 since 1980 can be compared with trends at Mirny (0.54 ± 0.64 K decade −1 for 1987 to 2011) and Casey (0.35 ± 0.97 K decade −1 from 1989 to mid-2012).Data from the nearby GF08 inland automatic weather station (68 • 29 36 S, 102 • 10 32 E, 2123 m elevation) for the period October 1986 to January 1998 also show a warming trend in both 4 m air temperature (1.12 ± 0.40 K decade −1 ) and 3 m subsurface firn temperature (0.38 ± 0.04 K decade −1 ).The time period chosen for Casey reflects data availability after repositioning of the station, while the time period for the Mirny analysis corresponds to the start of the GF08 record.For both Casey and Mirny the trends, while of comparable size, have large uncertainties.While the trends at GF08 are more significant, we note that both the air and subsurface senor height/depth were progressively lowering/deepening due to snow accumulation.The Casey and Mirny temperature changes are accompanied by increases in relative humidity.A piecewise linear data-fit allowing for a break in slope using the BREAKFIT algorithm (Mudelsee, 2009) shows a relative humidity trend at Casey for data from 1960 to mid-2012 of −0.38 ± 0.05 % yr −1 before 1989 and a trend of 0.51 ± 0.07 % yr −1 after.A corresponding analysis for Mirny over the period 1956-2011 shows a trend of 0.09 ± 0.09 % yr −1 before 1985 and a trend of 0.20 ± 0.11 % yr −1 after this date.
This picture of regional warming in East Antarctica is consistent with a recent temperature reconstruction (Steig et al., 2009) and also supports model simulations (Goosse et al., www.the-cryosphere.net/7/263/2013/The Cryosphere, 7, 263-273, 2013 2012, Fig. 3).Further, these results can be set in the context of wider Antarctic warming, with comparable Antarctic Peninsula warming (Mulvancy et al., 2012), a reconstructed recent temperature trend of up to 0.8 K decade −1 for the West Antarctic Ice Sheet divide (Orsi et al., 2012) and large but spatially variable temperature trends in inland Dronning Maud Land, East Antarctica (Muto et al., 2011).

Conclusions
The temperature profile down an ice borehole into the Mill Island icecap shows a distinct change of gradient around 49 m below the surface.Numerical modelling suggests that this temperature profile is due to a change in the surface energy budget that occurred around 1980/81 AD ±5 yr.This increase in the surface energy budget since circa 1980 has resulted in a warming of around 0.37 K decade −1 , although care must be taken in interpreting this effective ground surface temperature reconstruction.Factors that may have contributed to this trend of increasing ground surface temperature include an increase in air temperature from either local or regional influences, a change in the surface radiation budget and an increase in the amount (in either absolute or relative terms) of liquid water deposition.Comparison of the Mill Island temperature trend with instrumental data suggests that these trends are regional in nature and not dominated by local-scale effects.

Fig. 3 . 23 Fig. 3 .
Fig. 3.The influence of assumed vertical velocity profiles on the age distribution down the M borehole.

ρ
Fig. 4. (a) Density and (b) velocity profiles used in the ground surface temperature reconstruction.

Fig. 5 .Fig. 5 .
Fig. 5. Sensitivity of the reconstructed ground surface temperature reconstruction to the depth of the computation domain, and hence the position of the zero heat flux boundary condition.

Fig. 7 .
Fig. 7.The temperature profile down the Mill Island borehole by the forward model for the L median PSO surface temperature reconstructions. 27

Fig. 7 .
Fig. 7.The temperature profile down the Mill Island borehole by the forward model for the LSQR and median PSO surface temperature reconstructions.

22 Fig. 2. Sub
-surface observed temperature profile in the Mill Island borehole, also shown is the piecewise t to the data below 15 m.Open circles denote measurements in the seasonally varying zone iscarded for in this analysis.-surface observed temperature profile in the Mill Island borehole, also shown is the piecewise line data fit to the data below 15 m.Open circles denote measurements in the seasonally varying zone that were discarded in this analysis.
• C
simulates the down borehole temperature profile given a surface temperature history, taking into account advection, thermal diffusion and the heat associated with firn densification.

Table 3 .
Sensitivity of PSO median solution to the RMS threshold.