Skip to main content Accessibility help


  • Access
  • Cited by 4



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Millennially averaged accumulation rates for the Vostok Subglacial Lake region inferred from deep internal layers
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Millennially averaged accumulation rates for the Vostok Subglacial Lake region inferred from deep internal layers
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Millennially averaged accumulation rates for the Vostok Subglacial Lake region inferred from deep internal layers
        Available formats
Export citation


Accumulation rates and their spatio-temporal variability are important boundary conditions for ice-flow models. The depths of radar-detected internal layers can be used to infer the spatial variability of accumulation rates. Here we infer accumulation rates from three radar layers (26, 35 and 41 ka old) in the Vostok Subglacial Lake region using two methods: (1) the local-layer approximation (LLA) and (2) a combination of steady-state flowband modeling and formal inverse methods. The LLA assumes that the strain-rate history of a particle traveling through the ice sheet can be approximated by the vertical strain-rate profile at the current position of the particle, which we further assume is uniform. The flowband model, however, can account for upstream strain-rate gradients. We use the LLA to map accumulation rates over a 150 km × 350 km area, and we apply the flowband model along four flowbands. The LLA accumulation-rate map shows higher values in the northwestern corner of our study area and lower values near the downstream shoreline of the lake. These features are also present but less distinct in the flowband accumulation-rate profiles. The LLA-inferred accumulation-rate patterns over the three time periods are all similar, suggesting that the regional pattern did not change significantly between the start of the Holocene and the last ~20 ka of the last Glacial Period. However, the accumulation-rate profiles inferred from the flowband model suggest changes during that period of up to 1 cma–1 or ~50% of the inferred values.

1. Introduction

Vostok Subglacial Lake, East Antarctica, is the world’s largest known subglacial lake. Its interaction with the overlying East Antarctic ice sheet, through mass and energy exchange and reduction of basal drag, is of considerable glaciological interest (e.g. Bell and others, 2002; Pattyn and others, 2004). Radio-echo intensities of the ice–lake interface could potentially delineate regions of basal melting and accretion over the lake, but such analysis requires estimates of the radar attenuation over the lake, which is strongly temperature-dependent (e.g. MacGregor and others, 2007). Ice-flow models can predict ice temperatures, but such models require maps of accumulation rates and their temporal variations to make accurate predictions. Mapping accumulation rates over the Vostok Subglacial Lake region is therefore valuable for future glaciological studies there, as well as the paleoclimatic interpretation of the Vostok ice core (e.g. Parrenin and others, 2004) and possible future ice cores in this region.

Acquiring extensive field-based accumulation-rate measurements is a significant undertaking (e.g. Qin and others, 1994). However, such field data are inevitably sparse relative to the size of Antarctica and are subject to significant interannual and spatial variability, complicating the interpretation of trends in those data (e.g. Magand and others, 2007). Atmospheric modeling can estimate regional-scale accumulation rates (e.g. Rignot and others, 2008), but these studies have low spatial resolution and must be corrected for model biases using other observations (Magand and others, 2007). Satellite-microwave data (e.g. Arthern and others, 2006) can be used to interpolate field-based measurements, but the spatial resolution of the satellite data is also limited (25 km). This resolution is generally not fine enough to capture local accumulation-rate anomalies (e.g. an apparent accumulation-rate high along the western/upstream shoreline of Vostok Subglacial Lake (Leonard and others, 2004)).

Here we infer the spatial pattern of accumulation rates in the Vostok Subglacial Lake region averaged over three time periods dating back to 41 ka ago using layer-depth data from a gridded airborne radar dataset (Fig. 1). We use both a steady-state flowband model with inverse methods and a local one-dimensional strain-rate model to infer accumulation rates. The former approach considers the effects of both upstream accumulation-rate gradients and ice-thickness gradients upon the layer shapes, while the latter approach ignores them. We apply the flowband model along four flowbands and compare accumulation rates inferred from this method to an accumulation-rate map for the entire study area based on the simpler model that neglects upstream strain-rate gradients. We also discuss temporal changes of the accumulation-rate pattern over the three layer ages and compare our accumulation-rate map to previous estimates for this region.

Fig. 1. (a–c) Color maps of the depth to the three internal layers (A–C) used in this study. Layer A could be tracked over only about two-thirds of the study area. The center flowlines of the four flowbands (1–4) are shown as white lines and labeled at the head of the flowbands, and the flow directions are shown with white arrows. The locations of surface-velocity and surface accumulation-rate data are labeled. The edge of Vostok Subglacial Lake is outlined in black. The gray fill around each center flowline in (b) shows the width variation of the flowband. The start of each flowband has the same width (5 km). The black lines in (c) are the airborne-radar flight-lines.

2. Data and Methods

2.1. Radar data

