Groundwater dynamics beneath a marine ice sheet

Sedimentary basins beneath many Antarctic ice streams host substantial volumes of groundwater, which can be exchanged with a “shallow” subglacial hydrological system of till and channelised water. This exchange contributes substantially to basal water budgets, which in turn modulate the flow of ice streams. The geometry of these sedimentary basins is known to be complex, and the groundwater therein has been observed to vary in salinity due to historic seawater intrusion. However, little is known about the hydraulic properties of subglacial sedimentary basins and the factors controlling groundwater exfiltration and infiltration. We develop a mathematical model for two-dimensional groundwater flow beneath a marine-terminating ice stream on geological timescales, taking into account the effect of seawater intrusion. We find that seawater may become “trapped” in subglacial sedimentary basins, through cycles of grounding line advance and retreat or through “pockets” arising from basin geometry. In addition, we estimate the sedimentary basin permeability, which reproduces field observations of groundwater salinity profiles from beneath Whillans Ice Stream in West Antarctica. Exchange of groundwater with the shallow hydrological system is primarily controlled by basin geometry, with groundwater being exfiltrated where the basin becomes shallower and re-infiltrating where it becomes deeper. Our results highlight the potential importance of groundwater flow in sedimentary basins for the subglacial hydrology of ice streams.

1 Introduction

In order to understand the stability of the West Antarctic Ice Sheet and its potential contribution to future sea level rise, it is critical to understand the dynamics of marine-terminating glaciers (Seroussi et al., 2020). The flow of such glaciers is strongly modulated by a shallow hydrological system at the base of the ice, which includes pressurised water in films, channels, and lakes, along with a water-saturated layer of deformable sediments (Tulaczyk et al., 2000; Bougamont et al., 2011; Gustafson et al., 2022). In addition, the fast-flowing glaciers known as ice streams, which drain much of West Antarctica, typically overlie large sedimentary basins, which host substantial volumes of groundwater. The exchange of this deeper groundwater with the shallow hydrological system has emerged as an important potential contributor to the basal conditions, which dominate ice stream flow (Christoffersen et al., 2014; Siegert et al., 2018).

Geophysical data can provide constraints on the geometry and properties of Antarctic sedimentary basins (Aitken et al., 2023; Li and Aitken, 2024). Airborne magnetic investigations of the sedimentary basin beneath the Ross Ice Shelf have revealed a complex geometry, consisting of peaks, troughs, and faults (Tankersley et al., 2022b). Moreover, magnetotelluric (MT) measurements of electrical conductivity have shown an increase in salinity with depth in these aquifers, which has been attributed to seawater intrusion during a Holocene glacial minimum (Gustafson et al., 2022).

Such evidence of past seawater intrusions is useful for constraining the grounding line history of the West Antarctic Ice Sheet and thus gaining insight into its future stability (Kingslake et al., 2018; Venturelli et al., 2020, 2023). Previous studies have used evidence of historic seawater intrusion to infer the effects of continental ice sheets on groundwater systems during glaciations (Aquilina et al., 2015).

Several recent studies have considered the mathematical modelling of groundwater in subglacial sedimentary basins. These studies have focused on exchange of groundwater with the shallow hydrological system, driven by the poro-elastic response of sedimentary basins to ice overpressure changes on timescales of decades to centuries (Gooch et al., 2016; Li et al., 2022; Robel et al., 2023). However, these models mostly focus on vertical flow (although Li et al. (2022) consider two-dimensional flow) and do not consider the effect of salinity differences or sedimentary basin geometry.

In contrast, the mathematical theory of seawater intrusion in coastal aquifers is of longstanding interest due to its implications for groundwater management, beginning with the investigations of Badon-Ghyben and Herzberg (1896). A rich array of mathematical and numerical techniques have emerged to investigate the effects of overextraction and, more recently, sea level rise on coastal groundwater flow (Bear, 2013; Werner et al., 2013; Ketabchi et al., 2016; Mondal et al., 2019; Richardson et al., 2024). However, such models are yet to be applied to a subglacial context.

In this paper, we consider a mathematical model for groundwater flow in a sedimentary basin beneath a marine ice sheet on geological timescales, accounting for seawater intrusion, ice sheet growth and retreat, and basin geometry. We show that grounding line motion and basin geometry may give rise to “saltwater trapping”, where seawater remains a resident in portions of the aquifer that would be saturated with freshwater in the simplest steady state. The phenomenon of saltwater trapping depends heavily on hydraulic properties of the aquifer, which are summarised in a dimensionless parameter K, along with past grounding line motion.

We then apply this model to the sedimentary basin beneath the Ross Ice Shelf in West Antarctica, enabling us to predict how groundwater is exchanged with the shallow hydrological system. In addition, comparing modelled saltwater trapping to field observations enables us to predict appropriate values of the parameter K. Our results evidence the importance of accounting for horizontal flow and saltwater intrusion when modelling groundwater dynamics in a subglacial sedimentary basin and the wider implications of these sedimentary basins for subglacial hydrology.

2 Model design

2.1 Derivation of the governing equation

We investigate seawater intrusion beneath a marine ice sheet in a vertical cross-section of a permeable aquifer, given by $\mathrm{0}<x<x_{\mathrm{g}}\left(t\right)$, $b\left(x\right)<z<S\left(x\right)$. The setup is depicted in Fig. 1. Here, x_g represents the grounding line, where the grounded ice sheet transitions to a floating ice shelf, and x=0 represents the water divide within the aquifer. For x>x_g, the upper surface of the aquifer is therefore exposed to the ocean, which rapidly saturates this region of the aquifer with seawater. In addition, the floating ice shelf necessarily imposes the same overpressure as the seawater it displaces. We take the present-day sea level as the reference z=0.

Figure 1 Schematic of model. The aquifer is bounded by z=b(x) below and z=S(x) above and has thickness $H\left(x\right)=S\left(x\right)-b\left(x\right)$. The saltwater region occupies $b\left(x\right)<z<s\left(x\right)$ and has thickness $h\left(x\right)=s\left(x\right)-b\left(x\right)$, and the freshwater region occupies $s\left(x\right)<z<S\left(x\right)$ and has thickness H(x)−h(x). The “shallow” subglacial hydrological system is depicted as a layer of water between the ice and the aquifer at z=S(x). The grounding line x=x_g(t) is marked by a vertical dashed line.

The surfaces z=b(x) and z=S(x) represent the basement and upper surface of the aquifer, which therefore has a thickness $H\left(x\right)=S\left(x\right)-b\left(x\right)$. We assume that the aquifer is homogeneous, with porosity ϕ and permeability k, although a generalisation of the model to relax these assumptions is straightforward (as discussed in Appendix B). For x<x_g(t), the upper surface z=S(x) is overlain by an ice sheet above a lubricating layer of freshwater. The aquifer is assumed to be saturated at all times but exchanges water with this lubricating layer at a rate q_E(x,t), where q_E>0 corresponds to exfiltration of water from the sedimentary basin. This layer represents an idealised basal hydrological system. We treat this layer in the model as having zero thickness, on the basis that the corresponding “shallow” hydrological system has a depth on the order of metres to tens of metres, whereas the aquifer is hundreds to thousands of metres deep (Gustafson et al., 2022).

We also assume that the aquifer is rigid. This assumption is valid provided that any dynamic deformation of the sediments is small. For a rigid aquifer, q_E represents the rate of exfiltration driven by horizontal pressure gradients emerging from topography, ice thickness, and seawater intrusion on long timescales. However, in reality the poro-elastic response of the sediments produces an additional exfiltration over shorter timescales, which has been considered in previous models by Li et al. (2022) and Robel et al. (2023). These models find that this exfiltration may reach tens to hundreds of millimetres per year for a highly permeable aquifer. This distinct exfiltration, which we do not consider, requires further modelling to combine with the long-timescale exfiltration q_E(x,t). This assumption of rigidity does not extend to the deformable till comprising the upper few metres of the sedimentary basin, which we regard as part of the basal hydrological system.

We use the sharp-interface approximation, a frequent assumption in saltwater intrusion problems, which asserts that mixing between freshwater and saltwater through molecular diffusion and hydrodynamic dispersion is negligible (Bear, 2013; Mondal et al., 2019). Since the inversion results of Gustafson et al. (2022) suggest a smooth variation in salinity through the aquifer rather than a distinct sharp interface, we include in Sect. 6 a discussion of how a similar model could account for these mixing processes, although it is possible that this smoothing could be in part due to the inversion method used.

The sharp-interface approximation is warranted provided that the Péclet numbers associated with molecular diffusion and hydrodynamic dispersion are large (Bear, 2013; Dentz et al., 2006; Koussis and Mazi, 2018). For a molecular diffusivity $D={\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹ and a hydrodynamic dispersivity length $A={\mathrm{10}}^{-\mathrm{2}}$ m, which are appropriate for salt in groundwater in sedimentary rocks (Aquilina et al., 2015), these Péclet numbers are respectively

$$\begin{array}{}\text{(1)}& {\mathit{Pe}}_{\text{diff}}={\displaystyle \frac{[z{]}^{\mathrm{2}}}{D\left[t\right]}}\approx \mathrm{300},{\mathit{Pe}}_{\text{disp}}={\displaystyle \frac{[z{]}^{\mathrm{2}}}{A\left[x\right]}}\approx \mathrm{200},\end{array}$$

suggesting that the sharp-interface assumption is reasonable here. Although we assume a sharp interface between the two regions, the fluids are miscible, meaning that there is no surface tension between them. We can therefore ignore the possibility of residual trapping at the pore scale (Bear, 2013).

The aquifer therefore consists of a freshwater-saturated region $s(x,t)<z<S\left(x\right)$ and a saltwater-saturated region $b\left(x\right)<z<s(x,t)$, where $z=s(x,t)$ is the position of the “sharp” freshwater–saltwater interface. Freshwater and saltwater are physically distinguished by their densities ρ_f and ρ_s. We use $h(x,t)=s(x,t)-b\left(x\right)$ to denote the thickness of the saltwater-saturated region so that the thickness of the freshwater-saturated region is $H\left(x\right)-h(x,t)$. The following model is solved for h(x,t), along with the pore pressure $p(x,z,t)$ and the Darcy flux $\left(u\right(x,z,t),w(x,z,t\left)\right)$ of water in the aquifer.

The model consists of mass conservation for the pore water,

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial w}{\partial z}}=\mathrm{0},\end{array}$$

