Introduction

Vostok lake is the largest known subglacial lake in East Antarctica, about 200 km downstream from the Ridge B ice divide. The lake is approximately 260 km long and 80 km wide (Reference Tabacco, Bianchi, Zirizzotti, Zuccheretti, Forieri and Della VedovaTabacco and others, 2002), with an estimated volume of 5200 km³ (M. Studinger and others, unpublished data). Subglacial melting occurs where ice is thicker (the north), and lake water refreezes to the ice-sheet underside where the ice is thinner (the south). The refrozen or ‘accreted’ ice found at the bottom of the Vostok ice core (VIC) is effectively a sample of the lake water, modified by processes which occur during freezing.

Reference JouzelJouzel and others (1999) noted that the accreted ice δD and δ¹⁸O compositions do not lie on the ‘Vostok precipitation line’, which is produced by a well-obeyed linear relationship between δD and δ¹⁸O in meteoric ice that forms the top part of an ice core at Vostok station (Fig. 1). The apparent enrichment of heavy isotopes in the accreted ice (AI) was attributed to the incorporation of lake water into large frazil ice crystals during freezing. The amount of trapped water is large, accounting for some 30–58% of the total mass of ice crystals. It follows from simple isotope mass balance that differences between the lake’s inputs and outputs (Table 1) must be compensated by a change in the lake water composition. We can therefore consider whether or not the lake is in steady state with respect to δD and δ¹⁸O using box modelling.

Fig. 1. Vostok station ice-core δD and δ¹⁸O data for both the meteoric ice (0–3352 m) and accreted ice (3539–3610 m). Data were provided by J. Jouzel and F. Vimeux.

Table 1 Observed and modelled δD and δ¹⁸O values

Steady state is achieved in the δD and δ¹⁸O composition of Vostok lake if (i) the annual input and output of D, H, ¹⁸O and ¹⁶O, via melt and accreted ice formation respectively, is relatively constant compared to the mass of these species in the lake, and (ii) the volume of the lake is constant and well mixed. Reservoirs with simple inputs and outputs are able to reach a new steady state within three to four residence times after any perturbation has occurred (Reference LermanLerman, 1979). Estimates of the residence time of Vostok lake range from 36 to 300 kyr (after Reference SiegertSiegert and others, 2001; Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002; M. Studinger and others, unpublished data).

The notion of water pockets being incorporated into growing ice crystals (Reference JouzelJouzel and others, 1999) makes the amounts of gas and solute hard to account for. (The mass-balance calculations assume that the gas concentration and solute content of the VIC are comparable to those of ice melt feeding Vostok lake.) The average gas content of the VIC is ∼90cm³ kg^–1, whilst that of the AI is ∼1–2cm³ kg^–1, almost two orders of magnitude lower. Gas from the VIC must therefore be accumulating in the lake. The dissolved gas content of lake water is believed to build up to ∼30 times that of the VIC before clathrate hydrates become stable (Reference Lipenkov and IstominLipenkov and Istomin, 2001). Assuming that all gas is rejected from frazil ice during formation, any gas in the AI must be derived from trapped lake water. The maximum amount of water captured by growing ice crystals is ∼1–2% if the lake water has a gas content similar to the VIC, whereas a value of 0.03–0.07% is more realistic if sufficient gas has accumulated in the lake to form clathrate hydrates. These estimates are more than an order of magnitude less than that required by the steady-state isotope mass balance. Chloride (Cl^–), considered to be the most conservative major ion, has an average concentration of 4.4 μeq L^–1 in the VIC compared to 0.8–1.1 μeq L^–1 in the cleanest portions of AI (Reference Siegert, Tranter, Ellis-Evans, Priscu and LyonsSiegert and others, 2003). Assuming that all Cl^– is excluded from the growing frazil ice, such that Cl^– in the AI is derived from trapped water, the most lake water that can be trapped in the frazil ice is 18–25% (assuming that ice melt is the only source of Cl^– to the lake water). This estimate is a factor of approximately two less than that required by stable-isotope mass balance. Particulate sources of Cl^– to the lake water or to the clean AI will reduce this percentage.