We use the 60MHz airborne ice-penetrating radar data collected over the Vostok Subglacial Lake region in an approximately 150 km ×350km grid by the US Support Office for Aerogeophysical Research (SOAR) at the University of Texas Institute for Geophysics. Blankenship and others (2001) described the system characteristics. This dataset was used to derive both the surface and bed topographies (Studinger and others, 2003). The radar lines were flown along an orthogonal grid; lines with a roughly east–west orientation had a 7.5 km line spacing, and lines with a roughly north–south orientation had an 11.25 or 22.5 km line spacing (Fig. 1c). North–south profiles with the smaller line spacing were closer to the shoreline because the main goals of the aerogeophysical survey were to delineate the lake shoreline and to determine the geological controls on the lake. From these data, we use the two shallowest internal layers picked by Tikku and others (2004) (layers A and B) that were traceable over most of our study area (Fig. 1a and b), and another picked for this study (layer C, the deepest layer; Fig. 1c). Layer and bed depths were calculated using a radio-wave speed of 168.4mμs–1 averaged over the ice column to the ice–lake interface (Popov and others, 2003), and the layer depths varied between 334 and 1290 m. The layer, surface and bed elevations were linearly interpolated onto a 1.5 km grid. The location of the lake shoreline is determined as the edges of the region where the bed echo is bright and flat (Studinger and others, 2004).

We estimated the layer ages by first interpolating the gridded layer depths to the Vostok ice-core site and then interpolating the depth–age scale for the ice core (Parrenin and others, 2004) to the depth of each layer there. From shallowest to deepest, the layers labeled A, B and C are 26.2, 34.8 and 41.0 ka old, respectively. We assign an uncertainty of 2.0 ka to these layer ages (5–8% of their ages), because of the uncertainty in the layer depths, the uncertainty in the depth–age scale and the layer-depth interpolation that was used.

2.2. Inferring accumulation rates from layer depths

There are several approaches to inferring accumulation rates from internal-layer depths in ice sheets (Waddington and others, 2007). The simplest approach is the shallow-layer approximation (SLA), where the accumulation rate ḃ is equal to the ice-equivalent layer depth z divided by its age A. However, the SLA is valid only for very shallow layers (<1% of ice-sheet thickness) that have not yet experienced significant cumulative dynamic strain. For deeper layers, the local-layer approximation (LLA) may be valid. The LLA assumes that the actual strain-rate history of a particle traveling through the ice sheet can be approximated using the vertical strain-rate profile at its current location. If we further assume that the vertical strain rate is uniform with depth, the LLA can be used to infer accumulation rates ḃLLA from the layer depths and ice thickness Has


This approach has been applied previously to radar data from the Vostok Subglacial Lake region (e.g. Siegert, 2003). However, as also acknowledged by Siegert (2003), the LLA is not appropriate for deep internal layers composed of particles that have traveled significant horizontal distances and thus experienced non-steady and non-uniform strain-rate patterns. To quantify the suitability of the LLA, Waddington and others (2007) defined the non-dimensional depth number D. This number is based on the relationship between the horizontal distance (L path) that a particle has traveled from the surface to that layer, and the characteristic length scales of accumulation-rate (Lb) and ice-thickness variability (LH ) as


If D≪1, the layer is not ‘deep’, so that the spatial gradients in accumulation rates and ice thickness do not significantly affect the layer depths, and hence the LLA is sufficiently accurate to infer accumulation rates. We note that values of D help determine where the LLA is applicable, but there is no simple quantitative relationship between D and uncertainty in ḃLLA. L path is the product of the layer age and the depth-averaged horizontal velocity ū experienced by the particle as it traveled along its trajectory from the surface to the layer:


We use the measured surface speed at Vostok station (2ma–1; Wendt and others, 2006) to approximate u ~over the entire study area. At Vostok station, the measured surface speed should be equivalent to u ~because the vertical profile of horizontal velocity is uniform for floating ice. Characteristic lengths Lḃ and LH are derived from the along-flow (x-direction) gradients in ḃ and H, respectively:



To calculate dḃ/dx and dH/dx, we use ḃLLA and H values averaged over 60 km along-flow, which is close to the mean L path value for the three layers (68 km). Because our goal is to evaluate the suitability of the LLA, these db=dx values should be of the correct order of magnitude and thus adequate for the principal purpose of calculating D.

2.3. Flowband model and inverse solution procedure

A method that can provide more accurate accumulation-rate patterns from deeper layers is a combination of flowband modeling and formal inverse theory. We use the steady-state flowband model and inverse solution procedure described by Waddington and others (2007), and we refer the reader to that study for full details of that model. An inverse problem begins by solving a forward problem that is the set of governing equations, boundary conditions and parameter values. In our forward problem, a flowband model calculates layer depths using steady-state ice-surface elevations, ice velocities and temperatures. It is a 2.5-dimensional ice-flow model the flowband width of which can vary, but it assumes that all properties are uniform transverse to the direction of ice flow. The flowband model also requires initial estimates of the layer age and the input ice flux at the upstream end of the flowband. This ice flux is estimated kinematically using the width, ice thickness and estimated depth-averaged horizontal velocity at the start of the flowband. Uncertainty in the ice flux is set at 50% of its estimated value, because of the uncertainty in past accumulation rates (section 2.4). The ice flux and layer age are adjusted by the inverse solution procedure as necessary to match the data.

Flowbands were identified using the long-term flow directions determined by Tikku and others (2004) using structure tracking in the internal layers, not modern surface-slope gradients. At Vostok station, their flow direction matches the modern flow direction determined by Wendt and others (2006) to within 7 ˚. We do not use modern streamlines because they are difficult to constrain, due to the low surface slopes over the lake (Fig. 1a) and the limited spatial extent of reliable surface-velocity data. None of the airborne-radar profiles followed ice flowlines, so we calculated flowlines within the gridded flow field and linearly interpolated the layer, surface and bed elevations along them. We chose four flowbands that traversed Vostok Subglacial Lake, evenly divided the study area, and took advantage of available surface-velocity and accumulation-rate data. The southernmost flowband (labeled ‘1’) passes through Vostok station (Fig. 1a). Flowbands 2 and 4 pass through sites B37 and B78, respectively, for which modern accumulation-rate data are available. Flowband 3 crosses the lake roughly halfway between flowbands 2 and 4. We calculated the flowband widths by finding two flowlines that began 2.5 km away from the start of the center flowline and normal to the initial direction of flow, and then calculating the distance between those two flowlines (the edges of the flowband) along-flow. The flowband widths vary between 40% (flowband 1) and 620% (flowband 2) of their initial values.

