Research ArticleOCEANOGRAPHY

Dynamic flows create potentially habitable conditions in Antarctic subglacial lakes

See allHide authors and affiliations

Science Advances  17 Feb 2021:
Vol. 7, no. 8, eabc3972
DOI: 10.1126/sciadv.abc3972


Trapped beneath the Antarctic ice sheet lie over 400 subglacial lakes, which are considered to be extreme, isolated, yet viable habitats for microbial life. The physical conditions within subglacial lakes are critical to evaluating how and where life may best exist. Here, we propose that Earth’s geothermal flux provides efficient stirring of Antarctic subglacial lake water. We demonstrate that most lakes are in a regime of vigorous turbulent vertical convection, enabling suspension of spherical particulates with diameters up to 36 micrometers. Thus, dynamic conditions support efficient mixing of nutrient- and oxygen-enriched meltwater derived from the overlying ice, which is essential for biome support within the water column. We caution that accreted ice analysis cannot always be used as a proxy for water sampling of lakes beneath a thin (<3.166 kilometers) ice cover, because a stable layer isolates the well-mixed bulk water from the ice-water interface where freezing may occur.


The Antarctic continent is covered with ice, growing and shrinking over periods of tens to hundreds of thousands of years, since at least the last 14 million years (1). Over 250 hydrologically stable subglacial lakes (in which water inputs are constantly balanced by outputs) trapped between the bed and the ice are known to exist at and close to the ice sheet center (2). They comprise a wide variety of sizes and glaciological and topographic settings (3) and have been hypothesized as potential habitats for the in situ development of microbial organisms (4). Such remote, extreme, and isolated places qualify as analogs to extraterrestrial environments where life may occur, such as the subsurface oceans on Jovian and Saturnian moons (5). A further ~130 hydrologically active lakes, which experience rapid water discharges and large volume changes, exist toward the margin of the ice sheet (6, 7). While these may contain microbial life (8), they are not considered as isolated habitats where microbes can adapt independently over long periods due to the flushing of water in and out of their systems and their potentially ephemeral nature.

The Antarctic ice and bed material carry life’s building blocks, with oxygen and minerals held within dust in the former, and minerals trapped inside sediments and bedrock in the latter. Numerical models and radar observations have shown that the ice sheet base above subglacial lakes typically melts where the ice is thickest and freezes where it is thinnest (4, 9, 10). Thus, oxygen and minerals are released at the top of the water column. The rate at which this happens is key to assessing the possibility of having a biome, but remains largely uncertain. Although microbial life is anticipated at the floors of subglacial lakes, where sediments are known to exist (11), dynamic flows and mixing of bottom water within the water column are essential for life to be widespread and detectable, avoiding, for example, anoxic conditions if oxygen-rich surface water is unable to access deeper parts of the lake.

Subglacial lakes are isolated from winds and solar heating but can experience vertical convection flows due to the upward geothermal flux [at a background level of roughly 50 mW/m2; (12)], and horizontal convection flows due to the ubiquitous—albeit variable—tilt of their ice ceiling (about 10 times and in opposite direction to the ice surface slope). Previous work has estimated that velocities of few tenths of a millimeter/second are generally required to suspend sediments in the water column (13). This is of the same order of magnitude as that predicted by ocean modeling for a handful of subglacial lakes, including Lake Vostok (14), Lake Ellsworth (15), and Lake Concordia (9). However, uncertainties are large, and velocities remain unknown for most subglacial lakes, including Lake CECs, which might be the first stable lake to be drilled into in a clean way in the coming years (16). As a result, plans for direct sampling can be helped by establishing hydrological conditions in subglacial lakes, and their variation between lake settings, to recognize where microbial life is most likely to thrive.

Here, we predict the intensity of turbulence and large-scale water circulation for the entire range of stable subglacial lakes found in Antarctica, i.e., with ice cover thicknesses up to 5 km and lake water depths up to 1.5 km (Fig. 1). Thus, our work complements previous studies on convection in lakes at atmospheric pressure (i.e., open) or with thin ice covers (17), and, more specifically, previous efforts that aimed to predict the hydrological conditions of individual subglacial lakes, including Lake Vostok (18). We demonstrate that most subglacial lakes have large supercritical convective parameters, i.e., the geothermal flux is much larger than the minimum critical heat flux required to trigger convective flows, such that they are in a regime of vigorous turbulent convection. We show that vertical convection is as important as horizontal convection and that the convective dynamics vary considerably based on the ice thickness, water depth, and ceiling slope. For simplicity, we restrict our attention to freshwater, because salt concentration is typically low in isolated subglacial lakes (4).

Fig. 1 Problem schematic.

We provide predictions about the characteristic velocity of the large-scale circulation Ulsc, the characteristic velocity of turbulent plumes U, the thickness δ of the top stable conductive layer, and the anomalous temperature of the well-mixed bulk Tbulk (i.e., in excess of the freezing temperature Tf). The problem parameters are the water depth h, the ice thickness H (or ice overburden pressure pi), the Coriolis frequency f (due to Earth’s rotation), and the geothermal flux F.

We first calculate the minimum critical heat flux, Fc, required to trigger thermally forced vertical convection (Fig. 3) by solving an eigenvalue problem for the local stability of subglacial lakes with a nonlinear equation of state (19). We show that Fc is much smaller than 50 mW/m2, which is (approximately) the average geothermal flux, for a wide range of geophysical conditions and conclude that most Antarctic subglacial lakes are unstable to convection. We then demonstrate that most subglacial lakes (Figs. 4 and 5) subject to a geothermal flux of 50 mW/m2 experience dynamic flows by applying state-of-the-art scaling laws of classical thermal convection to convection in cold-temperature high-pressure lake environments.