along with the horizontal component of Darcy's law:

$$\begin{array}{}\text{(3)}& u=-{\displaystyle \frac{k}{\mathit{\mu }}}{\displaystyle \frac{\partial p}{\partial x}}.\end{array}$$

Here k is the permeability of the porous medium, assumed constant for simplicity, and μ is the fluid viscosity, which we assume is the same for fresh and salty water. These equations are combined with the Dupuit approximation, namely that the pressure is hydrostatic:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial p}{\partial z}}=-{\mathit{\rho }}_{j}g,j\in \mathit{\{}s,f\mathit{\}}.\end{array}$$

This approximation is standard in models of groundwater flow and is justified provided that the aspect ratio of the aquifer is small; i.e. the aquifer is long and thin, which is true for our purposes (see Tables 1 and 2) (Gustafson et al., 2022). Equation (4) can equivalently be found as the vertical component of Darcy's law in a shallow limit.

We next prescribe boundary conditions for the model. The basement z=b(x) is taken to be rigid and impermeable:

$$\begin{array}{}\text{(5)}& u{\displaystyle \frac{\partial b}{\partial x}}-w=\mathrm{0}\text{on}z=b\left(x\right).\end{array}$$

In addition, the sharp interface is taken to be kinematic:

$$\begin{array}{}\text{(6)}& \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial t}}+u_{\pm }{\displaystyle \frac{\partial s}{\partial x}}-w_{\pm }=\mathrm{0}\text{on}z=s(x,t),\end{array}$$

where the subscripts ± denote the Darcy flux components on either side of the interface, which are not continuous in general. However, the pressure $p(x,z,t)$ must be continuous. Equation (6) ensures that the velocity of the interface equals the Darcy velocity of the fluid on either side so that mass is conserved across the interface.

On the upper surface, we prescribe that the pressure of the pore fluid is continuous with that of the lubricating layer and moreover assume that this pressure is equal to the hydrostatic ice overburden pressure:

$$\begin{array}{}\text{(7)}& p=p_S(x,t)={\mathit{\rho }}_{\mathrm{i}}g{H}_{\mathrm{i}}(x,t)\text{on}z=S\left(x\right),\end{array}$$

where ρ_i and H_i(x,t) are the density and thickness of ice respectively. This boundary condition is standard in models of subglacial aquifers and is equivalent to stipulating that the “effective pressure”, defined as $p_{\mathrm{e}}={\mathit{\rho }}_{\mathrm{i}}g{H}_{\mathrm{i}}-p$, is zero (Lemieux et al., 2008; Gooch et al., 2016; Li et al., 2022; Robel et al., 2023).

The Darcy flux of water exfiltrated by the aquifer is then given by

$$\begin{array}{}\text{(8)}& q_{\mathrm{E}}={\left.\left(w-u{\displaystyle \frac{\partial S}{\partial x}}\right)\right|}_{z=S\left(x\right)},\end{array}$$

where q_E<0 corresponds with the aquifer being recharged by basal water from the lubricating layer.

In addition, horizontal boundary conditions are required. At the water divide x=0, we impose that there is no horizontal flux of water:

$$\begin{array}{}\text{(9)}& Hu=\mathrm{0}\text{at}x=\mathrm{0}.\end{array}$$

At the grounding line x=x_g(t), the weight of the ice sheet is that of the displaced ocean, and the aquifer is saturated with seawater. The pressure is therefore hydrostatic:

$$\begin{array}{}\text{(10)}& p=-{\mathit{\rho }}_{\mathrm{s}}gz\text{at}x=x_{\mathrm{g}},\end{array}$$

taking the atmospheric pressure as the reference pressure. Note that by definition

$$\begin{array}{}\text{(11)}& p_S={\mathit{\rho }}_{\mathrm{i}}g{H}_{\mathrm{i}}=-{\mathit{\rho }}_{\mathrm{s}}gS\text{at}x=x_{\mathrm{g}},\end{array}$$

since this is where the ice achieves flotation.

We now seek to reduce the above model to an evolution equation for h(x,t). Firstly, integrating Eq. (4) subject to the condition (7) gives the pressure as

$$\begin{array}{}\text{(12)}& p=\left\{\begin{array}{l}{p}_S+{\mathit{\rho }}_{\mathrm{f}}g(S-z)\text{for}s<z<S,\\ p_S+{\mathit{\rho }}_{\mathrm{f}}g\left(S-z+\mathit{\delta }(s-z)\right)\text{for}b<z<s,\end{array}\right.\end{array}$$

where

$$\begin{array}{}\text{(13)}& \mathit{\delta }={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{f}}}{{\mathit{\rho }}_{\mathrm{f}}}}\approx \mathrm{0.025}\end{array}$$

is the scaled density difference between fresh and saltwater. Darcy's law (Eq. 3) then gives the horizontal Darcy flux:

$$\begin{array}{}\text{(14)}& u=\left\{\begin{array}{l}-\frac{k{\mathit{\rho }}_{\mathrm{f}}g}{\mathit{\mu }}\frac{\partial }{\partial x}\left(\frac{p_S}{{\mathit{\rho }}_{\mathrm{f}}g}+S\right)\text{for}s<z<S,\\ -\frac{k{\mathit{\rho }}_{\mathrm{f}}g}{\mathit{\mu }}\frac{\partial }{\partial x}\left(\frac{p_S}{{\mathit{\rho }}_{\mathrm{f}}g}+S+\mathit{\delta }s\right)\text{for}b<z<s.\end{array}\right.\end{array}$$

Note that in order to be consistent with the condition (9), we need that

$$\begin{array}{}\text{(15)}& H{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{p_S}{{\mathit{\rho }}_{\mathrm{f}}g}}+S\right)=\mathrm{0}\text{at}x=\mathrm{0}.\end{array}$$

We now take the mass conservation statement Eq. (2) and integrate over the saltwater region $b\left(x\right)<z<s(x,t)$, subject to the kinematic conditions (5) and (6), yielding a statement of depth-integrated mass conservation:

$$\begin{array}{}\text{(16)}& \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\underset{b}{\overset{s}{\int }}u\mathrm{d}z=\mathrm{0}.\end{array}$$

Substituting Eq. (14) into Eq. (16) gives the desired equation for h(x,t):

$$\begin{array}{}\text{(17)}& \mathit{\varphi }{\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{k{\mathit{\rho }}_{\mathrm{f}}g}{\mathit{\mu }}}{\displaystyle \frac{\partial }{\partial x}}\left[h{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{p_S}{{\mathit{\rho }}_{\mathrm{f}}g}}+S+\mathit{\delta }s\right)\right],\end{array}$$

recalling that $s=b+h$. This equation requires two boundary conditions in x. The no-flux condition (9) implies, in addition to Eq. (15), that we need

$$\begin{array}{}\text{(18)}& h{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{p_S}{{\mathit{\rho }}_{\mathrm{f}}g}}+S+\mathit{\delta }s\right)=\mathrm{0}\text{at}x=\mathrm{0}.\end{array}$$

In addition, the condition (10) gives

$$\begin{array}{}\text{(19)}& h=H\text{at}x=x_{\mathrm{g}}.\end{array}$$

A similar depth integration of Eq. (2) over the entirety of the aquifer $b<z<S$, using the definition in Eq. (8), provides an equation for the exfiltration q_E:

$$\begin{array}{}\text{(20)}& q_{\mathrm{E}}={\displaystyle \frac{k{\mathit{\rho }}_{\mathrm{f}}g}{\mathit{\mu }}}{\displaystyle \frac{\partial }{\partial x}}\left[H{\displaystyle \frac{\partial }{\partial x}}\left({\displaystyle \frac{p_S}{{\mathit{\rho }}_{\mathrm{f}}g}}+S\right)+\mathit{\delta }h{\displaystyle \frac{\partial s}{\partial x}}\right].\end{array}$$

The first term, containing H, depends on the basement geometry and ice overpressure, and therefore it depends on time only through the instantaneous position of the ice sheet. The second, containing h, accounts for the effect of saltwater intrusion on the freshwater recharge or discharge.

2.2 Non-dimensionalisation

We assume that the present-day grounding line position, measured relative to the onset of the sedimentary basin (Peters et al., 2006), and the approximate subglacial aquifer thickness at the measurement sites by Gustafson et al. (2022) provide suitable horizontal and vertical length scales [x] and [z]. Since both the aquifer and the ice sheet are hundreds to thousands of metres deep, we may use the latter for both the aquifer and ice sheet. We also assume that grounding line motion imposes a natural timescale [t]. We then scale

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& x,x_{\mathrm{g}}\sim \left[x\right],z,h,H,b,s,S\sim \left[z\right],t\sim \left[t\right],\\ & p_S\sim {\mathit{\rho }}_{\mathrm{f}}g\left[z\right],q_{\mathrm{E}}\sim \left[z\right]/\left[t\right].\end{array}\end{array}$$

For example, we let $x=\left[x\right]\widehat{x}$, where $\widehat{x}$ is non-dimensional, and do likewise for each of the variables listed with a corresponding scale in Eq. (21) (i.e. $z=\left[z\right]\widehat{z}$, $h=\left[z\right]\widehat{h}$, etc.). We then work with these non-dimensional variables but drop all hats ($\widehat{}$) from the notation. From here onwards, variables are dimensionless unless otherwise stated, but dimensional values are used in figures, based on the scalings provided in Tables 1 and 2 unless otherwise stated.

Under these scalings, Eq. (17) becomes

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{\partial h}{\partial t}}=K{\displaystyle \frac{\partial }{\partial x}}\left[h{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S+\mathit{\delta }s\right)\right],\end{array}$$

where

$$\begin{array}{}\text{(23)}& K={\displaystyle \frac{k{\mathit{\rho }}_{\mathrm{f}}g\left[z\right]\left[t\right]}{\mathit{\varphi }\mathit{\mu }[x{]}^{\mathrm{2}}}}\end{array}$$

is the dimensionless hydraulic conductivity of the groundwater flow. A larger K corresponds to faster groundwater flow and vice versa. The conditions (18) and (19) become

$$\begin{array}{}\text{(24)}& h{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S+\mathit{\delta }s\right)=\mathrm{0}\text{at}x=\mathrm{0},\end{array}$$