Although ice temperatures and strain_rates can strongly influence each other, the temperature and mechanical components of our flowband model are not coupled. Depth-averaged horizontal velocities are determined kinematically, and the vertical profiles of horizontal velocity are calculated by multiplying temperature-dependent non-dimensional shape functions by the depth-averaged horizontal velocity. The vertical velocities are then derived from the horizontal velocities using local incompressibility. Because there is no basal drag over the lake, the shape functions should be uniform and equal to unity there. Those shape functions are adjusted so that the initial temperature-dependent shape functions are close to uniform values; they are not exactly uniform but are smoothed so that there is no velocity discontinuity at the edges of the lake.

To calculate ice temperatures along the flowlines for the temperature-dependent vertical shape functions of horizontal velocity, we use a one-dimensional steady-state temperature model that includes vertical, but not horizontal, advection and diffusion of heat (Paterson, 1994, p. 218). We horizontally smooth the two-dimensional temperature field formed from this set of one-dimensional profiles over a length scale of approximately 5 km (116–200% of ice thickness) to approximate the effects of horizontal advection and diffusion of heat, which are neglected in the temperature model. This temperature model neglects basal melting and freezing that occur over the lake because these rates are believed to be generally small relative to the accumulation rates (<0.5cm a–1; Wendt and others, 2006) and are not yet well constrained over the entire lake. However, one recent model suggested that there is significant spatial variability in the melting/freezing pattern, and predicted basal freezing rates >5cm a–1 near the western shoreline (Thoma and others, 2008), which would affect the layer depths and hence accumulation rates inferred using this model.

The layer depths generally shallow as the ice flows across the lake (see Fig. 4 below). However, the layers in the northern part of the lake (flowband 4 in Fig. 4k below), where basal melting is predicted (e.g. Thoma and others, 2008), shallow more abruptly than those layers in the southern part of the lake, where basal accretion has been observed (Bell and others, 2002; Tikku and others, 2004). This layer-depth pattern suggests that basal melting and accretion have less influence upon the layer depths used in this study than transverse strain rates, which generally increase, causing layer thinning, as the flowbands widen over the lake.

The ḃfb profiles are inherently smoother because the inverse solution procedure imposes a smoothness constraint upon ḃfb. The real accumulation-rate pattern may have more structure at shorter wavelengths than can be inferred from the layer, but ḃfb more accurately represents the strain-rate histories of the particles along their paths than ḃLLA. From the surface to a deep layer, particles have traveled a horizontal distance at least an order of magnitude greater than their depth, and they have likely traversed significant strain-rate gradients. Both the real ice sheet and the particle paths calculated by the flowband model integrate these upstream gradients, so it is difficult to discern the effect of small-scale spatial changes in accumulation rate on deep layers by any method that uses a deep layer.

For a given along-flow accumulation-rate profile and boundary conditions, the flowband model produces a velocity field, from which modeled isochronous layers are calculated. We then compare the depths of our modeled layers to those of each dated internal layer in our study area.

The initial values of the accumulation rates (b LLA), input ice flux and layer age are adjusted iteratively using an inverse solution procedure that seeks the ‘best’ accumulation-rate pattern ( fb), which is the pattern that is spatially smooth and that fits the layer-depth, surface-velocity and accumulation-rate data at an expected tolerance based on the data uncertainties. We must incorporate these uncertainties using an expected tolerance because overfitting the data can produce spurious variations in ḃfb. Our solution is one whose ḃ profile has a low curvature and the modeled input flux and layer age of which have small deviations from their prescribed values. The smoothness and data-fitting constraints are competing requirements upon the final ḃfb profile; these requirements are balanced by the inverse solution procedure using a trade-off parameter.

The upstream (western) edge of our study area is more than 200 km from Ridge B, which is an ice-flow divide (Parrenin and others, 2004). Consequently, all the ice in the Vostok Subglacial Lake region is undergoing flank flow and the relatively shallow particle paths that we are modeling all have similar near-linear trajectories. Particle paths that begin at the upstream edge of our study area reach the layer depths only after traversing one-third to one-half of the length of the flowband. The portion of the layer that is not intersected by particle paths beginning within our study area therefore has no influence on the inferred accumulation-rate pattern. This flow regime destabilizes the inverse-solution procedure from one iteration to the next and can result in spuriously large changes in accumulation rate. Such behavior is non-physical, so we truncate small non-zero singular values to stabilize the inverse-solution procedure, as described in appendix A of Waddington and others (2007). Although this approach decreases the ability of formal inverse methods to detect accumulation-rate changes at the upstream and downstream ends of the flowbands (as opposed to not truncating the singular values), it allows us to recover physically reasonable solutions.

2.4. Surface-velocity and accumulation-rate data

