Journal of Geophysical Research : Earth Surface Snow Densification and Recent Accumulation Along the iSTAR Traverse , Pine Island Glacier , Antarctica

Neutron probe measurements of snow density from 22 sites in the Pine Island Glacier basin have been used to determine mean annual accumulation using an automatic annual layer identification routine. A mean density profile which can be used to convert radar two-way travel times to depth has been derived, and the effect of annual fluctuations in density on estimates of the depth of radar reflectors is shown to be insignificant, except very near the surface. Vertical densification rates have been derived from the neutron probe density profiles and from deeper firn core density profiles available at 9 of the sites. These rates are consistent with the rates predicted by the Herron and Langway model for stage 1 densification (by grain-boundary sliding, grain growth and intracrystalline deformation) and stage 2 densification (predominantly by sintering), except in a transition zone extending from ≈8 to ≈13 m from the surface in which 10–14% of the compaction occurs. Profiles of volumetric strain rate at each site show that in this transition zone the rates are consistent with the Arthern densification model. Comparison of the vertical densification rates and volumetric strain rates indicates that the expected relation to mean annual accumulation breaks down at high accumulation rates even when corrections are made for horizontal ice velocity divergence.


Introduction
Better estimates of the contribution to sea level rise from Antarctica require an improvement in the meteorological and snow densification models used to interpret satellite measurements of elevation change. The UK Natural Environment Research Council (NERC) research program iSTAR provided an opportunity to collect data from the Pine Island Glacier (PIG) basin to verify these models and hence provide an improved estimate of sea level rise from this sector of Antarctica.
In this paper we focus on the analysis of snow density profiles collected using a neutron probe at 22 sites along an 890 km traverse of the basin and repeat measurements made a year later. These allow recent accumulation to be determined from the snow stratigraphy and profiles of volumetric strain rate to be calculated.
Accumulation has been successfully determined from neutron probe density profiles in the dry snow area of Greenland (Morris & Wingham, 2011, 2015, but in the PIG basin the method faces a new challenge: complex meteorological conditions which do not always produce a clear annual variation in snow density. New techniques of analysis have therefore had to be developed. The PIG strain rate profiles are of particular interest because they come from sites with a wide range of mean annual accumulation and relatively high surface density. The transition density between stage 1 and stage 2 densification in the Herron and Langway (1980) model (and many others, see section 3) is reached about halfway down the profile, so that the transition region can be examined in some detail. We are also able to consider the effect of horizontal velocity divergence on densification using data from mid-ice stream sites.
In this paper we compare direct measurements of the change in density with time with the rates inferred indirectly from density profiles, assuming steady state conditions. This indirect method is based on an early suggestion made by Robin (1958) that in a natural snow cover the increase in overburden stress is linearly related to the proportional decrease in air volume in the mean density profile. Much of the existing data Journal of Geophysical Research: Earth Surface 10.1002/2017JF004357 on strain rates in polar snow is derived from profiles using the Robin hypothesis. We show that at the high accumulation sites in the PIG basin this method does not produce good estimates in the transition region between stage 1 and stage 2 densification. This is important because, as we shall show in section 6, in the PIG basin ≈31% of compaction occurs in the surface layer and stage 1 region, 10-14% in the transition region, and 37-33% in the stage 2 region. Thus, compaction in the transition region is a significant proportion of the total compaction as the surface snow turns to ice.

The iSTAR Traverse
The iSTAR traverse route (Figure 1) crosses the Pine Island Glacier basin from site 1 at the northeast upstream boundary to site 22 near the grounding line in Pine Island Bay to the west. The positions and climatological characteristics of the sites are given in Table 1, which gives values for the mean annual accumulation derived from airborne radar survey data  over the period 1985-2009 and the regional atmospheric climate model RACMO 2.3 (van Wessem et al., 2014) over the period 1979-2013. The RACMO2.3 values for the iSTAR sites have been determined by interpolation between model grid points. For convenience the table also includes the mean annual temperatures (from 10 m temperatures measured during the traverses), accumulation rates for the period between traverses and mean annual accumulation rates derived from iSTAR traverse data in section 4.2. The latter cover periods of 8 to 29 years up to and including 2012.
Mean annual temperature only varies by a few degrees over the traverse, but the accumulation increases significantly in the southwest of the basin. Satellite altimeter observations show elevation change rates also vary considerably over the basin, with the greatest thinning rates in areas of high ice velocity (Shepherd et al., 2001). In order to distinguish between dynamic thinning and the effect of changes in accumulation on elevation, the traverse route was chosen to pass through sites with low accumulation and low thinning (sites 2 and 10), low accumulation and high thinning (sites 7 and 13), high accumulation and low thinning (site 22), and high accumulation and high thinning (site 18).
In the austral summer of 2013/2014 (T1) the route was traversed from site 1 to site 22; in 2014/2015 (T2) the traverse began at site 22. Snow density profiles were measured at all sites during T1 using a neutron probe deployed in 13 m deep access holes. During T2 repeat measurements were made in the same access holes, and new profiles were measured nearby in order to record the density of the snow accumulated since T1. At sites 7,8,13,15,17, and 18 density profiles were collected over a nested grid to establish local spatial variability on the 1 m to 1 km scale. During T1 ground penetrating radar (GPR) and very high bandwidth radar (pRES) surveys were conducted between sites, and during T2 firn cores to 50 m depth were collected at sites 1, 4, 6, 7, 8, 10, 13, 15, 18, and 20. Airborne surveys using the ASIRAS radar in December 2014 covered the traverse from site 7 to site 22. Other airborne survey tracks covered areas north of site 9, south of site 15 and passed between sites 4 and 5 on a north-south line.
The neutron probe profiles serve to supplement the firn core data and provide the density data needed to derive accumulation from the various radar measurements of snow layering. This work will be reported elsewhere; in this paper we focus on the derivation of densification and accumulation rates directly from the neutron probe profiles. The densification analysis requires an estimate of the horizontal velocity divergence at each site. Direct measurements of this divergence are available from the repeat measurement of strain networks during T1 and T2 at sites 6 to 9 and 13 to 19. Estimates can also be made from velocity maps for the PIG basin derived from satellite data.

The Steady State Approximation
Snow densification models need to provide a constitutive law to describe the dependence of volumetric strain rate,̇, on the applied stress, , and on the strength of the snow, characterized by its bulk density, temperature, and sometimes by other variables such as grain size, coordination number, and impurity concentration. The material-following densification rate is related tȯby the definitioṅ  (Bamber et al., 2009), the background color shows surface velocity (Rignot et al., 2011) and the basin boundary is shown as a blue solid line. Note. x is the distance along the traverse, T * m is the mean annual temperature in degrees Celsius (Mulvaney & Smith, 2017 is the accumulation rate between Traverse 1 and Traverse 2, andā is the mean annual accumulation rate estimated using (1) airborne radar data , (2) the regional atmospheric climate model RACMO 2.3 (van Wessem et al., 2014), and (3) the iSTAR density profiles (see section 4.2). N is the number of years over whichā (3) is averaged.

10.1002/2017JF004357
The variable is the snow density and t time. In polar firn there are density fluctuations on an annual or subannual scale about a vertically smoothed density 0 , which monotonically increases with depth. In the upper layer of the snow, densification is enhanced by annual and subannual temperature variations, but below ≈2-3 m depth the effect of these variations is sufficiently diminished for it to be possible to assume a constant temperature to a first approximation. The steady state approximation applies below the surface layer for constant mean annual temperature T m and constant accumulation rateā.
The Robin hypothesis (Robin, 1958) where i is the density of ice and c is a site-specific constant. Hidden within this formulation is an assumption that the overburden stress , which increases with water equivalent depth, is offset by a corresponding increase in snow strength because of some other factor that also increases with depth, such as the time since deposition, . Hence, at a given site, the volumetric strain rate depends only on the mean density and c can be interpreted as a constant density-corrected volumetric strain rate.
In the one-dimensional case when the horizontal velocity divergencėH is negligiblė=̇z z wherė We define a coordinate system in which vertical distance, z, velocity, w, and water equivalent height, q, are positive upward. In this system density decreases with increasing z and q. Given a constant accumulation ratē a = − 0 w, equations (2) and (3) give The value of c can be estimated from a single steady state profile 0 (z) by plotting the best fit straight line through ln( 0 ∕( i − 0 ) as a function of z over a section of the profile, say from the point z where the water equivalent height is q to the point z + Z where the water equivalent height is q + Q. Writing the (negative) gradient of this line as −k i we obtain where k is a positive constant. To emphasis the point that k is derived from the rate of change of density with the vertical coordinate z, we shall refer to it as the "vertical densification rate." Equation (5) shows that, if we accept the Robin hypothesis, in the 1-D, steady state case the density-corrected volumetric strain rate c may be calculated from the vertical densification rate if the (constant) accumulation rate is known. A similar calculation could be made given any other specified density function, of course, but the Robin hypothesis leads to a particularly simple transformation.

Time-Varying Conditions
Suppose that in a natural firn layer, built up under time-varying conditions, the effect of density on strain rate is given by a factor ( i − )∕ by analogy with equation (2), that is, the Robin hypothesis still holds. Then given measurements of density 1 (q) for a material element at time t and 2 (q) for the same element at time t + Δt, we may calculate a density-corrected volumetric strain rate as To minimize the effect of density fluctuations we calculate an average value of F over a section of the profile from q to q + Q where Q is several years of accumulation. Then Let −k i now be the gradient of the best straight line through a plot of ln( 1 ∕( i − 1 ) as a function of z over the section of the profile from q to q + Q and letā be the mean annual accumulation over this section. The indirect method of estimating strain rate from a single profile is based on the assumption thatF by analogy with the steady state case (equation (5)). In this paper we determineF,ā, and k separately to see how far this assumption holds. Herron and Langway (1980) proposed different densification rates for snow with <550 kg m −3 and 550 kg m −3 ≤ ≤ 800 kg m −3 . Other authors have followed suit. This implies a sharp transition between "stage 1 densification" for the lower densities, in which grain boundary sliding and grain growth are thought to be the dominant mechanisms and "stage 2 densification" for the higher densities in which sintering is the primary mechanism. The proposed transition density is the maximum packing density of uniform spheres of ice. However, in natural snow, with a range of grain sizes, high-resolution measurements of density profiles (Hörhold et al., 2011) suggest that the transition density varies between sites and that the transition may be gradual. For this reason it is useful to define a transition region between the stages 1 and 2 regions and to use the subscripts 0, t, and 1 to distinguish values of the vertical densification rate k. (This has the advantage of being a simple and consistent notation, but means that our symbol k 1 is not the same as the k 1 used by Herron & Langway, 1980 for stage 2 densification). We envisage the transition as a gradual change in the balance of densification processes from those dominant in stage 1 to those dominant in stage 2.

Stage 1 and Stage 2 Densification
In this paper we shall allocate densities greater than 550 kg m −3 from depths greater than 15 m below the surface to the stage 2 region. Densities greater than 550 kg m −3 from depths that are less than 15 m below the surface will be regarded as potentially transitional. Note that this choice of how to bin the data does not imply that there are "hard" boundaries to the transition region in reality; it is a convenience for the sake of data analysis. Herron and Langway (1980) suggest an Arrhenius expression for the local (site-specific) stage 1 vertical densification rate k 0 . Here R = 8.314 J mol −1 K −1 is the gas constant, E 0 = (10.16 ± 0.94) kJ mol −1 is the stage 1 activation energy, T m mean annual temperature in degrees Kelvin and k * 0 = 11(+6, −4) (metre water equivalent (mwe)) −1 is a constant global stage 1 vertical densification rate. Equation (9) does not necessarily apply in the near-surface layer with strongly varying temperature. For stage 2 densification, Herron and Langway (1980) propose where E 1 = 21.40 kJ mol −1 is the stage 2 activation energy and the constant k * 1 = 575 mwe −1∕2 a −1∕2 . IfF = F 0 in the stage 1 region and F 1 in the stage 2 region, the Robin hypothesis implies

Journal of Geophysical
and Other authors have suggested modifications to the Herron and Langway (1980) model. For example, Arthern et al. (2010) suggest that both k 0 and k 1 are independent ofā from consideration of the physics and find a good fit to Antarctic strain rate data from sites withā = 0.13-1.04 mwe a −1 using for stage 1 and for stage 2. All the iSTAR sites lie within this mean annual accumulation range. Ligtenberg et al. (2011) found it necessary to add empirical tuning parameters containingā to the Arthern model in order to match profiles from Antarctica over a wider range of accumulation. They suggest multiplying the right-hand side (RHS) of equation (13) by the ratio and the RHS of equation (14) by Since these tuning parameters have been derived using data from a wider range of climatological conditions than the data used to derive the Arthern model itself, there is an element of extrapolation involved.

Journal of Geophysical Research: Earth Surface
10.1002/2017JF004357 Figure 4. Estimates of mean annual accumulation rate given by the automatic layer identification routine StratiCounter as a function of the estimates derived by manual layer identification alone. Points lying below the 1:1 line arise when StratiCounter detects a layer missed by the manual method, which is used to provide initial estimates of the annual layering. The vertical bars show the uncertainty in the mean annual accumulation arising from the uncertainty in the number of layers determined by StratiCounter.

The Effect of Horizontal Velocity Divergence
We now consider the densification of snow above a horizontal ice surface subject to strain. At the ice/snow boundarẏ since the ice is incompressible. This horizontal velocity divergence is imposed upon the overlying snow and is constant with depth in the snow layer. If the main density peaks lie initially at water equivalent heights q 1 (i) and, after time interval Δt, at depths q 2 (i) conservation of mass in a snow column initially of unit area gives Given a suitable choice of origin for q 2 the gradient of the best straight line fitted through the points q 1 (i) − q 2 (i), q 2 (i)Δt can be used to estimate the horizontal velocity divergence from the neutron probe data albeit with an uncertainty which, for the iSTAR sites, is mostly of the same order aṡH. Fortunately, for the PIG basin more precise estimates are available, either from direct measurement of strain networks or from velocity maps derived by remote sensing. Horizontal velocity divergences range from −36.7 to +89 ⋅10 −4 a −1 .
In the presence of horizontal velocity divergence the accumulation in an annual layer years old is reduced by the factor 1∕(1 +̇H ). This factor is not significant for the time scale of the neutron probe data, but the deeper firn core annual accumulation rates do need to be corrected.
The density-corrected vertical strain rate, F z , may be estimated by subtracting a correction for the effect oḟH: where m is the mean density over time Δt. That is, for a given density change the vertical compression of a volume element increases in magnitude if there is horizontal loss of mass from the element.
In the steady state case the constant accumulation rateā≈− 0 w(1+̇H ) and substitution for w in equation (3) no longer leads to equation (4). We, therefore, do not necessarily expect thatF z can be estimated from −āk( 1 ) unless is sufficiently small. For example, given a horizontal divergence of 0.01 a −1 , foṙH to be an order of magnitude less than 1, must be less than 10 years.

Errors
Since neutron emission is a random process, there is a random error in density measured using the neutron probe, which in our case is ≈ 0.02 . From the neutron probe calibration equation, the relative systematic errors in 2 and q 2 , which arise if the nominal value R for the access hole radius has a systematic error, * R can be estimated as * If the radius of the access hole changes over time Δt because of horizontal velocity divergence * R R ≈̇HΔt (21) FoṙH ≈ 10 −3 a −1 and Δt ≈ 1 a, the relative systematic error in density and water equivalent height arising from horizontal velocity divergence is ≈ 4 ⋅ 10 −4 , that is, very much smaller than the random error. However, larger systematic errors can arise if the access hole becomes enlarged by careless drilling or if the probe does not remain resting against the sidewall of the hole.
F has a random error, F , which arises from the fluctuations in F, and a systematic error * which is ≈0.6̇H for a representative value of 2 of 550 kg m −3 . Using this value the potential systematic errors are of the same order as the random errors forF. The random error iṅH is unknown so we simply set Fz = F .

Methods
As an example of the data obtained using the neutron probe, Figure 2 shows profiles of density at site 21 as a function of water equivalent height. The techniques involved have been described in a series of papers (Morris & Cooper, 2003;Morris, 2008;Morris & Wingham, 2011, 2015 so will not be repeated here. Very near the surface both snow and atmosphere are included in the measurement volume, so the apparent snow density decreases. Missing near-surface data are represented in the profiles by a constant density (vertical line).

Accumulation Rates in 2014
At a given site, the water equivalent Δq of snow accumulated over the time between measurements can be determined from the short T2 record, using the T1 record to help identify the position of the T1 surface. An error of order ± 0.01 mwe may arise because of the missing near-surface data and will be more significant at low accumulations.

Annual Accumulation Series
The stratigraphic method of determining annual accumulation series depends on the identification of an annual variation in the physical properties of the surface snow. Benson (1962) singles out the "fall (autumn) sequence" of coarse-grained, low-density snow, often containing hoar crystals, overlain by a finer-grained, harder layer of higher density and uses this discontinuity to mark the boundary of annual layers in dry snow areas of the Greenland Ice Sheet. Lister (1959) also uses a "criterion" layer of depth hoar as the annual marker in Antarctic snow, but finds that this is not present in some years while in others there are hoar layers both above and below a hard summer layer. He stresses the importance of examining the cyclical structure of the data as a whole, rather than relying on a single marker. Alley et al. (1997) comment that Benson's simple fall sequence is not always found in Greenland snow. There may be a series of depth hoar/surface hoar complexes formed during the summer, separated by finer-grained, denser snow. These authors therefore suggest that the annual marker should be taken at the center of the hoar complexes, that is, at the midsummer rather than fall surface.
Where there is a clear annual variation in density, as at site 21 (Figure 2), the density peaks can be selected for use as annual markers without too much difficulty (Morris & Wingham, 2011). Here late summer/autumn low-density hoar layers alternate with winter snow which has densified under the influence of warm summer temperatures and the mass between peaks is an estimate of accumulation from mid-winter (June/July) to mid-winter. However, at some sites along the iSTAR traverse, for example, at site 7, the stratigraphy is quite complex and difficult to interpret, possibly because there are competing coastal influences from east and west. In an attempt to improve the manual estimates of annual layering we used an automatic method to help ensure the criteria for layer identification were applied consistently.  Note. The parameters k 0 and k t are stage 1 and transitional vertical densification rates derived from T1 neutron probe profiles; k 1 is the stage 2 vertical densification rate derived from T2 firn core density profiles.̇H is the horizontal divergence derived from velocity maps or strain networks (in bold). F 0 and F t are means of the density-corrected volumetric strain rate and F z0 , F zt are means of the density-corrected vertical strain rate, with their errors, and * F the systematic error arising froṁH.

StratiCounter
The layer counting algorithm StratiCounter (release date 30/4/2015) was developed for ice core analysis (Winstrup et al., 2012;Winstrup, 2016) and allows a variety of data records (e.g., chemical species, isotopic ratios, and stratigraphy) on a common depth scale to be used together to determine the annual layering in a core. Some preprocessing of the density data is required before it can be used. We fit a fourth-order polynomial 0n to each of the n available profiles from a given site and then normalize so that the density peaks are not damped with depth. This produces n data series of the form which are input to StratiCounter.
The manual peak picks are read in as a first estimate of the layering and have a strong influence on the outcome. This is because StratiCounter looks for repeated patterns, not single peaks; it is up to the user to decide whether a typical annual layer should have a single peak or a more complex density variation. There is then an iteration process to determine the best positions of the annual markers. The standard deviation of the log-normal distribution of annual accumulation, a , is held constant during iteration, usually at the value determined from the manual peak picks.
As an example of the output from StratiCounter, Figure 3 shows the annual layers inferred by StratiCounter for site 4. Two normalized density profiles are available for this site, from T1 and T2. The first estimate of annual layering was made by identifying peaks in the T1 profile. The standard deviation in layer thickness is a = 0.24, within the expected range for annual accumulation determined from chemical species (0.15-0.3). The upper part of the figure shows the same curves, with the annual layering determined automatically by giving equal weight to each density profile. One extra layer has been identified at the bottom of the profile, and the position of some of the annual markers has changed to reflect peaks in the T2 record.

Journal of Geophysical
Site 4 is moderately difficult to interpret because of the multiple density peaks in some winters; it requires experience to recognize these by eye. At sites where density peaks have not formed in some years, it is even harder to determine a good manual estimate of layering from the density profile alone. Given a range of possible initial estimates StratiCounter may converge on different solutions for the layering. In this case further information from neighboring sites must be used to fix the number of layers expected between the surface and a recognizable peak at a given depth. This is akin to layer tracing in an accumulation radar return or identification of volcanic horizons in ice cores. StratiCounter is designed to be able to include such tie points in automatic analysis of ice cores, but the neutron probe density records are much shorter, so the extra information is simply used to select from competing solutions. This was necessary only at sites 4 and 7, where we also used the ice core chemistry to guide the initial estimate of annual layering.
The great advantages of the automatic method are that missed or spurious layers can be identified using a rigorous statistical method and that the uncertainty in the number of annual layers is defined. Hence, errors can be attached to the mean annual accumulation rates derived using StratiCounter. Note that layers are "missed" in the sense that a consistent criterion for layer selection has not been applied at the manual stage and StratiCounter has been able to detect this and rectify it by adding a layer. The converse is the situation in which StratiCounter subtracts a "spurious" layer to rectify an inconsistency. Using StratiCounter improves the internal logical consistency of the stratigraphic method but does not, of course, guarantee that the annual layers have been correctly identified. For this ice core data, ideally from more than one chemical species, are required. Figure 4 shows that the automatic rates are generally slightly less than those derived by manual peak picking, an indication that missed layers have been found. An uncertainty of one layer is, of course, more important at higher accumulations where the number of annual layers in the 13 m profile is smaller. However, this is offset by the fact that the higher accumulation regions have less complex stratigraphy and are easier to interpret. Hence, there is no correlation between the percentage error (+ or −) derived using StratiCounter and mean annual accumulation. Figure 5 shows the large-scale variation of density with distance along the iSTAR traverse. Each density column (derived from a single profile) is centered over a site and covers the distance between the midpoints with neighboring sites. Note that there are loops in the traverse path so that, for example, sites 8 and 12 are quite close together. Some of the annual layers can be traced between nearby sites, but there is no obvious feature that can be traced along the traverse. The nominal density for the transition between stages 1 and 2 is reached at 7-9 m below the surface and the density at the bottom of the profile (z = −13 m) varies from 600 to 625 kg m −3 .

Accumulation
We derive annual accumulation series for each profile by retaining the manual peak positions but adding or subtracting peaks to conform with the number of layers determined using StratiCounter. This allows us to retain a common definition of the annual marker between sites. The mean value over the profile is shown in Table 1 as is the 2014 accumulation rate a (2014) ≈ Δq/Δt.

Vertical Densification Rates, k
At most sites neutron probe profiles were obtained to around 13 m depth. We also have some firn core density profiles to 50 m. Local vertical densification rates k 0 , k t , and k 1 are given by the gradients of linear fits to profiles of ln( ∕( i − )), where we write k t for the possibly transitional value obtained from the lower part of the neutron probe profiles and k 1 for the Stage 2 value obtained from the firn cores. As an example, Figure 6 shows the fits for the neutron probe profiles for two sites with similar mean annual temperature but contrasting mean annual accumulation. Data from the upper part of each profile have been excluded from the Stage 1 fit, because here densification has been enhanced by warm summer temperatures. Each section to be fitted runs from density peak to density peak, to avoid skewing the gradient. However, there is clearly a potential for a systematic error in k 0 and k t , since the annual fluctuations in density are not regular and are relatively large compared to the range of the ordinate. The fits for all sites are shown in the supporting information as Figures S1 to S3.
Values of k are given in Table 2 with random errors derived by minimizing the sum of the deviations from the best straight line. Stage 1 values range from (704 ± 32) ⋅10 −4 mwe −1 at site 17 to (1,234 ± 49) ⋅10 −4 mwe −1 at site 5. Transitional values range from (372 ± 12) ⋅10 −4 mwe −1 at site 5 to (541 ± 18) ⋅10 −4 mwe −1 at site 19. Stage 2 values are lower still, ranging from (257 ± 8) ⋅ 10 −4 mwe −1 at site 20 to (458 ± 14) ⋅10 −4 mwe −1 at site 10. Figure 7 shows an example of the variation of F with water equivalent height for sites 11 and 22. The mean values F 0 and F t are calculated from peak to peak, over exactly the same ranges of water equivalent as used in the calculation of k 0 and k t . There is no obvious step change between the two values. Table 2 shows that mean values of F 0 range from −0.080 a −1 at site 15 to −0.022 a −1 at site 10. Transitional vales F t range from −0.055 a −1 at site 19 to −0.015 a −1 at site 10. The low values at site 2 may be inaccurate because this was the only site where 2 does not come from the same access hole as 1 . F z0 is significantly different from F 0 at sites 4, 7, 9, 12, 14, 16, 18, and 19 (higher horizontal velocity divergence) and site 2 where there is a possible systematic error in R. F zt is significantly different from F t at more sites (1,4,5,6,7,8,9,11,12,13,14,15,16,17,18,19, and 2), because the random error F decreases with depth and the fluctuations in density are attenuated.

Spatial Variation of Density
In order to analyze radar returns from continuous tracks across a basin, a single global density profile is often used (see, for example, Richardson et al., 1997). Such a profile, derived from the neutron probe profiles from sites 6 to 22, has been used to calculate two-way travel (twt) times for iSTAR GPR data (H. Konrad, personal communication). Fitting a third-order polynomial to this profile gives the global density profile  Figure 8 shows the deviation from G at each site. Sites 1-5 are not noticeably different from the sites used to produce G . The annual and subannual fluctuations in density are up to ±8%. However, the cumulative effect of these fluctuations with depth is much smaller. Using the polynomial leads to a maximum underestimate of twt of 0.44 ns, that is, ≈5 cm in depth for any of the 22 sites. This is significant only for location of reflections very near the surface.
Taking G = 385 kg m −3 at z = 0 as an estimate of the density of the surface snow, we may now calculate the percentage of total compaction that takes place in the surface layer and stage 1 (i.e., to ≤ 550 kg m −3 ) as 31% and in stage 3 ( i.e., ≥ 800 kg m −3 ) as 22%. Hence, stage 2 and transitional compaction form 47% of the whole. For practical reasons we have allocated neutron probe data with z > −15 m to the transitional region and firn core data with z ≤ −15 m to stage 2; the division occurs at = 600-625 kg m −3 . This leads to an estimate of 10-14% compaction in the transitional zone and 37-33% in stage 2.
Throughout the traverse the spatial variation is complex and it is clear that mean annual temperature and accumulation are not the only factors controlling it. For example, sites 10 and 11 have similar values of T m andā (Table 1) but at site 10 (mid-ice stream valley) the density is greater than that at site 11 (interstream ridge). Similarly, 20 and 21 have similar T m andā and again the valley site 20 has denser snow than the ridge site 21. In a multiple linear regression analysis of 40 measurements of surface snow density, Kaspers et al. (2004) found that T m contributed ≈70% to the surface density,ā ≈ 2.5%, and the annual average wind speed at 10 m ≈7.5%. If the iSTAR valley sites have an increased surface density because of higher (katabatic) wind speeds, the difference in these profiles can be explained. Figure 11. Mean annual accumulation rate derived from RACMO 2.3 data (van Wessem et al., 2014) compared to mean annual accumulation rate from the density profiles. The gradient of the best linear fit through the origin for all sites is 1.08 ± 0.05 (r 2 = 0.96). Sites above 1,000 m above sea level (asl) are shown in red.
Although we do not have measurements of surface density from the probe, an approximate value of S can be obtained for each site from third-order polynomial fits to the density profiles (z). The regression equation S = (1620 ± 460 kg m −3 ) − (5 ± 1.8 kg m −3 K −1 ) T m + (68 ± 14 kg m −3 mwe −1 a)ā r 2 = 0.58 (25) shows a sensitivity toā similar to the 67 kg m −3 mwe −1 a found by Kaspers et al. (2004) who used measured values of S . However, the range of T m over the traverse is small, and the (unknown) wind speed cannot be explicitly included in the regression, so we do not put too much weight on this equation.

Accumulation
The accumulation values derived from the density profiles will in due course be used with radar and core data to produce an accumulation map of the PIG basin. At that stage it will be appropriate to discuss the spatial and temporal variability of accumulation over the basin and make a full comparison with meteorological model data. In this section we need only to establish that the point values of accumulation at the sites where we have measured densification rates are reliable. This will allow us to proceed to test models of densification that use mean annual accumulation rate as a parameter. Vertical densification rates as a function of mean annual temperature, T m : (a) stage 1 rates, k 0 , estimated from neutron probe data with ≤ 550 kg m −3 , (b) transitional rates, k t , estimated from neutron probe data with ≥ 550 kg m −3 and −13 m ≤ z ≤ −8 m and corrected for accumulation by a factorā 1∕2 , and (c) stage 2 rates, k 1 , estimated from firn core density data with 800 kg m −3 > > 550 kg m −3 and z ≥ −8 m, also corrected for accumulation by the factorā 1∕2 . Straight lines show the Herron and Langway (1980) densification rates for stages 1 (blue) and 2 (red).
Hence, over the traverse as a whole the accumulation rate in 2014 is not significantly different from the mean. Bearing in mind that the accumulation rates are point values, it is not surprising that there is a fair amount of scatter. However, no one point stands out as anomalous. This is a useful independent check of the stratigraphic estimates of mean annual accumulation.
A nonzero intercept suggests a constant (absolute) difference betweenā (1) andā (3) . If we discount this, the best linear fit through the origin is (Note that the coefficients of determination r 2 and r * 2 are not calculated in the same way for regression with nonzero and zero intercept and so are not strictly comparable). Equation (27) is preferable on purely statistical grounds. However, on physical grounds equation (28) is perhaps more realistic. In either case there are no obvious outliers, so that comparison with the radar data is again a useful check of the stratigraphic accumulation rates.
suggesting either that the model slightly overestimates accumulation over the whole traverse or that accumulation was higher over the more recent period. Recalculatingā (2) over N years for each site, that is, matching the averaging period for model and field data, produces essentially the same gradient ( Figure 11). This suggests that the model estimates may indeed be a little too high, at least at some sites. We note that Medley et al. (2013) have shown that for the neighboring Thwaites Glacier basin, radar and RACMO2 mean annual accumulations are very close for the elevation range 1,000-1,400 m above sea level (asl). Restricting the fit to this range givesā (2) = (0.95 ± 0.03)ā (3) r * 2 = 1.00.
so that the model values for the higher-elevation sites 1, 2, 3, 20, and 21 are slightly too low. Overall, the model values support the stratigraphic values and there are no obvious outliers.
Having established reliable point values of accumulation at the iSTAR sites, we can now use our data to test various densification models. We first compare the iSTAR vertical densification rates with the predictions of the Herron and Langway (1980) model. This model has parameters derived by fitting vertical densification rates derived from firn core profiles so we may think of it as a "k model." We then compare the iSTAR volumetric strain rates with the Arthern et al. (2010) model, which has parameters derived by fitting strain rates (an "F model") and the Ligtenberg et al. (2011) model which adds further empirical parameters to the Arthern model. Since these parameters are derived by fitting the depth to transition densities in firn core profiles, this becomes a "hybrid model."

Vertical Densification Rates
We can test whether the local vertical densification rates from the iSTAR sites given in Table 2 are consistent with the Herron and Langway (1980) model by plotting ln(k 0 ) as a function of 1∕T m (see equation (9) for stage 1 and ln(k 1 )ā 1∕2 as a function of 1∕T m for stage 2 (see equation (10)). In the Herron and Langway (1980) model our transitional region is treated as stage 2 so we plot ln (k t )ā 1∕2 as a function of 1∕T m . Figure 12a shows that iSTAR values of ln(k 0 ) are consistent with the Herron and Langway (1980) stage 1 curve. The temperature range is too narrow, and the scatter too great, for independent values of k * 0 and E 0 to be determined from these data alone. However, it is worth noting that site 4 and sites 15-22 do not have anomalous values of k 0 , suggesting that the model can be used for sites with mean annual accumulations outside the 0.02-0.50 mwe a −1 range for which it was derived. Figure 12c shows that the values of ln(k 1 )ā 1∕2 derived from the firn cores are consistent with the Herron and Langway (1980) stage 2 curve and are less scattered, as might be expected as the length of the profile from which the values are derived is nearly 10 times greater than for stage 1. Again, the high accumulation sites do not have anomalous values. The values of ln (k t )ā 1∕2 (Figure 12b) lie between the stage 1 and stage 2 curves; hence, our description of these as transitional.

Volumetric Strain Rates
The Arthern et al. (2010) model implies an activation energy of 17.6 kJ mol −1 for both stages 1 and 2 once T ≈ T m . We test whether the iSTAR strain rates are consistent with this model by plotting −F 0 exp and , with E 0 = E 1 = 17.6 kJ mol −1 , as a function ofā for all sites except site 2 ( Figure 13) for which the data are not reliable. The best linear fit through the origin for the PIG data in stage 1 has a gradient of 421 ± 12 mwe −1 with r 2 = 0.98, whereas the gradient expected from the model is 687 mwe −1 . That is, the Arthern model overpredicts the strain rate in stage 1. This is not unexpected as the Arthern model was tuned using data which included the surface snow layer, in which temperature varies and densification is enhanced, within the nominal stage 1. The best linear fit through the origin for the PIG data in the transition stage has a gradient of 273 ± 10 mwe −1 with r 2 = 0.97. This is close to the gradient of 294 mwe −1 expected from the Arthern model. This is nominally a stage 2 value, but was derived from strain rates in snow with 550 kg m −3 < ⪅ 700 kg m −3 , that is, from data covering the transitional densification region and only part of the stage 2 region. This may be why the model fits the PIG transitional data so well.
In order to simulate the two nominal transition depths at which = 550 and 800 kg m −3 in firn cores, Ligtenberg et al. (2011) found it necessary to supplement the Arthern model with functions depending on mean annual accumulation (equations (15) and (16)). In stage 1, at low accumulation rates, the Ligtenberg correction improves the fit to the PIG data. At higher accumulation rates in stage 1, and in the transition region, the strain rate is underpredicted. These differences can be explained by examining the data used to derive the Ligtenberg corrections. The depth at which = 550 kg m −3 depends on densification both in the surface layer and in stage 1 (as we define it), with the relative importance of enhanced densification in the surface layer dependent onā. Similarly the depth range between = 550 and 800 kg m −3 covers the whole of stage 2 and the transition region, with the relative importance of higher densification in the transition region dependent onā. It is therefore physically reasonable that there should be correction factors to the Arthern model, which depend on the accumulation rate. However, since the Arthern model has been derived for sites withā = 0.13-1.04 mwe a −1 , ice core data from sites withā outside this range cannot necessarily be used in the determination of these correction factors.

F/k
Finally, we test whether the negative ratio of the density-corrected volumetric strain rate F to the vertical densification rate k is approximately equal to the mean annual accumulation rateā, as we expect for steady state conditions when horizontal flow divergence can be neglected (equation (8)). Figure 14a shows this is the case for stage 1 (where the mean valueF = F 0 and k = k 0 ), but Figure 14b shows that transitional values of −F∕k (withF = F t and k = k t ) clearly deviate from the 1:1 curve for higher accumulations. Figures 14c and 14d show the ratio of the density-corrected vertical strain rate F z to the vertical densification rate k as a function ofā. At low accumulations in stage 1, −F z0 ∕k 0 is slightly closer than −F 0 ∕k 0 toā; at high accumulations in the transitional stage the opposite is true and −F t ∕k t is slightly closer than −F zt ∕k t toā. For some high accumulation sites (14, 15, and 16),̇H > 0.05 and we would not expect −F z ∕k to give an improved approximation toā (see section 3.4) . However, −F t ∕k ≠ā for higher accumulation rates even wheṅH is negligible. We must conclude that the Robin hypothesis, which leads to equation (8), is the source of the problem. If the effect of bulk density on the strain rate is not well described by the term ( i − )∕ over a certain range, thenF ≠ −āk for that range.
The consequence is, therefore, that it is only in the case of stage 1 densification with negligible horizontal divergence that the Herron and Langway (1980) vertical densification model can be transformed to a strain rate model usingā. We test whether the iSTAR data are consistent with such a transformed model by plotting −F 0 exp against the accumulation rateā for each site except site 2, 15, 16, 18, and 19 wherėH is high (Figure 15). E 0 is fixed at 10.16 kJ mol −1 . The gradient of the best linear fit through the origin is 12.0 ± 0.4 mwe −1 (r 2 = 0.98); the expected value is k * 0 = 11(+6, −4) mwe −1 .

Conclusion
Deriving accumulation rates from snow density profiles in the PIG basin has proved challenging, but by introducing the automatic layer identification routine Straticounter (Winstrup et al., 2012) we have been able to produce estimates of mean annual accumulation which are consistent with previous airborne measurements  and with direct measurements of accumulation over 1 year. Chemical analysis of ice core samples remains the most accurate way to determine annual layering, especially when information from several different species can be combined. However, the stratigraphic method does have some advantages: high-resolution density profiles can be collected relatively easily using the neutron probe and can be analyzed in situ. In the overall iSTAR project they have an important role to play in extending the accumulation data set to sites where ice core data were not collected. The density profiles also provide information needed to derive accumulation from radar returns. We have derived a mean density profile which can be used to convert two-way travel times to depth and have shown that the effect of annual fluctuations in density on depth estimates is insignificant, except very near the surface.
Other members of the iSTAR team will construct an accumulation map for the PIG basin using the radar data to extend the site measurements and will make a detailed comparison of measured climatological variables with those derived from meteorological models. In this paper we have made a limited comparison of our accumulation rates with surface mass balance given by the RACMO2.3 model (van Wessem et al., 2014). Although there is good agreement for sites above 1,000 m a.s.l., the model results for most of the iSTAR traverse sites have to be regarded with some caution. Hence, in this paper, we discuss densification in terms of the measured mean climatological parameters T m andā alone. A further study using time-varying meteorological input data will be undertaken when suitable meteorological models have been properly assessed against the limited in situ meteorological data available from the PIG basin.

10.1002/2017JF004357
Stress, Pa. Random error in density measurement, kg m −3 . * Systematic error in density measurement, kg m −3 . F Random error in density-corrected volumetric strain rate, a −1 . * F Systematic error in density-corrected volumetric strain rate, a −1 . * q Systematic error in mass measurement, mwe * R Systematic error in access hole radius, m. a Standard deviation of log-normal distribution of annual accumulation.
Time since deposition of snow, a.