$$\begin{array}{}\text{(25)}& h=H\text{at}x=x_{\mathrm{g}}.\end{array}$$

We can also obtain from Eq. (20) a dimensionless expression for the exfiltration flux:

$$\begin{array}{}\text{(26)}& q_{\mathrm{E}}=K{\displaystyle \frac{\partial }{\partial x}}\left[H{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S\right)+\mathit{\delta }h{\displaystyle \frac{\partial s}{\partial x}}\right].\end{array}$$

Values of the scalings in Eq. (21) and the corresponding values of dimensionless parameters are provided in Tables 1 and 2 respectively.

van der Wel et al. (2013); Gustafson et al. (2022); Li et al. (2022)Peters et al. (2006); Gustafson et al. (2022); Li et al. (2022)Bear (2013)Gooch et al. (2016); Li et al. (2022); Robel et al. (2023)Bear (2013)Bear (2013)Gooch et al. (2016); Michaud et al. (2017); Gustafson et al. (2022)Pollard and DeConto (2009); Willeit et al. (2019)Cuffey and Paterson (2010)Schoof (2007)Cuffey and Paterson (2010)

Table 1 Dimensional scalings used in numerical solutions.

Table 2 Values of dimensionless parameters used in numerical solutions.

In deriving Eq. (22), we have assumed that the top of the saltwater region $z=s(x,t)$ is a kinematic interface. However, it is possible for the aquifer to become saturated with saltwater throughout its depth so that $s(x,t)=S\left(x\right)$ (i.e. h=H). If this is the case, then saltwater, rather than freshwater, is discharged across the upper surface at a rate of

$$\begin{array}{}\text{(27)}& q_{\mathrm{E}}=K{\displaystyle \frac{\partial }{\partial x}}\left[H{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S+\mathit{\delta }{\displaystyle \frac{\partial S}{\partial x}}\right)\right].\end{array}$$

The state in which s=S can be maintained as long as the outflux given by Eq. (27) is non-negative; otherwise, freshwater influx will create a freshwater layer.

To allow for this case, we replace Eq. (22) by the following complementarity problem:

$$\begin{array}{}\text{(28)}& \begin{array}{r}{\displaystyle }\left(H-h\right)\left(-{\displaystyle \frac{\partial h}{\partial t}}+K{\displaystyle \frac{\partial }{\partial x}}\left[h{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S+\mathit{\delta }s\right)\right]\right)=\mathrm{0},\\ {\displaystyle }-{\displaystyle \frac{\partial h}{\partial t}}+K{\displaystyle \frac{\partial }{\partial x}}\left[h{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S+\mathit{\delta }s\right)\right]\ge \mathrm{0},\\ {\displaystyle }H-h\ge \mathrm{0}.\end{array}\end{array}$$

This formulation establishes that either Eq. (22) applies and h≤H or h=H and saltwater is discharged according to (27). We describe the numerical method used to solve this problem in Appendix A.

In Sect. 3 we consider the solution of Eq. (28) in a uniform aquifer geometry, firstly looking at steady-state solutions, before moving to a periodic time-dependent solution, which introduces the possibility of saltwater trapping via grounding line movement. In Sect. 4, we generalise to non-uniform aquifer geometries, which enables further potential saltwater trapping via so-called pockets. Finally in Sect. 5, we use Eq. (28) to model groundwater flow in the sedimentary basin beneath the Ross Ice Shelf.

2.3 Ice sheet forcing

The ice overpressure p_S is related to the ice thickness H_i by Eq. (7). To model the ice sheet, we use the shallow-shelf approximation (MacAyeal, 1989; Schoof, 2007), assuming a Weertman sliding law with exponent $\mathrm{1}/\mathrm{3}$ and that extensional stresses in the ice are negligible. We also assume that the accumulation is spatially uniform and that the dynamics of the ice sheet are quasi-steady for a given x_g. This leads to equations of mass and momentum conservation:

$$\begin{array}{}\text{(29)}& H_{\mathrm{i}}{u}_{\mathrm{i}}=ax,\mathit{\beta }{u}_{\mathrm{i}}^{\frac{\mathrm{1}}{\mathrm{3}}}={\mathit{\rho }}_{\mathrm{i}}g{H}_{\mathrm{i}}\left|{\displaystyle \frac{\partial }{\partial x}}\left(H_{\mathrm{i}}+S\right)\right|,\end{array}$$

where u_i is the (vertically uniform) ice velocity, a is the constant accumulation rate (measured in m s⁻¹), and β is a constant in the sliding law (measured in Pa m${}^{-\mathrm{1}/\mathrm{3}}$ s${}^{\mathrm{1}/\mathrm{3}}$). The latter equation represents a balance between the driving stress and the basal friction, which has a power-law dependence on the velocity under the Weertman sliding law. These together lead to the dimensional equation

$$\begin{array}{}\text{(30)}& {\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}g}{\mathit{\beta }}}\right)}^{\mathrm{3}}{H}_{\mathrm{i}}^{\mathrm{4}}{\left|{\displaystyle \frac{\partial }{\partial x}}\left(H_{\mathrm{i}}+S\right)\right|}^{\mathrm{3}}=ax.\end{array}$$

In dimensionless variables, after scaling H_i with [z], Eqs. (7) and (30) become

$$\begin{array}{}\text{(31)}& p_S=r_{\mathrm{i}}{H}_{\mathrm{i}},H_{\mathrm{i}}^{\mathrm{4}}{\left|{\displaystyle \frac{\partial }{\partial x}}\left(H_{\mathrm{i}}+S\right)\right|}^{\mathrm{3}}=\mathit{\alpha }x,\end{array}$$

where

$$\begin{array}{}\text{(32)}& r_{\mathrm{i}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{f}}}}\approx \mathrm{0.917},\mathit{\alpha }={\displaystyle \frac{a{\mathit{\beta }}^{\mathrm{3}}[x{]}^{\mathrm{4}}}{\left({\mathit{\rho }}_{\mathrm{i}}g{)}^{\mathrm{3}}\right[z{]}^{\mathrm{7}}}}.\end{array}$$

Equation (31) is consistent with the condition (15) provided that $\partial S/\partial x\left(\mathrm{0}\right)=\mathrm{0}$, and is solved by numerical integration subject to the condition (11). In dimensionless terms, Eq. (11) becomes

$$\begin{array}{}\text{(33)}& H_{\mathrm{i}}=-{\displaystyle \frac{\mathrm{1}+\mathit{\delta }}{r_{\mathrm{i}}}}S\text{at}x=x_{\mathrm{g}}.\end{array}$$

Since r_i is known, this model effectively contains one parameter α, which we treat as constant (although another version of this model could have α varying with time to reflect changes in accumulation). Sensitivity of p_S to α is low due to the high powers appearing in Eq. (31). Our chosen value of α in Sect. 5 is based on present-day measurements and historic reconstructions of ice thickness and is consistent with realistic values of β and a. In Sects. 3 and 4, a lower but plausible value of α is chosen for illustrative purposes.

We have chosen the above ice sheet model because it is simple and quick to solve for a given grounding line position and is physically appropriate for an ice stream context, where fast ice flow is dominated by basal slip (Cuffey and Paterson, 2010). This is because our main purpose is to model the dynamics of groundwater, rather than those of the ice sheet. A more general model could replace the Weertman sliding law with other sliding laws. However, such sliding laws generally involve coupling to the effective pressure in the shallow hydrological system, the modelling of which is beyond the scope of this paper. We include in Sect. 6 a discussion of the limitations of this model and of how it might be generalised.

3 Solutions in an idealised aquifer geometry

3.1 Steady states

Although, for realistic parameter values (see Table 1), the system is usually far from a steady state, it is useful to understand these states in order to gain insight into solutions of Eq. (28) more generally. The steady states of Eq. (22) satisfy

$$\begin{array}{}\text{(34)}& h{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S+\mathit{\delta }s\right)=\mathrm{0},\end{array}$$

using the no-flux condition (24). We therefore have

$$\begin{array}{}\text{(35)}& h=\mathrm{0}\text{or}{p}_S+S+\mathit{\delta }s=\text{const.}\end{array}$$

By the condition (25), the latter must hold near the grounding line x_g. There are therefore two basic forms of steady-state solution (excluding possible steady states of the complementarity problem (28) in which s=S over some region), a “lens” case and “nose” case (Bear, 2013), which we discuss below.

Figure 2 (a) Interface s(x) and (b) exfiltration flux q_E, for (i) “lens” and (ii) “nose” steady states arising from different grounding line positions x_g1 and x_g2, for a uniform depth aquifer with $b=-\mathrm{3000}$ m, $S=-\mathrm{1000}$ m. Parameters are as listed in Tables 1 and 2, with $k={\mathrm{10}}^{-\mathrm{12}}$ m² (K=0.41) and α=0.1.

3.1.1 Lens case

One possibility is that h>0 everywhere, resulting in a freshwater “lens” beneath the ice sheet. In this case, conditions (11) and (25) imply

$$\begin{array}{}\text{(36)}& p_S+S+\mathit{\delta }s=\mathrm{0}\text{for}\mathrm{0}<x<x_{\mathrm{g}}.\end{array}$$

Since we must have s>b, such a solution requires that

$$\begin{array}{}\text{(37)}& p_S+S+\mathit{\delta }b<\mathrm{0}\text{for}\mathrm{0}<x<x_{\mathrm{g}}.\end{array}$$

This generally occurs when the ice sheet is shallow relative to the depth of the sedimentary basin or when the grounding line has retreated far inland. The interface in the lens solution does not coincide with the upper surface z=S, provided that s<S, i.e.

$$\begin{array}{}\text{(38)}& S+\mathit{\delta }b\le -p_S\le (\mathrm{1}+\mathit{\delta })S.\end{array}$$

3.1.2 Nose case

If the condition (37) is violated for some x, then the lens solution must be replaced locally by the solution h=0, since h cannot become negative. The solution therefore displays a saltwater “nose” at x=x_n, where

$$\begin{array}{}\text{(39)}& x_n=max\mathit{\{}x<x_{\mathrm{g}}:p_S(x)+S(x)+\mathit{\delta }b(x)=\mathrm{0}\mathit{\}}.\end{array}$$

The basic nose solution is given by Eq. (36) for x>x_n and by h=0 for $\mathrm{0}<x<x_n$, though such a steady state is not always unique, as we shall discuss shortly.