Subglacial lakes can experience dynamic vertical convection flows [also known as Rayleigh-Bénard (RB) convection (20)], because the lake’s deepest waters, heated by Earth’s geothermal flux, are generally buoyant and will tend to rise through and mix with the rest of the water column. How far up and how quickly bottom water masses rise depend on the geothermal flux, F, the water depth, h, and the ice cover thickness, H (or ice overburden pressure pi). Convection in subglacial lakes is complex, because the thermal expansion coefficient of freshwater α(p, T), which indicates how fluid parcels contract or expand with changes in temperature (i.e., are buoyant), depends on water pressure p > pi and the temperature T itself (21). For relatively thick ice cover, i.e., HH* = 3166 m, the thermal expansion coefficient is always positive (i.e., the density decreases with temperature) and increases with pressure and temperature, such that convection becomes more vigorous as F, h, and H increase. For ice covers less than the critical ice depth H*, or ice pressure pi < p* = 2848 dbar (which we refer to as the critical ice pressure), however, the thermal expansion coefficient changes sign with temperature, such that density does not simply decrease with temperature but becomes a nonlinear and nonmonotonic function of T. Specifically, as is shown in Fig. 2, for H < H*, α increases with temperature but is first negative for Tf(pi) ≤ T < Td(p) (close to the ice ceiling), where Tf is the freezing temperature and Td is the temperature of maximum density, before becoming positive at higher temperature. Having α < 0 close to the ice ceiling for H < H* means that the density stratification is always stable at the top of the lake and that the bottom layer is buoyant only if the bottom temperature exceeds Td, i.e., such that the density stratification is top heavy near the bottom. Having α > 0 on the bottom boundary is a necessary condition for deep water masses to be buoyant but, however, is not sufficient for convection to set in. The geothermal flux must be also larger than the adiabatic heat flux and adequate to sustain a buoyancy force that can overcome viscous dissipation and thermal diffusion.

Fig. 2 Thermal expansion coefficient.

Plot of the thermal expansion coefficient α as a function of (T, p) superimposed with profiles of the temperature of maximum density Td (red solid line) and freezing temperature Tf (black solid line) with pressure. For small pressures p < p*, with p* the critical inversion pressure (blue dashed line), Td > Tf such that there exists a range of temperatures Tf < T < Td for which α is negative (area appearing with red colors) and water masses become anomalously denser with increasing temperatures. For p < p* and T > Td, or pp*, the water becomes monotonically lighter as temperature increases, which is the typical behavior of most fluids.

Here, we estimate the minimum critical heat flux that overcomes dissipation effects and permits convection in subglacial lakes from a stability analysis of the perturbation equations for a water column subject to geothermal heating and the Coriolis force due to Earth’s rotation. We consider a realistic nonlinear equation of state for freshwater using the Thermodynamic Equation of Seawater 2010 (TEOS-10) toolbox (19) and the Coriolis frequency at 80°S. We take into account the adiabatic temperature gradient by including compressibility effects in the energy equation. We perform the calculations for a wide range of ice thicknesses and water depths up to 20 m. For water depth, h > 20 m, the eigenvalue problem becomes too difficult to solve, so we use scaling laws that are either conservative or inferred from classical RB convection results in the limit of rapid rotation (22).

Figure 3A shows the minimum critical heat flux Fc that permits vertical convection for a wide range of ice pressures and water depths relevant to Antarctic subglacial lakes. Fc is large at small pressure and small water depth (top left corner) but then decreases with h and pi in most of the parameter space. We find Fc < 50 mW/m2 (as shown by the black isocontour labeled “50”), which is a typical background value for Earth’s geothermal flux around Antarctica, in most of the parameter space. We predict that Lake CECs and South Pole Lake (SPL) are unstable to vertical convection if subject to a 50 mW/m2 flux, i.e., their critical heat flux Fc is less than 50 mW/m2, despite having relatively thin ice covers H < H* (shown by the gray dashed line). Here, we have centered the vertical axis of Fig. 3A on the critical pressure p* by using the shifted ice pressure variable pip* and a symmetric logarithmic scale. As a result, the transition from a fully convective water column (for pi > p*) to a convective water column with a stably stratified upper layer (for pip*) is smooth even though Fc increases rapidly as pi decreases below p*. Figure 3B shows that subglacial lakes reported in the last inventory (2) have ice cover thicknesses almost equally distributed on either side of H*. Thus, the pi = p* isobar, which separates lakes that are fully convective from lakes that are only partially convecting, is important not only for the theoretical calculation of Fc but also in practice. Almost half of the subglacial lakes (with pi < p*) may be expected to have a top layer that is stable, although possibly modified by the dynamics near the ice ceiling and overshooting convection. We remark that Fc is constrained primarily by the condition of having α > 0 on the bottom boundary for small pi and h (top left corner) and by the condition that it must exceeds the adiabatic flux for large pi and h (bottom right corner). In between, viscous dissipation and thermal diffusion dominate the calculation of Fc. Note that our prediction of Fc for small pi and large h is conservative and may overestimate the true Fc. We provide further details about the calculation of Fc in Materials and Methods and in the Supplementary Materials.

Fig. 3 Critical heat flux.

(A) Minimum heat flux required to trigger vertical convection in subglacial lakes as a function of lake depth (bottom axis) and ice overburden pressure (left axis) or ice thickness (right axis). Solid lines are isocontours in mW/m2 of required heat flux, while filled circles highlight the positions of five well-known lakes in parameter space (see legend to the right). (B) Ice thickness distribution of isolated subglacial lakes from the last published inventory (2). The dashed gray lines highlight the critical thickness H*.

