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. Reference Blankenship, Alley and BindschadleBlankenship and others (2001) described the system characteristics. This dataset was used to derive both the surface and bed topographies (Reference StudingerStudinger 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 Reference Tikku, Bel, Studinge and ClarkTikku 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 (Reference Popov, Sheremet’ye, Masolo, Luki, Mirono and LuchininoPopov 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 (Reference Studinger, Bel and TikkStudinger 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 (Reference Parrenin, Rém, Rit, Sieger and JouzeParrenin 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 (Reference Waddington, Neuman, Koutni, Marshal and MorsWaddington 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. Reference SiegertSiegert, 2003). However, as also acknowledged by Reference SiegertSiegert (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, Reference Waddington, Neuman, Koutni, Marshal and MorsWaddington 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 (L_H ) 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; Reference WendtWendt 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 L_H 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 Reference Waddington, Neuman, Koutni, Marshal and MorsWaddington 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 Reference Tikku, Bel, Studinge and ClarkTikku 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 Reference WendtWendt 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 (Reference PatersonPaterson, 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; Reference WendtWendt 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 (Reference Thoma, Grosfel and MayeThoma 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. Reference Thoma, Grosfel and MayeThoma and others, 2008), shallow more abruptly than those layers in the southern part of the lake, where basal accretion has been observed (Reference Bell, Studinge, Tikk, Clark, Gutne and MeertenBell and others, 2002; Reference Tikku, Bel, Studinge and ClarkTikku 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 (Reference Parrenin, Rém, Rit, Sieger and JouzeParrenin 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 Reference Waddington, Neuman, Koutni, Marshal and MorsWaddington 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; Reference WendtWendt and others, 2006). Long-term accumulation-rate data are sparse (Reference MagandMagand 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; Reference Lipenkov, Yekayki, Barko and PurshLipenkov 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, Reference Parrenin, Rém, Rit, Sieger and JouzeParrenin 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 Reference Leonard, Bel, Studinge and TremblaLeonard 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.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 (Reference Arthern, Winebrenne and VaughaArthern 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 Reference Leonard, Bel, Studinge and TremblaLeonard 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 (Reference Leonard, Bel, Studinge and TremblaLeonard 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_ḃ/L_H varies between 2.7 and 5.2 for the different layers. The smaller value of either L_ḃ or L_H 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 (Reference Parrenin, Rém, Rit, Sieger and JouzeParrenin 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.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

Reference Arthern, Winebrenne and VaughaArthern 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 Reference Qin, Peti, Jouze and StievenarQin and others (1994) are in better agreement with our ḃ_LLA map for layer A than the map of Reference Arthern, Winebrenne and VaughaArthern and others (2006). Although the absolute accumulation-rate values differ between our map and that of Reference Arthern, Winebrenne and VaughaArthern and others (2006), there is a similar pattern of increasing values to the north in our study area. However, the map of Reference Arthern, Winebrenne and VaughaArthern 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.

Reference SiegertSiegert (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 Reference SiegertSiegert (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 Reference SiegertSiegert (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 Reference SiegertSiegert (2003) and Reference Leysinger Vieli, Sieger and PaynLeysinger 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).

Reference Leonard, Bel, Studinge and TremblaLeonard 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 (Reference Qin, Peti, Jouze and StievenarQin 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 (Reference Vaughan, Cor, Doak and WaddingtoVaughan and others, 1999; Reference King, Anderso, Vaugha, Man, Mobb and VospeKing 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 (Reference Tikku, Bel, Studinge and ClarkTikku 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. Reference Pattyn, de Smed and SouchePattyn 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. Reference Lipenkov, Yekayki, Barko and PurshLipenkov, 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. Reference Price, Waddingto and ConwaPrice 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. Reference Steen-Larsen, Koutni and WaddingtoSteen-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 (Reference Leysinger Vieli, Sieger and PaynLeysinger Vieli and others, 2004; Reference Salamatin, Tsyganov, Popo, Lipenko and HondohSalamatin 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. Reference Salamatin, Tsyganov, Popo, Lipenko and HondohSalamatin 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.