In addition to radar-layer depths, the inverse-solution procedure can be constrained by surface velocities and accumulation rates averaged over the ages of the layers. Modern surface-velocity data are only available for flowband 1 at Vostok station (2ma–1; Wendt and others, 2006). Long-term accumulation-rate data are sparse (Magand and others, 2007) and are available only for flowbands 1, 2 and 4 at the three locations shown in Figure 1a (2.2cm a–1 at Vostok, 4.0cma–1 at B37, 3.5 cm a–1 at B78; Lipenkov and others, 1998). All accumulation rates presented here and elsewhere in this paper are in ice equivalent.

To constrain our steady-state flowband model using these data, it is necessary to estimate constant accumulation rates and ice-flow speeds averaged over the past 41 ka. However, Parrenin and others (2004) found that the accumulation rates inferred from the Vostok ice core increased two-fold going from the Last Glacial Maximum (LGM) 18 ka ago into the Holocene. We account for this change in accumulation rates by multiplying the accumulation-rate data for the two oldest layers (B and C) by a factor of 0.75. This factor is calculated as the sum of the fractions of the layer age spent in each period (Holocene or glacial) multiplied by the ratios of the accumulation rate during each period to the modern value (b modern ≈ ḃHolocene ≈ 2 glacial). To maintain a steady state, smaller accumulation rates require lower ice-flow speeds, and Leonard and others (2004) also inferred 50–65% slower ice velocities between 26 and 41 ka ago from the hinge points in the depths of layers A and C. We therefore apply the same reduction fraction (50%) for the ice speed when applying the flowband model to layers B and C. We also assume an uncertainty of 20% in all of these data when attempting to match them using inverse methods.

3. Results

3.1. LLA-inferred accumulation-rate maps

Figure 2b shows ḃLLA inferred from layer C (41 ka old). ḃLLA values are generally higher in the western half of our study area, and there is a 2.5 cm a–1 difference between the lowest and highest values of ḃLLA inferred from layer C. The lowest values are found 10–20km downstream of the eastern/downstream shoreline of the lake. The highest values are focused in the northwestern corner of the lake and are consistent with the larger layer depths in that area (Fig. 1c).

Fig. 2. (a) Surface topography (contour interval is 10 m) over a color map of modern ice-equivalent accumulation rates (25 km grid spacing) derived from satellite-microwave emission and field-based data (Arthern and others, 2006). (b) LLA-inferred accumulation rates from layer C. (c) Difference between LLA-inferred accumulation rates from layer A and those for layer C. (d) Difference between LLA-inferred accumulation rates from layer B and those for layer C. The color bar to the right of (d) is for both (c) and (d).

The differences between ḃLLA values inferred from layers A and B in relation to layer C are shown in Figure 2c and d. The mean difference in ḃLLA between layers A and C is 0.34 cma–1; between layers B and C it is 0.09 cma–1. The larger mean difference between layers A and C than between B and C is consistent with the larger difference between their ages. Forty-six percent of the age of layer A is spent in the Holocene, whereas smaller fractions of the ages of layers B and C are spent in the Holocene (24% and 29%, respectively). Layer A is therefore composed of particles that experienced more of the Holocene accumulation-rate history, which has higher values than the Glacial Period that dominates layers B and C. The apparent accumulation-rate high near the upstream shoreline observed by Leonard and others (2004) is more prominent in layer A than it is in layer B or C. It is displaced east/downstream from the lake shoreline because the equivalent trough in the layer depths has flowed downstream during the intervening 26 ka (the age of layer A). However, it is only displaced about 10 km from the lake shoreline, not ±50 km (26 ka ×2ma–1), suggesting that either ū is lower there either now or in the past (Leonard and others, 2004), or that the cause of the trough in layer depths does not originate at the lake shoreline. Because the difference in age between layers B and C (6 ka) is small relative to their ages (<20%) and both layers are older than the age of the Holocene–LGM accumulation-rate change, there is little difference between their LLA values (Fig. 2c).

3.2. Suitability of the LLA (D values)

Figure 3 shows non-dimensional depth number D (Equation (2)) for our study area for the three different layers; Figure 4 shows interpolated values of D along the flowbands. The mean values of D for the study area are 0.28, 0.44 and 0.50 for layers A, B and C, respectively. D is generally larger for progressively deeper/older layers because D is proportional to L path, which increases with layer age. The mean ratio of L/LH varies between 2.7 and 5.2 for the different layers. The smaller value of either L or LH will tend to dominate their contribution to D, so ice-thickness gradients in our study area have a greater influence on ice flow than accumulation-rate gradients.

Fig. 3. (a–c) Along-flow values of D for all three internal layers (Equation (2)).

Although we are using the three shallowest spatially extensive internal layers in the radar data, nearly all of our study area (96–98%) has D values >0.1 that do not satisfy the D≪1 criterion, suggesting that the LLA may not accurately infer accumulation rates. Over the eastern (downstream) half of the lake, the LLA may be acceptable, but the large ice-thickness gradients near the lake shoreline suggest that the accumulation-rate pattern inferred with the LLA should be further investigated using a more sophisticated approach. We calculated D from values that were spatially averaged over 60 km along-flow, which is similar to L path for any of the layers, but L path is not well constrained for individual particle paths without using a flowband model. However, averaging over this distance produces a more meaningful value of D because it more accurately captures the length scales of changes that particles have experienced.

3.3. Accumulation-rate profiles inferred along flowbands