Examples of the lens and nose cases are shown in Fig. 2. (Note that here and in subsequent figures, the aspect ratio is exaggerated, and a schematic ice shelf is depicted for illustrative purposes only.) Freshwater is recharged into the aquifer far inland and discharged closer to the grounding line, with the total recharge and discharge balancing one another in the steady state. There is no flow in the saltwater region in the steady state. The infiltration q_E<0 is large in magnitude but finite at x=0, because the ice overpressure gradient $\partial P_S/\partial x$ is zero at x=0 whilst changing over a relatively short distance.

In the nose case, $\partial s/\partial x$ has a discontinuity at x_n, leading to a kink in the streamlines and a jump discontinuity in q_E at this point. This is a quirk of the steady-state solution and is only possible because the flux of saltwater is identically zero: as we shall see in Sect. 4.1, the transient solution never reaches the h=0 state from a nonzero initial condition, meaning that this kink in h and the corresponding jump in q_E are smoothed out in practice. However, the presence of this jump is an indication that saltwater intrusion may have non-negligible effects on freshwater exchange with the basal hydrological system by diverting the flow of fresh groundwater.

3.2 Periodic forcing

The above steady states can only be realised if the grounding line position x_g(t) and ice overpressure p_S(x,t) remain steady. However, Antarctic ice streams such as those in the Ross Embayment have undergone significant grounding line movement on timescales of tens to hundreds of thousands of years. For typical values of the dimensionless hydraulic conductivity K (see Table 2), the present-day groundwater dynamics are therefore unlikely to have reached a steady state since the most recent grounding line migration.

We therefore investigate dynamic solutions of Eq. (28), considering in particular periodic solutions. This is motivated by the fact that global temperatures during the Pleistocene display periodicity of around 100 kyr, associated with variations in the Earth's orbit of the Sun (Willeit et al., 2019), which suggests that a periodic form of the grounding line forcing x_g(t) is appropriate. In particular, we use a dimensionless forcing:

$$\begin{array}{}\text{(40)}& x_{\mathrm{g}}\left(t\right)={\overline{x}}_{\mathrm{g}}-\mathrm{\Delta }{x}_{\mathrm{g}}\cos \left(\mathrm{2}\mathit{\pi }t\right),\end{array}$$

wherein ${\overline{x}}_{\mathrm{g}}$ and Δx_g are constants.

Figure 3 Stable periodic solution for s(x,t) during (i) advance and (ii) retreat and q_E(x,t) during (iii) advance and (iv) retreat, with periodic grounding line forcing given by Eq. (40), for (a) K=10 (high permeability), (b) K=1 (moderate permeability), and (c) K=0.1 (low permeability). Instantaneous steady states are plotted for $x_{\mathrm{g}}={\overline{x}}_{\mathrm{g}},{\overline{x}}_{\mathrm{g}}\pm \mathrm{\Delta }{x}_{\mathrm{g}}$. Parameters are as listed in Tables 1 and 2 with α=0.1. For full solutions, see Supplementary Animations S1–S3. An example for K=0.01 is shown in Fig. S1 in the Supplement.

Figure 3, along with Supplementary Animations S1–S3, shows solutions to Eq. (28) for periodic grounding line forcing at various values of K. The solution $z=s(x,t)$ is initialised in a steady state and then calculated over several cycles of advance and retreat until it reaches a stable periodic solution. Only the final cycle of advance and retreat, which represents this periodic solution, is plotted.

For all values of K, the interface lies between the would-be steady states at the minimum and maximum grounding line positions within the cycle. Although the forcing is symmetric in time, the interface $z=s(x,t)$ takes different shapes during advance and retreat.

When K=10 (Fig. 3a(i)), groundwater flow is fast compared to grounding line movement. The interface therefore adjusts relatively rapidly and is consequently close to the quasi-steady state corresponding to the instantaneous grounding line position. We see that the steady-state solutions are most relevant in the limit of large K.

When K=1 (Fig. 3b(i)), groundwater flow takes place at a moderate speed. There is now a substantial region of saltwater, which lies below the interface z=s throughout the whole cycle but above the steady-state position of the interface at the most advanced grounding line position (black dashed line). We can regard this as dynamically “trapped” saltwater, which would be displaced by freshwater in the steady state for the most advanced grounding line position but remains present in the periodic solution due to the unsteadiness of x_g(t). For this particular value of K, the lens height $H(\mathrm{0},t)-h(\mathrm{0},t)$ achieves its maximum near the minimum grounding line position, so that the dynamics of the interface at x=0 are nearly in antiphase with those at the grounding line. This results in an additional region of trapped freshwater, which always lies above the interface s(x,t) but below its steady-state position at the most retreated grounding line position.

Conversely, when K=0.1 (Fig. 3c(i)), groundwater flow is slow compared to grounding line movement. In this case, s(x,t) is near-steady throughout most of the domain, apart from a thin freshwater region near the grounding line. The corresponding amount of saltwater trapped by the periodic motion is very large, while the corresponding amount of trapped freshwater is small.

Notably, the quasi-steady portion of the solution for small K does not correspond to the steady state at the mean grounding line position ${\overline{x}}_{\mathrm{g}}$. We hypothesise that as K→0, this state would approach the steady state for the most retreated grounding line position. (For an example with K=0.01, which supports this hypothesis, see Supplement Fig. S1.)

This is because, during grounding line retreat, seawater instantaneously displaces any freshwater from the newly exposed region x>x_g. However, during re-advance, the freshwater lens grows at a finite rate, determined by the infiltration q_E<0 and depending on K. This is the key cause of asymmetry between grounding line advance and retreat. We infer from this example that the furthest inland grounding line position is most important for determining the extent of saltwater intrusion whenever grounding line movement is fast compared to the natural timescale of groundwater flow.

The exfiltration flux q_E also depends strongly on K, as shown in Fig. 3a–c(ii). Equation (27) shows that q_E is multiplied through by a factor of K, but it also depends on the position of the interface s, which is itself affected by K. When K is large (Fig. 3a(iii) and (iv)), q_E therefore follows s in being close to its instantaneous quasi-steady-state value, resembling that shown in Fig. 2b(i). However, when K is small (Fig. 3c(iii) and (iv)), exfiltration is small apart from a prominent peak near x_g during grounding line retreat, where freshwater is rapidly flushed out of the aquifer as seawater saturates the region beyond the grounding line. Seawater intrusion could therefore significantly influence the shallow hydrological system via its effect on q_E, as we shall see again in Sect. 5.

4 Solutions in a non-uniform aquifer geometry

4.1 “Pocket” steady states

For certain basement geometries b and S and overpressures p_S, Eq. (22) also admits steady states that, in addition to a nose, contain saltwater “pockets”, where h>0 on a discrete region that is not connected to the saltwater nose. An illustrative example of these “pocket” steady states is shown in Fig. 4. Such a pocket may be formed, for example, if the ice has advanced sufficiently rapidly to its current position that saltwater is still present relatively far inland once the ice has halted (see later Figs. 5 and 6). These pockets represent another form of saltwater trapping, which can occur in the steady-state solution rather than requiring ongoing grounding line movement. For a single pocket in $x_q<x<x_p$, where x_p<x_n, we have

$$\begin{array}{}\text{(41)}& {\displaystyle }\begin{array}{rl}{p}_S\left(x\right)& +S\left(x\right)+\mathit{\delta }s\left(x\right)=p_S\left(x_p\right)+S\left(x_p\right)+\mathit{\delta }s\left(x_p\right)\\ & \text{for}{x}_q<x<x_p,\end{array}\text{(42)}& {\displaystyle }h=\mathrm{0},{\displaystyle \frac{\partial h}{\partial x}}\le \mathrm{0}\text{at}x=x_{p-},\end{array}$$

where the “−” subscript denotes a limiting value taken from the left-hand side, and

$$\begin{array}{}\text{(43)}& h=\mathrm{0}\text{at}x=x_q.\end{array}$$

Such a pocket is therefore possible if and only if

$$\begin{array}{}\text{(44)}& {\displaystyle \frac{\partial }{\partial x}}\left(p_S+S+\mathit{\delta }b\right)\ge \mathrm{0}\text{at}x=x_p.\end{array}$$

The pocket solution is then given by

$$\begin{array}{}\text{(45)}& \begin{array}{rl}h& ={\displaystyle \frac{\mathrm{1}}{\mathit{\delta }}}\left[p_S\right(x_p)-p_S(x)+S(x_p)-S(x\left)\right]+b\left(x_p\right)\\ & -b\left(x\right)\text{for}{x}_q\le x\le x_p,\end{array}\end{array}$$

where x_q is the largest x<x_p for which this expression is zero. We impose that h=0 for x<x_q to avoid the above solution becoming negative, making x_q effectively a second nose (although in general, further pockets could exist in the region x<x_q). (It is also possible for a pocket to form a saltwater lens if the above expression for h is positive for $\mathrm{0}\le x\le x_p$.)

Since $\partial p_S/\partial x<\mathrm{0}$ in a realistic setting, and $|\partial S/\partial x|$ is usually small compared to $|\partial b/\partial x|$ (see Sect. 5), Eq. (44) effectively requires that $\partial b/\partial x$ is positive and relatively large (of the order of $\mathrm{1}/\mathit{\delta }\approx \mathrm{40}$). However, the Dupuit approximation (Eq. 4) remains valid in such a geometry, since the aspect ratio of the aquifer is around $\mathrm{1}/\mathrm{500}$ (see Table 2).

Figure 4 (a) A case with an infinite number of possible “pocket” steady states, with the criteria in Eq. (44) illustrated. The geometry is given by $S=-\mathrm{1000}$ m, $b=(-\mathrm{2500}+\mathrm{1000}\exp (-\left(\right(x-\mathrm{125}\mathrm{km})/\mathrm{12.5}\mathrm{km}{)}^{\mathrm{2}}))$ m. (b) Freshwater exfiltration and (c) streamlines, with no pocket and the largest possible pocket. Parameters are as listed in Tables 1 and 2, with $k={\mathrm{10}}^{-\mathrm{12}}$ m², K=0.41, and α=0.05.

In the example of Fig. 4, the sedimentary basin contains an idealised bottleneck due to a local rise in the basement rock. Behind this obstacle, the basement b(x) is sufficiently steep for pockets to form via the criterion (44).