This brief analysis leads us to believe that trapping of water in frazil ice crystals is an unlikely explanation for the difference between VIC and AI δD and δ¹⁸O. Instead, in this paper, we adopt a simple box modelling approach to synthesize the AI composition from a lake input reconstructed using the known VIC chemistry.

Model Description

Vostok lake can be conceptualized as a simple box model (Fig. 2). The body of the lake is the reservoir which receives ice melt (flow-in) and loses water when it freezes to the ice ceiling (flow-out). The model assumes that flow-in balances flow-out, so that lake volume remains constant. The model is driven by the estimated δD and δ¹⁸O of ice melting into the lake. This ice originates close to Dome B which is approximately 200km north of Vostok station. We are able to reconstruct an ice core for Dome B by applying the north–south isotope gradient observed in surface snows by Reference Qin, Petit, Jouzel and StievenardQin and others (1994) to the VIC. Using the difference between surface snow values at Vostok station and the average value for two stations, 217 and 267 km north on a transect to Komsomolskaya, leads to enrichments of 18.2‰ and 2.7‰ for δD and δ¹⁸O respectively (Table 1). This north–south correction is applied to each pair of δD and δ¹⁸O values for the VIC (supplied by J. Jouzel and F. Vimeux) to produce a hypothetical Ridge B ice core which flows (and melts) into Vostok lake. The isotope gradient for contemporary surface snows and ice deposited during the Last Glacial Maximum is similar, suggesting a conservative gradient with respect to climatic changes (J. Jouzel and others, unpublished information). We therefore assume that the gradient has not changed significantly during the 460 kyr record represented by the VIC. The sensitivity of the north–south correction is examined below, since the precise location of the melting area above Vostok lake is not known.

Fig. 2. Box model conceptualization.

The east–west isotope gradient in this area is believed to be small (personal communication from J.-R. Petit, 2003), and no correction is made in the first instance to account for the east–west offset between the Vostok-station–Komsomolskaya transect and Ridge B. However, an isotopic gradient might be expected between Vostok station and Ridge B, since there is an elevation difference of >200 m. We explore the possible effects of different east–west isotope gradients below.

The initial lake-water composition is set to the average input over the last four glacial cycles from the reconstructed Dome B record (Table 1). At each model time-step a portion of ice melts into the lake and an equivalent portion of ice accretes to the ice ceiling. There is no fractionation of H and O isotopes during melting, but there is fractionation during freezing, with the heavier isotopes being preferentially incorporated into the ice (Reference JouzelJouzel and others, 1999). This is described by the equilibrium fractionation factor, αx , such that (1000 + δ _ice X) = αx(1000 + δ _water X), where X denotes ¹⁸O or D. In reality, melting and freezing are spatially discrete, and it is envisaged that melting occurs over one-third of the lake’s surface area, with freezing over the other two-thirds (Reference StudingerStudinger and others, 2003). We use these values to notionally set the relative rates of melting and freezing.

The lake area is taken as 14 000km² (Reference Kapitsa, Ridley, Robin, Siegert and ZotikovKapitsa and others, 1996), and the volume is estimated as 5200 km³ (M. Studinger and others, unpublished data). We assume a spatially averaged melt rate of 4 cm a^–1 (after Reference Mayer, Grosfeld and SiegertMayer and others, 2003), which forces a freezing rate of 2 cm a^–1. Hence, the implied water flux is 0.19 km³ a^–1, giving a residence time of ∼28 kyr assuming the lake is fully mixed (homogeneous). This is a possibility on time-scales of this order, since both vertical and horizontal water motion is relatively fast (Reference Wüest and CarmackWüest and Carmack, 2000). Theoretical and experimental fractionation factors vary from 1.0171 (± 0.005) to 1.0235 and from 1.00048 to 1.0031 for δD and δ¹⁸O respectively (Reference Moser, Stichler, Fritz and FonteMoser and Stichler, 1980). We have chosen the higher values because the accretion of ice is slow enough for isotopic equilibrium to be reached. These values were used by Reference JouzelJouzel and others (1999), facilitating comparison between studies. Standard model variables are listed in Table 2.