For a geothermal flux F greater than the critical heat flux Fc, it is of interest to know whether the convective instability results in high or low flow velocities. In general, estimates of flow velocities require dedicated simulations or laboratory experiments. For the case of turbulent vertical convection, however, various scientific communities have proposed predictive laws of hydrodynamic variables based on control parameters, which can be used as leading-order estimates for turbulence intensity in convective subglacial lakes [see (20) and references therein]. The canonical problem of natural convection relevant to our work is known as rotating RB convection and has applications in many different fields, including geophysics (23), astrophysics, and engineering (20). Hydrodynamic variables such as turbulent flow velocities, large-scale flow velocities, and temperature fluctuations are predicted on the basis of the value of the Rayleigh number Ra of the system, which is a dimensionless measure of the available convective energy; the Prandtl number Pr, which compares viscous dissipation to thermal diffusion; and the Ekman number Ek, which compares viscous dissipation to the Coriolis force. The Rayleigh number, Prandtl number, and Ekman number for subglacial lakes in Antarctica can be written asRaF=gαeffFheff4νκk,Pr=νκ,Ek=νfheff2(1)where g is the gravitational acceleration, αeff is the characteristic thermal expansion coefficient, heff is the effective water depth where convection occurs, ν is the kinematic viscosity, κ is the thermal diffusivity, k is the thermal conductivity, and f is the Coriolis frequency (note that we use ∣f∣ in our definition of Ek > 0, since f < 0 in the Southern Hemisphere). The subscript F of RaF means that the Rayleigh number of subglacial lakes is a flux-based Rayleigh number, since it is based on a prescribed geothermal flux rather than a prescribed temperature difference, which is more common in idealized studies of natural convection (20). Note that we neglect compressibility effects for the prediction of variables in the turbulent regime because Earth’s geothermal flux is several orders of magnitude larger than the adiabatic heat flux.

In the context of subglacial lakes, the geothermal flux is sufficiently large that the lake water is in a fully turbulent state that is almost not affected by rotation, i.e., FFc, and the effect of rotation is weak [see the Supplementary Materials and (24)]. Thus, here we use scaling laws derived for fully turbulent nonrotating convection to make predictions about hydrodynamic variables. The variables of interest are the thickness of the conductive layer near the ice ceiling δ, the anomalous temperature of the well-mixed bulk Tbulk (in excess of Tf), the characteristic turbulent flow velocity U, and the length scale of turbulence 𝓁, which represents the typical distance between thermal plumes (Fig. 1). We assume a geothermal flux of F = 50 mW/m2 throughout and use scaling laws derived from numerical simulations (25) as well as scaling laws inferred from the Grossmann-Lohse (GL) unifying theory of RB convection, which is based on theoretical arguments (26).

We first show in Fig. 4 (A to D) the results for δ and Tbulk based on the scaling laws derived in (25) for a wide range of ice thicknesses and water depths. The thickness of the conductive layer near the ice ceiling is almost independent of lake depth but varies substantially with ice pressure (Fig. 4, A and B). For thin ice cover, the top conductive layer consists of a layer with a stable density stratification attached to the ice ceiling (where α < 0), which can be several meters thick (thickness δS), and a turbulent transition layer (just above the convective bulk), which is typically on the order of a few centimeters or smaller (thickness δt), i.e., δ = δS + δt. For ice thickness, H < 2000 m, δS is between 10 and 40 m. The stable layer thickness δS decreases with H and vanishes for H > H* (since α > 0 everywhere in this case), such that the full conductive layer is limited—for a thick ice cover—to a turbulent boundary layer attached to the ice ceiling, which is small. In all cases, the top stable layer transfers heat by conduction only. Thus, the temperature increases linearly from Tf(pi) near the ice to Tf(pi) + δF/k at the bottom of the stratified layer. The anomalous temperature of the well-mixed convective bulk (above freezing) can then be approximated as Tbulk = δF/k (Fig. 4, C and D), and hence shows similar trends as δ. The white area on the top left corners of Fig. 4 (A and C) highlights subglacial lakes that are stable because the thermal expansion coefficient is negative everywhere in the water column. The filled squares in Fig. 4 (B and D) show δ and Tbulk based on scaling laws derived in (26) (labeled “GL”). There is a good agreement between predictions based on the scaling laws in (25) (shown by circles) and (26).

Fig. 4 Conductive layer thickness and anomalous bulk temperature.

(A and B) Thickness δ of the conductive stratified layer at the top of subglacial lakes assuming a geothermal flux of 50 mW/m2. (A) δ as a function of lake depth (bottom axis) and ice pressure (left axis). (B) δ as a function of ice pressure only for selected lake depths of 32, 156, and 1067 m [shown by vertical lines in (A)]. GL refers to results obtained with the GL theory and h = 1067 m. (C and D) Same as (A) and (B) but for the anomalous bulk temperature Tbulk (above Tf) of the well-mixed convective layer.

Figure 5 shows the predicted lake velocity U and horizontal length scale 𝓁 based on previously derived scaling laws (25). Figure 5 (A and B) shows that U is almost independent of ice pressure but increases with lake depth, up to about 1 cm/s for h = 1500 m. Figure 5 (C and D) shows that 𝓁 increases slowly with lake depth and remains on the order of 1 m for all ice pressures. Both U and 𝓁 appear discontinuous at pi = p* because the effective thermal expansion coefficient αeff, which we estimate conservatively (cf. Materials and Methods) and use in Eq. 1 for RaF, decreases rapidly across the p* isobar for small water depths. For instance, αeff decreases from 3 × 10−6 ° C−1 to 3 × 10−7 ° C−1 between pi = p* + 100 dbar and pi = p* for h = 10 m. We expect that the discontinuity would become less sharp but would not completely disappear upon relaxing our conservative approximation for αeff, because α will always be (overall) much smaller in lakes with a thin ice cover (pip*) than in lakes with a thick ice cover (pi > p*). We also show in Fig. 5B the prediction for the characteristic velocity Ulsc based on the scaling laws derived in (26) (labeled GL) for an ice thickness H = 3945 m (filled squares), which is the ice thickness above Lake Vostok. Ulsc is smaller than U (shown by the filled circles) by up to a factor 5 because the GL theory focuses on the velocity of the large-scale circulation, while U is the characteristic root-mean-square velocity of turbulent plumes (25), which is likely to be faster than the mean large-scale flow. Figure 5B also shows a prediction for the horizontal velocity Vhc of the baroclinic flow expected along a sloped ice-water interface, using scaling laws inferred from recent results on horizontal convection (27). We show the prediction for Vhc assuming two different ice-water interface slopes, i.e., s = 10−3 and s = 10−2, the steepest slope resulting in the largest horizontal velocity due to the increased temperature gradient along the ice ceiling. The horizontal velocity of the baroclinic flow is of the same order (for a steep slope, s = 10−2) or one order of magnitude smaller (for a moderate slope, s = 10−3) than the large-scale velocity of vertical convection (Ulsc).