We see that the bottleneck in the sediment geometry also produces local exfiltration and subsequent re-infiltration of freshwater, as the flux of freshwater (set by the ice overpressure) is forced through a shallowing and subsequently deepening aquifer. However, the presence of a trapped saltwater pocket extends the distance over which the freshwater region grows shallower, leading to smaller exfiltration over a larger area and substantially modifying the spatial dependence of q_E. We infer that potential saltwater pockets, along with saltwater intrusion more generally, may significantly impact the role of sedimentary basins in basal hydrology.

When pocket steady states exist, they are non-unique, since we may choose any x_p from the interval on which Eq. (44) holds, or simply have h=0 for all x<x_n. Figure 4 also illustrates this infinite family of steady states and their relation to the criterion in Eq. (44). Choosing the smallest such x_p leads to the “maximal” pocket, which contains the largest possible volume of trapped saltwater. Depending on the geometry, it may also be possible to satisfy Eq. (44) on multiple disjoint intervals so that steady states exist with several discrete saltwater pockets. However, it is not yet clear how, if at all, the existence of pocket steady states affects the transient solution of Eq. (28). It is natural to ask, firstly, whether pocket steady states are stable and, if so, whether a realistic combination of initial conditions and forcings can lead to such states. We address these questions in the following subsection.

4.2 Stability of pocket steady states

For the case of a steady ice sheet p_S and grounding line x_g, the transient solution of Eq. (28) relaxes towards a steady state. However, in the case where there are multiple possible steady states, it is unclear which will be achieved. In this and the following subsection, we demonstrate that a range of final steady states with or without pockets can be obtained, depending on initial conditions and on the dimensionless hydraulic conductivity K.

Figure 5 (a) Relaxation from an initial state h=H towards a pocket steady state. (b) Relaxation towards a steady state from an initially perturbed pocket-free state, with the direction of 𝒰(x) indicated. Note that $x_p^{*}$ is not at the maximum of b(x) but a point where $\partial b/\partial x>\mathrm{0}$. Parameters are as in Tables 1 and 2, with K=10 and α=0.05. For full solutions, see the Supplementary Animations S4 and S5.

Figure 5a and the Supplementary Animation S4 show a transient solution for the geometry and ice sheet shown in Fig. 4, in which the aquifer is initially saturated with saltwater throughout (h=H). As freshwater enters the aquifer from above, a significant volume of saltwater is trapped behind the obstacle, resulting in an eventual pocket steady state. However, as becomes evident after a very long period of time, the final state does not have the largest possible pocket, as might be expected. This overshoot occurs because the slope of the steady interface $\partial s/\partial x$ fails to be continuous at x_p and x_q, where the solution switches between h=0 and h>0. On the other hand, the transient solution s(x,t) has a continuous derivative due to the regularity imposed by Eq. (28). As a result, this solution is smoothed compared to the steady state, as seen for instance at 10² kyr. The smoothed region near x_p in turn allows a small amount of trapped saltwater to escape the pocket.

We next consider the opposite case of a small initial perturbation to the pocket-free state, as shown in Fig. 5b and the Supplementary Animation S5. Near h=0, the linearised version of Eq. (22) takes the form

$$\begin{array}{}\text{(46)}& {\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left[\mathcal{U}\left(x\right)h\right]=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(47)}& \mathcal{U}\left(x\right)=-K{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S+\mathit{\delta }b\right).\end{array}$$

This is a hyperbolic equation representing advection with velocity 𝒰. Equation (47) indicates that, at leading order, the flow is driven by the overburden and topographic pressure gradients, and the second-order buoyancy term is negligible.

The direction in which a perturbation (i.e. a small volume of saltwater) travels is determined by the sign of 𝒰, which is shown in Fig. 5b. The condition (44) for a possible pocket in $x_q<x<x_p$ is equivalent to 𝒰(x_p)<0, and the largest possible pocket occupies $x_q<x<x_p^{*}$, where $\mathcal{U}\left(x_p^{*}\right)=\mathrm{0}$. Therefore, a small volume of saltwater in $x<x_p^{*}$ close to $x_p^{*}$ will move in the negative x direction and become trapped in the region $x<x_p^{*}$. This trapped saltwater will in turn eventually form a pocket in the final steady state. Meanwhile, a small volume of saltwater in $x_p^{*}<x<x_n$ will travel in the positive x direction and reach the saltwater nose.

To illustrate this, Fig. 5b and the Supplementary Animation S5 show the result of initially perturbing a pocket-free state. The eventual steady state towards which the solution relaxes includes a pocket of trapped saltwater. The volume of saltwater in $\mathrm{0}<x<x_p^{*}$, given by the integral of h(x,t) over this region, changes by less than 2 % from the initial to final state, strongly suggesting that this pocket has been formed by trapping the saltwater initially located in this region, via the linearised flow 𝒰. Conversely, the saltwater in $x>x_p^{*}$ has been evacuated towards the ocean.

4.3 Pocket creation by grounding line advance

If grounding line movement brings about a change from a system with a unique steady state to one with multiple steady states, it is unclear towards which of the new available steady states the solution of Eq. (28) will relax. Here we show that it is possible for the system to reach either a pocket state or one without, depending on the intermediate states and on the parameter K.

For the purposes of plotting our results, we also define the total volume of freshwater:

$$\begin{array}{}\text{(48)}& V_{\mathrm{f}}=\underset{\mathrm{0}}{\overset{x_{\mathrm{g}}}{\int }}H-h\mathrm{d}x.\end{array}$$

Figure 6 and the Supplementary Animation S6 include an example wherein the grounding line advances from an inland position past a bottleneck similar to that shown in Fig. 4. The final grounding line position accommodates multiple steady states, although only the maximal and minimal such states are shown. Different final steady states are reached depending on the relative speed of groundwater flow, described by the parameter K. In this case we select values of K that illustrate a difference in final states whilst being relatively fast to approach these states, for the sake of readability of the figure.

Figure 6 (a–d) Solutions $z=s(x,t)$ plotted at various points in time for two different K. (e) Grounding line position x_g. (f) Difference in freshwater volumes $\mathrm{\Delta }{V}_{\mathrm{f}}\left(t\right)=V_{f,K=\mathrm{10}}-V_{f,K=\mathrm{0.5}}$ between the two cases. Parameters are as listed in Tables 1 and 2, with α=0.1. For full solutions, see Supplementary Animation S6.

When groundwater flow is sufficiently slow compared to grounding line movement (e.g. for K=0.5), the interface slowly relaxes from above to the maximal pocket state, resulting in trapped saltwater. On the other hand, when groundwater flow is fast (e.g. for K=10), the interface rapidly relaxes to the quasi-steady state for the instantaneous grounding line position. If one of these intermediate quasi-steady states inundates the region of the would-be pocket with freshwater, the final steady state has no pocket and no saltwater is trapped.

Similarly to the case of periodic forcing in Sect. 3.2, we see from this example that the trapping of saltwater (in this case via a pocket steady state) is dependent on the dimensionless hydraulic conductivity K as well as previous grounding line forcing. Indeed, this model demonstrates the possibility of “hysteresis”, where grounding line retreat and re-advance induces a transition between two possible steady states. This hysteresis would affect the exfiltration q_E into the shallow hydrological system, because the presence of a saltwater pocket influences q_E (e.g. as seen in Fig. 4). This in turn would modify the distribution of water in the shallow hydrological system beneath an ice sheet and could therefore feed back on the dynamics of the ice via basal sliding.

5 Application: Ross Sea, West Antarctica

5.1 Design

Having explored the above model of seawater intrusion in idealised settings, we now seek to apply it to a more realistic scenario, using the sedimentary basin beneath the Ross Ice Shelf in West Antarctica as a case study. As well as being the location of observed paleo-seawater intrusions by Gustafson et al. (2022), the thickness of this sedimentary basin has been modelled in detail using airborne magnetic data by Tankersley et al. (2022b), and the grounding line history of this region has been subject to extensive investigation and modelling (Venturelli et al., 2023, 2020; Neuhaus et al., 2021; Lowry et al., 2019; Kingslake et al., 2018). Since the data are three-dimensional, we investigate a transect of the data corresponding roughly to a flowline of Whillans Ice Stream (more precisely, we average between a set of adjacent transects to reduce uncertainty due to lateral variation). A summary of the basement geometry and grounding line data used, along with our choices of transect, and the sites of magnetotelluric measurements of salinity by Gustafson et al. (2022), is plotted in Fig. 7.

Figure 7 Measurements of (a) upper surface S, (b) basement b, and (c) thickness H for the Ross Sea sedimentary basin, which we sample along transects (Fretwell et al., 2013; Rignot et al., 2013; Greene et al., 2017; Mouginot et al., 2017; Kingslake et al., 2018; Gustafson et al., 2022; Tankersley et al., 2022a). (d) Averaged cross-section along transects.

5.2 Aquifer geometry

Our estimates for H(x) and b(x) along the selected flowline come from the results of Tankersley et al. (2022b), which are derived from airborne magnetic data. We combine this with bedrock topography and bathymetry data for S(x) from Fretwell et al. (2013). Since the measurements of b and H only extend a short distance beyond the present-day grounding line, we linearly extrapolate b in this region, assuming that H(0)=0 based on localised seismic measurements (Peters et al., 2006). The onset of the sedimentary basin x=0 corresponds to the onset of ice streaming, which can be inferred from surface ice flow velocity. We also select a point on the flowline, marked with a red cross, corresponding to the Subglacial Lake Whillans (SLW) observation site, where the salinity of the sedimentary basin has been modelled by Gustafson et al. (2022).

5.3 Grounding line and ice sheet model

Our reconstruction of the grounding line position x_g(t) is based on paleo-ice-sheet modelling supported by geophysical observations, which suggests the following (Kingslake et al., 2018):

1. slow advance to a maximum, located at the edge of the Antarctic continental shelf, for 100–20 kyr before present;

2. rapid retreat to a minimum, several hundred kilometres inland of the present-day grounding line position, for 20–10 kyr before present;

3. re-advance from this minimum to the present-day grounding line for 10–0 kyr before present.

Since there is substantial uncertainty in the precise speed and timing of grounding line advance and retreat, we simply linearly interpolate x_g(t) between the maximum, minimum, and present-day position, as shown in Figs. 10a and 11a.