Table 2 Standard model variables

The location of subglacial melting over the lake determines the back trajectory of the inflowing ice, and hence the origin and isotopic composition of this ice on Ridge B. The sensitivity of model results to the location of the ice in the vicinity of Ridge B was assessed by running the model with ice cores reconstructed from the VIC, but with north–south corrections derived from locations 217 km (δD=+13.2 and δ¹⁸O=+2.0) and 267 km (δD=+23.4 and δ¹⁸O=+3.4) north of Vostok station using surface snow data from Reference Qin, Petit, Jouzel and StievenardQin and others (1994), rather than the average of these two sites as used above (Fig. 3a and b).

Fig. 3. The effect of a north–south isotopic gradient on accreted ice compositions. Corrections derived from (a) 217 km north (δD=+13.2 and δ¹⁸O=+2.0) and (b) 267 km north (δD=+23.4 and δ¹⁸O=+3.4). The model assumes that no east–west isotopic gradient exists.

In the absence of published data, we show the likely impact of an east–west isotope gradient using the deuterium content as predicted by the NASA/GISS model of Reference Jouzel, Russell, Koster, Suozzo, White and BroeckerJouzel and others (1987) (from Reference Qin, Petit, Jouzel and StievenardQin and others, 1994, Fig. 4b). The δD east–west gradient in the vicinity of Vostok station is approximately half that in the north–south direction (7.7‰). To examine the sensitivity of this gradient, the model was run with δD gradients of 7.7‰, 5.8‰, 3.9‰ and 1.9‰ per 100 km (Fig. 4a–d). δ¹⁸O corrections were derived from the δD/δ¹⁸O relationship in surface snow data in Reference Qin, Petit, Jouzel and StievenardQin and others (1994) (δD =–22.1140 + 7.3138 δ¹⁸O; R ²=0.94). Since Dome B is approximately 200 km west of Vostok station, this gave maximum corrections of 15.4, 11.6, 7.8 and 3.8‰ for δD. Each depth of the VIC was first corrected for the north–south gradient (using the average correction from sites 217 and 267 km north), and then for the east–west gradient by adding a proportion of the maximum correction based on the age of the ice, a proxy for the distance that ice has travelled to Vostok station. In this way, no east–west correction is made to the bottom part of the core, which is assumed to have originated as precipitation near Ridge B, but the maximum correction is made at the top of the core.

Fig. 4. The effect of east–west isotopic gradients on accreted ice compositions. δD gradients of 7.7‰ (a), 5.8‰ (b), 3.9‰ (c) and 1.9‰ (d) per 100 km. δ¹⁸O corrections were derived from the δD/δ¹⁸O relationship in surface snow data in Reference Qin, Petit, Jouzel and StievenardQin and others (1994) (see main text for equation).

Results

The δD and δ¹⁸O values of the modelled accreted ice (MAI) using the standard model are compared with the 71 samples of accreted ice in Figure 5. The MAI becomes progressively depleted in the heavy isotopes, moving from the top right towards the bottom left of Figure 5. The ‘jagged’ nature of the lines is an artefact of the variable incoming ice composition. The standard model (Fig. 5) eventually reaches steady state, where the modelled ice composition ‘bounces’ at the bottom left end of the line around the average value of the Ridge B input (δD=–442.2 and δ¹⁸O=–56.7). The data point for 3539m is a clear outlier beyond the range of modelled ice compositions. This depth marks the transition from deformed meteoric to accreted ice, so this ice may not be refrozen lake water and can be ignored when considering whether the lake has reached steady state. Figures 3 and 4 present the results of the sensitivity of the model to variations in north–south and east–west isotope correction.

Fig. 5. Real and modelled accreted ice compositions, assuming that the δD and δ¹⁸O of ice flowing and melting into Vostok lake is 18.2‰ and 2.7‰ heavier than the VIC respectively, and that no east–west isotopic gradient exists.

Discussion

