Layered seawater intrusion and melt under grounded ice
- 1School of Earth and Atmospheric Sciences, Georgia Institute of Technology, Atlanta, GA, USA
- 2Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA, USA
- 3Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA
- 4Thayer School of Engineering, Dartmouth College, Hanover, NH, USA
Correspondence: Alexander A. Robel (firstname.lastname@example.org)
Increasing melt of ice sheets at their floating or vertical interfaces with the ocean is a major driver of marine ice sheet retreat and sea level rise. However, the extent to which warm, salty seawater may drive melting under the grounded portions of ice sheets is still not well understood. Previous work has explored the possibility that dense seawater intrudes beneath relatively light subglacial freshwater discharge, similar to the “salt wedge” observed in many estuarine systems. In this study, we develop a generalized theory of layered seawater intrusion under grounded ice, including where subglacial hydrology occurs as a macroporous water sheet over impermeable beds or as microporous Darcy flow through permeable till. Using predictions from this theory, we show that seawater intrusion over flat or reverse-sloping impermeable beds may feasibly occur up to tens of kilometers upstream of a glacier terminus or grounding line. On the other hand, seawater is unlikely to intrude more than tens of meters through permeable till. Simulations using the Ice-sheet and Sea-level System Model (ISSM) show that even just a few hundred meters of basal melt caused by seawater intrusion upstream of marine ice sheet grounding lines can cause projections of marine ice sheet volume loss to be 10 %–50 % higher. Kilometers of intrusion-induced basal melt can cause projected ice sheet volume loss to more than double. These results suggest that further observational, experimental and numerical investigations are needed to determine the conditions under which seawater intrusion occurs and whether it will indeed drive rapid marine ice sheet retreat and sea level rise in the future.
Where ice sheets come into contact with seawater, ice may be lost through melting and dissolution. Increasing ocean-induced ice sheet melt has led to glacier retreat and upstream glacier thinning in Greenland and Antarctica, contributing significantly to global mean sea level rise (Straneo and Heimbach, 2013; Shepherd et al., 2018). The grounding line (the location where glacier ice loses contact with the bed) has traditionally been considered a hydraulic barrier to the intrusion of seawater beneath grounded ice due to the horizontal hydropotential gradient imposed by the increasing weight of a glacier that thickens upstream. Previous studies have hypothesized that seawater intrusion may only be possible when transient tidal variations in the hydrostatic pressure of seawater at the grounding line either overcome this hydraulic barrier (Walker et al., 2013) or cause it to migrate upstream from the grounding line (Sayag and Worster, 2013). In such models, the compositional difference of seawater and subglacial discharge is not considered, leaving the horizontal variation in overburden pressure as the primary control on the subglacial intrusion of seawater. However, the interface between freshwater discharge and seawater can have a vertical structure in which dense seawater intrudes horizontally under the lighter freshwater forming a “salt wedge”, as has been observed in estuaries (Geyer and Ralston, 2011), enclosed wastewater outfalls (Adams et al., 1994) and coastal karst channels (Dermissis, 1993).
Wilson et al. (2020) first showed in theory and experiments that layered seawater intrusion is possible in laterally confined subglacial channels (i.e., Rothlisberger or Nye channels) and may extend several kilometers upstream of glacier termini under realistic conditions. A diverse range of observations have also indicated the possibility for seawater intrusion upstream of grounding lines, including the absence of a transition in bed reflectivity across the grounding line at Whillans and Kamb ice streams (MacGregor et al., 2011), a large subglacial channel crossing the grounding line of Whillans ice stream imaged with active-source seismic methods (Horgan et al., 2013), a similarly large radar-imaged channel crossing the grounding line at Roi Baudouin Ice Shelf (Drews et al., 2017) and elevated subglacial melt upstream of the grounding line measured by satellite-based interferometric synthetic aperture radar at Thwaites Glacier (Milillo et al., 2019). Such observations provide, at best, plausible arguments for the presence of seawater upstream of the grounding line but provide few constraints on the composition, temperature and vertical structure of these potential intrusions.
Understanding the extent of possible seawater intrusion is also directly relevant to projections of mass loss of marine ice sheets, which are highly sensitive to the intensity of ocean melt occurring near the grounding line (Arthern and Williams, 2017; Reese et al., 2018; Goldberg et al., 2019) and the manner in which ocean melt is applied at the grounding line in models (Seroussi and Morlighem, 2018). Ice sheet models which unintentionally include ocean melt kilometers upstream of the grounding line (due to numerical inaccuracies) simulate up to 2 times more Antarctic ice sheet mass loss in response to climate change over the next century (Seroussi et al., 2019). Thus, investigating the possibility that seawater intrusion may drive melt under grounded ice is of first-order importance to the problem of accurately simulating the response of marine ice sheets to climate change.
In this study, we generalize the channelized intrusion theory of Wilson et al. (2020) to a broader range of subglacial hydrology types which encompass most glacier grounding lines and termini (Sect. 2). This generalized theory predicts that for the wide range of parameters describing subglacial hydrology, seawater intrusion can either be suppressed entirely or extend tens of kilometers inland (Sect. 3). In a state-of-the-art ice sheet model, we show that the transient rate of ice sheet mass loss depends sensitively on the extent of warm seawater intrusion upstream of marine ice sheet grounding lines (Sect. 4). Finally, we suggest observations which may be useful to better constrain seawater intrusion and subglacial hydrology near the ice–ocean interface and discuss the implications for simulations of future ice sheet mass loss in response to ocean warming (Sects. 5 and 6).
The goal of this study is to consider the possibility that seawater will intrude under grounded portions of ice sheets. In this section, we consider a generalized theory for the horizontal distance over which this intrusion may occur in the cases of hard (i.e., impermeable) and soft (i.e., permeable) beds. In both cases, the intrusion distance depends sensitively on the characteristics of the subglacial hydrological system. In Sect. 3, we discuss the range of predictions from this theory and the relationship to observations of subglacial hydrology.
2.1 Hard beds
On hard beds that are impermeable to vertical drainage of water, observations and theory indicate that subglacial hydrology organizes into systems that are either macroporous sheets (i.e., “inefficient” or linked-cavity drainage; Creyts and Schoof, 2009) or channelized (i.e., “efficient” drainage; Rothlisberger, 1972). If such a system is sufficiently permeable in horizontal directions, the flow can be described by the shallow water equations. If the system is less permeable in the horizontal (i.e., microporous), inertial terms become unimportant, and the flow of water is better described by Darcy's law, as we discuss further in Sect. 2.2.
Wilson et al. (2020) originally considered the dynamics of a two-layer shallow water flow for the case of a laterally and vertically confined subglacial channel on an impermeable bed, with free-flowing, fresh subglacial discharge entering a large saltwater body at rest. However, where efficient channelization does not occur, subglacial water flow likely occurs through macroporous sheets (i.e., not laterally confined) which are kept open by water pressure and the influence of large clasts protruding from the bed (Creyts and Schoof, 2009; Hewitt, 2011). We generalize the seawater intrusion theory of Wilson et al. (2020) by considering two-layer macroporous water flow through a vertically confined subglacial sheet over an impermeable bed (illustrated in Fig. 1a). We will show later that the channelized system of Wilson et al. (2020) is a special intermediate case of the generalized theory derived in this study.
We consider a vertically confined, two-layer (fresh and saline) shallow water system, in which the grounding line or terminus is located at X=0 and the layers are confined under ice where X<0 (Fig. 1). Mass conservation in both layers is given by
where H1 is the freshwater layer thickness, H2 is the seawater layer thickness, is the thickness of the full water sheet (i.e., the confined thickness between the bed and ice), and Q1=H1U1 and Q2=H2U2 are the area fluxes of water of each layer, with fluid velocities U1 and U2.
In the two-layer macroporous system, flow may be influenced by the relative buoyancy of the two layers, external barotropic pressure gradients, drag from overlying ice and underlying bed, drag at the interface of the two layers, and drag from the embedded obstacles (i.e., clasts). The momentum balance for both layers is then
where X is the distance from the grounding line (with X<0 upstream of the grounding line), t is time, is the along-flow barotropic pressure gradient, Cice, Cbed, Cobs and Ci are the drag coefficients for flow past the ice, bed, obstacles and layer interface, respectively, ϕ1,2 are the bulk porosities of the sheet in each layer (defined as the fraction of the water sheet volume occupied by obstacles), d1,2 are the mean obstacle diameters in the layers, is the reduced gravity for density difference Δρ between the two layers, and tan θ is the bed slope. The main difference between this two-layer flow and that described in Wilson et al. (2020) is that the final term in Eqs. (3) and (4) is given by the bulk drag due to flow through a field of obstacles (Rominger and Nepf, 2011) instead of drag due to friction at side walls. is related to a length scale of the rough surfaces of the obstacles that are encountered by water flow. We assume hereafter that the drag coefficient is the same for the ice, bed and obstacles () within the macroporous water sheet. We also assume that the obstacles' diameters and porosities are the same in both layers of the water sheet (, ). Though these are simplifying assumptions, the uncertainty of drag coefficients and obstacle properties is sufficiently large that it encompasses any variations that exist within the macroporous water sheet.
As can be seen in Eqs. (3) and (4), the externally imposed barotropic pressure gradient (), which is known as the subglacial hydropotential gradient among glaciologists, acts equally on both layers. A horizontal variation in this barotropic pressure may result from ice surface slope or a difference in slope between the ice base and bed (which produces a gradient in the total subglacial layer thickness, H). However, because the fresh and saline layers are vertically layered, barotropic pressure gradients act equally on both layers. Consequently, when the momentum equations are ultimately combined (below), the pressure gradient falls out of the solution for intrusion distance. This model of subglacial seawater intrusion is focused on the horizontal extent of the saline layer, to which the hydropotential gradient is not directly relevant, and is fundamentally different than previous efforts which focus on the horizontal interface between subglacial discharge and seawater without considering their compositional differences as we do here (through the inclusion of buoyancy and consideration of the potential for vertical layering). Nonetheless, subglacial barotropic pressure far upstream of the grounding line does ultimately set the inflowing freshwater flow velocity and thus does have an indirect influence on the water velocity scale (as we explore further below).
The velocity scale is set by the local dimensionless densimetric Froude number, , which changes with the local freshwater velocity. Where freshwater discharge is sufficiently slow, the flow of water may be considered subcritical (i.e., Fr <1), and mixing between the two layers is negligible (which would otherwise enter as additional terms in the above equations) as in slow estuarine and subglacial flows, as shown in experiments (Wilson et al., 2020) and observations (Ralston et al., 2010).
We proceed by considering the steady-state solution to Eqs. (1)–(4) where the saline layer is at rest (U2=0). We also non-dimensionalize Eqs. (1)–(4) by making the following substitutions: , , , , and . C0 is a characteristic scale for drag coefficients, and γ is a dimensionless parameter capturing the resistance to flow by the macroporous matrix. Combining Eqs. (1)–(4) to eliminate the barotropic pressure gradient and using the relationship , we arrive at the following differential equation for the non-dimensional freshwater layer thickness (h):
The Froude number can be written in terms of a Froude number scale:
where is the Froude number of the freshwater inflow upstream where it occupies the entire water sheet (). In the above equation, Fr changes as the thickness and also velocity of the freshwater layer both change. We note that when Θ=0, it must be the case that the sign of is set by Fr2−1. Since we assume in this study that Fr≤1, h always increases upstream of the grounding line, requiring that h=1 as .
At the grounding line or terminus (x=0), the freshwater flow becomes unconfined and hence supercritical (i.e., Fr=1 in Eq. 6), so we set the boundary condition:
The dimensionless intrusion distance, ℓ, is where the freshwater layer first occupies the entire subglacial water layer thickness () upstream from the grounding line, which is also where the saline layer no longer occupies any of the subglacial water layer. In Appendix A, we describe how Eq. (5) can be numerically integrated from x=0 until in order to find the intrusion distance, ℓ. We can also derive analytical approximations for ℓ in two limiting cases to describe how the intrusion distance depends on certain parameters describing the subglacial hydrology.
We make three further simplifying assumptions to derive analytical approximations for the intrusion distance. First, we only consider flat beds (Θ=0), though we will numerically explore how sloped beds modify the intrusion distance in Sect. 3.3. Second, we assume that interfacial drag is negligible (Ci=0), as most estimates of interfacial drag between saline and fresh water layers are O(10−4) (MacCready and Geyer, 2010; Geyer and Ralston, 2011), whereas even limiting cases of drag on relatively smooth ice walls are of the order of 10−3, and drag is greater yet for realistic subglacial surfaces. Third, we can assume that subglacial freshwater flow is highly subcritical (Fr0≪1), which is appropriate for places like Antarctica, where there is little injection of surface meltwater to the bed and subglacial water flow velocities are low (centimeters per second or less). With these three assumptions, Eq. (5) reduces to
where we have re-written the dimensionless parameter groups involving the local Froude number in terms of h and Fr0. This differential equation can be analytically integrated for a solution, but the resulting expression is unwieldy and difficult to interpret, so we consider two end-member cases which provide reasonable lower and upper bounds on expected seawater intrusion distance.
The subglacial clasts which obstruct flow in the water sheet may be spaced far apart and still provide sufficient support to maintain a separation between the ice and bed (or water pressure may maintain this separation on its own, as in Hewitt, 2011). In this scenario, the porosity of subglacial obstacles (ϕ) is low and so γh≪1. Equation (8) can then be integrated exactly to yield an approximation for the intrusion distance for an unobstructed water sheet:
In this end-member, the dominant drag on the intrusion of seawater into the subglacial water system is the drag of the water against the bed and against the ice. This is also the limiting case of intrusion in which the width of channels is much greater than their height, considered previously in numerical calculations by Wilson et al. (2020).
The contrasting limiting case is when clasts which obstruct flow within the water sheet are closely spaced, and so the macroporosity of the water sheet is high (though not sufficiently so to produce a non-inertial Darcy flow regime, which we consider separately in Sect. 2.2). In such a case, we expect that γh≫1, in which Eq. (8) is integrated to yield an approximation for the intrusion distance for a macroporous water sheet:
Re-dimensionalizing ℓu and ℓp yields upper and lower bounds on Lh (), the dimensional seawater intrusion distance on hard beds, which can be written as
We define an “intrusion length scale” for hard beds:
which may also be derived from Eq. (8) using dimensional analysis. As discussed later in Sect. (3), we generally expect γ≲2, meaning that the intrusion distance is likely to be a fraction of the intrusion length scale between and .
We briefly note that for the case of water flow through a subglacial channel (described in Wilson et al., 2020), the shallow water equations have identical structure except that , capturing the effect of drag from channel side walls (not obstacle drag as in this model). Observations indicate that the geometry of subglacial channels spans a range from circular (Rothlisberger, 1972), corresponding to γ=2 in Eq. (10), to flat (Nye, 1976), which would correspond to an unobstructed water sheet (Eq. 9). Thus, although we have focused our discussion on the macroporous water sheet, the two limiting cases of ℓu and ℓp also encompass intrusion distances expected in subglacial channels. The approach of calculating two end-member solutions for intrusion distance then provides a general theory that we expect to apply in a wide range of settings with hard beds, regardless of the organization of the subglacial drainage.
In Fig. 2, we evaluate the validity of the assumptions that were made on the way to derive these analytical approximations for the intrusion distance (purple circles) by comparing them to numerically calculated solutions to Eq. (5) (black lines). For large γ (heavily obstructed seawater intrusion), Fig. 2a shows that ℓp (Eq. 10) is an excellent approximation of the intrusion distance if Fr0 is small. Figure 2b shows that the ℓp approximation breaks down as Fr0 becomes O(1). At intermediate γ (Fig. 2c), ℓp is off by O(1) as the γ≫1 approximation becomes less appropriate. Finally, at small γ (unobstructed sheet-like flow), Fig. 2d shows that ℓu (Eq. 9) is a good approximation if Ci≪Cd. We note that, in general, seawater intrusion distance increases as γ and Fr0 decrease (these dependencies are discussed in much more detail in Sect. 3), as predicted by the analytical approximations. Ultimately, where these approximations break down, they are typically too large by factors of O(1). As we discuss later, a lack of constraints on physical parameters which enter into these analytical approximations produces uncertainties in intrusion distance that range over 1–4 orders of magnitude. Thus, these approximations are good enough to provide insights into the physical controls on intrusion distance. However, where accuracy is required, numerically solving Eq. (5) is preferable (as we do in Sect. 3).
2.2 Soft beds
Darcy's law is a non-inertial, creeping, incompressible simplification of the Navier–Stokes equations which describes steady-state fluid flow in microporous media. Subglacial water flow can be described solely by Darcy's law where there is a microporous till layer beneath grounded ice (without significant sheet flow at the ice–till interface). This is in contrast to the layered shallow water flow considered in the previous section where inertial and other terms are important.
The particular problem of seawater intrusion into microporous coastal aquifers where Darcy's law is a valid description of water flow has been extensively studied due to the important implications for drinking water in coastal areas (for review see Werner et al., 2013). The canonical case of seawater intrusion into an aquifer considers Darcy flow where freshwater discharge flows from inland at a prescribed rate towards the ocean (illustrated in Fig. 1b). The most well-known method, due jointly to Dupuit (1863) and Forchheimer (1886), starts with the approximation that freshwater flow within an aquifer (Q1) is set by horizontal gradients of the hydraulic head:
where K is the hydraulic permeability of the microporous medium, and ψ is the elevation of the hydraulic head (which is related to the pore pressure within the till). When the aquifer is confined, the right-hand side is multiplied by to account for H1, the thickness of the freshwater layer, and H, the thickness of the till layer. Ghyben (1888) and Herzberg (1901) later showed that when seawater intrudes into an aquifer, the thickness of the buoyant freshwater layer in an aquifer must be hydrostatically compensated for by the local thickness of seawater leading to an expression for the freshwater layer thickness:
where is a dimensionless parameter determined by the density difference of freshwater (ρf) and saltwater (ρs), and Hs is the sea level relative to the depth of the bottom of the till at the grounding line. The expression above is derived under the assumption that variations in the overburden pressure on the surface of the till layer are small relative to variations in pore pressure due to hydrostatic loading from the ocean. From these two approximations, Strack (1976) showed that by rearranging Eq. (14) for ψ, taking the spatial derivative of Eq. (14), and inserting it into Eq. (13) we arrive at
Integrating this result with respect to X and H1 to find the distance at which the freshwater layer occupies the entire aquifer (H1=H), we finally arrive at the horizontal extent of seawater intrusion inland in a flat, confined aquifer:
where Ls is the dimensional seawater intrusion distance through soft permeable beds. This solution is widely used in the hydrology community and has been shown to predict observed and laboratory-measured seawater intrusion well (Werner et al., 2013).
On sloped beds, the geometry of the saltwater–freshwater interface relative to the bed slope becomes a potentially important factor in determining the intrusion distance. Where the bed deepens away from the ocean (referred to as a retrograde or reverse-slope bed in glaciology), the distance of seawater intrusion will be extended as the seawater intrusion flows downhill. For sloped, confined aquifers, Eq. (14) is
where θ1 and θ2 are the slopes of the lower ice surface and the bed, respectively. For subglacial till layers, we do not expect the lower ice surface slope and bed slope to be significantly different over large areas (which would require the thickness of the till layer to change significantly), so we can safely assume that for our case . This allows us to simplify the above implicit equation and solve explicitly for the intrusion distance, yielding
where seawater intrusion occurs if Ls>0. In the next section we will explore the range of predicted seawater intrusion distances from this soft-bed theory and the hard-bed theory of the previous section.
The theories described in the preceding section can be used to make predictions for the seawater intrusion distance expected over hard and soft beds. In this section we explore the range of seawater intrusion distances that would be predicted to occur, using the wide range of parameter values measured or indirectly inferred from observations of subglacial hydrology near grounding lines and glacier termini. We summarize the range of parameter values we discuss and the corresponding predicted intrusion distances in Table 1.
3.1 Hard beds
There are just a few parameters that play a role in determining the horizontal seawater intrusion distance, as can be deduced by examining either the differential equation that exactly specifies the intrusion distance in the model we consider (Eq. 5) or the analytical approximations of the intrusion distance. In particular the intrusion length scale (Eq. 12) provides an excellent starting point for understanding which parameters play a role in setting the seawater intrusion distance. We start with the reduced gravity of the two layer system, g′, which is generally near 0.27 m s−2 if the fresh and saline layers remain distinct and unmixed (since it depends on the density difference between layers). Compositional variability in either layer may lead to variations in g′, but these are likely to be less than 10 %. As previously discussed, we assume that entrainment between layers is negligible, though we note here that g′ will decrease if entrainment mixes the layers. The drag coefficient, Cd, for water flow past ice or obstacles within the subglacial hydrology is typically taken to be between 10−3 and 10−2 for a range of observations typically under sea ice (Lu et al., 2011) and lab experiments for idealized walls and objects (Ezhova et al., 2018). We note though that these prior observations are experiments for turbulent flow past ice, and laminar flows are likely to have even lower drag coefficients, though exact values are not well-constrained in the literature on water flow past ice.
For realistic clast sizes, Hewitt (2011) estimated, using a mathematical model, that the natural thickness scale for a macroporous subglacial water sheet would be 1–5 cm. Creyts and Schoof (2009) found that the thickness of steady-state subglacial water sheets is set by the size of clasts within the sheet and the water pressure maintained within the system. Antarctic tills are characterized by a wide range of particle sizes (Clarke, 2005), though tills sampled directly from Antarctic ice streams indicate the widespread presence of ploughing boulders up to tens of centimeters in size. Consistent with this, MacGregor et al. (2011) estimated that a subglacial water sheet of at least 20 cm thickness exists in the grounding zones of Whillans and Kamb ice streams, with the possibility that seawater intrusion explains the persistence of high bed reflectivity across the putative grounding line (discussed further in Sect. 5).
In the same mathematical model discussed above, Hewitt (2011) estimated water velocities in a centimeter-scale water sheet of 0.5–1 cm s−1 near the ice sheet margin. Carter and Fricker (2012) modeled water velocities inferred for the Siple Coast from subglacial melt production and subglacial hydropotentials, also finding velocities of approximately 1 cm s−1. Other studies have estimated subglacial water fluxes from basal melt (Joughin et al., 2009; Pattyn, 2010), with fluxes that are consistent with ∼ 1 cm s−1 water velocities through an approx. centimeters-thick water sheet or ∼ 0.1 cm s−1 through a ∼ 10 cm thick water sheet (Re ∼ 50 for both cases, indicating laminar flow). Though a centimeter-scale water sheet may seem too thin to maintain stratification between layers, a highly stratified salt wedge with low mixing rates (Kh ∼ 10−9 m2 s−1 is the molecular diffusivity of salt in water; Kunze, 2003) will mix on a timescale of a day, by which time the stratification will have been renewed by freshwater discharge of the order of centimeters per second from upstream over a kilometer-scale seawater intrusion. For a ∼ 10 cm water sheet, mixing of highly stratified layers would take months. Thus, we consider it completely plausible to maintain distinct layers within centimeter-scale water sheets if mixing between the layers is slow. We also note that even at more moderate levels of mixing (Kh ∼ 10−7 m2 s−1 due to, for example, double diffusive convection; van der Boog et al., 2021), it is still possible to maintain two distinct layers in water sheets ∼ 10 cm thick.
For subglacial macroporous water sheets, is dependent on the density and size of obstacles in the macroporous subglacial water sheet. We expect that the thickness of the sheet (H) will, at most, be the typical diameter of clastic obstacles (d) which maintain the ice–bed opening. Spherical objects of nearly similar size have a maximum three-dimensional packing density of (from the “Kepler conjecture”; Hales, 2005). Consequently, we expect that the maximum value of γ is likely to be near . This argument provides a reasonable upper bound on γ, though it should be kept in mind that nature does not necessarily arrange clasts in an optimal packing configuration and that clasts may in reality vary in size.
Using this range of values for the thinnest water sheet (H=1 cm), the highest drag for a wall with less than centimeter-scale roughness (Cd=0.01) and the highest water velocity (Uin=1 cm s−1), we calculate (from Eq. 12) a hard-bed intrusion length scale of m, giving an intrusion distance of 450 m for a densely obstructed water sheet or 700 m for an unobstructed water sheet (from the lower and upper bounds in Eq. 11, respectively). A more moderate range of parameters (H=5 cm, Cd=0.005, Uin=5 mm s−1) gives intrusion distances in the range of tens of kilometers. At the other end of the spectrum, a thick water sheet suggested by radar to exist at Siple Coast grounding zones (H=20 cm; MacGregor et al., 2011), with drag coefficients appropriate for large roughness or dense clasts (Cd=0.01) and low water velocity (Uin=1 mm s−1), would suggest an intrusion length scale easily spanning the entire ice sheet interior. Maintaining steady water sheet thickness in this range would require very high water pressures which may lead to turbulent flow, violating one of the base assumptions of our theory (Fr<1). We also would not expect such conditions to persist all the way into the ice sheet interior (where subglacial water discharge and water sheet thickness are both likely to be much lower, in addition to potential variations in bed slope). Therefore, we note this possibility but remain skeptical that our theory would provide a reliable prediction for seawater intrusion distance in this case. We should also note there have been indirect observations of very thin water sheets approx. millimeters thick (Engelhardt and Kamb, 1997) at ice–till interfaces in Antarctica, but we do not envision that maintaining distinct layering would be possible in such a circumstance.
This extreme range in potential intrusion distance reflects the fact that subglacial hydrological parameters are uncertain over orders of magnitude, resulting in many orders of magnitude uncertainty in the intrusion distance. Nonetheless, we conclude that on hard beds, seawater intrusion distances of at least hundreds of meters are possible even under very conservative assumptions, and intrusion distances of tens of kilometers are plausible. In Sect. 4, we explore the implications of intrusion-induced basal melt over such length scales upstream of marine ice sheet grounding lines.
3.2 Soft beds
When subglacial drainage occurs by Darcy flow through confined layers of till (i.e., a “soft bed”), there is a distinct (from the hard-bed case discussed above) horizontal length scale governing seawater intrusion, given by Eq. (16). The α parameter (analogous to g′) is fairly well known for freshwater–seawater systems to be 40. H is the thickness of the confined till layer, which can range anywhere from centimeters to tens of meters. Similarly, the inflow rate of freshwater through till layers (Uin) has been estimated to be on the order of 0.1–1 mm s−1 in ice streams where frictional heating at the ice–till interface generates high melt rates (Joughin et al., 2009). However, in regions of lower basal melt rates, we might expect inflow rates closer to 0.001–0.01 mm s−1 (Pattyn, 2010).
Though H and Uin are uncertain over 1–4 orders of magnitude, by far the most uncertain property of subglacial till relevant to seawater intrusion is K, the hydraulic conductivity of till. Laboratory measurement of K in till from formerly and currently glaciated field sites ranges from 10−12 to m s−1 (Freeze and Cherry, 1979; Cuffey and Paterson, 2010). Other methods have also been used to measure till hydraulic conductivity including in situ measurement of borehole fluids far upstream of the grounding line at West Antarctic ice streams (Engelhardt et al., 1990) and inference from measurements of grounding line migration in different West Antarctic ice streams (Warburton et al., 2020). These other approaches to constraining till properties give a range of 10−9 to 10−4 m s−1 for till hydraulic conductivity.
Even though the range on the till hydraulic conductivity, thickness and inflow discharge velocities is large, we can confidently conclude that, based on the scaling given in Eq. (16), seawater intrusion into flat, confined till layers will extend, at most, meters upstream of the grounding line. For most of the range of till properties discussed above, the length scale of seawater intrusion on flat beds will be negligible. This can be compared to the range of seawater intrusion into coastal aquifers of meters to kilometers. Hydraulic conductivity of consolidated subglacial till is well known to be low compared to a typical aquifer, which is largely the reason for such a small intrusion distance through flat soft beds. It is, however, important to note that previous field-based studies of till properties (particularly those reporting higher hydraulic conductivities) suggest that saturation of the subglacial till layer may prevent infiltration of basal melt below the ice–till interface, leading to sheet-like flow between ice and till, more similar to the hard-bed sheet-like hydrology discussed in Sects. 2.1 and 3.1 (though the thickness of such a supra-till water sheet is not entirely clear). Thus, as observations from West Antarctica suggest (MacGregor et al., 2011; Horgan et al., 2013), even in places where there is substantial evidence for the existence of thick subglacial till layers, seawater intrusion may still occur through shallow water flow between ice and till.
3.3 Sloping beds
The intrusion of seawater under grounded ice is ultimately driven by gravity. The salt wedge is the thickest near the ocean and thins towards the glacier interior (Fig. 1), producing a pressure gradient (proportional to the salt wedge thickness gradient, , in Eq. 5) that maintains the salt wedge against resistance from drag at walls, obstacles and the layer interface. On flat beds, the layer interface becomes flatter towards the interior (e.g., Fig. 2) such that eventually the drag exceeds gravitational driving, causing the termination of seawater intrusion. A bed that deepens into the glacier interior (θ>0 in Eq. 4 or Eq. 19) slopes in the same direction as the salt wedge interface and thus drives further intrusion over what would be expected on a flat bed. A sufficiently reverse-sloping bed will drive additional and potentially unbounded seawater intrusion distance (i.e., L→∞). Here we explore the criteria for the “critical bed slope” beyond which seawater intrusion extends unbounded for both hard and soft beds. We will also numerically solve for the intrusion distance over a range of bed slopes and other parameter values to determine generally when bed slope becomes an important driver of intrusion.
Following Wilson et al. (2020), we note that where the right-hand side of Eq. (5) becomes less than or equal to zero, the drag terms will never be sufficient to overcome the gravitational driving of the seawater intrusion, implying a seawater layer with constant or growing thickness (towards the glacier interior). This implies a “critical bed slope” criterion:
where η is a dimensionless number that is the minimum of the function , and C0 is the drag coefficient scale (which is set by the obstacle/wall drag, Cd). Analytical approximation of η generally yields unwieldy expressions in terms of , and γ. However, we can make two general observations about the scales of η before turning to the numerical solution. For an unobstructed water sheet (γ≪1), η is generally O(1). For a macroporous water sheet (γ∼ O(1)), η is generally O(γ)). Thus, for most realistic configurations of the subglacial water sheet (γ≲), we expect the critical bed slope to be .
To show how bed slope affects seawater intrusion distance generally, we numerically solve Eq. (5) and in Fig. 3 plot the seawater intrusion distance (L) over a wide range of bed slopes and Froude numbers. What we see largely reflects the inequality in Eq. (20). For realistically low flow rates of subglacial discharge (Fr0 ∼ O(0.1)) and drag coefficients (C0 ∼ ), seawater intrusion becomes unbounded when the bed has a reverse slope that is steeper than 1–. Additionally, prograde bed slopes that are common () effectively rule out seawater intrusion of more than 10 m. We generally see that, for subcritical water flow (Fr0<1) over reverse bed slopes, the critical bed slope is O(10−2) or flatter. These bed slopes are well within the range of reverse bed slopes found at many grounding lines at glacier termini in Greenland and Antarctica (Ross et al., 2012; Morlighem et al., 2017, 2020).
For a sloped subglacial porous aquifer, there is a similar limit at which the intrusion length becomes unbounded. This occurs when the slope of the seawater intrusion interface is less than the bed slope, thus meaning that they will never make contact (unless the bed slope or freshwater inflow changes). This occurs when the logarithm in Eq. (19) becomes unbounded, which is the case for
assuming θ≈tan (θ). is a dimensionless parameter that quantifies the relative importance of freshwater discharge and hydraulic conductivity of the till, which together determine the “resistivity” of the till to seawater intrusion. When this parameter is large, seawater will have a difficult time intruding into the till due to large freshwater discharge pushing back and/or low hydraulic conductivity preventing easy flow. At the upper bound of the range of observed till hydraulic conductivity – m s−1 and the lowest feasible freshwater flow rates in till from the interior of the catchment () – this “resistivity” is O(0.1), and the critical bed slope is in the range of the absolute steepest reverse bed slopes in either Greenland or Antarctica. Indeed, calculating seawater intrusion distance from Eq. (19) for a wide range of θ values (Fig. 4) shows that while there are some circumstances under which seawater intrusion into till aquifers might be non-negligible (i.e., greater than meters), such conditions are seemingly rarely attained. However, it is important to emphasize that if the hydraulic conductivity of till was observed to be any higher than 10−4 m s−1 or freshwater discharge was very low (perhaps under a slowly flowing glacier), seawater intrusion into till would become a distinctly non-negligible possibility on reverse-sloping beds.
Where water comes into contact with ice sheets, heat and salt are exchanged across the ice–water boundary layer (Jenkins, 1991), which may lead to dissolution or melting of ice. To estimate the effect of such intrusion-induced basal melt of grounded ice on the evolution of marine ice sheets, we incorporate a simple parameterization of intrusion melt into the Ice-sheet and Sea-level System Model (ISSM; Larour et al., 2012). The parameterization assumes that basal melt rates () decrease linearly from the grounding line upstream, reaching zero at a specified intrusion distance, L:
where is the melt rate at the grounding line (either parameterized or determined via coupling to an ocean model). In this parameterization in ISSM, x is the horizontal distance to the nearest grounding line. In an ice sheet model with explicitly simulated subglacial hydrology, x could be taken as a horizontal distance along hydropotential minima. We also specify that this parameterization is only applied upstream of the grounding line, where . Figure 5 shows an illustration of this parameterized melt rate as a function of distance from the terminus or grounding line.
Observations show that there are many places where a cold, fresh layer exists between floating ice and the warm, salty seawater. In such circumstances, basal melt still occurs through double-diffusive convection (Kimura et al., 2015; Begeman et al., 2018), though it is lower than it would be if a fully mixed and turbulent boundary layer existed at the ice–ocean interface (Rosevear et al., 2021). Where there is intrusion of a warm, saline layer under fresh, cold subglacial discharge beneath grounded ice, with subcritical water flow and limited entrainment (as may exist for seawater intrusion over hard beds), we argue that double-diffusive convection may occur and maintain distinct water layers while driving some ice melt that is dependent on properties of the saline layer, as has been shown in experiments and direct numerical simulations of the sub-ice boundary layers (Martin and Kauffman, 1977; Turner and Veronis, 2004). Correspondingly, in our first model configuration, we consider two sets of simulations. In the first set of simulations, grounding line melt rates are tens of meters per year, in line with previous benchmark simulations (described in the next section). In a second set of simulations, grounding line melt rates are meters per year, similar to those observed in a thin (i.e., tens of centimeters) double-diffusive staircase beneath the George VI ice shelf (Kimura et al., 2015). It is important to note two other possible outcomes from such a strongly stratified subglacial water layer: (1) where the sub-ice fresh layer is very thick (i.e., meters), double-diffusive convection drives even lower basal melt rates (Begeman et al., 2018), and (2) where seawater temperature is below the freshwater freezing point, there will be freeze-on at the ice base and not melting (Martin and Kauffman, 1977; Pedley et al., 1988). Rather, the goal of this section is to explore the consequences of intrusion-induced melt where the freshwater layer is sufficiently thin and the seawater layer is sufficiently warm, which currently occurs at some locations in Greenland and Antarctica and may occur at more locations in the future.
Intrusion-induced basal melt fluxes would need to be balanced by an additional influx of seawater to maintain steady-state seawater intrusion distance (and layer thicknesses). For intrusion-induced basal melt rates of the form of Eq. (22), the integrated melt flux would be . From the grounding line boundary condition for freshwater flow speed over a hard bed (Eq. 7), we can deduce the dimensional saline layer thickness at the grounding line: . Ocean currents can cause an incoming flux of , where Uo is the incoming current velocity from the ocean. For intrusion distances of hundreds to thousands of meters, subglacial water layers thicknesses of 1–10 cm, and freshwater discharge velocities of 0.1–1 cm s−1, Uo needs to be 0.1–10 cm s−1 to balance basal melt rates of the order of 1–100 m yr−1. This is well within the range of current speeds measured near the grounding lines of ice streams in Antarctica (Begeman et al., 2018). Thus, we consider these melt rates to be plausible purely from a consideration of the conditions necessary to maintain a steady state. However, we do note that such considerations do limit the possible range of steady-state intrusion distances and basal melt rates (i.e., intrusive melt rates greater than 100 m yr−1 over tens of kilometers would only be possible to maintain in a steady state with very rapid ocean inflow to compensate). Future work should explore the internal circulation of steady-state seawater intrusion layers in more detail.
Our goal in this section is to implement a simple parameterization of intrusive melt to estimate the effect of intrusion-driven melt under grounded ice on marine ice sheet evolution in two well-studied model configurations. However, we do recognize that the parameterized linear decrease in melt rates in Eq. (22) is clearly a simplification of the shape of the warm salt wedges shown in Fig. 2 and does not consider the potential complexities of sub-ice boundary layer processes. In reality, a more realistic fluid dynamic model including turbulence and heat fluxes would be needed to fully understand melting in seawater intrusions. In Sect. 5, we discuss the myriad ways that the representation of intrusion-driven melt can be made more realistic in numerical models.
In this study, we have adapted the “sub-element melt 2” (SEM2) parameterization described by Seroussi and Morlighem (2018) to incorporate melt from seawater intrusion on grounded ice by calculating the “level set” horizontal distance from the grounding line everywhere in the ISSM domain (x in Eq. 22). Intrusion distance L is specified as a model parameter, and sub-shelf melt rate is prescribed as a constant over floating ice (though it could equally well come from a more sophisticated model for floating melt rates). In this section, we will explore the effect of different intrusion distances on the transient evolution of a marine-terminating glacier in two benchmark model configurations: the MISMIP+ idealized bed topography (Asay-Davis et al., 2016) and future evolution of Thwaites Glacier in West Antarctica (based on a well-tested model configuration; Seroussi et al., 2017; Robel et al., 2019). All simulations in this section are conducted at 125 m horizontal resolution, where errors due to the precise form of the numerical implementation of basal melt near the grounding line are negligible, as shown in Seroussi and Morlighem (2018). A convergence study (not plotted) indicates that all results presented here using the intrusion melt parameterization are less than 3 % different from equivalent simulations conducted at 250 m horizontal resolution, confirming that these results are indeed converged in horizontal resolution.
4.1 ISSM intrusion melt: MISMIP bed topography
To demonstrate the effect of melt from seawater intrusion on ice sheet mass loss, we start by using a common benchmark of marine ice sheet models, the MISMIP+ model configuration (Asay-Davis et al., 2016; Cornford et al., 2020). In this idealized configuration, a marine-terminating glacier with a buttressing ice shelf begins in a steady state with the grounding line on a reverse-sloping bed, with no basal melt applied to the floating ice, and with constant snowfall accumulation rate in space. For simplicity, in this study we consider a variant on the “Ice1r” transient experiment, in which we apply a constant melt rate to all floating ice and then permit melt to occur with a linearly decreasing profile upstream of the grounding line according to our intrusion melt parameterization (Eq. 22 and Fig. 5). This constant-melt variant ensures that any dynamic feedback in response to melt forcing is due to ice sheet dynamics alone and not complicated feedbacks due to melt rate dependencies on ice sheet and bed geometry. This simplification will also allow us to straightforwardly compare the effect of adding intrusion melt over certain distances to cases without intrusion melt (in the discussion below).
In Fig. 6, we show the transient ice loss (in terms of volume above floatation, VAF) in the MISMIP+ configuration in response to various seawater intrusion distances (L in Eq. 22) and high prescribed floating melt rates ( m yr−1 in Eq. 22). The solid colored lines (in panel a) and colored circles (in panel b) show ice loss over 100 years. In this “high melt” scenario, we find a robust and significant effect of melt from seawater intrusion, with a 10 %–50 % increase in the rate of ice volume loss for seawater intrusion over hundreds of meters and a 50 %–105 % increase in loss rate for intrusion over kilometers. The direct increase in ice loss from basal melt in the intrusion region generally amounts to less than 5 % of the increased volume loss across experiments. Rather, it is the dynamic marine ice sheet response to melting upstream of the grounding line which is responsible for the significant increase in ice loss in the experiments with intrusion-induced basal melt. Intrusion melt turns grounded ice into floating ice, which is then subject to higher melt rates and induces a greater flux of ice from upstream. The effect of 500 m of seawater intrusion is equivalent to doubling the floating melt rate without intrusion (dashed line), and the effect of 1 km of seawater intrusion is equivalent to tripling the floating melt rate without intrusion (dotted line).
As noted previously, it may be difficult to maintain high basal melt rates under grounded ice in steady state. So, we also conduct comparable “low melt” simulations with an order of magnitude lower baseline basal melt rates ( m yr−1), which is closer to the melt rates expected in locations with double-diffusive convection driving melt (Kimura et al., 2015). In Fig. 7, we find a 5 %–20 % increase in the rate of ice volume loss for seawater intrusion over hundreds of meters and a 20 %–55 % increase in loss rate for intrusion over kilometers. The effect of 500–1000 m of seawater intrusion is equivalent to increasing the floating melt rate without intrusion by (dashed line), and the effect of 4 km of seawater intrusion is equivalent to doubling the floating melt rate without intrusion (dotted line). Thus, while these intrusion-induced increases in ice loss rate for the “low melt” scenario are less than in the “high melt” scenario, intrusion does still have a first-order effect on ice loss even for basal melt rates more typical of double-diffusive convection.
Though these simulations are instructive as to the large effect of intrusion-induced melt, we do emphasize that the exact sensitivity of ice loss to prescribed floating ice melt rate and intrusion distance is dependent on the particular model configuration (i.e., bed topography, extent of buttressing, initial glacier state). In the next section, we test another model configuration of interest to provide a non-idealized configuration as a point of comparison.
4.2 ISSM intrusion melt: Thwaites Glacier, West Antarctica
We also test effect of melt from seawater intrusion on ice loss from Thwaites Glacier (TG), a marine-terminating glacier in West Antarctica that is the focus of intense interest due to its recent acceleration and contribution to global sea level rise. The extent and intensity of submarine melting under the floating portions of TG have been studied extensively using field and remote sensing observations (e.g., Bevan et al., 2021; Wåhlin et al., 2021) and ocean modeling (e.g., Seroussi et al., 2017; Nakayama et al., 2019). In this section, we adapt the model configuration of ISSM from Seroussi et al. (2017), with a domain encompassing the TG catchment and a fine horizontal resolution. Submarine melt rates are again set to be constant at 60 m yr−1, which promotes a rapid retreat of the TG grounding line and is consistent with observations of Thwaites sub-shelf melt rates (Milillo et al., 2019). Surface mass balance is held constant in time and is based on regional simulations with RACMO2 averaged over the 1979–2010 time period (Lenaerts et al., 2012). The characteristics of subglacial hydrology of TG have been the focus of many recent studies, with evidence for a soft bed of variable strength and channelized hydrology near the grounding line (Schroeder et al., 2013, 2016) which sits on a steeply reverse-sloping bed (Morlighem et al., 2020) that, in places, exceeds the critical bed slope for unbounded seawater intrusion discussed in Sect. 3.3. These characteristics of the subglacial environment are generally favorable for the intrusion of seawater, even in particular channelized regions, which is consistent with the observations of localized high basal melt upstream of the TG grounding line by Milillo et al. (2019). However, observations of TG subglacial hydrology are still sufficiently uncertain that we simply explore the potential effect of intrusion-induced melt upstream of the grounding line for the evolution of TG (similarly to the MISMIP+ case).
In Fig. 8, we show the transient ice loss over 500-year simulations in which TG retreats completely through its catchment (in all simulations). Seawater intrusion distances of hundreds of meters accelerate the onset of the most rapid retreat by decades, and seawater intrusion over kilometers accelerates the onset of the most rapid retreat by 100–200 years. In the first 100 years of TG retreat (a common simulation period), hundreds of meters of seawater intrusion increase the rate of ice loss at TG by up to 20 %, while kilometers of seawater intrusion increase the rate of ice loss by 20 %–100 %. The effect of 2–4 km of seawater intrusion is equivalent to increasing the submarine melt rate by 50 % without intrusion (dashed line). Though there are some quantitative differences between TG and MISMIP+ simulations in the sensitivity of ice loss rate to seawater intrusion, it is generally the case that seawater intrusion of hundreds to thousands of meters (well within the range of feasibility as shown in Sect. 3) leads to substantially more ice loss from retreating marine-terminating glaciers. These results suggest that melt from seawater intrusion could be a very important process in driving ice loss from marine ice sheets. In the next section, we discuss the implications of these results for observing the ice–ocean interface and simulating intrusion in ice sheet models.
Future evolution of marine ice sheets is strongly sensitive to the potential for seawater intrusion under grounded ice. However, as we have shown, predictions of the horizontal extent of seawater intrusion remain uncertain over many orders of magnitude primarily due to the lack of constraints on subglacial hydrology near the ice–ocean interface. To make better predictions of how seawater intrusion may affect future marine ice sheet evolution will require more sophisticated numerical models of the ice–ocean interface, as well as more observational and experimental constraints on the properties of subglacial hydrology near the ice–ocean interface. In particular, future work can focus on answering three questions: (1) can we definitively confirm that seawater intrusion and melt are happening under grounded ice at some glaciers? (2) What are the properties of the subglacial hydrologic system under grounded ice near the ice–ocean interface? (3) Which processes play a role in the fluid dynamics of the freshwater–seawater interface?
5.1 “Effective” intrusion melt in large-scale ice sheet models
We have shown that the extent of potential seawater intrusion and melt under grounded ice is highly consequential for the future evolution of marine ice sheets. This is consistent with the conclusion of previous studies that have found that grounding line migration is most sensitive to basal melt at the grounding line (Parizek et al., 2013; Arthern and Williams, 2017; Reese et al., 2018; Goldberg et al., 2019). Seroussi and Morlighem (2018) further showed that, in some ice sheet models, basal melt intended for floating ice is “erroneously” applied to areas just upstream of the grounding line. Such misapplication of basal melt will lead to more rapid marine ice sheet retreat and ice loss (e.g., by up to 100 % compared to cases without such “errors”; Seroussi et al., 2019). Other studies of past and future ice sheet evolution due to ocean warming (e.g., Golledge et al., 2017) have used such a coarse-resolution ice sheet model in which ocean melt is applied to partially grounded elements, effectively including many kilometers of seawater intrusion in their simulations. Such models simulate the collapse of the West Antarctic Ice Sheet and portions of the East Antarctic Ice Sheet during the Pliocene warm period, explaining the much higher global sea level during this period that had been a challenging target for other ice sheet models (Dutton et al., 2015). Such models have also been used to simulate the future contribution of the Antarctic Ice Sheet to sea level rise (Golledge et al., 2019), falling near the high end of the range of future model projections (Seroussi et al., 2020). Indeed, ice sheet intercomparisons have demonstrated that marine ice sheet models which include such “effective” intrusion-induced basal melt on grounded ice are consistently more sensitive to ocean forcing than models in which no melt is applied to these elements (increasing the Antarctic contribution to sea level rise by 50 %–100 %; Seroussi et al., 2019). Other ice sheet dynamical mechanisms (particularly those controlling the rate of ice loss from the ice–ocean interface) have also been incorporated into models as a way of explaining the higher marine ice sheet sensitivity to climate change implied by Pliocene sea level (e.g., DeConto and Pollard, 2016). Seawater intrusion would appear to have a similar degree of support as these other sensitivity-boosting mechanisms both from fundamental theory and contemporary glaciological observations. The potential role of seawater intrusion in explaining past high sea levels and substantially increasing future projections of sea level rise suggests that this process deserves more attention through observational and modeling efforts to constrain and predict subglacial hydrology near the ice–ocean interface.
5.2 Models of seawater intrusion
In this study, we showed that intrusion distance is sensitively dependent on the bed slope and subglacial hydrology. A more sophisticated treatment of seawater intrusion in an ice sheet model would dynamically calculate L at all grounding lines or termini, depending on the local bed type (hard or soft), bed slope and state of subglacial hydrology. Unfortunately, the extent of soft and hard beds is still not well-known below the Greenland and Antarctic ice sheets, most marine ice sheet models do not dynamically simulate subglacial hydrology, and those that do differ widely in the assumptions and equations used to describe evolving drainage (de Fleurian et al., 2018).
The theory in this study also makes many assumptions in order to explore the generic importance and extent of seawater intrusion. One assumption is that the two laminar layers experience drag at their interface and with the surrounding solid boundaries (ice, bed, clasts or till grains). Another important assumption is that subglacial discharge is subcritical (i.e., Fr>1), which permits the existence of layered subglacial flow. Though we posit that seawater intrusion could produce basal melt via a process such as double-diffusive convection, we do not model the sub-ice boundary layer in detail here. Additionally, compositional heterogeneities and flow from the ocean could further complicate our simplified assumptions of a uniformly salty body of seawater at rest. Transient evolution of the salt wedge from tidal pumping (as in Walker et al., 2013; Sayag and Worster, 2013) or evolution of subglacial hydrology through melting by seawater may also complicate matters. Lateral heterogeneities in subglacial hydrology and other bed properties (e.g., channelization, bed slope) may lead to intrusive melt at some parts of the ice–ocean interface but not others. As we show, the sensitivity of ice sheet evolution to intrusion-induced melt is strongly sensitive to the extent of seawater intrusion and heat fluxes from the intrusion into grounded basal ice. The potentially substantial importance of seawater intrusion in determining the response of marine ice sheets to ocean warming necessitates detailed considerations of all the complexities listed above. Future studies should utilize high-resolution two- and three-dimensional numerical modeling of the fluid and thermal transition between subglacial water flow and the ocean circulation, potentially coupled to ice sheet models. This study should provide the motivation for future model developments in this direction.
5.3 Observations and experiments on seawater intrusion
Prior studies have found intriguing evidence of seawater intrusion in observations of grounding lines. MacGregor et al. (2011) observe no discernable change in the bed reflectivity in radar transects across the grounding lines at Whillans and Kamb ice streams, indicating a continuous water layer without a strong change in composition across this transition region (up to tens of kilometers upstream). Though they ultimately attribute these observations to other factors (leaching from subglacial sediments), they do consider the possibility that seawater intrusion may explain these observations (though ultimately dismiss this possibility due to the subglacial hydropotential gradient argument, discussed in Sect. 2.1). Studies by Horgan et al. (2013) and Drews et al. (2017) both observe subglacial channels across grounding lines in West and East Antarctica, and both argue that weak hydraulic potentials would suggest the potential for tidal pumping of seawater upstream of the grounding line as in an estuary. Such an observation is fully consistent with the independent observations of MacGregor et al. (2011), as well as the “estuary-like” salt wedge described here and previously for channels in Wilson et al. (2020), though the theory does not require tides for seawater to intrude past the grounding line. Milillo et al. (2019) observed subglacial melt from satellite-based interferometric synthetic aperture radar kilometers upstream of the grounding line along channelized pathways. Though catchment-spanning seawater intrusion is likely ruled out by the lack of seawater observed in boreholes drilled to the bed far upstream (hundreds of kilometers) of the grounding line (Tulaczyk et al., 2014), there are indications of high-salinity water far upstream of the grounding line in recent electromagnetic observations (Gustafson, 2020).
Interferometry, sounding radar and seismic methods provide complementary observations of the ice–bed interface (e.g., bed reflectivity, melt rate, acoustic impedance). More recent advances in radar technology and processing (e.g., ApRES, full-waveform inversion) and electromagnetic methods (Key and Siegfried, 2017) hold promise to better constrain the structure of subglacial water system (e.g., H), pore-water properties, the salinity of subglacial water, and basal melt rates at high resolution across the grounding line. Providing better constraints on the velocity of subglacial water flow from upstream (Uin) in either till or subglacial water sheets would likely require either direct access through borehole drilling near the grounding line (as in Engelhardt and Kamb, 1997, but closer to the grounding line) or model-based estimates of meltwater production upstream (Joughin et al., 2009; Pattyn, 2010; Carter and Fricker, 2012). Additionally, drag coefficients for a wider range of subglacial conditions and obstacle types can be better constrained with experimental techniques, though only limited efforts have been made in this direction (Prohaska, 2017). Ultimately, targeted efforts to measure seawater intrusion in the field and in laboratory experiments would provide important constraints on the theory developed in this study.
The interface of ice, ocean and the solid Earth plays host to a complex array of processes that drive much of the changes observed at marine ice sheets. This study considers one such process, the intrusion of seawater under grounded ice. We extend the theory of layered seawater intrusion under grounded ice developed for channels by Wilson et al. (2020) to a generalized theory for seawater intrusion under grounded ice in many possible subglacial hydrological systems. We find that there is a wide range of seawater intrusion distances predicted by this generalized theory, which includes cases with effectively no intrusion but also up to tens of kilometers of intrusion (and potentially even further) in macroporous water sheets (or channels) over impermeable beds. Seawater intrusion is generally negligible through microporous till layers, but under the right circumstances (steep reverse bed slopes, high hydraulic conductivity) till can support substantial intrusion, analogous to seawater intrusion in porous coastal aquifers. Critically, if seawater intrusion melts grounded ice, there can be a substantial increase in the rate of transient ice loss from marine ice sheets. This finding is in line with previous studies which have found strong sensitivity of the projected contribution of marine ice sheets to future sea level rise to the numerical treatment of basal melt just upstream of grounding lines in ice sheet models (Seroussi et al., 2019).
This study demonstrates that seawater intrusion under grounded ice is theoretically possible (even expected) and has substantial implications for projections of future sea level rise from marine ice sheet retreat under climate change. To determine whether seawater intrusion actually occurs in reality, targeted observational campaigns and experiments are needed to investigate the subglacial hydrology and melt rates near grounding lines and glacier termini. Additionally, the sensitivity of sea level projections to uncertainties in intrusion-induced melt should be tested in a deliberate fashion (rather than as a numerical artifact). High-resolution experimental and numerical approaches are also needed to understand the potential role of turbulence, tides and other important processes in either supporting or suppressing seawater intrusion. Though seawater intrusion has thus far been understudied as a potential driver of ice sheet changes, we hope that synthesizing theoretical and modeling approaches to intrusion will catalyze future efforts to better understand this elusive process.
In this study, we present various analytic approximations in certain limits for the seawater intrusion distance. The seawater intrusion distance is defined as the location where h(x)=1 governed by the ordinary differential equation
where and for which we have a single boundary condition . This problem is an initial value problem in which the unknown is the value of x at which a certain criterion is satisfied. Therefore, it is simple enough to integrate the above equation numerically, starting from x=0 and continuing to march towards lowers value of x until the criterion is satisfied.
However, there is one aspect of this problem which requires care in the marching scheme, namely that introduces a singularity in . Though there are sophisticated ways to handle such an issue analytically, to facilitate a straightforward numerical solution, we simply set where is a machine-precision perturbation to the Froude number and use a variable-resolution marching scheme, in which
where dx0 is a larger step size (we use 0.1). This approach will ensure this method accurately captures the boundary layer near x=0, where h changes rapidly and the step size needs to be quite small, while using a more computationally efficient large step size dx0 away from the boundary layer where h is changing more slowly.
All Python and MATLAB scripts used to produce the figures in this paper are freely available in a public GitHub repository: https://github.com/aarobel/IntrusionUnderIceSheets (last access: 2 February 2022) or as a persistent Zenodo repository: https://doi.org/10.5281/zenodo.5949310 (Robel, 2022). The results of Sects. 2–3 can also be reproduced interactively in the cloud using the binderized Jupyter notebook: https://mybinder.org/v2/gh/aarobel/IntrusionUnderIceSheets/master?filepath=SeawaterIntrusionUnderIce.ipynb (last access: 2 February 2022). The ISSM software package is publicly available for download from https://issm.jpl.nasa.gov/ (Larour et al., 2012). The subglacial intrusion melt parameterization will be included in future releases.
AAR and EW conceived the study based on EW's prior work. AAR and HS modified ISSM, and AAR conducted the ISSM simulations. AAR wrote the manuscript with contributions from all authors.
The contact author has declared that neither they nor their co-authors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Thanks to Winnie Chu, Adrian Jenkins, Matthew Siegfried, Samer Naif, Joyce Sim, Jeremy Bassis, Peter Washam, Chris Chungkei Lai and Vincent Verjans for helpful discussions during the completion of this work. The final product was improved by helpful reviews from Carolyn Begeman and an anonymous reviewer, in addition to editing from Nicolas Jourdain and Nanna Karlsson. Computing resources were provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology, Atlanta. We are thankful for PACE research scientist Fang (Cherry) Liu's assistance on HPC challenges.
Alexander A. Robel was supported by startup funding from the Georgia Institute of Technology and the University System of Georgia.
This paper was edited by Nicolas Jourdain and reviewed by Carolyn Begeman and one anonymous referee.
Adams, E. E., Sahoo, D., Liro, C. R., and Zhang, X.: Hydraulics of seawater purging in tunneled wastewater outfall, J. Hydraul. Eng., 120, 209–226, 1994. a
Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone, R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland, P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.: Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1), Geosci. Model Dev., 9, 2471–2497, https://doi.org/10.5194/gmd-9-2471-2016, 2016. a, b
Begeman, C. B., Tulaczyk, S. M., Marsh, O. J., Mikucki, J. A., Stanton, T. P., Hodson, T. O., Siegfried, M. R., Powell, R. D., Christianson, K., and King, M. A.: Ocean stratification and low melt rates at the Ross Ice Shelf grounding zone, J. Geophys. Res.-Oceans, 123, 7438–7452, 2018. a, b, c
Bevan, S. L., Luckman, A. J., Benn, D. I., Adusumilli, S., and Crawford, A.: Brief communication: Thwaites Glacier cavity evolution, The Cryosphere, 15, 3317–3328, https://doi.org/10.5194/tc-15-3317-2021, 2021. a
Clarke, G. K.: Subglacial processes, Annu. Rev. Earth Pl. Sc., 33, 247–276, 2005. a
Cornford, S. L., Seroussi, H., Asay-Davis, X. S., Gudmundsson, G. H., Arthern, R., Borstad, C., Christmann, J., Dias dos Santos, T., Feldmann, J., Goldberg, D., Hoffman, M. J., Humbert, A., Kleiner, T., Leguy, G., Lipscomb, W. H., Merino, N., Durand, G., Morlighem, M., Pollard, D., Rückamp, M., Williams, C. R., and Yu, H.: Results of the third Marine Ice Sheet Model Intercomparison Project (MISMIP+), The Cryosphere, 14, 2283–2301, https://doi.org/10.5194/tc-14-2283-2020, 2020. a
Cuffey, K. and Paterson, W.: The Physics of Glaciers, Pergamon, 3rd Edn., ISBN: 9780123694614 0123694612, 2010. a
DeConto, R. M. and Pollard, D.: Contribution of Antarctica to past and future sea-level rise, Nature, 531, 591–597, 2016. a
de Fleurian, B., Werder, M. A., Beyer, S., Brinkerhoff, D. J., Delaney, I. A. N., Dow, C.F., Downs, J., Gagliardini, O., Hoffman, M. J., Hooke, R. L., and Seguinot, J.: SHMIP The subglacial hydrology model intercomparison Project, J. Glaciol., 64, 897–916, 2018. a
Drews, R., Pattyn, F., Hewitt, I., Ng, F., Berger, S., Matsuoka, K., Helm, V., Bergeot, N., Favier, L., and Neckel, N.: Actively evolving subglacial conduits and eskers initiate ice shelf channels at an Antarctic grounding line, Nat. Commun., 8, 1–10, 2017. a, b
Dupuit, J.: Études théoriques et pratiques sur le mouvement des eaux dans les canaux découverts et à travers les terrains perméabls: avec des considérations relatives au régime des grandes eaux, au débouché à leur donner, et à la marche des alluvions dans les rivières à fond mobile, 2nd Edn., Dunod, 1863. a
Dutton, A., Carlson, A., Long, A., Milne, G., Clark, P., DeConto, R., Horton, B., Rahmstorf, S., and Raymo, M.: Sea-level rise due to polar ice-sheet mass loss during past warm periods, Science, 349, aaa4019, https://doi.org/10.1126/science.aaa4019, 2015. a
Engelhardt, H., Humphrey, N., Kamb, B., and Fahnestock, M.: Physical conditions at the base of a fast moving Antarctic ice stream, Science, 248, 57–59, 1990. a
Ezhova, E., Cenedese, C., and Brandt, L.: Dynamics of three-dimensional turbulent wall plumes and implications for estimates of submarine glacier melting, J. Phys. Oceanogr., 48, 1941–1950, 2018. a
Forchheimer, P.: Über die Ergiebigkeit von Brunnen-Anlagen und Sickerschlitzen [On the yield of well fields and drainage trenches], Zeitschrift des Architekten-und Ingenieur-Vereins zu Hannover, 32, 540–564, 1886. a
Freeze, R. A. and Cherry, J. A.: Groundwater, Englewood Cliffs, N.J. Prentice-Hall, ISBN: 0133653129 9780133653120, 604 pp., 1979. a
Ghyben, B. W.: Nota in verband met de voorgenomen putboring nabij, Amsterdam, The Hague, 21, 1888 (in Dutch). a
Golledge, N. R., Thomas, Z. A., Levy, R. H., Gasson, E. G. W., Naish, T. R., McKay, R. M., Kowalewski, D. E., and Fogwill, C. J.: Antarctic climate and ice-sheet configuration during the early Pliocene interglacial at 4.23 Ma, Clim. Past, 13, 959–975, https://doi.org/10.5194/cp-13-959-2017, 2017. a
Golledge, N. R., Keller, E. D., Gomez, N., Naughten, K. A., Bernales, J., Trusel, L. D., and Edwards, T. L.: Global environmental consequences of twenty-first-century ice-sheet melt, Nature, 566, 65–72, 2019. a
Gustafson, C. D.: Electromagnetic Investigations of Submarine and Subglacial Hydrogeologic Systems, PhD thesis, Columbia University, Palisades, NY, 2020. a
Hales, T. C.: A proof of the Kepler conjecture, Ann. Math., 162, 1065–1185, 2005. a
Herzberg, A.: Die wasserversorgung einiger Nordseebader, J. Gasbeleucht. Wasserversorg., 44, 842–844, 1901. a
Jenkins, A.: A one-dimensional model of ice shelf-ocean interaction, J. Geophys. Res.-Oceans, 96, 20671–20677, 1991. a
Joughin, I., Tulaczyk, S., Bamber, J. L., Blankenship, D., Holt, J. W., Scambos, T., and Vaughan, D. G.: Basal conditions for Pine Island and Thwaites Glaciers, West Antarctica, determined using satellite and airborne data, J. Glaciol., 55, 245–257, 2009. a, b, c
Key, K. and Siegfried, M. R.: The feasibility of imaging subglacial hydrology beneath ice streams with ground-based electromagnetics, J. Glaciol., 63, 755–771, 2017. a
Kunze, E.: A review of oceanic salt-fingering theory, Prog. Oceanogr., 56, 399–417, 2003. a
Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res.-Earth, 117, F01022, https://doi.org/10.1029/2011JF002140, 2012. a, b
Lenaerts, J., Van den Broeke, M., van de Berg, W. J., van Meijgaard, E., and Kuipers Munneke, P.: A new, high-resolution surface mass balance map of Antarctica (1979–2010) based on regional atmospheric climate modeling, Geophys. Res. Lett., 39, L04501, https://doi.org/10.1029/2011GL050713, 2012. a
Lu, C., Xin, P., Kong, J., Li, L., and Luo, J.: Analytical solutions of seawater intrusion in sloping confined and unconfined coastal aquifers, Water Resour. Res., 52, 6989–7004, 2016. a
MacCready, P. and Geyer, W. R.: Advances in estuarine physics, Annu. Rev. Mar. Sci., 2, 35–58, 2010. a
MacGregor, J. A., Anandakrishnan, S., Catania, G. A., and Winebrenner, D. P.: The grounding zone of the Ross Ice Shelf, West Antarctica, from ice-penetrating radar, J. Glaciol., 57, 917–928, 2011. a, b, c, d, e, f
Martin, S. and Kauffman, P.: An experimental and theoretical study of the turbulent and laminar convection generated under a horizontal ice sheet floating on warm salty water, J. Phys. Oceanogr., 7, 272–283, 1977. a, b
Milillo, P., Rignot, E., Rizzoli, P., Scheuchl, B., Mouginot, J., Bueso-Bello, J., and Prats-Iraola, P.: Heterogeneous retreat and ice melt of Thwaites Glacier, West Antarctica, Sci. Adv., 5, eaau3433, https://doi.org/10.1126/sciadv.aau3433, 2019. a, b, c, d
Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., and Fenty, I.: BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multibeam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, 11–51, 2017. a
Morlighem, M., Rignot, E., Binder, T., Blankenship, D., Drews, R., Eagles, G., Eisen, O., Ferraccioli, F., Forsberg, R., Fretwell, P., and Goel, V.: Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet, Nat. Geosci., 13, 132–137, 2020. a, b
Nakayama, Y., Manucharyan, G., Zhang, H., Dutrieux, P., Torres, H. S., Klein, P., Seroussi, H., Schodlok, M., Rignot, E., and Menemenlis, D.: Pathways of ocean heat towards Pine Island and Thwaites grounding lines, Sci. Rep., 9, 1–9, 2019. a
Nye, J.: Water flow in glaciers; jökulhlaups, tunnels and veins, J. Glaciol., 17, 181–207, 1976. a
Parizek, B., Christianson, K., Anandakrishnan, S., Alley, R., Walker, R., Edwards, R., Wolfe, D., Bertini, G., Rinehart, S., Bindschadler, R., and Nowicki, S. M. J.: Dynamic (in) stability of Thwaites Glacier, West Antarctica, J. Geophys. Res.-Earth, 118, 638–655, 2013. a
Pedley, M., Paren, J., and Potter, J.: Localized basal freezing within George VI Ice Shelf, Antarctica, J. Glaciol., 34, 71–77, 1988. a
Prohaska, Y.: Laboratory Experiment of Subglacial Drainage: Prototype Development, PhD thesis, ETH Zurich, 2017. a
Ralston, D. K., Geyer, W. R., and Lerczak, J. A.: Structure, variability, and salt flux in a strongly forced salt wedge estuary, J. Geophys. Res.-Oceans, 115, C06005, https://doi.org/10.1029/2009JC005806, 2010. a
Robel, A. A., Seroussi, H., and Roe, G. H.: Marine ice sheet instability amplifies and skews uncertainty in projections of future sea-level rise, P. Natl. Acad. Sci. USA, 116, 14887–14892, 2019. a
Rominger, J. T. and Nepf, H. M.: Flow adjustment and interior flow associated with a rectangular porous obstruction, J. Fluid Mech., 680, 636–659, 2011. a
Rosevear, M. G., Gayen, B., and Galton-Fenzi, B. K.: The role of double-diffusive convection in basal melting of Antarctic ice shelves, P. Natl. Acad. Sci. USA, 118, e2007541118, https://doi.org/10.1073/pnas.2007541118, 2021. a
Ross, N., Bingham, R. G., Corr, H. F., Ferraccioli, F., Jordan, T. A., Le Brocq, A., Rippin, D. M., Young, D., Blankenship, D. D., and Siegert, M. J.: Steep reverse bed slope at the grounding line of the Weddell Sea sector in West Antarctica, Nat. Geosci., 5, 393–396, 2012. a
Schroeder, D. M., Blankenship, D. D., and Young, D. A.: Evidence for a water system transition beneath Thwaites Glacier, West Antarctica, P. Natl. Acad. Sci. USA, 110, 12225–12228, 2013. a
Schroeder, D. M., Grima, C., and Blankenship, D. D.: Evidence for variable grounding-zone and shear-margin basal conditions across Thwaites Glacier, West Antarctica, Geophysics, 81, WA35–WA43, 2016. a
Seroussi, H., Nakayama, Y., Larour, E., Menemenlis, D., Morlighem, M., Rignot, E., and Khazendar, A.: Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation, Geophys. Res. Lett., 44, 6191–6199, https://doi.org/10.1002/2017GL072910, 2017. a, b, c
Seroussi, H., Nowicki, S., Simon, E., Abe-Ouchi, A., Albrecht, T., Brondex, J., Cornford, S., Dumas, C., Gillet-Chaulet, F., Goelzer, H., Golledge, N. R., Gregory, J. M., Greve, R., Hoffman, M. J., Humbert, A., Huybrechts, P., Kleiner, T., Larour, E., Leguy, G., Lipscomb, W. H., Lowry, D., Mengel, M., Morlighem, M., Pattyn, F., Payne, A. J., Pollard, D., Price, S. F., Quiquet, A., Reerink, T. J., Reese, R., Rodehacke, C. B., Schlegel, N.-J., Shepherd, A., Sun, S., Sutter, J., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., and Zhang, T.: initMIP-Antarctica: an ice sheet model initialization experiment of ISMIP6, The Cryosphere, 13, 1441–1471, https://doi.org/10.5194/tc-13-1441-2019, 2019. a, b, c, d
Seroussi, H., Nowicki, S., Payne, A. J., Goelzer, H., Lipscomb, W. H., Abe-Ouchi, A., Agosta, C., Albrecht, T., Asay-Davis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., Galton-Fenzi, B. K., Gladstone, R., Golledge, N. R., Gregory, J. M., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huybrechts, P., Jourdain, N. C., Kleiner, T., Larour, E., Leguy, G. R., Lowry, D. P., Little, C. M., Morlighem, M., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Reese, R., Schlegel, N.-J., Shepherd, A., Simon, E., Smith, R. S., Straneo, F., Sun, S., Trusel, L. D., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., Zhao, C., Zhang, T., and Zwinger, T.: ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070, https://doi.org/10.5194/tc-14-3033-2020, 2020. a
Shepherd, A., Ivins, E., Rignot, E., Smith, B., van den Broeke, M., Velicogna, I., Whitehouse, P., Briggs, K., Joughin, I., Krinner, G., and Nowicki, S.: Mass balance of the Antarctic Ice Sheet from 1992 to 2017, Nature, 556, 219–222, 2018. a
Straneo, F. and Heimbach, P.: North Atlantic warming and the retreat of Greenland's outlet glaciers, Nature, 504, 36–43, 2013. a
Tulaczyk, S., Mikucki, J. A., Siegfried, M. R., Priscu, J. C., Barcheck, C. G., Beem, L. H., Behar, A., Burnett, J., Christner, B. C., Fisher, A. T., and Fricker, H. A.: WISSARD at Subglacial Lake Whillans, West Antarctica: scientific operations and initial observations, Ann. Glaciol., 55, 51–58, 2014. a
Turner, J. and Veronis, G.: The influence of double-diffusive processes on the melting of ice in the Arctic Ocean: laboratory analogue experiments and their interpretation, J. Marine Syst., 45, 21–37, 2004. a
van der Boog, C. G., Dijkstra, H. A., Pietrzak, J. D., and Katsman, C. A.: Double-diffusive mixing makes a small contribution to the global ocean circulation, Commun. Earth Environ., 2, 1–9, 2021. a
Wåhlin, A., Graham, A., Hogan, K., Queste, B., Boehme, L., Larter, R., Pettit, E., Wellner, J., and Heywood, K.: Pathways and modification of warm water flowing beneath Thwaites Ice Shelf, West Antarctica, Sci. Adv., 7, eabd7254, https://doi.org/10.1126/sciadv.abd7254, 2021. a
Walker, R. T., Parizek, B. R., Alley, R. B., Anandakrishnan, S., Riverman, K. L., and Christianson, K.: Ice-shelf tidal flexure and subglacial pressure variations, Earth Planet. Sc. Lett., 361, 422–428, 2013. a, b
Werner, A. D., Bakker, M., Post, V. E., Vandenbohede, A., Lu, C., Ataie-Ashtiani, B., Simmons, C. T., and Barry, D. A.: Seawater intrusion processes, investigation and management: recent advances and future challenges, Adv. Water Resour., 51, 3–26, 2013. a, b
Wilson, E. A., Wells, A. J., Hewitt, I. J., and Cenedese, C.: The dynamics of a subglacial salt wedge, J. Fluid Mech., 895, A20, https://doi.org/10.1017/jfm.2020.308, 2020. a, b, c, d, e, f, g, h, i, j, k, l