We assume that the grounding line behaviour obeys 100 kyr periodicity, in accordance with the periodic behaviour of global temperatures during the Pleistocene (Pollard and DeConto, 2009; McKay et al., 2012). We therefore repeat our approach from Sect. 3.2 wherein the model is run for a series of forcing cycles until a stable periodic solution is reached.

To model p_S, we employ the steady shallow ice sheet model described in Eq. (30), with the constant α=5 based on the results of paleo-ice sheet reconstructions using more sophisticated models (Huybrechts, 2002; Kingslake et al., 2018). We choose to use these simplified models on the basis that the exact forms of p_S and x_g are a relatively small source of model uncertainty compared to, for example, lateral variations in H and b or heterogeneity in aquifer properties such as k and ϕ. We have therefore prioritised capturing appropriate values for the grounding line maximum and minimum and their timing. A more sophisticated approach could use ice core data for historic accumulation and temperature to force a model of both p_S and x_g, but such a model may require excessive fine-tuning to reproduce existing estimates of the grounding line extrema.

5.4 Aquifer properties and timescales

The largest source of uncertainty is given by the hydraulic properties of the aquifer, namely k and ϕ. This includes model uncertainty (e.g. possible spatial heterogeneity in ϕ and k, dependence of k on ϕ, and possible anisotropy in k) and parameter uncertainty (assuming that k is constant, what should its value be?). Previous studies have accounted for this by solving models for a wide range of possible k, typically 10⁻¹⁷–10⁻¹³ m² (Gooch et al., 2016; Li et al., 2022; Robel et al., 2023). In the dimensionless model, the values of k and ϕ are both accounted for in the single parameter K, which takes small to moderate values for the above range of k (see Table 3). Note that we disregard values larger than $k=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{14}}$ m², which gives K≈0.01, on the basis we do not expect solutions for very small values of K to differ substantially, and that we include an artificially high value $k=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ m² (K≈40) for comparison. In Appendix B, we discuss how Eq. (22) can be generalised to the case of variable $k(x,z,t)$ and $\mathit{\varphi }(x,z,t)$.

Figure 8 Potential steady-state pockets for the basement geometry z=b(x) and z=S(x) (black curves) estimated along a historic flowline of the ice in the Ross Sea, West Antarctica, as shown in Fig. 7. The freshwater–saltwater interface z=s(x) and ice surface are plotted for various values of the ice sheet parameter α. (Any pockets in the region (x<400 km) where b(x) has been linearly interpolated are ignored.) Parameters are as listed in Tables 1 and 2.

Figure 9 Interface s(x,t) for varying k, plotted (a) during grounding line advance (b) at grounding line maximum (c) at grounding line minimum (d) at present day, with the SLW site marked by a black dashed line. Parameters are as listed in Tables 1 and 2, with α=5. For full solutions see Supplementary Animations S7–S12.

Figure 10 (a) Grounding line position x_g(t). (b) Dimensionless spatially averaged exfiltration ${\overline{q}}_{\mathrm{E}}$ after scaling with K. (c) Dimensional ${\overline{q}}_{\mathrm{E}}$. Parameters are as listed in Tables 1 and 2, with α=5.

Figure 11 (a) Basement geometry b(x) and S(x), coloured according to $-\partial H/\partial x$. (b) Space–time plot of q_E(x,t), with red corresponding to $q_{\mathrm{E}}(x,t)>\mathrm{0}$ (exfiltration) and blue to q_E<0 (infiltration). Values with $\left|q_{\mathrm{E}}\right|>\mathrm{500}$ mm yr⁻¹ are plotted using the same colours as ±500. Parameters are as listed in Tables 1 and 2, with α=5 and $k=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{12}}$ m².

Table 3 Depth of the freshwater lens at the SLW site for various values of k.

5.5 Results

Figure 8 shows where potential steady-state saltwater pockets can exist in the above geometry for a range of ice thicknesses, which in our model are characterised by the parameter α. For our choice of α=5, which is based on historic ice thickness modelled by Kingslake et al. (2018) and Huybrechts (2002), the ice surface is too steep to permit saltwater pockets from existing in the steady state. However, such pockets are possible for shallower ice sheets or in the transient case.

The results of our time-dependent simulations are shown in Figs. 9, 10, and 11 and Table 3, as well as Supplementary Animations S7–S12. For each value of k, we calculate the present depth of the freshwater lens at the SLW site:

$$\begin{array}{}\text{(49)}& d_{\mathrm{f},\mathrm{SLW}}=S\left(x_{\text{SLW}}\right)-s(x_{\text{SLW}},\mathrm{0}),\end{array}$$

along with the space-averaged exfiltration

$$\begin{array}{}\text{(50)}& {\overline{q}}_{\mathrm{E}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{x_{\mathrm{g}}}}\underset{\mathrm{0}}{\overset{x_{\mathrm{g}}\left(t\right)}{\int }}{q}_{\mathrm{E}}(x,t)\mathrm{d}x.\end{array}$$

Figure 9 shows the interface s(x,t) at selected points in time for various values of k, and Table 3 shows the corresponding values of d_f,SLW. Figure 10 shows a time series of the space-averaged exfiltration ${\overline{q}}_{\mathrm{E}}$, both dimensionally and dimensionlessly after scaling with K, and Fig. 11 shows a space–time plot of q_E(x,t) for the case $k=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{12}}$ m² (K=1.2), alongside the basement geometry for comparison.

For the very highly permeable test case ($k=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ m²) in Fig. 9, the interface is close to the quasi-steady state for the current grounding line position, as expected. The resemblance to the steady state is also observed for $k=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{12}}$ m², in addition to a thin layer of trapped saltwater. For the successively lower permeabilities, the depth of the transient freshwater lens becomes successively less, and correspondingly more saltwater is trapped due to grounding line movement. For the less permeable examples, we see that both possible solution forms within the complementarity problem (28) are required, as the solution exhibits small but significant regions of saltwater exfiltration where h=H (e.g. near x=600 m and x=1200 m in Fig. 9b).

Table 3 shows the present-day depth of the transient freshwater lens d_f,SLW at the location of SLW in the above solutions. This increases with the permeability k, reaching the full depth H of the layer for the very high value of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ m², where the solution is nearly quasi-steady. Modelled salinity profiles found by inversion of magnetotelluric data by Gustafson et al. (2022) suggest that the groundwater salinity reaches that of seawater around 600 m below the ice bed and 50 % of this salinity about 400 m below the ice bed. These profiles show a smooth increase in salinity as a result of either mixing of freshwater and saltwater or due to the inversion process itself, as opposed to our results obtained assuming a sharp freshwater–saltwater interface. This introduces some uncertainty in how best to infer the “freshwater lens depth” from these results. However, it is clear that the high permeability value $k=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{12}}$ m², for which $d_{\mathrm{f},\mathrm{SLW}}=$ 400 m, gives the results that can be considered to agree best with those of Gustafson et al. (2022). We discuss how this model might be adapted to allow smoothly varying salinity in Sect. 6.

The exfiltration q_E scales roughly linearly with K (see Eq. 27), but it also depends on K through s(x,t). Figure 10b allows us to see the importance of the latter dependence: if q_E were determined by basement geometry and ice overpressure alone (i.e. if seawater intrusion were ignored), the curves would collapse onto one another. However, the magnitude of ${\overline{q}}_{\mathrm{E}}/K$ is at points several times larger for $k=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{14}}$ m², which includes the most trapped saltwater, than for $k=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{10}}$ m², where effectively no saltwater is trapped. We therefore note that saltwater intrusion and trapping may have a substantial effect on groundwater exchange with the shallow hydrological system. We also note that ${\overline{q}}_{\mathrm{E}}<\mathrm{0}$ throughout most of the cycle, meaning that, on average across the basin, more water infiltrates than is exfiltrated (ultimately driven by the downstream hydraulic potential gradient due to the slope of the ice surface). However, around the glacial minimum exfiltration dominates. Over the duration of a periodic cycle, the net flux of freshwater into the domain must be zero, which is achieved by localised discharge of freshwater at the grounding line x_g during retreat.

Figure 10c allows us to compare the dimensional values of ${\overline{q}}_{\mathrm{E}}$. For the large values $k=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{12}}$ m² and $k=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{12}}$ m², the magnitude of ${\overline{q}}_{\mathrm{E}}$ is in the low tens of millimetres per year, dropping to a few millimetres per year for lower permeabilities. However, local values of q_E may be several times larger than the average ${\overline{q}}_{\mathrm{E}}$, as seen in Fig. 11. Previous modelling has found that groundwater contributes significantly to subglacial groundwater budgets (Christoffersen et al., 2014). This suggests that a realistic value of q_E should be similar in size to other sources of basal water, such as the basal melt rate, which may reach 20–100 mm yr⁻¹ beneath an ice stream (Joughin et al., 2009). To achieve a mean exfiltration of this order of magnitude, we require high permeabilities $k\gtrsim {\mathrm{10}}^{-\mathrm{12}}$ m², which is consistent with our previous estimate based on freshwater lens depth. For smaller permeabilities, such values of q_E are locally achievable but are counterbalanced by infiltration (q_E<0) elsewhere.

Moreover, any estimate obtained under the assumptions of our model will fail to consider short-term exfiltration due to poro-elastic deformation of the sediments, which previous modelling suggests may reach tens to hundreds of mm yr⁻¹ for $k=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{13}}$ m² (Li et al., 2022; Robel et al., 2023). This short-term exfiltration may enable sedimentary basins to act as a net source of basal water, as is suggested by the modelling of Christoffersen et al. (2014), rather than a sink. However, further modelling is required to understand the combined effects of short-term and long-term exfiltration.

In addition, q_E is highly heterogeneous in space, as seen in Fig. 11, displaying localised regions of significant exfiltration and infiltration. The sign of q_E is closely correlated to that of $-\partial H/\partial x$, with exfiltration generally occurring where the aquifer is shallowing in the direction of flow and infiltration where it is deepening. This correlation occurs because the terms $\partial S/\partial x$ and $\mathit{\delta }\partial s/\partial x$ are relatively small in Eq. (27) so that we have

$$\begin{array}{}\text{(51)}& q_{\mathrm{E}}\approx K{\displaystyle \frac{\partial }{\partial x}}\left(H{\displaystyle \frac{\partial p_S}{\partial x}}\right)\approx -K{\displaystyle \frac{\partial H}{\partial x}}\left|{\displaystyle \frac{\partial p_S}{\partial x}}\right|,\end{array}$$