Figure 4 shows the LLA profiles for all four flowbands, which were used to initialize the inverse solution procedure, and the final fb profiles inferred from that procedure. Figure 4 also shows the surface, layer and bed elevations, and D values along each flowband. For all four flowbands, accumulation rates are higher at their upstream end (0.2–1.0 cma–1) than downstream end, although the detailed structure of the accumulation-rate profiles varies substantially between each flowband. They also decrease smoothly across the lake, with a typical gradient of approximately –0.01cm a–1 km–1. Layers B and C are close in age (5 ka) relative to the age of layer C (12% of 41 ka), and both their LLA and fb profiles are similar (mean difference of 8% of layer B’s values for flowbands 1, 2 and 4). Despite being 9 and 15 ka younger than layers B and C, respectively, the large-scale structure (>50 km) of layer A’s accumulation-rate profiles are often similar to those inferred from the older layers, especially for flowband 1. Accumulation rates along the flowband that intersects Vostok station (flowband 1) were twice as high during the Holocene as during the Glacial Period (Parrenin and others, 2004), which is consistent with the difference in magnitude between LLA and fb inferred from layer A versus those profiles inferred from layers B and C. Based on these results, we argue that the spatial pattern of relative accumulation rate has changed in our study area over the last 41 ka, although our ability to resolve these changes is limited by the steady-state models that we used, and the changes are greater for some flowbands (e.g. flowband 4) than for others (e.g. flowband 1). We note that there is no initial expectation built into the flowband model that the accumulation-rate patterns over this period would be similar.

4. Discussion

4.1. Comparison between LLA and fb profiles

Our initial evaluation of the non-dimensional number D (Fig. 3) suggested that the LLA generally could be inaccurate for our study area. The mean D values along each flowband range between 0.2 and 0.5, with larger values for deeper layers, and fb is generally lower than LLA where D is low along the flowband. The LLA and fb profiles often have different shapes, except for layer A for flowband 1, where they are nearly identical. These patterns suggest that our D≪1 criterion is sufficient for the suitability of the LLA but that it is not always necessary.

Structures in the LLA profiles are often translated further downstream than the fb profiles for progressively deeper layers (e.g. the upstream accumulation-rate high in flowband 4). This pattern is expected, as the inverse solution procedure should correctly assign these structures to their upstream origin on the surface, i.e. show less erroneous downstream translation of the structures in the accumulation-rate profiles. Flowband 1 has a relatively smooth fb profile and generally lower values than the other flowbands. The accumulation-rate high at the upstream shoreline visible in the b LLA profile is less prominent in the b fb profile. Flowband 2 has larger differences between LLA and fb, and those differences are consistent for all layers. For example, there is no significant accumulation-rate high along the upstream shoreline in fb, and fb is also higher over the lake. Upstream of the lake, fb is consistently lower than ḃLLA, whereas east/downstream of the lake the opposite is true. Flowband 2 is the longest of the four flowbands and it extends further downstream of the lake than any of the other flowbands, which may explain why it captures a larger difference between LLA and fb east/downstream of the lake. Flowband 4 is similar to flowband 1 in that it has a smoother decrease in both LLA and ḃfb over the lake. It has a larger difference between ḃLLA for layers B and C than for the other flowbands, and the same is true for fb. The ḃfb profile for layer C in flowband 4 deviates from the other ḃfbprofiles at the upstream and downstream ends of the lake. There, ḃfb is about 20% higher for layer C than would be expected based on the ḃfb profiles of layers A and B. This difference suggests either a change in accumulation rates there between 35 and 41 ka ago or that the flowband model has difficulty reproducing layer C.

Flowband 3 shows large differences between ḃLLA and ḃfb that are not consistent between the three layers. These ḃfb profiles are also very smooth and have no structures in common with the ḃLLA profiles. Although it is able to converge to a solution, we suspect that the inverse solution procedure has failed to accurately model flowband 3 because of the unusual ḃfb profiles and for two additional reasons. First, flowband 3 is not constrained by any surface-velocity or accumulation-rate data, thus it is the least constrained of all the flowbands we used. Second, the flowband model likely does not adequately represent the more complicated ice dynamics along flowband 3, which crosses over a shallow embayment before crossing the main body of Vostok Subglacial Lake. Flowband 1 also crosses over a shallow embayment, but it is better constrained by other data. The inverse solution procedure cannot resolve features the wavelength of which is close to L path, i.e. such features are in its nullspace, which may cause the features in the ḃfb profiles of flowband 3.

For all three layers along flowbands 1, 2 and 4, the mean differences between ḃfb and ḃLLA are 4%, 16% and 12% of ḃfb, respectively. The differences between ḃfb and ḃLLA are greater along the lake shoreline and away from the lake, and are generally larger for deeper layers. These small relative differences, along with the above comparisons of the LLA and ḃfb profiles, suggest that the ḃLLA map is an acceptable proxy for the real accumulation-rate map, particularly over the northern and southern ends of the lake, and less so in its middle. If numerous along-flow radar profiles are not available for a region of an ice sheet (as is generally the case), then the LLA is the best way to predict the regional accumulation-rate pattern from radar layers. However, we have not investigated some features of this ḃLLA map, such as the accumulation-rate high in the northwestern corner of the lake (Fig. 2), and we note that features of the ḃLLA map on spatial scales smaller than the radar-line spacing are less reliable.