Fig. 5 Characteristic turbulent flow velocity and length scale.

(A and C) Same as Fig. 4 (A and C), but for (A) the turbulent flow (plume) velocity U and (C) the characteristic length scale 𝓁 in the convective layer. (B and D) Turbulent flow velocity U and length scale 𝓁 as functions of lake depth only, for selected ice thicknesses H = 3945, 2653, and 1000 m [shown by horizontal lines in (A) and (C)]. GL refers to the GL predictions for the large-scale velocity Ulsc of vertical convection for H = 3945 m (shown by filled black squares). In (B), we also show a prediction for the horizontal velocity Vhc of the baroclinic flow along a tilted ice-water interface, assuming either a steep slope s = 10−2 (green stars) or a moderate slope s = 10−3 (tilted blue crosses).


We have demonstrated that the critical heat flux leading to vertical convection in subglacial lakes is much less than 50 mW/m2 for a broad range of ice overburden pressures and water depths (Fig. 3). Thus, it should be considered that most—if not all—Antarctic subglacial lakes are dynamic hydrologic environments. We expect that the same conclusion holds for isolated subglacial lakes in Greenland and elsewhere in the solar system (5, 28). We note that our prediction of the critical heat flux Fc is conservative for large water depth and small ice pressure (see the Supplementary Materials). Also, estimates of Fc exceeding 50 mW/m2 (as is the case for a lake 20 m deep and under 1 km of ice), hence suggesting stable subglacial lakes, could be verified qualitatively through direct sampling and revised if measurements demonstrate a dynamic environment.

Vertical convection in subglacial lakes is different from vertical convection in the canonical RB problem, mainly because the thermal expansion coefficient (α) of freshwater in subglacial lakes changes with pressure and temperature (Fig. 2), while it is typically constant in RB studies. α is negative at low pressures and low temperatures, such that a layer of stable density stratification exists at the top of subglacial lakes beneath a thin ice cover (H < H*). Here, we have used state-of-the-art scaling laws of RB convection and took into account the variability of the thermal expansion coefficient to make predictions about the thickness of the stable layer near the ice ceiling (δ), the anomalous temperature of the well-mixed turbulent bulk (Tbulk), the characteristic velocity of plumes (U), the characteristic velocity of the large-scale circulation (Ulsc), and the characteristic distance between plumes (𝓁). For completeness, we have also used state-of-the-art scaling laws of horizontal convection to make predictions about the typical horizontal velocity (Vhc) of the baroclinic flow that develops along a tilted ice-water interface.

The predictions for the different hydrodynamic variables are shown in Figs. 4 and 5. The key result of Fig. 4 is that subglacial lakes with a thin ice cover have an upper conductive layer several meters thick and a warm turbulent bulk (up to 1 K above freezing), whereas subglacial lakes with a thick ice cover have a thin conductive layer (centimeter scale) and a cold turbulent bulk beneath (0.01 K or less above freezing). The key result of Fig. 5 is that subglacial lakes deeper than 100 m experience substantial flow velocities, specifically U ≈ 4 mm/s and Ulsc ≈ 1 mm/s for a 1-km deep lake. These vertical velocities are larger than the horizontal velocity associated with the baroclinic flow along a tilted ice-water interface, even if the ice slope is as large as s = 10−2. The ratio Ulsc/Vhc is not much larger than 1 for steep slopes. However, if we assume that the vertical velocity of the baroclinic flow scales like ∼Vhch/L with Lh the horizontal length of the lake, then UlscL/(Vhch) ≫ 1. Thus, geothermal heating is a key factor—if not the dominant one—controlling hydrological conditions in Antarctic subglacial lakes. We provide predictions of flow properties for five well-studied subglacial lakes in Table 1.

Table 1 Properties and expected characteristics of five Antarctic subglacial lakes.

The last column is the predicted maximum diameter of particulates maintained in suspension in the mixed bulk by the large-scale circulation of vertical convection (see Discussion section). Geophysical characteristics are obtained from (2, 9, 10, 16, 57, 58), while flow conditions are derived from scaling laws discussed in the Results section of the main text and described in detail in the sections, “Scaling laws for nonrotating vertical convection” and “Scaling laws for rotating horizontal convection,” in the Materials and Methods. Ice drop refers to the difference in ice thickness above the lake due to the mean slope of the ice-water interface.

View this table:

Our analysis assumes that vertical convection and horizontal convection are decoupled. This limitation comes from the fact that numerical simulations and laboratory experiments tackling both dynamics simultaneously (either in a realistic or idealized setting) are lacking. A handful of coarse numerical models have provided some insights into the large-scale circulation of select subglacial lakes (9, 10, 14), but they rely on approximations (including parameterized turbulent processes) and are too expensive to run to allow the derivation of scaling laws of combined vertical and horizontal convection. There has been only one attempt so far at a laboratory analog of subglacial lake dynamics (29) in which the combined convective dynamics, dominated by rotation and taking the form of columnar vortex structures, was observed. The possibility to have vortices extending throughout the entire water column at full scale is an open question, which would be worth exploring.