with the latter approximation holding because $\partial p_S/\partial x$ varies relatively little compared to H throughout the aquifer. The striped pattern of Fig. 11 indicates that this term dominates throughout most of the domain over time. However, there are also visible deviations from this pattern due to the other terms, such as the role of saltwater intrusion as discussed above. These are most prominent near the grounding line, where the interface s(x,t) and ice thickness H_i(x,t) are steepest.

Overall, we therefore find that a high permeability $k\gtrsim {\mathrm{10}}^{-\mathrm{12}}$ m² is required to reproduce the observed depth of the freshwater lens at the SLW site and is consistent with the modelled contribution of groundwater to the basal water budget. This is an unusually high value for a sedimentary rock, suggesting that geological features such as faults and fractures, which our model has not considered explicitly, but which may macroscopically be captured by a large permeability, may play an important role in subglacial groundwater transport. Our chosen value of k represents an idealisation of a sedimentary aquifer, which may in reality be highly heterogeneous and anisotropic.

6 Discussion

In this paper, we have developed a mathematical model for groundwater flow and seawater intrusion in the sedimentary basin beneath a marine ice sheet, considering in particular the effects of grounding line advance and retreat and of sedimentary basin geometry. Both of these lead to the phenomenon of saltwater trapping, where saltwater permanently resides in regions where a simplistic steady-state solution would predict freshwater. Saltwater trapping is also highly dependent on hydraulic properties, which we have accounted for with the parameter K. The concept of saltwater trapping, via pocket steady states, should be generalised from a subglacial context to a wide range of seawater intrusion problems provided that the aquifer geometry is sufficiently nonuniform. We discuss this generalisation in Appendix B.

In order to understand the role of groundwater in the subglacial hydrological system, it is crucial to constrain the dimensionless hydraulic conductivity K, since the flux q_E of water from the aquifer to the shallow hydrological system scales with this value. However, the appropriate value of K is highly uncertain, due to the observational challenges in measuring the porosity and permeability of such an aquifer, along with the fundamental limitations of modelling such an aquifer as homogeneous and isotropic (Appendix B discusses how this assumption may be relaxed). However, by modelling groundwater flow beneath the Ross Ice Shelf, using our knowledge of grounding line history, we are able to gain insight into the effect of K on saltwater trapping. By combining this modelling with present-day observations of trapped saltwater, we are able to estimate a physically appropriate range of values for K. We find that high values of the permeability ($k\sim {\mathrm{10}}^{-\mathrm{12}}$ m²) are required to reproduce observations of the freshwater lens depth and are also consistent with modelled contributions of groundwater to the basal water budget via exfiltration.

Our modelling also suggests that the spatial exfiltration q_E of groundwater into the shallow hydrological system is dominated by sedimentary basin geometry. Groundwater is generally exfiltrated where the aquifer is thinning in the direction of flow and infiltrates where the aquifer is thickening. However, the ice overpressure p_S(x,t) and the freshwater–saltwater interface s(x,t) also affect q_E(x,t). In the latter case, the presence of trapped or intruding saltwater changes the shape of the freshwater-saturated region of the aquifer and in some cases may entirely displace freshwater so that saltwater is exfiltrated into the groundwater system.

The latter possibility highlights a limitation of our model, namely our idealised treatment of the shallow hydrological system in assuming that only freshwater re-infiltrates the aquifer and that the effective pressure p_e=0. A more complex model that relates q_E and p_e by separately accounting for the dynamics of the shallow hydrological layer and that considers both freshwater and saltwater would provide a more subtle understanding of the role of groundwater exfiltration and infiltration in subglacial hydrology.

Such a model would also allow for coupling between the ice sheet model and subglacial hydrology. This could be achieved, for instance, by replacing the Weertman sliding law used in the ice sheet model above with a sliding law coupled to the effective pressure p_e (e.g. a regularised Coulomb sliding law Schoof (2005)). The use of such an ice sheet model would likely lead to ice that is thinner and shallower near the grounding line, leading to more saltwater intrusion in the steady state and a higher likelihood of trapping seawater in pockets, as well as leading to lower exfiltration/infiltration. However, the same modelled ice sheet may be steeper inland, resulting in a larger $\left|q_{\mathrm{E}}\right|$ further inland. It is also possible that feedbacks could exist between groundwater flow and the ice sheet. For example, groundwater exfiltration could lower the effective pressure in the shallow hydrological system, which would affect the sliding of the ice and hence the profile of the ice sheet, whose gradient feeds back into q_E.

Since the SLW site is relatively near the grounding line at the present day, we expect that an ice sheet model which predicts shallower ice here would lead to a thinner freshwater lens and hence a higher predicted permeability. However, the results of Sect. 5 suggest that, when periodic grounding line movement prevents groundwater dynamics from reaching a steady state, the shape of the ice sheet is less important than the permeability and geometry of the sedimentary basin. The resulting model uncertainty is therefore small compared to that resulting from, for example, aquifer heterogeneity and anisotropy or lateral variations in basement geometry.

Since basement geometry is important for the results of the model, particularly q_E, it should be noted that the use of a cross-sectional model introduces uncertainty even after transverse averaging. Moreover, the data of Tankersley et al. (2022a) include some model uncertainty introduced when inverting magnetic measurements. Future work is therefore required to explore the dynamics of subglacial groundwater flow in all three dimensions, which we have discussed in Appendix C. We have assumed, following Gustafson et al. (2022) and Tankersley et al. (2022b), that the basement may be treated as impermeable, although an extension of this model could include a basement which is weakly permeable. The inclusion of basement permeability would weaken the effects of basement geometry on q_E by providing an alternative route for groundwater to leave or re-enter the sedimentary basin.

Moreover, another possible development of this model could consider the effect of exfiltration and infiltration on subglacial heat transport. Previous modelling by Gooch et al. (2016) suggests that exchange of groundwater with the subglacial hydrological system may substantially affect geothermal heat flux to the ice bed and subsequently basal melting. Tankersley et al. (2022b) note that sedimentary basin geometry, in particular the presence of highly permeable faults, may affect ice sheet dynamics via the effect of groundwater on heat flux and on basal melting. Such a model would require an equation of heat transport via advection and diffusion and could include freshwater generation via basal melting of the ice. For a groundwater flux scale of [q_E]=10 mm yr⁻¹, a porous medium conductivity of k_c=3 W m⁻¹ K⁻¹, and a specific heat capacity of water of c_w=4000 J kg⁻¹ K⁻¹ (Gooch et al., 2016), the Péclet number associated with vertical heat transport is

$$\begin{array}{}\text{(52)}& {\mathit{Pe}}_{\text{heat}}={\displaystyle \frac{\left[q_{\mathrm{e}}\right]\left[z\right]{\mathit{\rho }}_{\mathrm{f}}{c}_{\mathrm{w}}}{k_{\mathrm{c}}}}\approx \mathrm{0.4},\end{array}$$

which is O(1), suggesting that both advection by groundwater and diffusion are important to subglacial vertical heat transport.

In addition, we have assumed that the aquifer is rigid and neglected the dynamic response of the sediments to loading. This enables us to model background rates of exfiltration and infiltration over the long timescales associated with grounding line movement and horizontal groundwater flow. However, a full account of the exchange of groundwater with the shallow hydrological system at a given point in time should also include exfiltration due to dynamic effects, namely the poro-elastic response of sedimentary basins to ice sheet loading, which has been separately modelled by Gooch et al. (2016), Li et al. (2022), and Robel et al. (2023). Indeed, short-term exfiltration may substantially exceed the long-term value: for example, Li et al. (2022) and Robel et al. (2023) obtain short-term exfiltration rates of 10–100 mm yr⁻¹ for $k={\mathrm{10}}^{-\mathrm{13}}$ m², compared to our modelled average ${\overline{q}}_{\mathrm{E}}\approx -\mathrm{2}$ mm yr⁻¹ for the same k. We have also neglected basal freeze-on, which could drive substantial exfiltration of water from sediments over the timescales of glacial advance and retreat, considered for instance by Christoffersen and Tulaczyk (2003). However, the dynamic effects of this exfiltration are typically confined to the upper 5–50 m of sediment, meaning that the corresponding pressure gradients do not substantially affect the overall flow of groundwater throughout the depth of the approximately 1 km thick sedimentary basin. Such exfiltration could be considered in a possible extension of this model that includes heat transport, as discussed above.