4.2. Comparison with previous studies of Vostok Subglacial Lake accumulation rates

Arthern and others (2006) presented a 25 km resolution map of modern accumulation rates across Antarctica, inferred from surface-measured accumulation-rate data and microwave emission recorded by satellites (Fig. 2a). It shows a broad low in accumulation rates (about 4 cma–1) over the southern half of the lake and higher values (>5 cma–1) northeast of the lake. It also shows lower accumulation rates directly over the lake, which may be an artifact due to lower decimeter-scale surface roughness over the lake. Field measurements from Qin and others (1994) are in better agreement with our ḃLLA map for layer A than the map of Arthern and others (2006). Although the absolute accumulation-rate values differ between our map and that of Arthern and others (2006), there is a similar pattern of increasing values to the north in our study area. However, the map of Arthern and others (2006) does not resolve the large accumulation-rate high in the northwestern corner of the lake that we infer from the internal layers (Fig. 2b). It also does not resolve the ḃLLA high along the upstream lake shoreline because of its coarser resolution.

Siegert (2003) inferred accumulation rates using radar-layer depths along a radar flight-line that crossed over Ridge B and Vostok station and that is close to our flowband. He included a layer of a similar age (46 ka) to the oldest layer in this study (layer C, 41 ka). Our work improves upon that study because Siegert (2003) used only the LLA for layers of ages similar to, or older than, those used in this study, and because our flowband 1 follows the ice flow more accurately. The values of D shown in Figures 3c and 4b suggest that the LLA is generally not suitable for flowband 1, but the small relative difference between ḃLLA and ḃfb for flowband 1 (5%) suggests that the LLA is acceptable. The change in accumulation rate across Vostok Subglacial Lake that Siegert (2003) found using the 46 ka old layer is close to that which we infer using layer C (<0.5 cma–1; Fig. 4c). Our results and those of Siegert (2003) and Leysinger Vieli and others (2004) infer higher accumulation rates along this flowband upstream from the western edge of the lake (Fig. 4c).

Fig. 4. Along-flow characteristics of flowbands 1 (a–c; top left), 2 (d–f; top right), 3 (g–i; bottom left) and 4 (j–l; bottom right). The vertical gray bands represent the portion of each flowband that overlies Vostok Subglacial Lake. (a, d, g, j) Surface, layer and bed elevations along the flowbands. The black (blue) lines represent the elevations of the surface and bed (three internal layers), and their vertical scale is in black (blue) and shown on the left (right). The vertical scale for the internal layers has a smaller range to better show their structure. Blue dots along the deepest layer (C) represent the points at which the radar lines cross the flowband, which shows where the two-dimensional grid interpolation may have introduced spurious structure into the internal layer shapes. Black circles along the surface-elevation profile show the location of field data used to constrain the inverse solution procedure. (b, e, h, k) Smoothed D shown in logarithmic scale for all three layers (A: black; B: blue; C: red). (c, f, i, l) Ice-equivalent accumulation-rate profiles inferred with the LLA ( LLA; dashed) and with formal inverse methods and a flowband model ( fb; solid) using the three layers (A: black; B: blue; C: red). The horizontal range for all panels is the length of the longest flowband (2).

Leonard and others (2004) identified a stationary accumulation-rate high along the upstream shoreline of Vostok Subglacial Lake that is likely due to the relatively large changes in surface slopes there. Higher accumulation rates have also been measured there using 1 m snow pits (Qin and others, 1994). Accumulation rates are often higher slightly downstream of steep surface slopes because katabatic winds are stronger (weaker) on steeper (shallower) slopes, so snow will be transported and then accumulate slightly downstream of the slope inflection points (Vaughan and others, 1999; King and others, 2004). Our ḃfb profiles do not show higher values near the upstream lake shoreline, but the expected accumulation-rate high near there may not be resolvable using the flowband model because the localized high can be averaged out by ice flow. These results suggest that ice-flow changes due to flotation over Vostok Subglacial Lake, as well as accumulation-rate changes along the upstream shoreline, cause the observed trough in the layer depths near the upstream shoreline from which an accumulation-rate high was inferred.

4.3. Uncertainty in fb and improvements to the flowband model and inverse solution procedure

Figure 4 shows that uncertainties in ḃfb can be large (>50%) given limited data to constrain them, as was the case for flowband 3. A problem that affects all flowbands is the limited information on the upstream surface and bed topographies. Additional radar data that both follow flowlines and survey the area upstream of Vostok Subglacial Lake to the ice divide at Ridge B would be valuable for constraining the inverse solution procedure. More field-measured surface-velocity and accumulation-rate data would provide further constraints (section 2.4). The flow field around the perimeter of Vostok Subglacial Lake is poorly constrained due to the low surface-slope gradients and few internal structures that can be tracked (Tikku and others, 2004). Additional layer picking may resolve this issue. If we had further confidence in the flow field, then we could model a spatially denser set of flowbands and interpolate the differences between ḃLLA and ḃfb along those flowbands across the entire study area to produce a more accurate accumulation-rate map.