Future work should also consider investigating the importance of the coupling dynamics between the lake circulation and melting/freezing processes along the ice-water interface. For a flat ice ceiling, we may expect that melting and freezing patterns emerge where vertical convection drives upwelling and downwelling, respectively. Such patterns would be separated by a distance equal to the lake water depth, which is the characteristic length scale of the large-scale circulation of vertical convection. For a tilted ice roof, state-of-the-art numerical models of subglacial lake dynamics typically predict that melting occurs where the ice is thickest (9, 10). However, vertical convection is not well represented in these models such that uncertainties are large regarding melting patterns and the back reaction of melting and freezing processes on the underlying lake dynamics. For instance, melting induced by the baroclinic flow may intensify local vertical convection if the melt water is dense. The possibility that topographical features emerge because of variable melt rates along the ice-water interface and influence the long-term flow dynamics is another interesting point that has yet to be addressed.

Melting and freezing occur as a result of heterogeneous heat fluxes along the ice-water interface driven by the lake water circulation. Melting of the ice ceiling into the lake releases oxygen and minerals trapped in dust particulates, and sediments can be incorporated into the lake from upstream (30). An important question is: What happens to particulates released from the lake roof and how are they dispersed by the lake circulation? Particulates in subglacial lakes can most likely be considered as passive tracers because (i) their characteristic spherical radii, which are in the range 1 μm < r < 100 μm, are much smaller than the Kolmogorov length scale, which is η = h/Re3/4 ≈ 1 cm (20) for a typical water depth of h = 100 m; (ii) their density ρs is larger but of the same order as the density of water, i.e., ρs ≈ 3ρ0 with ρ0 ≈ 999 kg/m3 the mean density of water; and (iii) particulates’ loading is expected to be dilute (13). In a quiescent fluid, sub-Kolmogorov particulates settle by gravity with speed W = 2gr2s − ρ0)/9η, assuming a linear Stokes drag, with η = 0.0017 m2 s−1 the dynamic viscosity of water. In a convecting fluid, particulates can either settle with a similar velocity or stay suspended provided that the large-scale circulation is upward and has velocity Ulsc > W, where, here, Ulsc is estimated from the GL theory. We report in the last column of Table 1 twice the maximum radius rmax of particulates maintained in suspension by the large-scale flow, i.e., such that W(rmax) = Ulsc. We find 2rmax > 7.8 μm in all cases (2rmax > 20 μm for all lakes but SPL), which means that a broad range of particulates as observed in Vostok’s accreted ice, and qualifying as “fine silt,” may be suspended in all five lakes. For Lake Vostok, it may be noted that we predict a maximum radius (36 μm) larger than that (23 μm) reported by (13) (and for nonspherical particles, such as micas, the longest axis may be even larger). This difference arises because we calculate larger flow velocities in the lake. The flow velocities and suspended particulates, which we predict for Lake Vostok, would certainly be observable by direct measurements.

In addition to the large-scale circulation, subglacial lakes experience fast and turbulent motions that can lift sediments from the bed and oppose particulates’ settling by dispersing them. The mean vertical distribution of small particulates with low inertia can be approximated by an advection-diffusion equation. The steady-state distribution in such a model is an upward-decaying exponential for the concentration of particulates n(z) ∼ ezW/D in the bulk, which shows that increased turbulence increases the particulates’ concentration by raising the background effective diffusivity, which we denote by D, in the water column. The background effective diffusivity in the open ocean is well documented and typically ranges from D = 10−5 m2 s−1 to D = 10−4 m2 s−1 (31). For a particulate with radius 4 μm and settling velocity W = 0.04 mm s−1, the corresponding e-folding decay length scale ranges from 25 cm to 2.5 m. The effective diffusivity in subglacial lakes is unknown but may be estimated from our predictions for the characteristic velocity U and length scale 𝓁 of plumes as DU𝓁. For most subglacial lakes, we predict 0.1 mm/s < U < 1 cm/s and 𝓁 ∼ 1 m (Fig. 5), such that DU𝓁 ∼ 10−4 − 10−2 m2 s−1 and the e-folding decay length scale goes from 2.5 to 250 m. Whether an effective diffusivity based on the velocity and length scale of plumes is more applicable than an effective diffusivity typical of the open ocean is an open question. The effective diffusivity based on the properties of plumes is most likely an upper bound, since plumes are intermittent. Thus, we might expect that the mean concentration of particulates, i.e., uniform in space and time, is controlled by a weak diffusivity (∼10−4) and decays by at least one order of magnitude every 10 m. This means that future explorations limited to sampling in the bulk of the lake would have to rely on intermittent plumes and local upwelling of the large-scale circulation to bring particulates upward. The mean number N of particulates in the water column, or turbidity, is key to fully assessing the habitability of subglacial lakes, in addition to the concentration of oxygen molecules derived from the ice above (11). Our calculations demonstrate that mixing of subglacial lake water is highly likely and would encourage dispersion of oxygen-rich water throughout the water column and down to the lake floor sediments, where microbial life is likely to be most abundant. A comparison of predictions for N based on advection-diffusion models as well as inferred from particulates’ concentration in basal and accreted ice, as already done for Lake Vostok (32, 33), will be key to assessing the robustness of the hydrological conditions predicted in this paper and of future particulate distribution models.

This paper provides predictions for flow velocities (0.01 to 1 cm/s), turbulent length scales (1 m), top stable layer thickness (0.01 to 10 m), temperature fluctuations (0.001 to 1 K), and the radius of particulates suspended in the water column (1 to 40 μm) due to vertical convection in Antarctic subglacial lakes. Those predictions will be verifiable by future explorations sampling lake waters and sediments using, e.g., conductivity, temperature, and depth (CTD) profilers, such as envisioned for Lake CECs and as was initially planned for Lake Ellsworth (34). To date, planning for the exploration of Lake Vostok has hinged on the analysis of accreted ice from the lake’s water in ice cores (32, 33). Our work shows that such an approach might prove inappropriate for lakes with ice covers thinner than H* = 3166 m, such as Lake CECs, since, in this case, a thick meterscale stable layer at the top of the water column prevents the upwelling of deep water and its freezing at the ice-water interface. It also means that sampling from Lake CECs should not take place at and close to the ice-water interface; instead, we predict essential measurements are required at least 1 m below the surface of the lake and likely along the entire water column.