Our model has made the assumption of a sharp interface dividing freshwater and saltwater, on the basis that the Péclet numbers associated with diffusion and hydrodynamic dispersion are large. This greatly simplifies the task of solving the problem, at the cost of being unable to account for the smoothly varying salinity modelled by Gustafson et al. (2022), although this smoothness may in part a result of the model used to invert the magnetotelluric data. To account for the effect of mixing, it is possible (though more complicated) to solve the full problem of density-coupled salt transport, known as the Henry problem (e.g. Croucher and O'Sullivan, 1995), or to introduce a “mixing layer” around the sharp interface where these effects are accounted for (Van Duijn and Peletier, 1992; Paster and Dagan, 2007). In steady-state seawater intrusion problems, the inclusion of saltwater–freshwater mixing reduces the extent of seawater intrusion by inducing a circulatory flow in the saltwater-saturated region (Cooper, 1959; Koussis and Mazi, 2018). We therefore expect that including mixing would result in a more retreated “interface” (which would need to be newly defined as e.g. the 50 % salinity contour) than that achieved with the sharp-interface assumption at the same permeability. This would lead to a lower permeability estimate than that obtained above. The process of saltwater mixing would also reduce the trapping of saltwater in pockets via entrainment into the freshwater flow.

7 Conclusions

In conclusion, we have used an idealised model to demonstrate several previously unidentified features of a classical seawater intrusion problem that emerge when attempting to model groundwater flow beneath a marine ice sheet, namely the trapping of saltwater via grounding line movement and aquifer geometry. Our work provides new insights into the modelling of groundwater flow in subglacial sedimentary basins on geological timescales and the contribution of these sedimentary basins to Antarctic subglacial hydrology. Using this model, we are able to make quantitative estimates of groundwater exfiltration and infiltration throughout space and time and to infer hydraulic properties of the aquifer based on local measurements of salinity. Our work provides a starting point for more developed models of groundwater flow in sedimentary basins, which could include a more sophisticated treatment of hydraulic properties, three-dimensional flow, coupling groundwater flow to shallow hydrology and the ice sheet itself, or quantification of the uncertainty in modelling the above. Further modelling is required to advance our understanding of sedimentary basin hydrology and its role in determining the stability of the West Antarctic Ice Sheet.

Appendix A: Numerical method

In order to numerically solve Eq. (22) on the changing domain $\mathrm{0}<x<x_{\mathrm{g}}\left(t\right)$, we apply the coordinate transformation

$$\begin{array}{}\text{(A1)}& x=x_{\mathrm{g}}\left(t\right)X.\end{array}$$

The rescaled equation (after rearrangement into a more canonical conservation form) becomes

$$\begin{array}{}\text{(A2)}& {\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{\partial }{\partial X}}\left[K{\displaystyle \frac{h}{x_{\mathrm{g}}^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial X}}\left(p_S+S+\mathit{\delta }s\right)+{\displaystyle \frac{{\dot{x}}_{\mathrm{g}}}{x_{\mathrm{g}}}}Xh\right]-{\displaystyle \frac{{\dot{x}}_{\mathrm{g}}}{x_{\mathrm{g}}}}h,\end{array}$$

which may then be solved on the fixed domain $\mathrm{0}<X<\mathrm{1}$ subject to the conditions (18) and (19). To solve Eq. (A2) as part of the more general complementarity problem (28), we formulate the discrete complementarity problem:

$$\begin{array}{}\text{(A3)}& {\displaystyle }\begin{array}{rl}{U}^{\mathrm{2}}& =-{\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial }{\partial X}}\left[K{\displaystyle \frac{h}{x_{\mathrm{g}}^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial X}}\left(p_S+S+\mathit{\delta }s\right)+{\displaystyle \frac{{\dot{x}}_{\mathrm{g}}}{x_{\mathrm{g}}}}Xh\right]\\ & -{\displaystyle \frac{{\dot{x}}_{\mathrm{g}}}{x_{\mathrm{g}}}}h,\end{array}\text{(A4)}& {\displaystyle }{V}^{\mathrm{2}}=H-h,\text{(A5)}& {\displaystyle }UV=\mathrm{0}.\end{array}$$

This formulation replaces Eq. (22) with the more general “complementarity condition” (A5), and Eqs. (A3) and (A4) enforce the constraints that their respective right-hand sides are non-negative. This formulation ensures that either Eq. (A2) is obeyed and h≤H or h=H and the value $\partial h/\partial t$ implied by Eq. (A2) is non-negative (equivalent to q_E≥0 by Eq. 26). If this implied value of $\partial h/\partial t$ is negative (equivalent to q_E<0), the solution must detach from h=H.

We solve the system in Eqs. (A3), (A4), and (A5) using a finite volume method, which retains the conservative properties of Eq. (A2). The “diffusive” buoyancy term multiplied by δ is discretised using a second-order midpoint arithmetic mean, whereas the remaining “advective” flux terms are discretised using a first-order upwind method. The calculated flux terms may then be used to find q_E via Eq. (26). The time derivative is discretised using a first-order fully implicit scheme. The discrete equation at each time step is solved using MATLAB's inbuilt function fsolve. Between 100 and 200 spatial mesh points are used, with time steps between Δt=0.001 and 0.01, both depending on the specific problem in question. (Larger time steps up to Δt=1 are used for the late stages of the long-term simulation shown in Fig. 5a.) Smaller time steps are required during periods of faster flow (e.g. rapid grounding line movement) to ensure stability of the first-order upwind method.

Appendix B: Developments of the model

We have presented the above model with a number of simplifying assumptions, in order to illustrate the key features introduced by aquifer geometry and grounding line movement. However, several generalisations of this model are possible, a few of which we outline below.

B1 Prescribed recharge/unconfined aquifer

Our derivation of the above model assumes that p_S, the pressure at z=S, is known and given by the ice sheet overburden, which is consistent with the existence of a shallow layer of basal water in which the effective pressure is zero. However, it is also possible to write a model in which the exfiltration (or infiltration) q_E is prescribed, and p_S must be found. In this case, Eq. (26) becomes an equation for p_S:

$$\begin{array}{}\text{(B1)}& H{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S\right)+h\mathit{\delta }{\displaystyle \frac{\partial s}{\partial x}}=\underset{\mathrm{0}}{\overset{x}{\int }}{q}_{\mathrm{E}}\mathrm{d}x=:-Q,\end{array}$$

using the no-flux condition at x=0. This can then be substituted into Eq. (22) to obtain

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\partial h}{\partial t}}=K{\displaystyle \frac{\partial }{\partial x}}\left[h\left(-{\displaystyle \frac{Q}{H}}+\mathit{\delta }\left(\mathrm{1}-{\displaystyle \frac{h}{H}}\right){\displaystyle \frac{\partial s}{\partial x}}\right)\right].\end{array}$$

This equation displays the same features as Eq. (22). The ice overburden becomes entirely decoupled in this limit, only appearing implicitly through the basal hydrological system, which exchanges freshwater with the aquifer in $\mathrm{0}<x<x_{\mathrm{g}}$. Between the cases of prescribing p_S or q_E on z=S, there are also potential intermediate cases, where a relationship between p_S and q_E is prescribed in order to encapsulate the dynamics of the shallow basal hydrological system.

When considering a glacier that is not marine-terminating (e.g. a continental ice sheet during a glacial period), we would also encounter the possibility of an unconfined and unsaturated region of the aquifer between the ice sheet margin and the ocean. This region could also include permafrost, which substantially reduces aquifer permeability.

The region beneath the ice sheet is modelled as in this paper, with a confined aquifer in $b<z<S$ and the pressure p=p_S prescribed on the upper surface. To model the unsaturated region, we introduce a free boundary $z=\widehat{S}(x,t)=b+\widehat{H}$, representing the water table, which lies below the ground surface z=S. On $z=\widehat{S}$, both the exfiltration q_E<0 (representing rainfall or snow melt) and the pressure p_S=p_atm are prescribed. The resulting model consists of two equations:

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\partial h}{\partial t}}=K{\displaystyle \frac{\partial }{\partial x}}\left[h{\displaystyle \frac{\partial }{\partial x}}\left(\widehat{S}+\mathit{\delta }s\right)\right],\text{(B4)}& {\displaystyle \frac{\partial \widehat{H}}{\partial t}}=K{\displaystyle \frac{\partial }{\partial x}}\left[\widehat{H}{\displaystyle \frac{\partial \widehat{S}}{\partial x}}+h\mathit{\delta }{\displaystyle \frac{\partial s}{\partial x}}\right]-q_{\mathrm{E}},\end{array}$$

with the fluxes of freshwater and saltwater both zero at x=0, and $h=\widehat{H}=H$ at the coastline x=x_c. Moreover, we must have $\widehat{H}=H$ at the ice margin x=x_m (where $H_{\mathrm{i}}\left(x_m\right(t),t)=\mathrm{0}$) so that the water table is continuous with the subglacial region. This problem is more similar to classical seawater intrusion problems, yet it still permits multiple steady states for variable b. However, since x_c changes only through sea level changes and not through the ice sheet response to climatic forcing, the variation in x_c is less pronounced than that of x_g.

B2 Aquifer heterogeneity

When the porosity ϕ(x,z) and permeability k(x,z) are heterogeneous (but k remains isotropic), it is straightforward to derive the generalised form of Eq. (22):

$$\begin{array}{}\text{(B5)}& {\displaystyle \frac{\partial }{\partial t}}\left(h{\overline{\mathit{\varphi }}}_h\right)=K{\displaystyle \frac{\partial }{\partial x}}\left[h{\overline{k}}_h{\displaystyle \frac{\partial }{\partial x}}\left(p_S+S+\mathit{\delta }s\right)\right],\end{array}$$

where

$$\begin{array}{}\text{(B6)}& {\overline{\mathit{\varphi }}}_h={\displaystyle \frac{\mathrm{1}}{h}}\underset{b}{\overset{s}{\int }}\mathit{\varphi }\mathrm{d}z,{\overline{k}}_h={\displaystyle \frac{\mathrm{1}}{h}}\underset{b}{\overset{s}{\int }}k\mathrm{d}z.\end{array}$$

For example, ϕ may display an “Athy profile”, in which (Gustafson et al., 2022; Gooch et al., 2016)

$$\begin{array}{}\text{(B7)}& \mathit{\varphi }={\mathit{\varphi }}_{\mathrm{0}}{e}^{-\mathit{\beta }(S-z)},{\overline{\mathit{\varphi }}}_h={\mathit{\varphi }}_{\mathrm{0}}{e}^{-\mathit{\beta }H}{\displaystyle \frac{e^{\mathit{\beta }h}-\mathrm{1}}{\mathit{\beta }h}}.\end{array}$$

The permeability k could be a piecewise constant, representing different sedimentary strata (Li et al., 2022) or obey a relation k(ϕ) (Freeze and Cherry, 1979). The porosity ϕ may also be time-dependent, for example due to sediment compaction.

B3 Three-dimensional model

In reality, coastal aquifers are three-dimensional structures. The three-dimensional analogue of Eq. (22) is

$$\begin{array}{}\text{(B8)}& {\displaystyle \frac{\partial h}{\partial t}}=K\mathrm{\nabla }\cdot \left[h\mathrm{\nabla }\left(p_S+S+\mathit{\delta }s\right)\right],\end{array}$$

where $\mathrm{\nabla }=\left({\partial }_x,{\partial }_y\right)$ is the horizontal gradient. This equation should be solved on a region R, whose boundary comprises an inland ice divide ∂R_i, on which $\mathrm{\nabla }h\cdot n=\mathrm{0}$, and a grounding line ∂R_g, on which h=H. As before, the steady states are given by

$$\begin{array}{}\text{(B9)}& h=\mathrm{0}\text{or}{p}_S+S+\mathit{\delta }s=\text{const.}\end{array}$$

However, it is not as straightforward to establish a criterion for saltwater pockets analogous to Eq. (44), as this requires knowledge of the global topology of steady-state solutions. We must first find the “minimal” steady state, consider the subset of R₀⊂R on which this state has h=0, and determine whether there exists a constant C such that $p_S+S-C>\mathit{\delta }b$ on a nonempty compact subset of R₀.