Dynamic flows create potentially habitable conditions in Antarctic subglacial lakes
CCouston and Siegert,
Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
O C E A N O G R A P H Y
Dynamic flows create potentially habitable conditions in Antarctic subglacial lakes
Louis-Alexandre Couston * and Martin Siegert 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.
INTRODUCTION
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 ( ). 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 ( ). They comprise a wide variety of sizes and glaciological and topographic settings ( ) and have been hypothe-sized as potential habitats for the in situ development of microbial organisms ( ). 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 ( ). A fur-ther ~130 hydrologically active lakes, which experience rapid water discharges and large volume changes, exist toward the margin of the ice sheet ( , ). While these may contain microbial life ( ), they are not considered as isolated habitats where microbes can adapt inde-pendently 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 min-erals 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 ( , , ). 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 ( ), dynamic flows and mixing of bottom water within the water column are essential for life to be widespread and detectable, avoid-ing, 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/m ; ( )], and hori-zontal 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 sedi-ments in the water column ( ). This is of the same order of magni-tude as that predicted by ocean modeling for a handful of subglacial lakes, including Lake Vostok ( ), Lake Ellsworth ( ), and Lake Concordia ( ). 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 ( ). 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 micro-bial 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 ( ), and, more specifically, previous efforts that aimed to predict the hydrological conditions of individual sub-glacial lakes, including Lake Vostok ( ). 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 atten-tion to freshwater, because salt concentration is typically low in iso-lated subglacial lakes ( ).We first calculate the minimum critical heat flux, F c , 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 ( ). We show that F c is much smaller than 50 mW/m , 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 British Antarctic Survey, Cambridge CB3 0ET, UK. Department of Applied Mathe-matics and Theoretical Physics, University of Cambridge, Cambridge CB3 0WA, UK. Univ Lyon, ENS de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique, F-69342 Lyon, France. Grantham Institute and Department of Earth Science and Engineering, Imperial College London, London, SW7 2AZ, UK.*Corresponding author. Email: [email protected]
Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY). on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC BY). on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E geothermal flux of 50 mW/m experience dynamic flows by applying state-of-the-art scaling laws of classical thermal convection to con-vection in cold-temperature high-pressure lake environments. RESULTS
Subglacial lakes can experience dynamic vertical convection flows [also known as Rayleigh-Bénard (RB) convection ( )], 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 p i ). Convection in subglacial lakes is complex, because the thermal expansion coefficient of fresh-water ( p , T ), which indicates how fluid parcels contract or expand with changes in temperature (i.e., are buoyant), depends on water pressure p > p i and the temperature T itself ( ). For relatively thick ice cover, i.e., H ≥ H * = 3166 m, the thermal expansion coefficient is always positive (i.e., the density decreases with temperature) and in-creases 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 p i < 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 non-monotonic function of T . Specifically, as is shown in Fig. 2, for H < H * , increases with temperature but is first negative for T f ( p i ) ≤ T < T d ( p ) (close to the ice ceiling), where T f is the freezing temperature and T d 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 tem-perature exceeds T d , 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 buoy-ancy force that can overcome viscous dissipation and thermal diffusion.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 fresh-water using the Thermodynamic Equation of Seawater 2010 (TEOS-10) toolbox ( ) and the Coriolis frequency at 80°S. We take into account the adiabatic temperature gradient by including compressibility ef-fects 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 ( ).Figure 3A shows the minimum critical heat flux F c that permits vertical convection for a wide range of ice pressures and water depths relevant to Antarctic subglacial lakes. F c is large at small pressure and small water depth (top left corner) but then decreases with h and p i in most of the parameter space. We find F c < 50 mW/m (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 parame-ter space. We predict that Lake CECs and South Pole Lake (SPL) are unstable to vertical convection if subject to a 50 mW/m flux, i.e., their critical heat flux F c is less than 50 mW/m , despite having rel-atively 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 p i − p * and a symmetric logarithmic scale. As a result, the transition from a fully convective water column (for p i > p * ) to a convective water column with a stably stratified upper layer (for p i ≤ p * ) is smooth even though F c increases rapidly as p i decreases below p * . Figure 3B shows that subglacial lakes reported in the last inventory ( ) have ice cover thicknesses almost equally distributed on either side of H * . Thus, the p i = 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 F c but also in practice. Almost half of the subglacial lakes (with p i < 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 F c is constrained pri-marily by the condition of having > 0 on the bottom boundary Fig. 1. Problem schematic.
We provide predictions about the characteristic velocity of the large-scale circulation U lsc , the characteristic velocity of turbulent plumes U , the thickness of the top stable conductive layer, and the anomalous temperature of the well-mixed bulk T bulk (i.e., in excess of the freezing temperature T f ). The problem parameters are the water depth h , the ice thickness H (or ice overburden pressure p i ), the Coriolis frequency f (due to Earth’s rotation), and the geothermal flux F . 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 T d (red solid line) and freezing temperature T f (black solid line) with pressure. For small pressures p < p *, with p * the critical inversion pressure (blue dashed line), T d > T f such that there exists a range of temperatures T f < T < T d for which is nega-tive (area appearing with red colors) and water masses become anomalously denser with increasing temperatures. For p < p * and T > T d , or p ≥ p *, the water becomes mono-tonically lighter as temperature increases, which is the typical behavior of most fluids. on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
Plot of the thermal expansion coefficient as a function of ( T , p ) superimposed with profiles of the temperature of maximum density T d (red solid line) and freezing temperature T f (black solid line) with pressure. For small pressures p < p *, with p * the critical inversion pressure (blue dashed line), T d > T f such that there exists a range of temperatures T f < T < T d for which is nega-tive (area appearing with red colors) and water masses become anomalously denser with increasing temperatures. For p < p * and T > T d , or p ≥ p *, the water becomes mono-tonically lighter as temperature increases, which is the typical behavior of most fluids. on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E for small p i and h (top left corner) and by the condition that it must exceeds the adiabatic flux for large p i and h (bottom right corner). In between, viscous dissipation and thermal diffusion dominate the calculation of F c . Note that our prediction of F c for small p i and large h is conservative and may overestimate the true F c . We provide fur-ther details about the calculation of F c in Materials and Methods and in the Supplementary Materials.For a geothermal flux F greater than the critical heat flux F c , 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 tur-bulent 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 ( ) 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 ( ), astrophysics, and engineering ( ). 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 as Ra F = g eff Fh eff ─ k , Pr = ─ , Ek = ─ ∣ f ∣ h eff (1) where g is the gravitational acceleration, eff is the characteristic thermal expansion coefficient, h eff is the effective water depth where convection occurs, is the kinematic viscosity, is the thermal dif-fusivity, 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 Ra F 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 pre-scribed temperature difference, which is more common in idealized studies of natural convection ( ). Note that we neglect compress-ibility effects for the prediction of variables in the turbulent regime because Earth’s geothermal flux is several orders of magnitude larg-er than the adiabatic heat flux.In the context of subglacial lakes, the geothermal flux is suffi-ciently large that the lake water is in a fully turbulent state that is almost not affected by rotation, i.e., F ≫ F c , and the effect of rota-tion is weak [see the Supplementary Materials and ( )]. Thus, here we use scaling laws derived for fully turbulent nonrotating convec-tion to make predictions about hydrodynamic variables. The vari-ables of interest are the thickness of the conductive layer near the ice ceiling , the anomalous temperature of the well-mixed bulk T bulk (in excess of T f ), 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/m throughout and use scaling laws derived from nu-merical simulations ( ) as well as scaling laws inferred from the Grossmann-Lohse (GL) unifying theory of RB convection, which is based on theoretical arguments ( ).We first show in Fig. 4 (A to D) the results for and T bulk based on the scaling laws derived in ( ) for a wide range of ice thicknesses and water depths. The thickness of the conductive layer near the ice 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/m 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 ( ). The dashed gray lines highlight the critical thickness H *. on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
Plot of the thermal expansion coefficient as a function of ( T , p ) superimposed with profiles of the temperature of maximum density T d (red solid line) and freezing temperature T f (black solid line) with pressure. For small pressures p < p *, with p * the critical inversion pressure (blue dashed line), T d > T f such that there exists a range of temperatures T f < T < T d for which is nega-tive (area appearing with red colors) and water masses become anomalously denser with increasing temperatures. For p < p * and T > T d , or p ≥ p *, the water becomes mono-tonically lighter as temperature increases, which is the typical behavior of most fluids. on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E for small p i and h (top left corner) and by the condition that it must exceeds the adiabatic flux for large p i and h (bottom right corner). In between, viscous dissipation and thermal diffusion dominate the calculation of F c . Note that our prediction of F c for small p i and large h is conservative and may overestimate the true F c . We provide fur-ther details about the calculation of F c in Materials and Methods and in the Supplementary Materials.For a geothermal flux F greater than the critical heat flux F c , 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 tur-bulent 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 ( ) 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 ( ), astrophysics, and engineering ( ). 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 as Ra F = g eff Fh eff ─ k , Pr = ─ , Ek = ─ ∣ f ∣ h eff (1) where g is the gravitational acceleration, eff is the characteristic thermal expansion coefficient, h eff is the effective water depth where convection occurs, is the kinematic viscosity, is the thermal dif-fusivity, 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 Ra F 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 pre-scribed temperature difference, which is more common in idealized studies of natural convection ( ). Note that we neglect compress-ibility effects for the prediction of variables in the turbulent regime because Earth’s geothermal flux is several orders of magnitude larg-er than the adiabatic heat flux.In the context of subglacial lakes, the geothermal flux is suffi-ciently large that the lake water is in a fully turbulent state that is almost not affected by rotation, i.e., F ≫ F c , and the effect of rota-tion is weak [see the Supplementary Materials and ( )]. Thus, here we use scaling laws derived for fully turbulent nonrotating convec-tion to make predictions about hydrodynamic variables. The vari-ables of interest are the thickness of the conductive layer near the ice ceiling , the anomalous temperature of the well-mixed bulk T bulk (in excess of T f ), 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/m throughout and use scaling laws derived from nu-merical simulations ( ) as well as scaling laws inferred from the Grossmann-Lohse (GL) unifying theory of RB convection, which is based on theoretical arguments ( ).We first show in Fig. 4 (A to D) the results for and T bulk based on the scaling laws derived in ( ) for a wide range of ice thicknesses and water depths. The thickness of the conductive layer near the ice 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/m 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 ( ). The dashed gray lines highlight the critical thickness H *. on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E ceiling is almost independent of lake depth but varies substantially with ice pressure (Fig. 4, A and B). For thin ice cover, the top con-ductive 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 centime-ters 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 conduc-tion only. Thus, the temperature increases linearly from T f ( p i ) near the ice to T f ( p i ) + 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 T bulk = 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 T bulk based on scaling laws derived in ( ) (labeled “GL”). There is a good agree-ment between predictions based on the scaling laws in ( ) (shown by circles) and ( ).Figure 5 shows the predicted lake velocity U and horizontal length scale 𝓁 based on previously derived scaling laws ( ). 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 p i = p * because the effective thermal expan-sion coefficient eff , which we estimate conservatively (cf. Materials and Methods) and use in Eq. 1 for Ra F , 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 p i = p * + 100 dbar and p i = p * for h = 10 m. We expect that the discontinuity would become less sharp but would not completely disappear upon relaxing our 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/m . (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 T bulk (above T f ) of the well-mixed convective layer. on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
Plot of the thermal expansion coefficient as a function of ( T , p ) superimposed with profiles of the temperature of maximum density T d (red solid line) and freezing temperature T f (black solid line) with pressure. For small pressures p < p *, with p * the critical inversion pressure (blue dashed line), T d > T f such that there exists a range of temperatures T f < T < T d for which is nega-tive (area appearing with red colors) and water masses become anomalously denser with increasing temperatures. For p < p * and T > T d , or p ≥ p *, the water becomes mono-tonically lighter as temperature increases, which is the typical behavior of most fluids. on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E for small p i and h (top left corner) and by the condition that it must exceeds the adiabatic flux for large p i and h (bottom right corner). In between, viscous dissipation and thermal diffusion dominate the calculation of F c . Note that our prediction of F c for small p i and large h is conservative and may overestimate the true F c . We provide fur-ther details about the calculation of F c in Materials and Methods and in the Supplementary Materials.For a geothermal flux F greater than the critical heat flux F c , 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 tur-bulent 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 ( ) 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 ( ), astrophysics, and engineering ( ). 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 as Ra F = g eff Fh eff ─ k , Pr = ─ , Ek = ─ ∣ f ∣ h eff (1) where g is the gravitational acceleration, eff is the characteristic thermal expansion coefficient, h eff is the effective water depth where convection occurs, is the kinematic viscosity, is the thermal dif-fusivity, 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 Ra F 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 pre-scribed temperature difference, which is more common in idealized studies of natural convection ( ). Note that we neglect compress-ibility effects for the prediction of variables in the turbulent regime because Earth’s geothermal flux is several orders of magnitude larg-er than the adiabatic heat flux.In the context of subglacial lakes, the geothermal flux is suffi-ciently large that the lake water is in a fully turbulent state that is almost not affected by rotation, i.e., F ≫ F c , and the effect of rota-tion is weak [see the Supplementary Materials and ( )]. Thus, here we use scaling laws derived for fully turbulent nonrotating convec-tion to make predictions about hydrodynamic variables. The vari-ables of interest are the thickness of the conductive layer near the ice ceiling , the anomalous temperature of the well-mixed bulk T bulk (in excess of T f ), 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/m throughout and use scaling laws derived from nu-merical simulations ( ) as well as scaling laws inferred from the Grossmann-Lohse (GL) unifying theory of RB convection, which is based on theoretical arguments ( ).We first show in Fig. 4 (A to D) the results for and T bulk based on the scaling laws derived in ( ) for a wide range of ice thicknesses and water depths. The thickness of the conductive layer near the ice 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/m 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 ( ). The dashed gray lines highlight the critical thickness H *. on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E ceiling is almost independent of lake depth but varies substantially with ice pressure (Fig. 4, A and B). For thin ice cover, the top con-ductive 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 centime-ters 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 conduc-tion only. Thus, the temperature increases linearly from T f ( p i ) near the ice to T f ( p i ) + 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 T bulk = 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 T bulk based on scaling laws derived in ( ) (labeled “GL”). There is a good agree-ment between predictions based on the scaling laws in ( ) (shown by circles) and ( ).Figure 5 shows the predicted lake velocity U and horizontal length scale 𝓁 based on previously derived scaling laws ( ). 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 p i = p * because the effective thermal expan-sion coefficient eff , which we estimate conservatively (cf. Materials and Methods) and use in Eq. 1 for Ra F , 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 p i = p * + 100 dbar and p i = p * for h = 10 m. We expect that the discontinuity would become less sharp but would not completely disappear upon relaxing our 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/m . (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 T bulk (above T f ) of the well-mixed convective layer. on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E conservative approximation for eff , because will always be (overall) much smaller in lakes with a thin ice cover ( p i ≤ p * ) than in lakes with a thick ice cover ( p i > p * ). We also show in Fig. 5B the predic-tion for the characteristic velocity U lsc based on the scaling laws de-rived in ( ) (labeled GL) for an ice thickness H = 3945 m (filled squares), which is the ice thickness above Lake Vostok. U lsc 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 ( ), which is likely to be faster than the mean large-scale flow. Figure 5B also shows a prediction for the horizontal velocity V hc of the baroclinic flow expected along a sloped ice-water interface, using scaling laws inferred from recent results on horizontal convection ( ). We show the prediction for V hc 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 baro-clinic 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 ( U lsc ). DISCUSSION
We have demonstrated that the critical heat flux leading to vertical convection in subglacial lakes is much less than 50 mW/m for a broad range of ice overburden pressures and water depths (Fig. 3). Thus, it should be considered that most—if not all—Antarctic sub-glacial lakes are dynamic hydrologic environments. We expect that the same conclusion holds for isolated subglacial lakes in Greenland and elsewhere in the solar system ( , ). We note that our prediction of the critical heat flux F c is conservative for large water depth and small ice pressure (see the Supplementary Materials). Also, estimates 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 charac-teristic 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 U lsc of vertical convection for H = 3945 m (shown by filled black squares). In (B), we also show a prediction for the horizontal velocity V hc 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). on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
We have demonstrated that the critical heat flux leading to vertical convection in subglacial lakes is much less than 50 mW/m for a broad range of ice overburden pressures and water depths (Fig. 3). Thus, it should be considered that most—if not all—Antarctic sub-glacial lakes are dynamic hydrologic environments. We expect that the same conclusion holds for isolated subglacial lakes in Greenland and elsewhere in the solar system ( , ). We note that our prediction of the critical heat flux F c is conservative for large water depth and small ice pressure (see the Supplementary Materials). Also, estimates 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 charac-teristic 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 U lsc of vertical convection for H = 3945 m (shown by filled black squares). In (B), we also show a prediction for the horizontal velocity V hc 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). on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E of F c exceeding 50 mW/m (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 mea-surements 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 predic-tions about the thickness of the stable layer near the ice ceiling ( ), the anomalous temperature of the well-mixed turbulent bulk ( T bulk ), the characteristic velocity of plumes ( U ), the characteristic velocity of the large-scale circulation ( U lsc ), 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 ( V hc ) 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 U lsc ≈ 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 U lsc / V hc is not much larger than 1 for steep slopes. However, if we assume that the vertical velocity of the baroclinic flow scales like ∼ V hc h / L with L ≫ h the horizontal length of the lake, then U lsc L /( V hc h ) ≫ 1. Thus, geothermal heating is a key factor—if not the dominant one—controlling hydrological con-ditions in Antarctic subglacial lakes. We provide predictions of flow properties for five well-studied subglacial lakes in Table 1.Our analysis assumes that vertical convection and horizontal convec-tion are decoupled. This limitation comes from the fact that numerical simulations and laboratory experiments tackling both dynamics simul-taneously (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 ( , , ), 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 ( ) 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 cir-culation 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 ( , ). However, vertical convection is not well represented in these models such that uncer-tainties 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 dy-namics 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 ( ). 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 / Re ≈ 1 cm ( ) 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 with ≈ 999 kg/m the mean density of water; and (iii) particulates’ loading is expected to be dilute ( ). In a quiescent fluid, sub- Kolmogorov particulates settle by gravity with speed W = 2 gr ( s − )/9 , 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 ( , , , , , ), 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. Ice thickness (m) Ice drop (m) Lake length (km) Water depth (m) (m) T bulk (K) ℓ (m) U (mm/s) U lsc (mm/s) V hc (mm/s) 2 r max ( m) CECs 2653 159 10.35 300 7.7 0.69 1.6 0.97 0.32 0.041 22SPL 2857 30 10 32 4.7 0.42 0.8 0.10 0.04 0.010 7.8Ellsworth 3400 300 10 156 0.077 0.0069 1.2 0.69 0.26 0.066 20Vostok 3945 600 280 1067 0.066 0.0059 2.3 3.80 0.85 0.066 36Concordia 4055 168 45 126 0.063 0.0056 1.0 0.83 0.31 0.044 22 on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
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 ( , , , , , ), 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. Ice thickness (m) Ice drop (m) Lake length (km) Water depth (m) (m) T bulk (K) ℓ (m) U (mm/s) U lsc (mm/s) V hc (mm/s) 2 r max ( m) CECs 2653 159 10.35 300 7.7 0.69 1.6 0.97 0.32 0.041 22SPL 2857 30 10 32 4.7 0.42 0.8 0.10 0.04 0.010 7.8Ellsworth 3400 300 10 156 0.077 0.0069 1.2 0.69 0.26 0.066 20Vostok 3945 600 280 1067 0.066 0.0059 2.3 3.80 0.85 0.066 36Concordia 4055 168 45 126 0.063 0.0056 1.0 0.83 0.31 0.044 22 on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E assuming a linear Stokes drag, with = 0.0017 m 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 U lsc > W , where, here, U lsc is estimated from the GL theory. We report in the last column of Table 1 twice the maximum radius r max of particulates maintained in suspension by the large-scale flow, i.e., such that W ( r max ) = U lsc . We find 2 r max > 7.8 m in all cases (2 r max > 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 ( ) (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 particu-lates, which we predict for Lake Vostok, would certainly be observ-able by direct measurements.In addition to the large-scale circulation, subglacial lakes experi-ence 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 approxi-mated by an advection-diffusion equation. The steady-state distri-bution in such a model is an upward-decaying exponential for the concentration of particulates n ( z ) ∼ e − zW / 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 m s −1 to D = 10 −4 m s −1 ( ). For a particulate with radius 4 m and set-tling 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 predic-tions for the characteristic velocity U and length scale 𝓁 of plumes as D ∼ U 𝓁 . For most subglacial lakes, we predict 0.1 mm/s < U < 1 cm/s and 𝓁 ∼ D ∼ U 𝓁 ∼ −4 − 10 −2 m 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 prop-erties 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 ( ∼ −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 particu-lates upward. The mean number N of particulates in the water column, or turbidity, is key to fully assessing the habitability of sub-glacial lakes, in addition to the concentration of oxygen molecules derived from the ice above ( ). 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 ( , ), will be key to assessing the robustness of the hydro-logical 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 particu-lates 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 ( ). 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 ( , ). 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 meter-scale 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 environ-ment adjacent to the ice ceiling (even if flat). Internal gravity waves ( ) generated by penetrative convection ( ) can propagate within the stable layer and affect particulate settling ( ). How much energy is transferred from convective motions to internal waves depends on the ratio of the buoyancy frequency of the stable layer, f S , to the convective frequency, f c . For a subglacial lake, such as Lake CECs, we have f S = √ ___________ − g ρ −1 d ρ / dz ≅ 1 min −1 and f c ∼ U lsc / h ∼ −1 . Thus, f S ≫ f c , 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 con-vection) ( ). 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 devel-oped 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 ( ). MATERIALS AND METHODS
Equation of state
We use the TEOS-10 toolbox in MATLAB ( ) to estimate values for (i) the density of water, e ( T , p , S ), as a function of in situ tem-perature T , water pressure p , and absolute salinity S ; (ii) the freezing temperature, T f e ( p i , S ) , as a function of p i (the pressure at the ice-water interface) and S ; and (iii) the temperature of maximum density, T d e ( 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 p a = 10.1325 dbar. Note that p is the full water pressure, which includes pressure con-tributions from the ice cover, such that p > p i , with p i the ice over-burden pressure. The ice pressure is related to the ice thickness through p i = i gH /10 , with p i in decibars and H in meters, and we assume a mean ice density i = 917 kg/m . Unless stated otherwise, all variables use SI units except temperature variables, which are in on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
We use the TEOS-10 toolbox in MATLAB ( ) to estimate values for (i) the density of water, e ( T , p , S ), as a function of in situ tem-perature T , water pressure p , and absolute salinity S ; (ii) the freezing temperature, T f e ( p i , S ) , as a function of p i (the pressure at the ice-water interface) and S ; and (iii) the temperature of maximum density, T d e ( 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 p a = 10.1325 dbar. Note that p is the full water pressure, which includes pressure con-tributions from the ice cover, such that p > p i , with p i the ice over-burden pressure. The ice pressure is related to the ice thickness through p i = i gH /10 , with p i in decibars and H in meters, and we assume a mean ice density i = 917 kg/m . Unless stated otherwise, all variables use SI units except temperature variables, which are in on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E 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 , T f , and T d from the TEOS-10 exact values. The freezing temperature and the temperature of maximum density can be well approximated by quadratic polynomials in p i and p , respectively. For 0 < p , p i < 10 dbar, we find that the best-fit polynomials (with p , p i in dbar) T f = 4.7184 × 10 −3 − 7.4584 × 10 −4 p i − 1.4999 × 10 −8 p i (2) T d = 3.9795 − 2.0059 × 10 −3 p − 6.2511 × 10 −8 p (3)approximate T f e and T d e to within 0.002 K. We approximate the density of water by a quadratic bivariate polynomial, which is maximum at T = T d ( p ). For 0 < p < 10 dbar and T f < T < T f + 15 K, we find that the best-fit bivariate polynomial (with p in dbar) = + ( p ) + C ( p ) [ T − T d ( p ) ] (4)with = 9.9999 × 10 , = 4.9195 × 10 −3 p − 1.4372 × 10 −8 p C = − 7.0785 × 10 −3 + 1.8217 × 10 −7 p + 4.2679 × 10 −12 p (5)approximates e to within less than 0.01% relative error. This implies density errors less than 0.1 kg/m , 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 ─ ∂ ─ ∂ T ∣ p = − 2 C ( p ) [ T − T d ( p ) ] ─────────── (6)which is shown in Fig. 2 and changes sign at T = T d > T f for pres-sures 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 atmo-spheric pressure ( ) 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 com-pressibility effects in the energy equation because we are interested in the calculation of the (small) critical heat flux at the onset of con-vection, but compressibility effects can otherwise be neglected when considering the (large) geothermal heat flux. In a Cartesian coordi-nates system ( x , y , and z ) centred on a lake’s top boundary, the equa-tions for the velocity vector u , pressure p , density , and temperature T (in °C) read ( – ) D u ─ Dt + f e z × u = − ∇ p + ∇ u − g e z (7) ∇ ⋅ u = 0 (8) c p DT ─ Dt − ( T + T ) Dp ─ Dt = k ∇ T (9) where f is the Coriolis frequency, is the dynamic viscosity, k is the thermal conductivity, c p is the isobaric specific heat capacity, g is the gravitational acceleration, and T = 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 in-compressible fluid; and Eq. 9 is the energy equation including pres-sure 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 10 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 c p = 4.2174 × 10 J kg −1 K −1 ; and we recall that g = 9.81 m/s 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 pres-sure and temperature conditions of our interest such that we take them constants ( , ). We denote by = / = 1.7 × 10 −6 m /s the constant kinematic viscosity and by = k / c p = 1.3 × 10 −7 m /s the constant thermal diffusivity.Subglacial lake water must be at the freezing temperature at the upper lake boundary, i.e., T = T f ( p i ) at z = 0 m, while at the base of the lake, it is the heat flux that is enforced, i.e., k∂ z T = − F at z = − h , with F > 0 the (geothermal) heat flux and h > 0 the lake depth; also, p = p i 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 = T f ( p i ) − zF ─ k , d z ¯ p = − ( ¯ T , ¯ p ) g (10)i.e., the temperature increases linearly with depth, and the pressure is hydrostatic. For simplicity, we assume ≈ in the hydrostatic base-state equation, such that ¯ p = p i − gz at leading order (assum-ing 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 = p i − gz (overbar dropped) to be locally (statically) stable is ( , ) 1 ─ ∂ ─ ∂ T ∣ p ds ─ dz = − ( c p ─ T + T dT ─ dz − ─ dp ─ dz ) = − c p ─ T + T ( dT ─ dz − dT ad ─ dz ) < 0 (11)where s is the entropy, T = 273.15 K (we recall that we express tem-perature variables in °C), and we approximate ≈ at leading order. dT ad / dz is known as the adiabatic temperature gradient and reads dT ad ─ dz = − ( T + T ) g ─ c p (12)such that dT ad / 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 < dT ad / dz < 0, i.e., if the heat flux exceeds (in ab-solute value) the adiabatic heat flux. Thus, the two conditions for on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
The criterion for an ideal (dissipationless) compressible fluid with hydrostatic base-state pressure p = p i − gz (overbar dropped) to be locally (statically) stable is ( , ) 1 ─ ∂ ─ ∂ T ∣ p ds ─ dz = − ( c p ─ T + T dT ─ dz − ─ dp ─ dz ) = − c p ─ T + T ( dT ─ dz − dT ad ─ dz ) < 0 (11)where s is the entropy, T = 273.15 K (we recall that we express tem-perature variables in °C), and we approximate ≈ at leading order. dT ad / dz is known as the adiabatic temperature gradient and reads dT ad ─ dz = − ( T + T ) g ─ c p (12)such that dT ad / 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 < dT ad / dz < 0, i.e., if the heat flux exceeds (in ab-solute value) the adiabatic heat flux. Thus, the two conditions for on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E subglacial lakes experiencing geothermal heating (i.e., such that dT / dz < 0) to be locally unstable are > 0 (13) dT ─ dz < dT ad ─ dz (14)Note that both and dT ad / dz are functions of z when considering the base state of a subglacial lake heated from below. Specifically, increases with depth, while dT ad / 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 gradi-ent with a buoyancy anomaly that is also large enough to overcome viscosity and diffusivity effects. The calculation of the exact, mini-mum 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 investi-gating how small initial perturbations evolve over time. We expand the variables (generically represented by X ) as X = ¯ 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 ∂ u ′ ─ ∂ t + f e z × u ′= − ∇ p ′+ ∇ u ′+ ¯ T ′ g e z − ¯ p p ′ g e z (16) ∇ ⋅ u ′= 0 (17) c p ( ∂ t T ′− w ′ F / k ) − ¯ ( ¯ T + T ) ( ∂ t p ′− w ′ g ) = k ∇ T ′ (18)where ¯ p is related to the small compressibility of the background state and is derived from Eq. 4 as ¯ p = [ + 2 ¯ p + ( ¯ T − ¯ T d ) ( C + 2 C ¯ p ) + 2 ¯ C ( ¯ T − ¯ T d ) (− T d − 2 T d ¯ p ) ] (19)with subscripts 1,2 denoting the linear and leading coefficients of the polynomial expressions for , T d , and C , e.g., T d1 = − 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 form X ′( x , y , z , t ) = ˆ X ( z ) e t + i ( k x x + k y y ) + c . c . , (20) with as the growth rate, k x and k y as the wave numbers in the x and y directions, and c.c. as the complex conjugate. Assuming horizontal isotropy, i.e., k x = k y , 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 ( ). 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 ⊥ = √ _ k x + k y and heat flux F for a range of input parameters ( p i , h ). The critical minimum heat flux F c 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 F c in Fig. 3A. Note that the eigenvalue problem becomes challenging at large h , such that we limit the cal-culations 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 F c . 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 con-vection than geostrophic convection. As a result, we use scaling laws obtained in the idealized limit of nonrotating turbulent convection to make predictions about , T bulk , U , and 𝓁 for subglacial lakes sub-ject to a geothermal flux F = 50 mW/m .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 p i < p * . In other words, we write = t + S with t as the thick-ness of the turbulent boundary layer and S as the thickness of the stable layer. We obtain S as the positive solution of the equation T d ( p i + g S /10 ) = T f ( p i ) + S F / k , i.e., S is found as the location z = − S where the temperature of maximum density equals the tem-perature of the conductive base-state profile. The temperature at the base of the stable layer is correspondingly T S = T f + S F / k . The definition of the control parameter Ra F (Eq. 1) of subglacial lakes includes the effec-tive water depth h eff 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 h eff = h − S . Estimating eff inside a convective lake is difficult since the vertical pro-files of temperature and are unknown. Here, for simplicity, we use eff = ( p i + gh /10 , T S ), i.e., we take on the bottom boundary as the effective thermal expansion coefficient but assume that the tempera-ture of the lake does not exceed T S . This assumption is likely to under-estimate 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 ( , ) to predict t , T bulk , U , and 𝓁 . We re-late our flux-based Rayleigh number Ra F to Ra following previous works ( , ), i.e., such that Ra F = RaNu (21)The scaling laws for Nu in ( ) and Re in ( ) are Nu = 0.16 Ra (22) Re = 0.18 ( Ra − Ra ─ Nu ) Pr −1 (23) on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
We show in the Supplementary Materials that vertical convection in Antarctic subglacial lakes is better represented by nonrotating con-vection than geostrophic convection. As a result, we use scaling laws obtained in the idealized limit of nonrotating turbulent convection to make predictions about , T bulk , U , and 𝓁 for subglacial lakes sub-ject to a geothermal flux F = 50 mW/m .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 p i < p * . In other words, we write = t + S with t as the thick-ness of the turbulent boundary layer and S as the thickness of the stable layer. We obtain S as the positive solution of the equation T d ( p i + g S /10 ) = T f ( p i ) + S F / k , i.e., S is found as the location z = − S where the temperature of maximum density equals the tem-perature of the conductive base-state profile. The temperature at the base of the stable layer is correspondingly T S = T f + S F / k . The definition of the control parameter Ra F (Eq. 1) of subglacial lakes includes the effec-tive water depth h eff 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 h eff = h − S . Estimating eff inside a convective lake is difficult since the vertical pro-files of temperature and are unknown. Here, for simplicity, we use eff = ( p i + gh /10 , T S ), i.e., we take on the bottom boundary as the effective thermal expansion coefficient but assume that the tempera-ture of the lake does not exceed T S . This assumption is likely to under-estimate 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 ( , ) to predict t , T bulk , U , and 𝓁 . We re-late our flux-based Rayleigh number Ra F to Ra following previous works ( , ), i.e., such that Ra F = RaNu (21)The scaling laws for Nu in ( ) and Re in ( ) are Nu = 0.16 Ra (22) Re = 0.18 ( Ra − Ra ─ Nu ) Pr −1 (23) on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
10 of 11
Combining Eq. 1 with Eqs.21 and 22, we can estimate Ra = ( Ra F /0.16) and predict t = 0.5 h eff / Nu = 0.5 h eff /(0.16 Ra ), such that = 0.5 h eff ─ Ra + S (24)Assuming a conductive temperature profile in the stable and tur-bulent boundary layers, we then find T bulk = F ─ k (25)The turbulent flow velocity is inferred from Re = Uh eff / and Eqs. 22 and 23 as U = √ _____________ Ra − 6.25 Ra ───────────── Pr h eff (26)The estimate of the characteristic turbulent length scale, or dis-tance between plumes, is finally obtained from equation (6.3) of reference ( ) as ℓ = 0.8 √ _ RePr (3.125 Ra −2/7 ) h eff (27)We provide a second set of predictions for t , T bulk , and U lsc based on scaling laws inferred from the GL unifying theory of RB convec-tion ( ). 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 ( ). The procedure for de-riving t , T bulk , and U lsc 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.5 h eff / Nu , T bulk = F / k ( S is unchanged) and U lsc = Re / h eff . The scaling laws for Nu and Re combine the expressions derived in subregions I u , III u , and IV u of the GL theory [see ( ) and the Sup-plementary 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 Ek L and the Rayleigh number Ra L based on the lake’s horizontal length L ( , ), i.e. Ek L = ─ ∣ f ∣ L , Ra L = g i i L ─ (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 T f ( p i ) with the ice pressure (Eq. 2). The ice pressure drop along the ice-water interface is pi = sL i g /10 dbar, with s the slope and L in meters, such that i = T f ( p i ) − T f ( p i + p i ) > 0 (29)For i , we take the maximum of along the ice-water interface, i.e. i = max [ ∣ ( p i , T f ( p i ) ) ∣, ∣ ( p i + p i , T f ( p i + p i ) ) ∣] (30) We show in the Supplementary Materials that horizontal con-vection in Antarctic subglacial lakes is constrained by rotation. Thus, here we use a recent scaling law for the Reynolds number Re hc of rotation-constrained horizontal convection ( ) to estimate the char-acteristic velocity of the large-scale horizontal flow V hc . 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 vari-able 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 ( ) is Re hc = ( Ra L Ek L ) ─ Pr (31)such that V hc = Re hc ─ L (32)Note that V hc ∼ L , such that estimates for the horizontal veloc-ity would be only weakly affected (weakly decreasing) by decreasing the aspect ratio. SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/7/8/eabc3972/DC1
REFERENCES AND NOTES
1. M. J. Siegert, F. Florindo, Antarctic climate evolution.
Dev. Earth Environ. Sci. , 1–11 (2008). 2. A. Wright, M. Siegert, A fourth inventory of Antarctic subglacial lakes. Antarct. Sci. , 659–664 (2012). 3. M. J. Siegert, RESEARCH FOCUS: A wide variety of unique environments beneath the Antarctic ice sheet. Geology , 399–400 (2016). 4. M. J. Siegert, J. C. Ellis-Evans, M. Tranter, C. Mayer, J.-R. Petit, A. Salamatin, J. C. Priscu, Physical, chemical and biological processes in Lake Vostok and other Antarctic subglacial lakes. Nature , 603–609 (2001). 5. C. S. Cockell, E. Bagshaw, M. Balme, P. Doran, C. P. Mckay, K. Miljkovic, D. Pearce, M. J. Siegert, M. Tranter, M. Voytek, J. Wadham, Subglacial environments and the search for life beyond the Earth.
Antarct. Subglacial Aquat. Environ. , 129–148 (2011). 6. B. E. Smith, H. A. Fricker, I. R. Joughin, S. Tulaczyk, An inventory of active subglacial lakes in Antarctica detected by ICESat (2003–2008).
J. Glaciol. , 573–595 (2009). 7. M. R. Siegfried, H. A. Fricker, Thirteen years of subglacial lake activity in Antarctica from multi-mission satellite altimetry. Ann. Glaciol. , 42–55 (2018). 8. B. C. Christner, J. C. Priscu, A. M. Achberger, C. Barbante, S. P. Carter, K. Christianson, A. B. Michaud, J. A. Mikucki, A. C. Mitchell, M. L. Skidmore, T. J. Vick-Majors, W. P. Adkins, S. Anandakrishnan, G. Barcheck, L. Beem, A. Behar, M. Beitch, R. Bolsey, C. Branecky, R. Edwards, A. Fisher, H. A. Fricker, N. Foley, B. Guthrie, T. Hodson, H. Horgan, R. Jacobel, S. Kelley, K. D. Mankoff, E. McBryan, R. Powell, A. Purcell, D. Sampson, R. Scherer, J. Sherve, M. Siegfried, S. Tulaczyk; the WISSARD Science Team, A microbial ecosystem beneath the West Antarctic ice sheet. Nature , 310–313 (2014). 9. M. Thoma, K. Grosfeld, I. Filina, C. Mayer, Modelling flow and accreted ice in subglacial Lake Concordia, Antarctica.
Earth Planet. Sci. Lett. , 278–284 (2009). 10. J. Woodward, A. M. Smith, N. Ross, M. Thoma, H. F. J. Corr, E. C. King, M. A. King, K. Grosfeld, M. Tranter, M. J. Siegert, Location for direct access to subglacial Lake Ellsworth : An assessment of geophysical data and modeling.
Geophys. Res. Lett. , L11501 (2010). 11. A. M. Smith, J. Woodward, N. Ross, M. J. Bentley, D. A. Hodgson, M. J. Siegert, E. C. King, Evidence for the long-term sedimentary environment in an Antarctic subglacial lake. Earth Planet. Sci. Lett. , 139–151 (2018). 12. M. J. Siegert, J. A. Dowdeswell, Spatial variations in heat at the base of the Antarctic ice sheet from analysis of the thermal regime above subglacial lakes.
J. Glaciol. , 501–509 (1996). 13. G. Royston-Bishop, J. C. Priscu, M. Tranter, B. Christner, M. J. Siegert, V. Lee, Incorporation of particulates into accreted ice above subglacial Vostok lake, Antarctica. Ann. Glaciol. , 145–150 (2005). on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
J. Glaciol. , 501–509 (1996). 13. G. Royston-Bishop, J. C. Priscu, M. Tranter, B. Christner, M. J. Siegert, V. Lee, Incorporation of particulates into accreted ice above subglacial Vostok lake, Antarctica. Ann. Glaciol. , 145–150 (2005). on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ouston and Siegert, Sci. Adv. : eabc3972 17 February 2021 S C I E N C E A D V A N C E S | R E S E A R C H A R T I C L E
11 of 11
14. M. Thoma, K. Grosfeld, C. Mayer, Modelling mixing and circulation in subglacial Lake Vostok, Antarctica.
Ocean Dyn. , 531–540 (2007). 15. M. Thoma, K. Grosfeld, C. Mayer, A. M. Smith, J. Woodward, N. Ross, The “tipping” temperature within Subglacial Lake Ellsworth, West Antarctica and its implications for lake access. Cryosph. , 561–567 (2011). 16. A. Rivera, J. Uribe, R. Zamora, J. Oberreuter, Subglacial Lake CECs: Discovery and in situ survey of a privileged research site in West Antarctica. Geophys. Res. Lett. , 3944–3953 (2015). 17. D. Bouffard, A. Wüest, Convection in Lakes. Annu. Rev. Fluid Mech. , 189–215 (2019). 18. A. Wüest, E. Carmack, A priori estimates of mixing and circulation in the hard-to-reach water body of Lake Vostok. Ocean Model. , 29–43 (2000). 19. T. J. McDougall, P. M. Barker, Getting started with TEOS-10 and the Gibbs Seawater (GSW) (Oceanographic Toolbox, 2011), p. 28. 20. G. Ahlers, S. Grossmann, D. Lohse, Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection.
Rev. Mod. Phys. , 503–537 (2009). 21. M. Thoma, K. Grosfeld, A. M. Smith, C. Mayer, A comment on the Equation of State and the freezing point equation with respect to subglacial lake modelling. Earth Planet. Sci. Lett. , 80–84 (2010). 22. S. Chandrasekhar,
Hydrodynamic and Hydromagnetic Stability (Dover, 1961). 23. J. M. Aurnou, M. A. Calkins, J. S. Cheng, K. Julien, E. M. King, D. Nieves, K. M. Soderlund, S. Stellmach, Rotating convective turbulence in Earth and planetary cores.
Phys. Earth Planet. Inter. , 52–71 (2015). 24. R. E. Ecke, J. J. Niemela, Heat transport in the geostrophic regime of rotating rayleigh-bénard convection.
Phys. Rev. Lett. , 114301 (2014). 25. E. M. King, S. Stellmach, B. Buffett, Scaling behaviour in Rayleigh-Bénard convection with and without rotation.
J. Fluid Mech. , 449–471 (2013). 26. S. Grossmann, D. Lohse, Scaling in thermal convection: A unifying theory.
J. Fluid Mech. , 27–56 (2000). 27. C. A. Vreugdenhil, B. Gayen, R. W. Griffiths, Transport by deep convection in basin-scale geostrophic circulation: Turbulence-resolving simulations.
J. Fluid Mech. , 681–719 (2019). 28. J. S. Bowling, S. J. Livingstone, A. J. Sole, W. Chu, Distribution and dynamics of Greenland subglacial lakes.
Nat. Commun. , 2810 (2019). 29. M. G. Wells, J. S. Wettlaufer, Circulation in Lake Vostok: A laboratory analogue study. Geophys. Res. Lett. , L03501 (2008). 30. K. Winter, J. Woodward, N. Ross, S. A. Dunning, A. S. Hein, M. J. Westoby, R. Culberg, S. M. Marrero, D. M. Schroeder, D. E. Sugden, M. J. Siegert, Radar-detected englacial debris in the West Antarctic ice sheet. Geophys. Res. Lett. , 10454–10462 (2019). 31. W. Munk, C. Wunsch, Abyssal recipes II: Energetics of tidal and wind mixing. Deep Sea Res. Part 1 Oceanogr. Res. Pap. , 1977–2010 (1998). 32. J. C. Priscu, E. E. Adams, W. B. Lyons, M. A. Voytek, D. W. Mogk, R. L. Brown, C. P. Mckay, C. D. Takacs, K. A. Welch, C. F. Wolf, J. D. Kirshtein, R. Avci, Geomicrobiology of subglacial ice above Lake Vostok, Antarctica. Science , 2141–2144 (1999). 33. C. Gura, S. O. Rogers, Metatranscriptomic and metagenomic analysis of biological diversity in Subglacial Lake Vostok (Antarctica).
Biology , 55 (2020). 34. M. J. Siegert, R. J. Clarke, M. Mowlem, N. Ross, C. S. Hill, A. Tait, D. Hodgson, J. Parnell, M. Tranter, D. Pearce, M. J. Bentley, C. Cockell, M. N. Tsaloglou, A. Smith, J. Woodward, M. P. Brito, E. Waugh, Clean access, measurement, and sampling of Ellsworth subglacial lake: A method for exploring deep antarctic subglacial lake environments. Rev. Geophys. , RG1003 (2012). 35. B. R. Sutherland, Internal Gravity Waves (Cambridge Univ. Press, 2010). 36. L.-A. Couston, D. Lecoanet, B. Favier, M. Le Bars, The energy flux spectrum of internal waves generated by turbulent convection.
J. Fluid Mech. , (2018). 37. J. Magnaudet, M. J. Mercier, Particles, drops, and bubbles moving across sharp interfaces and stratified layers.
Annu. Rev. Fluid Mech. , 61–91 (2020). 38. S. D. Vance, The habitability of icy ocean worlds in the solar system. Handb. Exopl. , 2855–2877 (2018). 39. T. R. Osborn, P. H. LeBlond, Static stability in freshwater lakes.
Limnol. Oceanogr. , 544–545 (1974). 40. E. A. Spiegel, G. Veronis, On the Boussinesq approximation for a compressible fluid. Astrophys. J. , 442–447 (1960). 41. G. K. Vallis,
Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation (Cambridge Univ. Press, ed. 2, 2017). 42. T. Alboussière, Y. Ricard, Rayleigh-Bénard stability and the validity of quasi-Boussinesq or quasi-anelastic liquid approximations.
J. Fluid Mech. , 264–305 (2017). 43. P. Först, F. Werner, A. Delgado, The viscosity of water at high pressures – especially at subzero degrees centigrade.
Rheol. Acta. , 566–573 (2000). 44. M. L. Huber, R. A. Perkins, D. G. Friend, J. V. Sengers, M. J. Assael, I. N. Metaxa, K. Miyagawa, R. Hellmann, E. Vogel, New international formulation for the thermal conductivity of H O. J. Phys. Chem. Ref. Data Monogr. , 033102 (2012). 45. L. D. Landau, L. M. Lifschitz, Fluid Mechanics (Pergamon Press, 1959). 46. K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods.
Phys. Rev. Res. , 23068 (2020). 47. E. M. King, S. Stellmach, J. M. Aurnou, Heat transfer by rapidly rotating Rayleigh-Bénard convection. J. Fluid Mech. , 568–582 (2012). 48. H. Johnston, C. R. Doering, Comparison of turbulent thermal convection between conditions of constant temperature and constant flux.
Phys. Rev. Lett. , 64501 (2009). 49. E. H. Anders, B. P. Brown, J. S. Oishi, Accelerated evolution of convective simulations.
Phys. Rev. Fluids , 083502 (2018). 50. C. A. Vreugdenhil, B. Gayen, R. W. Griffiths, Mixing and dissipation in a geostrophic buoyancy-driven circulation. J. Geophys. Res. Ocean , 3372–3380 (2016). 51. H. J. S. Fernando, D. C. Smith IV, Vortex structures in geophysical convection.
Eur. J. Mech. B Fluids , 437–470 (2001). 52. R. P. J. Kunnen, R. Ostilla-Mónico, E. P. Van Der Poel, R. Verzicco, D. Lohse, Transition to geostrophic convection: The role of the boundary conditions. J. Fluid Mech. , 413–432 (2016). 53. J. S. Cheng, J. M. Aurnou, K. Julien, R. P. J. Kunnen, A heuristic framework for next-generation models of geostrophic convective turbulence.
Geophys. Astrophys. Fluid Dyn. , 277–300 (2018). 54. M. Plumley, K. Julien, Scaling laws in Rayleigh-Bénard convection.
Earth Sp. Sci. , 1580–1592 (2019). 55. E. H. Anders, C. M. Manduca, B. P. Brown, J. S. Oishi, G. M. Vasil, Predicting the Rossby number in convective experiments. Astrophys. J. , 138 (2019). 56. K. Julien, S. Legg, J. McWilliams, J. Werne, Rapidly rotating turbulent Rayleigh-Bénard convection.
J. Fluid Mech. , 243–273 (1996). 57. L. E. Peters, S. Anandakrishnan, C. W. Holland, H. J. Horgan, D. D. Blankenship, D. E. Voigt, Seismic detection of a subglacial lake near the South Pole, Antarica.
Geophys. Res. Lett. , L23501 (2008). 58. M. Studinger, R. E. Bell, G. D. Karner, A. A. Tikku, J. W. Holt, D. L. Morse, T. G. Richter, S. D. Kempf, M. E. Peters, D. D. Blankenship, R. E. Sweeney, V. L. Rystrom, Ice cover, landscape setting, and geological framework of Lake Vostok, East Antarctica. Earth Planet. Sci. Lett. , 195–210 (2003).
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.Submitted 22 April 2020Accepted 4 January 2021Published 17 February 202110.1126/sciadv.abc3972
Citation:
L.-A. Couston, M. Siegert, Dynamic flows create potentially habitable conditions in Antarctic subglacial lakes.
Sci. Adv. , eabc3972 (2021). on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm
Sci. Adv. , eabc3972 (2021). on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm ynamic flows create potentially habitable conditions in Antarctic subglacial lakes Louis-Alexandre Couston and Martin SiegertDOI: 10.1126/sciadv.abc3972 (8), eabc3972. Sci Adv
ARTICLE TOOLS http://advances.sciencemag.org/content/7/8/eabc3972
MATERIALSSUPPLEMENTARY http://advances.sciencemag.org/content/suppl/2021/02/12/7.8.eabc3972.DC1
REFERENCES http://advances.sciencemag.org/content/7/8/eabc3972
PERMISSIONS
Science Advances
York Avenue NW, Washington, DC 20005. The title (ISSN 2375-2548) is published by the American Association for the Advancement of Science, 1200 New
Science Advances
BY).Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution License 4.0 (CC Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of on F eb r ua r y , tt p :// ad v an c e s . sc i en c e m ag . o r g / D o w n l oaded f r o mm