We remark that having a stable density stratification at the top of the water column does not imply a completely quiescent environment adjacent to the ice ceiling (even if flat). Internal gravity waves (35) generated by penetrative convection (36) can propagate within the stable layer and affect particulate settling (37). How much energy is transferred from convective motions to internal waves depends on the ratio of the buoyancy frequency of the stable layer, fS, to the convective frequency, fc. For a subglacial lake, such as Lake CECs, we have fS=gρ01dρ/dz1 min−1 and fcUlsc/h ∼ 0.1 day−1. Thus, fSfc, such that convection is unlikely to penetrate far into the stratified layer and internal wave generation is weak (although this prediction neglects the possible influence of horizontal convection) (36). Nevertheless, it would be worth investigating whether internal waves in subglacial lakes can promote melting or freezing at the ice-water interface. Last, an analysis similar to the one developed in this work could be implemented for predicting dynamic conditions in icy moons in the solar system, where deep subsurface oceans exist and have attracted attention as potential habitats for extraterrestrial life (38).


Equation of state

We use the TEOS-10 toolbox in MATLAB (19) to estimate values for (i) the density of water, ρe(T, p, S), as a function of in situ temperature T, water pressure p, and absolute salinity S; (ii) the freezing temperature, Tfe(pi,S), as a function of pi (the pressure at the ice-water interface) and S; and (iii) the temperature of maximum density, Tde(p,S), as a function of p and S. Here, we restrict our attention to freshwater as opposed to seawater, i.e., we set S = 0, such that variables do not depend on S. We use superscript e to denote exact quantities, and we call the water pressure p the pressure for short, which is the absolute pressure minus atmospheric pressure pa = 10.1325 dbar. Note that p is the full water pressure, which includes pressure contributions from the ice cover, such that p > pi, with pi the ice overburden pressure. The ice pressure is related to the ice thickness through pi = ρigH/104, with pi in decibars and H in meters, and we assume a mean ice density ρi = 917 kg/m3. Unless stated otherwise, all variables use SI units except temperature variables, which are in degrees Celsius (°C) and pressure variables, which are in decibars (dbar), since °C and dbar are standard units in physical oceanography.

For simplicity, we derive explicit, approximate expressions for ρ, Tf, and Td from the TEOS-10 exact values. The freezing temperature and the temperature of maximum density can be well approximated by quadratic polynomials in pi and p, respectively. For 0 < p, pi < 104 dbar, we find that the best-fit polynomials (with p, pi in dbar)Tf=4.7184×1037.4584×104pi1.4999×108pi2(2)Td=3.97952.0059×103p6.2511×108p2(3)approximate Tfe and Tde to within 0.002 K. We approximate the density of water by a quadratic bivariate polynomial, which is maximum at T = Td(p). For 0 < p < 104 dbar and Tf < T < Tf + 15 K, we find that the best-fit bivariate polynomial (with p in dbar)ρ=ρ0+ρ1(p)+C(p)[TTd(p)]2(4)withρ0=9.9999×102,ρ1=4.9195×103p1.4372×108p2C=7.0785×103+1.8217×107p+4.2679×1012p2(5)approximates ρe to within less than 0.01% relative error. This implies density errors less than 0.1 kg/m3, which is an order of magnitude less than density variations expected with temperature alone. From Eq. 4, we can derive the approximate thermal expansion coefficientα=1ρ0ρTp=2C(p)[TTd(p)]ρ0(6)which is shown in Fig. 2 and changes sign at T = Td > Tf for pressures lower than p* = 2848 dbar.

We note that the nonmonotonic, anomalous behavior of water density at low pressure is well known for freshwater lakes at atmospheric pressure (39) and disappears progressively with increasing salt concentration. With typical salinities of S ≈ 35 g/kg, the density of Earth’s oceans decreases monotonically with increasing temperatures.

Evolution equations

The evolution equations for subglacial lakes are the Navier-Stokes equations in the Boussinesq approximation. Here, we include compressibility effects in the energy equation because we are interested in the calculation of the (small) critical heat flux at the onset of convection, but compressibility effects can otherwise be neglected when considering the (large) geothermal heat flux. In a Cartesian coordinates system (x, y, and z) centred on a lake’s top boundary, the equations for the velocity vector u, pressure p, density ρ, and temperature T (in °C) read (4042)ρ0DuDt+ρ0fez×u=p+μ2uρgez(7)·u=0(8)ρ0cpDTDtα(T+T0)DpDt=k2T(9)where f is the Coriolis frequency, μ is the dynamic viscosity, k is the thermal conductivity, cp is the isobaric specific heat capacity, g is the gravitational acceleration, and T0 = 273.15 K; D/Dt ≡ ∂t + u · ∇ denotes material derivative, with t the time derivative and ∇ is the gradient operator. Equation 7 is the momentum equation in the Boussinesq approximation; Eq. 8 is mass conservation for an incompressible fluid; and Eq. 9 is the energy equation including pressure effects, which are relevant in the limit of small temperature variations. Note that p is in pascal (Pa) in the above equations but is converted to dbar by dividing by 104 when used in Eqs. 3 to 6.

We consider the Coriolis frequency at 80°S, i.e., we take f = − 1.4363 × 10−4 rad/s; we use μ = 1.7 × 10−3 kg m−1 s−1 and k = 0.56 W m−1 K−1, which are the dynamic viscosity and thermal conductivity values at reference pressure p = 0 dbar and temperature T = 0.01°C; we use cp = 4.2174 × 103 J kg−1 K−1; and we recall that g = 9.81 m/s2 at and near Earth’s surface. Note that μ and k may be expected to vary with p and T. However, to the best of our knowledge, only few studies have investigated their dependence, in particular in the cold-temperature and high-pressure regimes, and reported little variations for the pressure and temperature conditions of our interest such that we take them constants (43, 44). We denote by ν = μ/ρ0 = 1.7 × 10−6 m2/s the constant kinematic viscosity and by κ = k0cp = 1.3 × 10−7 m2/s the constant thermal diffusivity.