Our results imply that in the absence of an east–west isotopic gradient, Vostok lake has reached isotopic steady state (Fig. 5). This suggests that the lake has not been subject to any recent major perturbations, such as volume changes. A volume increase caused by net melting would be expected to enrich the lake with δD and δ¹⁸O, while lake contraction caused by net freezing would leave the lake depleted in δD and δ¹⁸O. Hydrothermal input to the lake has been suggested by several authors (e.g. Reference StudingerStudinger and others, 2003). However, Figure 5 argues against a significant hydrothermal input. This would enrich the lake in δD and δ¹⁸O, depending on the crustal composition and reaction types in vents. The results also argue against the possibility of a significantly warmer or cooler climate before ∼460 kyr ago, since this would produce isotopically enriched and depleted precipitation respectively. This different input would be reflected in the accreted ice as the lake evolved towards steady state.

Varying the north–south isotope correction has a significant effect on the modelled steady-state composition of the AI and, thus, determining whether or not the lake has reached steady state. An input originating 217km north suggests that the AI should be lighter than observed at steady state (Fig. 3a), whereas an input 267 km north implies that the observed AI is too light to represent a system at steady state (Fig. 3b). Although the melting area, and hence the origin of this ice along Ridge B, are unlikely to be confined to a small geographical area, it is possible that extensive areas of less rapid freezing are balanced by small areas with high melt rates (Reference StudingerStudinger and others, 2003). In this way, the lake could be supplied by heavier or lighter ice than we assume in the standard model.

Introducing an east–west isotopic gradient has a noticeable effect on the modelled evolution of accreted ice towards steady state. Providing that this gradient is less than half the north–south gradient, the model continues to predict that the lake is in steady state (Fig. 4c and d). When the east– west gradient is three-quarters, or the same as the north– south gradient, the model predicts that the lake is not in steady state. Specifically, the accreted ice is isotopically lighter than the modelled steady-state composition. This anomaly could be explained by contraction of the lake or a colder climate before ∼460 kyr ago.

Contraction of Vostok lake by net freezing would cause the lake to become depleted in heavy isotopes, since these are preferentially incorporated into accreted ice during its formation. Changes in the lake’s volume require an external forcing process. The most likely is from the overriding ice sheet, which is known to change subtly over glacial–interglacial cycles in two ways. First, the ice-sheet surface (and therefore thickness) reduces during full glacial periods and increases during interglacials, as a consequence of the decrease and increase in ice accumulation rates, respectively. Reference Ritz, Rommelaere and DumasRitz and others (2001) calculate the variation in elevation at Vostok station over the last glacial cycle to be about 150 m. Second, the ice-sheet flow direction also changes slightly over glacial–interglacial cycles. Evidence in support of such changes across Vostok lake comes from internal radar layer structures that have been mapped out across the southern end of Vostok lake (Reference Bell, Studinger, Tikku, Clarke, Gutner and MeertensBell and others, 2002), which depict the ice-flow direction at the time it crossed the lake’s western margin, but which are at a distinct angle to the current flow direction measured by interferometric synthetic aperture radar (Reference Kwok, Siegert and CarseyKwok and others, 2000). Ice-sheet elevation and flow-direction changes have consequences for the lake volume, as they will affect the slope of the ice sheet crossing the lake’s western edge, the ice surface slope over the lake itself and, importantly, the slope of the ice/water interface above the lake. As the slope of the ceiling of Vostok lake has to be inclined at ten times the ice surface slope (for the water to be in hydrostatic equilibrium), even small changes to the ice surface will have a significant effect on the slope of the ice/water interface and, thus, the lake volume.

Conclusion

Simple box modelling of δD and δ¹⁸O in Vostok lake accretion ice has shown that it is in isotopic steady state if the current conceptualization of the system is correct. However, varying the north–south location of ice input to the lake in our model is critical to the expected steady-state composition of the lake. This representation of Vostok lake is only relevant if subglacial melting is highly concentrated in small areas, or if the ice sheet above the northern third is not the dominant input. If an east-to-west isotopic gradient exists, with a magnitude similar to that observed in a north– south transect, the lake may not be in steady state. In both cases, the most likely explanation for non-steady state is a periodic volume change induced by ice-thickness or ice-flow-direction changes.