There are several simplifications in our flowband model that could be improved upon. Our one-dimensional temperature model is unsophisticated relative to the temperature calculations in many existing thermomechanical ice-flow models (e.g. Pattyn and others, 2004). Our modeled temperature profile is consistently higher than the observed temperature profile (mean difference of 2.9 K) at Vostok station (personal communication from V.Ya. Lipenkov, 2006). However, this difference is largest in the 1500–2500m depth range, which is greater than the depth of the deepest layer used here, and ice temperatures determine only the normalized horizontal velocity field (its shape functions) and not its magnitudes. Our simplified temperature model is therefore acceptable for the purposes of the flowband model. A flowband model that includes longitudinal strain-rate gradients and a coupled temperature model (e.g. Price and others, 2007) would better represent the ice dynamics over the lake, although initial and boundary conditions would be more difficult to set. Finally, an alternate inverse solution procedure using Monte Carlo methods could better evaluate the sensitivity of the modeled accumulation rates to our initial guesses of the model parameters (e.g. Steen-Larsen and others, 2007).

We used a steady-state flowband model to infer accumulation rates, although previous modeling work using radar layers suggest that ice flow was non-steady in the Vostok Subglacial Lake region during the last 41 ka (Leysinger Vieli and others, 2004; Salamatin and others, in press). Surface elevations, ice velocities and the location of the Vostok flowline may have changed in the last 41 ka, but here we have implicitly assumed that such changes did not significantly affect the inferred accumulation rates. If the flowline that passes through Vostok has not changed in the last 41 ka, then our steady-state model should still recover the correct mean accumulation rate during this period. Salamatin and others (in press) tuned a non-steady thermomechanical flowband model along the Vostok flowline that produced layers that matched well with the observed layers. They did not use formal inverse methods to match their data, but their results and ours suggest that future accumulation-rate studies should combine the increased accuracy of non-steady flowband models with the computational efficiency of formal inverse methods.

5. Conclusions

We have presented an accumulation-rate map for the Vostok Subglacial Lake region inferred from internal-layer depths observed in radar data. By comparing this map with results from a formal inverse method that incorporates a flowband model, we find that the regional accumulation-rate pattern inferred from internal layers using the LLA is an acceptable estimate of accumulation rates in this study area. In terms of its spatial resolution and evaluation of its uncertainties, our regional map is a significant improvement upon previous studies of accumulation rates in this region. It reproduces some features of the spatial pattern that have been observed previously, including a broad low in the southern half of the lake and a high near the upstream shoreline of the lake. We also infer a broad accumulation-rate high in the northwestern corner of the lake that has not been identified previously. For at least two of our four modeled flowbands, this accumulation-rate pattern has changed significantly (up to 50%) over the past 41 ka, possibly during the transition from the last Glacial Period to the Holocene.

Airborne radar surveys are often designed as orthogonal grids that may not follow ice flowlines. Without extensive efforts to reconstruct layer depths along flowlines from sparsely crossing radar profiles, the LLA is a reasonable method for inferring the regional accumulation-rate pattern. However, our results show that the LLA occasionally predicts questionable abrupt accumulation-rate variations in response to upstream ice-thickness gradients, such as along the shoreline of Vostok Subglacial Lake. The steady-state flowband procedure accounts for these gradients at the cost of a smoother accumulation-rate profile that cannot capture some short-scale anomalies, such as the possible high along the upstream lake shoreline. Because the ḃLLA profiles might have differed substantially from the ḃfb profiles (e.g. Waddington and others, 2007), we emphasize the need to apply formal inverse methods for inferring past accumulation rates from deep layers to improve upon accumulation-rate patterns inferred from simpler methods.


The US National Science Foundation supported this work at the University of Washington (ANT 0538674 and 0440666) and the Lamont–Doherty Earth Observatory (ANT 0537752). The Applied Physics Laboratory at the University of Washington provided a research fellowship for J.A.M. We thank SOAR at the University of Texas for collection of the radar data, V.Ya. Lipenkov for providing the Vostok temperature data, A.A. Ekaykin for providing accumulation-rate data, F. Pattyn for discussions and motivation for this study, and T.A. Neumann for developing temperature-dependent components of the flowband model. We thank the scientific editor F. Parrenin and two anonymous referees for many suggestions that improved the clarity of the manuscript.