Subglacial lake water must be at the freezing temperature at the upper lake boundary, i.e., T = Tf(pi) at z = 0 m, while at the base of the lake, it is the heat flux that is enforced, i.e., k∂zT = − F at z = − h, with F > 0 the (geothermal) heat flux and h > 0 the lake depth; also, p = pi at z = 0. Equations 7 to 9, along with the equation of state (Eq. 4), have the stationary base-state solution (denoted by overbars)u¯=0, T¯=Tf(pi)zFk, dzp¯=ρ(T¯,p¯)g(10)i.e., the temperature increases linearly with depth, and the pressure is hydrostatic. For simplicity, we assume ρ ≈ ρ0 in the hydrostatic base-state equation, such that p¯=piρ0gz at leading order (assuming a pressure variable expressed in Pascals).

Static stability of an ideal compressible fluid

The criterion for an ideal (dissipationless) compressible fluid with hydrostatic base-state pressure p = pi − ρ0gz (overbar dropped) to be locally (statically) stable is (39, 45)1ρ0ρTpdsdz=α(cpT+T0dTdzαρ0dpdz)=αcpT+T0(dTdzdTaddz)<0(11)where s is the entropy, T0 = 273.15 K (we recall that we express temperature variables in °C), and we approximate ρ ≈ ρ0 at leading order. dTad/dz is known as the adiabatic temperature gradient and readsdTaddz=α(T+T0)gcp(12)such that dTad/dz > 0 if α < 0. Here, heating is provided from the bottom of the lake such that we always have dT/dz < 0. As a result, when α < 0, equation 11 is always satisfied and the lake is stable. For α > 0, equation 11 is not satisfied, and the lake is unstable to vertical convection if dT/dz < dTad/dz < 0, i.e., if the heat flux exceeds (in absolute value) the adiabatic heat flux. Thus, the two conditions for subglacial lakes experiencing geothermal heating (i.e., such that dT/dz < 0) to be locally unstable areα>0(13)dTdz<dTaddz(14)

Note that both α and dTad/dz are functions of z when considering the base state of a subglacial lake heated from below. Specifically, α increases with depth, while dTad/dz decreases with depth (note that it is negative and so increases in absolute value), such that equation 13 (resp. Eq. 14) is more readily satisfied at the bottom (resp. top) of the lake. As a result, it is possible to find cases where a subglacial lake is globally unstable but remains statically stable in some places. When this happens, convection is expected to occur in subregions of the water column where equations 13 and 14 are satisfied.

Equations 13 and 14 are necessary but not sufficient conditions for flow instability. The heat flux must sustain a temperature gradient with a buoyancy anomaly that is also large enough to overcome viscosity and diffusivity effects. The calculation of the exact, minimum critical heat flux leading to convection in subglacial lakes, with dissipation effects taken into account, is the result of the stability analysis described in the next section.

Linear stability analysis

We study the stability of the base-state solution (Eq. 10) by investigating how small initial perturbations evolve over time. We expand the variables (generically represented by X) asX=X¯+X(15)with primes denoting the perturbed variables. Substituting expanded variables in Eqs. 7 to 9, using Eq. 4, and linearizing, we obtain the perturbation equationsρ0ut+ρ0fez×u=p+μ2u+ρ0α¯Tgezρ¯ppgez(16)·u=0(17)ρ0cp(tTwF/k)α¯(T¯+T0)(tpwρ0g)=k2T(18)where ρ¯p is related to the small compressibility of the background state and is derived from Eq. 4 asρ¯p=[ρ11+2ρ12p¯+(T¯T¯d)2(C1+2C2p¯)+2C¯(T¯T¯d)(Td12Td2p¯)](19)with subscripts 1,2 denoting the linear and leading coefficients of the polynomial expressions for ρ1, Td, and C, e.g., Td1 = − 2.0059 × 10−3 K/dbar (Eqs. 3 and 5).

Since the perturbation equations do not depend explicitly on x, y, and t, the stability criterion can be inferred from the temporal evolution of plane waves of the formX(x,y,z,t)=X̂(z)eσt+i(kxx+kyy)+c.c.,(20)with σ as the growth rate, kx and ky as the wave numbers in the x and y directions, and c.c. as the complex conjugate. Assuming horizontal isotropy, i.e., kx = ky, and substituting variables of the form given by Eq. 20 in Eqs. 16 to 18, we derive a one-dimensional linear eigenvalue problem for the growth rate σ, which we solve numerically with the open-source pseudospectral Dedalus code and the Eigentools package (46). We expand variables in the z direction using Chebyshev modes and compute the largest growth rate σ(k, F) as a function of wave number k=kx2+ky2 and heat flux F for a range of input parameters (pi, h). The critical minimum heat flux Fc that destabilizes the base state is the minimum of F for which there exists a wave number k such that σ(k, F) > 0. We report Fc in Fig. 3A. Note that the eigenvalue problem becomes challenging at large h, such that we limit the calculations to h ≤ 20 m. We extrapolate to larger h using scaling laws that are either asymptotically valid or conservative, i.e., such that they may overestimate Fc. We describe the extrapolation procedure in detail in the Supplementary Materials.

Scaling laws for nonrotating vertical convection

We show in the Supplementary Materials that vertical convection in Antarctic subglacial lakes is better represented by nonrotating convection than geostrophic convection. As a result, we use scaling laws obtained in the idealized limit of nonrotating turbulent convection to make predictions about δ, Tbulk, U, and 𝓁 for subglacial lakes subject to a geothermal flux F = 50 mW/m2.

First, we remark that the conductive layer at the top of the lake includes the turbulent boundary layer and a stable layer where α < 0 for pi < p*. In other words, we write δ = δt + δS with δt as the thickness of the turbulent boundary layer and δS as the thickness of the stable layer. We obtain δS as the positive solution of the equation Td(pi + ρ0gδS/104) = Tf(pi) + δSF/k, i.e., δS is found as the location z = − δS where the temperature of maximum density equals the temperature of the conductive base-state profile. The temperature at the base of the stable layer is correspondingly TS = Tf + δSF/k. The definition of the control parameter RaF (Eq. 1) of subglacial lakes includes the effective water depth heff and the characteristic thermal expansion coefficient αeff. The effective water depth is simply the region of the water column where convection occurs (α > 0), such that we assume heff = h − δS. Estimating αeff inside a convective lake is difficult since the vertical profiles of temperature and α are unknown. Here, for simplicity, we use αeff = α(pi + ρ0gh/104, TS), i.e., we take α on the bottom boundary as the effective thermal expansion coefficient but assume that the temperature of the lake does not exceed TS. This assumption is likely to underestimate αeff but is the best conservative assumption possible without prior knowledge of the vertical temperature profile.

We use scaling laws for the Nusselt number Nu and the Reynolds number Re as functions of Ra inferred from recent state-of-the-art numerical simulations (25, 47) to predict δt, Tbulk, U, and 𝓁. We relate our flux-based Rayleigh number RaF to Ra following previous works (48, 49), i.e., such thatRaF=RaNu(21)

The scaling laws for Nu in (47) and Re in (25) areNu=0.16Ra2/7(22)Re=0.18(RaRaNu)1/2Pr1(23)

Combining Eq. 1 with Eqs.21 and 22, we can estimate Ra = (RaF/0.16)7/9 and predict δt = 0.5heff/Nu = 0.5heff/(0.16Ra2/7), such thatδ=0.5heff0.16Ra2/7+δS(24)

Assuming a conductive temperature profile in the stable and turbulent boundary layers, we then findTbulk=Fδk(25)

The turbulent flow velocity is inferred from Re = Uheff/ν and Eqs. 22 and 23 asU=ν0.18Ra6.25Ra5/7Prheff(26)

The estimate of the characteristic turbulent length scale, or distance between plumes, is finally obtained from equation (6.3) of reference (25) as=0.8RePr(3.125Ra2/7)3/2heff(27)

We provide a second set of predictions for δt, Tbulk, and Ulsc based on scaling laws inferred from the GL unifying theory of RB convection (26). We use subscript lsc for the velocity as the GL theory applies to the velocity of the large-scale circulation rather than the velocity of the plumes, as is the case in (25). The procedure for deriving δt, Tbulk, and Ulsc from the GL theory is the same as above, i.e., we combine Eqs. 1 and 21 with scaling laws for Nu and Re to derive δt = 0.5heff/Nu, Tbulk = Fδ/kS is unchanged) and Ulsc = νRe/heff. The scaling laws for Nu and Re combine the expressions derived in subregions Iu, IIIu, and IVu of the GL theory [see (26) and the Supplementary Materials].

Scaling laws for rotating horizontal convection

The ice-water interface of subglacial lakes is often sloped such that a baroclinic horizontal convection flow develops along the ice-water interface. Here, we provide approximate estimates of the horizontal velocity of the baroclinic flow to compare the dynamical importance of vertical convection to horizontal convection. Besides the Prandtl number Pr, the control parameters for horizontal convection are the Ekman number EkL and the Rayleigh number RaL based on the lake’s horizontal length L (27, 50), i.e.EkL=νfL2,RaL=gαiΔiL3νκ(28)where αi is the effective thermal expansion coefficient near the ice ceiling, and Δi is the temperature difference along the tilted ice-water interface due to the changing freezing temperature Tf(pi) with the ice pressure (Eq. 2). The ice pressure drop along the ice-water interface is δpi = sLρig/104 dbar, with s the slope and L in meters, such thatΔi=Tf(pi)Tf(pi+δpi)>0(29)

For αi, we take the maximum of α along the ice-water interface, i.e.αi=max[α(pi,Tf(pi)),α(pi+δpi,Tf(pi+δpi))](30)

We show in the Supplementary Materials that horizontal convection in Antarctic subglacial lakes is constrained by rotation. Thus, here we use a recent scaling law for the Reynolds number Rehc of rotation-constrained horizontal convection (27) to estimate the characteristic velocity of the large-scale horizontal flow Vhc. For simplicity, we assume a fixed aspect ratio of L/h = 250, i.e., the horizontal length is 250 times the water depth, even though subglacial lakes have variable aspect ratio. We note that a ratio of 250 is on the higher end of observed aspect ratios (see Table 1), such that our estimates of the lakes’ lengths may be closer to an upper than a lower bound. The scaling law for rotating horizontal convection inferred from (27) isRehc=(RaLEkL)2/3Pr(31)such thatVhc=νRehcL(32)

Note that VhcL1/3, such that estimates for the horizontal velocity would be only weakly affected (weakly decreasing) by decreasing the aspect ratio.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: We gratefully acknowledge fruitful discussions with K. Makinson at the British Antarctic Survey and C. Vreugdenhil at the University of Cambridge. L.-A.C. thanks T. Alboussière from ENS Lyon for the helpful discussions on the thermodynamics and stability of compressible fluids. Funding: This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement 793450. We acknowledge PRACE for awarding us access to Marconi at CINECA, Italy. Author contributions: All authors conceptualized and supervised the work. L.-A.C. performed the calculations and led the writing of the original draft and revision. All authors discussed the results and reviewed and edited the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. All data, code, and materials will be made available on a per request basis. The simulation code Dedalus used in this work is open source.

Stay Connected to Science Advances

Navigate This Article