Arthern, R.J., Winebrenne, D.P.. and Vaugha, D.G... 2006. Antarctic snow accumulation mapped using polarization of 4.3 cm wavelength microwave emissio.. J. Geophys. Res., 111(D6), D06107. (10.1029/2004JD005667.)
Bell, R.E., Studinge, M.., Tikk, A.A.., Clark, G.K.C.., Gutne, M.M.. and Meerten, C... 2002. Origin and fate of Vostok Subglacial Lake water frozen to the base of the East Antarctic ice shee.. Nature, 416(6878), 307–310.
Blankenship, D.D. and 9 others. 2001. Geologic controls on the initiation of rapid basal motion for West Antarctic ice streams: a geophysical perspective including new airborne radar sounding and laser altimetry results. In Alley, R.B. and Bindschadle, R.A.., eds. The West Antarctic ice sheet: behavior and environmen.. Washington, DC, American Geophysical Union, 105–121. (Antarctic Research Series 77.)
King, J.C., Anderso, P.S.., Vaugha, D.G.., Man, G.W.., Mobb, S.D.. and Vospe, S.B... 2004. Wind-borne redistribution of snow across an Antarctic ice ris.. J. Geophys. Res., 109(D11), D11104. (10.1029/2003JD004361.)
Leonard, K., Bel, R.E.., Studinge, M.. and Trembla, B... 2004. Anomalous accumulation rates in the Vostok ice-core resulting from ice flow over Vostok Subglacial Lak.. Geophys. Res. Lett., 31(24), L24401. (10.1029/2004GL021102.)
Leysinger Vieli, G.J.M., Sieger, M.J.. and Payn, A.J... 2004. Reconstructing ice sheet accumulation rates at ridge B, East Antarctic.. Ann. Glaciol., 39, 326–330.
Lipenkov, V.Y., Yekayki, A.A.., Barko, N.I.. and Pursh, M... 1998. O svyazi plotnosti poverkhnostnogo sloya snega v Antarktide so skorost’yu vetra [On the relationship between surface snow density in Antarctica with wind velocity.. Mater. Glyatsiol. Issled. 85, 148–158. [In Russian.]
MacGregor, J.A., Winebrenne, D.P.., Conwa, H.., Matsuok, K.., Mayewsk, P.A.. and Clo, G.D... 2007. Modeling englacial radar attenuation at Siple Dome, West Antarctica, using ice chemistry and temperature dat.. J. Geophys. Res., 112(F3), F03008. (10.1029/2006JF000717.)
Magand, O. and 6 others. 2007. An up-to-date quality-controlled surface mass balance data set for the 90˚ –180˚ E Antarctica sector and 1950–2005 perio.. J. Geophys. Res., 112(D12), D12106. (10.1029/2006JD007691.)
Parrenin, F., Rém, F.., Rit, C.., Sieger, M.J.. and Jouze, J.. 2004. New modeling of the Vostok ice flow line and implication for the glaciological chronology of the Vostok ice cor. J. Geophys. Res., 109(D2), D20102. (110.1029/2004JD004561.)
Paterson, W.S.B. 1994. The physics of glaciers. Third edition. Oxford, etc., Elsevier.
Pattyn, F., de Smed, B.. and Souche, R... 2004. Influence of subglacial Vostok lake on the regional ice dynamics of the Antarctic ice sheet: a model stud.. J. Glaciol., 50(171), 583–589.
Popov, S.V., Sheremet’ye, A.N., Masolo, V.N.., Luki, V.V.., Mirono, A.V.. and Luchinino, V.S.. 2003. Velocity of radio-wave propagation in ice at Vostok station, Antarctic. J. Glaciol., 49(165), 179–183.
Price, S.F., Waddingto, E.D.. and Conwa, H... 2007. A full-stress, thermomechanical flow band model using the finite volume metho.. J. Geophys. Res., 112(F3), F03020. (10.1029/2006JF000724.)
Qin, D., Peti, J.R.., Jouze, J.. and Stievenar, M... 1994. Distribution of stable isotopes in surface snow along the route of the 1990 International Trans-Antarctica Expeditio.. J. Glaciol., 40(134), 107–118.
Rignot, E. and 6 others. 2008. Recent Antarctic ice mass loss from radar interferometry and regional climate modellin.. Nature Geosci., 1(2), 106–110.
Salamatin, A.N., Tsyganov, E.A.., Popo, S.V.. and Lipenko, V.Ya... In press. Ice flow line modeling in ice core data interpretation: Vostok Station (East Antarctica). In Hondoh, T., ed. Physics of ice core recordsv, Vol. 2. Sapporo, Hokkaido University Press.
Siegert, M.J. 2003. Glacial–interglacial variations in central East Antarctic ice accumulation rate.. Quat. Sci. Rev., 22(5–7), 741–750.
Steen-Larsen, H.C., Koutni, M.R.. and Waddingto, E.D... 2007. Formulating an inverse problem to determine the accumulation rate pattern from deep internal layering in an ice shee.. Geophys. Res. Abstr., 9, 01181 (1607-7962/gra/EGU2007-A-01181.).
Studinger, M. and 11 others. 2003. Ice cover, landscape setting, and geological framework of Vostok Subglacial Lake, East Antarctic.. Earth Planet. Sci. Lett., 205(3–4), 195–210.
Studinger, M., Bel, R.E.. and Tikk, A.A... 2004. Estimating the depth and shape of Vostok Subglacial Lake’s water cavity from aerogravity dat.. Geophys. Res. Lett., 31(12) L12401. (10.1029/2004GL019801.)
Thoma, M., Grosfel, K.. and Maye, C... 2008. Modelling accreted ice in subglacial Vostok Subglacial Lake, Antarctic.. Geophys. Res. Lett., 35(11) L11504. (10.1029/2008GL033607.)
Tikku, A.A., Bel, R.E.., Studinge, M.. and Clark, G.K.C... 2004. Ice flow field over Vostok Subglacial Lake, East Antarctica, inferred by structure trackin.. Earth Planet. Sci. Lett., 227(3–4), 249–261.
Vaughan, D.G., Cor, H.F.J.., Doak, C.S.M.. and Waddingto, E.D... 1999. Distortion of isochronous layers in ice revealed by ground-penetrating rada.. Nature, 398(6725), 323–326.
Waddington, E.D., Neuman, T.A.., Koutni, M.R.., Marshal, H.-P.. and Mors, D.L... 2007. Inference of accumulation-rate patterns from deep layers in glaciers and ice sheet.. J. Glaciol., 53(183), 694–712.
Wendt, J. and 9 others. 2006. Geodetic observations of ice flow velocities over the southern part of subglacial Vostok Subglacial Lake, Antarctica, and their glaciological implication.. Geophys. J. Int., 166(3), 991–998.