Modeling mixing and melting in laminar seawater intrusions under grounded ice

Small-scale ice–ocean interactions near and within grounding zones play an important role in determining the current and future contribution of marine ice sheets to sea level rise. However, the processes mediating these interactions are simplified in large-scale models due to limited observations and computational resources, contributing to uncertainty in future projections. Previous modeling studies have demonstrated that seawater can interact with subglacial discharge upstream of the grounding zone, and recent observations appear to support this possibility. In this study, we investigate turbulent mixing of quasi-laminar intruded seawater and glacial meltwater under grounded ice using a computational fluid dynamics solver. In agreement with previous work, we demonstrate the strongest control on intrusion distance is the speed of subglacial discharge and the geometry of the subglacial environment. We show that, in the fluid regimes simulated here, and expected at ice shelf grounding zones, turbulent mixing plays a negligible role in setting intrusion distance. Basal melting from seawater intrusion produces buoyant meltwater, which may create a negative feedback by chilling and freshening near-ice water, therefore reducing further melting; however, this remains unquantified. The magnitude of modeled basal melt rates from seawater intrusion can be replicated by existing sub-ice-shelf melt parameterizations by modifying the traditionally used transfer coefficients. We conclude that, in times or places when subglacial discharge is slow, seawater intrusion can be an important mechanism of ocean-forced basal melting of marine ice sheets when considering added geometric complexities and ocean conditions.

1 Introduction

Marine-terminating glaciers in Greenland and Antarctica have experienced accelerating ice loss over the past several decades (Otosaka et al., 2023). Ocean melting of marine-terminating glaciers has driven a considerable amount of this mass loss (Depoorter et al., 2013; Rignot et al., 2013), but it is not yet accurately represented within current coupled ice–ocean models. Efforts to improve coupled models focus on improving current parameterizations of melt rates at the ice–ocean interface (Kimura et al., 2015; Middleton et al., 2022; Zhao et al., 2024; Washam et al., 2023), collecting in situ data (Stanton et al., 2013; Christianson et al., 2016; Jackson et al., 2017; Washam et al., 2020; Davis et al., 2023; Schmidt et al., 2023), and uncovering novel melt mechanisms not included in current coupled models (Rosevear et al., 2021, 2022). In this study, we focus on fluid processes that mediate one such novel melt mechanism under grounded ice.

Grounding lines are junctions between the grounded and floating portions of ice sheets. The grounding line has historically been considered to be a hydraulic barrier between the cold, fresh subglacial hydrologic system and the relatively warm, saline ocean. Models have suggested that tidal forcing may push seawater upstream of the grounding line (Sayag and Worster, 2013; Walker et al., 2013), causing tidally asymmetric melt in this “grounding zone” region (Gadi et al., 2023). Such asymmetry results in stronger melting during the ascent of high tide and weaker melting during the transition to low tide. Recent field observations have suggested that this zone is hydraulically active, with mixing occurring between the ocean and subglacial hydrology upstream of the grounding line (Macgregor et al., 2011; Horgan et al., 2013; Whiteford et al., 2022; Kim et al., 2024). Satellite observations find evidence for elevated rates of basal melt at and beyond the grounding line relative to very low values typically expected for grounded ice, contributing to retreat (Milillo et al., 2019; Ciracì et al., 2023).

More recently, seawater intrusion within and beyond grounding zones has been hypothesized to behave similarly to flow in estuaries, developing a wedge-shaped density front at which fresh glacier melt discharge flows over saline seawater. Wilson et al. (2020) and Robel et al. (2022) adapted a theoretical model for layered shallow water flows in estuaries (e.g., Krvavica et al., 2016) to demonstrate mathematically that freshwater velocity, the geometry of the subglacial environment, and the wall drag acting on the fluid all potentially exert important controls on the extent of seawater intrusion in subglacial hydrological systems. Due to the stratified nature of a salt wedge and lower fluxes, Robel et al. (2022) have hypothesized that ice loss from seawater intrusion is driven by double-diffusive convection. More recent work suggests that shear-driven melting beneath grounded ice can enlarge the subglacial cavity, enhancing seawater intrusion and potentially triggering a runaway positive feedback (Bradley and Hewitt, 2024). This feedback arises when melting exceeds ice advection, preventing upstream ice from replenishing the ablated region and allowing the conduit to grow unchecked. Bradley and Hewitt (2024) identified a regime in which seawater intrusions could become unbounded, effectively “intruding infinitely”.

Prior theories of ocean-driven melt have emphasized the importance of turbulent mixing in driving heat and salt transport towards the ice–ocean interface (Mcphee, 2008) and creating fully mixed boundary layers adjacent to the ice. Traditional melt parameterizations rest on this assumption of an ice-adjacent fully mixed boundary layer (Jenkins, 2011). However, in regions of seawater intrusion such assumptions may not hold, as a layer of buoyant freshwater lies above a dense salt wedge in a stratified environment (Wilson et al., 2020), with a boundary layer that evolves along the length of the seawater intrusion. This means the near-ice salinity goes to zero, and the equation of state changes and is therefore a distinctly different regime than the assumptions inherent in current melt parameterizations. Turbulent mixing can be a mechanism that reduces stratification by enhancing interfacial mixing between the salt wedge and freshwater layer, transporting heat and salt upward into the ice-adjacent freshwater layer. Conversely, excess buoyant forcing from basal melting induced by seawater intrusion may inhibit intrusion by increasing total freshwater discharge and strengthening stratification. Prior studies of seawater intrusion have omitted turbulent mixing in the interest of obtaining simple mathematical theories and have not considered feedbacks between intrusion-induced melting and intrusion persistence. In this study, we investigate these effects with the aid of a computational fluid dynamics solver, with three aims: (1) to test previously proposed controls on seawater intrusion distance, (2) to determine the effects of turbulent mixing on seawater intrusion, and (3) to investigate the dynamics of intrusion-induced basal melting.

2 Methods

To study the dynamics of seawater intrusion beneath grounded ice, we utilize ANSYS Fluent (ANSYS, 2022), a computational fluid dynamics (CFD) solver of the Reynolds-averaged Navier–Stokes (RANS) equations using the finite volume method. More details on the RANS equations and ANSYS Fluent are provided in Appendix B1. Using a CFD solver allows the simulations to have model resolutions on the scale of millimeters. Unlike previous studies, by using this CFD solver at such fine resolution, we are able to resolve heat and mass transfer through the entire water column, and appropriate turbulence closure schemes allow for the boundary layer to be resolved. Interfacial shear instabilities that might be expected of highly stratified flows such as the ones simulated in this study are not explicitly resolved in RANS models; however, the averaged mixing effect of such instabilities should be captured. Furthermore, a recent study using direct numerical simulations has identified that such instabilities do not arise at the lower Reynolds numbers we are simulating here (Zhu et al., 2023). ANSYS Fluent has been extensively validated across a wide range of flow geometries and conditions. The κ–ϵ turbulence closure model, widely used in computational fluid mechanics and implemented in ANSYS Fluent, has proven effective for simulating flow structures around complex bathymetry (Zangiabadi et al., 2015; Al-Zubaidy and Hilo, 2022) and sediment-laden plumes (Nguyen et al., 2020). Its applicability to multiphase modeling in geometries similar to those studied here has also been validated against experimental results (Sultan et al., 2019). Most relevant to this study, Chalá et al. (2024) demonstrated ANSYS Fluent's capability to simulate seawater intrusion in porous aquifers, showing strong agreement between experimental data and model predictions for intrusion length and shape. Additionally, ANSYS Fluent has been successfully applied to freeze desalination processes, using the volume of fluid method to model saltwater mixtures cooled at the base (Jayakody et al., 2017).

In this study, we test the effects of freshwater velocity (u_f), turbulent mixing, and subglacial geometry on seawater intrusion distance and vertical structure. The freshwater velocity is varied over 3 orders of magnitude (0.05–5 cm s⁻¹) to mimic a range of likely subglacial discharge velocities in Antarctica and Greenland in non-summer months (Carter et al., 2017; Davis et al., 2023; Washam et al., 2020) and previous experimentally tested speeds (Wilson et al., 2020). The corresponding Reynolds numbers for the given geometry and freshwater flux are 25, 250, and 2500 for the low, medium, and fastest freshwater cases presented. In the experimental setup in this study, we only consider quasi-laminar flow in order to facilitate comparison to previous studies (Wilson et al., 2020; Robel et al., 2022), which find that subcriticality (Fr<1) is required for intrusion development. To test the role of turbulent mixing, we varied a mixing parameter across three values along with simulating laminar flow cases. The simulations in which we range over subglacial discharge velocities with medium turbulent mixing (e.g., C_μ=0.09) are our reference simulations to which we will compare all other results. The purpose of the reference simulations is to compare with previous theory that considers quasi-laminar flow and no melting effects. Further testing on potential intrusion control variables is also considered and explored in detail. More discussion on turbulence modeling is given in Sect. 2.3.

In evaluating the geometric effect on intrusion dynamics, we tested both a retrograde slope with an angle of 0.5° relative to the horizontal and a thicker subglacial environment. We also simulated melting induced by seawater intrusion to investigate how this secondary source of buoyancy affects intrusion persistence and structure. Each simulation is initialized with a warm, salty ocean basin (S=30 ppt and T=0.5 °C) and a fresh, cold subglacial environment (S=0 ppt and $T=-\mathrm{0.73}$ °C) as shown in Fig. B1. The thermal driving associated with the initial conditions is approximately 1.2°C, and the reduced gravity is 0.23 m s⁻². The transient solver is then run for 12 h at 5 s time steps, with data exported every 20 time steps (100 s). A quasi-steady state is reached by all simulations, as evaluated by the time change of the average density of the subglacial space being less than 10⁻⁴ kg m⁻³. All results presented are time-averaged values from when the quasi-steady state is reached. A summary of the experimental setup, parameters, and key results is given in Tables 1, 2, and 3.

2.1 Domain and boundary conditions

We consider a two-dimensional subglacial domain, encompassing one vertical and one horizontal (orthogonal to the local grounding line) dimension. The domain is akin to an unbounded freshwater sheet that meets the ocean at a specified discharge point, representing the grounding line. Since we do not resolve ice dynamics, we prescribe the grounding line as a vertical ice face in the domain boundary geometry, instead of including an ice shelf to reduce the domain size needed in these simulations and limit geometric constraints to an idealized subglacial water sheet. The geometry of the bounding surfaces in this configuration does not change in time, with the vertical ice front chosen to limit the geometric influences on intrusion distance to only within the subglacial environment. Tides are not considered in this study, which would temporally alter the geometry of the subglacial environment and therefore be another factor influencing intrusion distance. The underlying bedrock is impermeable, and the subglacial environment between the ice and the bedrock is not obstructed by obstacles, both of which would introduce additional controls on the intrusion distance (Robel et al., 2022).

Figure 1 depicts the standard model configuration we use in this study. There are two velocity inlets common to all simulations: a seawater source at the inlet boundary of the tall ocean basin (dark-blue arrow in Fig. 1) and a freshwater source at the inlet boundary of the subglacial environment (light-blue arrow in Fig. 1). A pressure outlet boundary (red arrow in Fig. 1) is prescribed in the ocean basin employing a zero-gradient flux at the boundary and ensuring mass conservation in the model. In describing the results, we utilize the convention of downstream being towards the ocean basin and the upstream being towards the freshwater inlet. The seawater inlet velocity is prescribed as u_o=0.5 cm s⁻¹ across all cases, acting as a sustaining source of saline warm water into the model domain. Simulations with varied u_o (Table B1) indicate that the seawater inlet speed does not have a qualitative influence on the seawater intrusion distance or vertical structure over a range of relatively weak ocean current speeds that we consider appropriate for the constrained ocean cavity near the grounding line (0.05–5 cm s⁻¹), so we set u_o=0.5 cm s⁻¹ for all simulations going forward.

Finally, vertical and horizontal ice wall boundaries (where the ice is in contact with the fluid domain) are defined with characteristics that mimic a grounding line environment. The ice wall boundaries have a pressure- and salinity-dependent thermal boundary condition represented by the liquidus condition:

$$\begin{array}{}\text{(1)}& T_{\mathrm{b}}=S_{\mathrm{b}}{\mathit{\lambda }}_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{2}}+z_{\mathrm{b}}{\mathit{\lambda }}_{\mathrm{3}}.\end{array}$$

Here, λ₁, λ₂, and λ₃ are constants, and the boundary salinity is S_b. The depth of the ice is equal to z_b; in these simulations we set this to be 1000 m. The boundary salinity, S_b, is the salinity of the cell filled with water nearest to the ice face and is permitted to evolve during the simulation via the evolution in Eq. (B12). For the reference simulations, the boundary salinity evolves in time and space due to advection and diffusion of the subglacial discharge and seawater intrusion. In the simulations with melting, an additional source of freshwater is injected into the near-boundary grid cells, actively freshening the boundary salinity. Both the vertical and horizontal ice boundaries have a no-slip kinematic condition in the non-melting cases, forcing the freestream fluid velocity to be zero at the ice wall. The vertical ice front in this configuration resembles a tidewater glacier and is utilized to restrict the geometric controls on intrusion distance to those within the subglacial environment. Having a low-sloping ice-shelf bottom would introduce further constraints on the ability for seawater to intrude beyond the grounding line since intrusion distance is a function of the height of its environment (Robel et al., 2022; Wilson et al., 2020). Such a configuration would also require computational domains much larger than are feasible to simulate at the resolution needed to properly resolve the ice–water boundary layer. Despite the tidewater-like geometric configuration, the low subglacial discharge fluxes make the fluid domain appropriate for simulating conditions expected in Antarctica year-round or Greenland in non-summer months.

The “subglacial environment” is the domain upstream of the grounding line (x>0 in Fig. 1). Since the bedrock and ice-adjacent boundary layer thicknesses depend on the freestream fluid velocities, the mesh resolution changes for each freshwater velocity to accurately model near-wall processes with the chosen turbulence closure scheme. For the tested freshwater velocities, the vertical domain size hinders the development of a full boundary layer. Instead, everywhere in the domain, the fluid feels the effects of the wall boundary. Further discussion on domain and meshing is included in Sects. B4 and B5.

Figure 1 Schematic of the domain used in ANSYS Fluent. Ice block is not simulated and used to graphically depict the vertical and horizontal (red solid line) ice–ocean interfaces. The length and height of the subglacial environment vary between simulations. The ocean basin is 2 m wide by 5 m tall for all simulations. The dark-blue arrow on the left represents seawater input (u_o), while the light-blue arrow on the right represents freshwater input (u_f). The solid red arrow represents the zero-gradient flux boundary outlet. The red dashed arrows represent the meltwater input ($\dot{m}$) that is turned on for the simulations with melting as a function of near-wall temperature (Eq. 4).

2.2 Salt and heat transport

In addition to solving the RANS equation for fluid velocities, we configure the CFD solver to calculate the concentration of salt with a “species transport model” (advection–diffusion equations) (ANSYS, 2009). Salt is therefore transported as an active tracer within the fluid domain. To ensure mass transport within the computation domain is realistic and physical, there are two velocity (hence mass) inlets (i.e., the subglacial discharge and ocean inflow) in the non-melting case and three velocity inlets (subglacial discharge, the melting horizontal ice face, and the ocean inflow) in the melting case. Mass is conserved via a pressure outlet (red arrow in Fig. 1), which employs a zero-gradient flux boundary condition. At this boundary, the mass outflow rate is not specified and is determined as part of the numerical solution based on the requirement that all flow variables have zero gradients in the direction normal to the boundary. This kind of arrangement is typically used to emulate fluid flows in an infinite domain as in our case where subglacial channel discharge is released at the grounding line into an ocean with infinite extent. Energy, and therefore fluid temperature, is evolved via an energy conservation equation employed by the CFD solver resolving advection, conduction, and salt diffusion (ANSYS, 2009). Conduction represents heat transfer due to horizontal thermal gradients. Since salt is tracked as an active tracer, it transports enthalpy associated with its specific heat and concentration gradients. This must be accounted for in the energy equation Fluent solves. Further details on heat and salt transport in the model can be found in Appendix B3 alongside a discussion of the RANS formulation in Appendix B1.

The seawater inlet is prescribed with 30 ppt salinity and T_o=0.5 °C, while the freshwater inlet is prescribed with zero salinity and the corresponding pressure-dependent freezing point from Eq. (1). Seawater temperature is chosen to represent warm cavity Antarctic conditions (Middleton et al., 2022; Kimura et al., 2015). Seawater salinity is set to generate a water mass with density characteristic of the Southern Ocean. The species transport model also calculates water density (ρ_w) with a prescribed linear equation of state from Roquet et al. (2015):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{\rho }}_{\mathrm{w}}& =\mathrm{1000}+S_{\mathrm{w}}\cdot (\mathrm{0.7718}\cdot \mathrm{1000})\\ & +(T_{\mathrm{w}}-\mathrm{273.15})\cdot (-\mathrm{0.1775}).\end{array}\end{array}$$

The salinity, S_w, and temperature, T_w, of the fluid are found at the center of each cell. This simplified linear equation of state was chosen for ease of implementation as a user-defined function within ANSYS Fluent and with the assumption that cabeling or thermobaricity does not play a role in this environment. A higher-order equation of state was tested to compare the linear formulation used here. While the fluid density is different, the salinity and temperature of the fluid remain unchanged since they have their respective transport equations. Further discussion on choice and justification of the linear equation of state is presented in Appendix B8.

2.3 Turbulence closure scheme

Turbulent mixing can increase the exchange of heat and salt between seawater and fresh subglacial discharge, homogenizing the water column and potentially reducing the ability for seawater to intrude. The theory of Wilson et al. (2020) assumes no mixing between the two water masses, treating the interfacial drag (which dictates shear) as a free parameter dependent on fluid flow. Robel et al. (2022) treat interfacial drag as negligible to derive a closed-form prediction for intrusion distance. However, under scenarios of fast subglacial discharge, the experiments described in Wilson et al. (2020) exhibited interfacial mixing with the formation of wave crests at the top of the salt wedge. To understand the role turbulent mixing has in a seawater intrusion regime, we enable the commonly used κ–ϵ two-equation turbulence closure scheme (Mansour et al., 1989; Launder and Spalding, 1983). A closure scheme is necessary because averaging the Navier–Stokes equations introduces Reynolds stresses due to turbulent motion within the fluid. Reynolds stresses take the form $\left(\stackrel{\mathrm{‾}}{u_{\mathrm{i}}^{\prime }{u}_j^{\prime }}\right)$, the averaged product of turbulent velocity fluctuations. The closure scheme used here employs the turbulent–viscosity hypothesis, which relates the deviatoric Reynolds stresses to the mean strain rate via a positive scalar eddy viscosity (Pope, 2000). The κ–ϵ closure scheme solves for the eddy viscosity with

$$\begin{array}{}\text{(3)}& {\mathit{\mu }}_t={\mathit{\rho }}_{\mathrm{w}}{C}_{\mathit{\mu }}{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{\mathit{ϵ}}},\end{array}$$

where κ represents the turbulent kinetic energy, ϵ represents dissipation, and ρ_w is the fluid density. This form of turbulence closure provides an avenue to manipulate the degree of turbulent mixing by modulating the parameter C_μ. This is a model parameter that dictates the amount of turbulent transport (mixing) given some ${\mathit{\kappa }}^{\mathrm{2}}/\mathit{ϵ}$. Larger values of C_μ imply that for a given level of turbulent kinetic energy and dissipation, turbulent stresses will be increased due to increased eddy viscosity. We vary this parameter across a factor of 4 (0.045, 0.09, 0.18), which allows us to explore a range of possible turbulent mixing while maintaining model stability. The middle value, 0.09, is commonly adopted as the “standard value” and is derived from experiments with equilibrium shear flows, e.g., the log–law region and above in pipe flows. We employ the Yang–Shih low-Reynolds formulation of the κ–ϵ closure scheme (Yang and Shih, 1993), which uses damping equations near wall boundaries to adequately resolve the viscous sublayer, allowing for finer resolution near the ice boundaries and complete boundary layer resolution (Hrenya et al., 1995; ANSYS, 2009). Further discussion on the ANSYS Fluent model and formulation of the turbulence closure model can be found in Appendix B2.

Table 1 Simulations with varied freshwater discharge velocity and degree of turbulent mixing for flat geometry with H=5 cm.

2.4 Basal melting

In some simulations, we also simulate the added buoyancy flux resulting from intrusion-forced melting. Using the last time step from the corresponding non-melting simulation as the initial condition, we enable melting at the horizontal ice wall by changing it from a wall boundary condition to a velocity inlet (red dashed arrows in Fig. 1). In ANSYS Fluent, a wall boundary condition holds a constant temperature while providing a shearing force on the fluid in accordance with the no-slip wall boundary condition. When changed to a velocity inlet, the inflowing melt can have a prescribed temperature and salinity but does not provide any source of shear. Therefore, the melting cases have a free shear kinematic boundary condition. To address this model limitation, we also conducted experiments without melting and with a free-slip kinematic condition, for comparison with the intrusions in the melting experiments. The temperature of the inflowing meltwater is set by Eq. (1). Melting is only turned on for the horizontal ice face (red boundary in Fig. 1) and not the vertical ice face to isolate the added buoyancy effects from seawater intrusion forced melt only. The vertical ice face is susceptible to plume-driven melt from the buoyant subglacial discharge and is distinctly different than the seawater intrusion melt domain. The downward fluid velocity prescribed at the horizontal ice face is set by the melt rate, $\dot{m}$, and is a function of the difference between the near-wall cell's centroid temperature T_w and ice–ocean interfacial temperature T_b, thermal conductivity κ_T, and density of the ice ρ_i.

$$\begin{array}{}\text{(4)}& \dot{m}={\displaystyle \frac{{\mathit{\kappa }}_{\mathrm{T}}}{{\mathit{\rho }}_{\mathrm{i}}{L}_{\mathrm{i}}}}{\displaystyle \frac{T_{\mathrm{w}}-T_{\mathrm{b}}}{\mathrm{0.5}{H}_{\mathrm{c}}}}.\end{array}$$

The thermal forcing is divided by half of the near-wall cell height, H_c, to obtain the near-wall thermal gradient, $\partial T/\partial y$. L_i is the latent heat of ice. The values of all constants are presented in Tables A1 and A2. This framework represents the conservation of heat at the ice–ocean interface, which varies along the horizontal ice face as the near-wall thermal gradient changes due to seawater intrusion and vertical mixing. We can use Eq. (4) to calculate melt rates instead of a parameterization because we resolve the boundary layer directly and do not need to make assumptions about how heat and salt are transported through the boundary layer. Inherent in this is the assumption ANSYS Fluent accurately simulates all appropriate boundary layer transport processes necessary to get heat from the freestream fluid flow to the cell grid next to the ice face, and any melting experienced is due to the conservation of heat at the ice-adjacent cell. In setting a vertical velocity to mimic a moving interface, we introduce additional sources of momentum to the fluid. The vertical velocity arising from the meltwater inlet is small relative to the main flow and is therefore negligible, but the buoyancy that the meltwater brings into the fluid domain is an integral part of the heat-driven melting process. The input of fresh, cold water due to melting introduces a buoyancy flux into the domain, due to density differences between the meltwater and intrusion. An alternative to this formulation would be to replace κ_T with the product of thermal diffusivity (K_T), seawater density (ρ_w), and seawater heat capacity (c_w), allowing for varying seawater density to affect heat transfer to the boundary. However, back-of-envelope calculations show small variances in density (∼20 kg m⁻³) lead to small changes (<5 %) in the melt rate. Therefore, using a constant thermal conductivity is appropriate. Note the boundary does not move over time (as in Bradley and Hewitt, 2024), and an evolving geometry of the subglacial space is not tested in this work. This choice greatly simplifies the computational domain and considerations of meshing with turbulent closures.

Table 2 Simulations with a melting ice interface and simulations with no melting and a free-slip kinematic ice boundary condition.

3 Results

3.1 Characteristics of the seawater intrusion

In all simulations, warm seawater intrudes some distance beyond the defined grounding line (x=0 m in Fig. 2). A strong control on seawater intrusion distance is the freshwater discharge velocity, in line with previous work (Wilson et al., 2020; Robel et al., 2022; Krvavica et al., 2016). The simulation with the lowest flux of freshwater (Fig. 2c) has approximately a 25 m intrusion, while the fastest flux experiences only about 0.5 m of intrusion (Fig. 2b). This range of intrusion distances demonstrates a weaker dependence on freshwater velocity than suggested by Robel et al. (2022) where intrusion distance has an inverse quadratic dependence on freshwater velocity and therefore should vary by a factor of 1000 in response to the range of input velocities tested here. Turbulent mixing, as modulated by C_μ, affects intrusion distance to a lesser degree than freshwater discharge velocity when varied over a wide range encompassing likely values on the lower end for realistic estuarine-like mixing rates (Geyer et al., 2000, 2008). For the middle freshwater velocity (Fig. 2a), increased turbulent mixing slightly reduces the intrusion distance. To contrast the effects of turbulent mixing, we tested laminar flow cases with no turbulent mixing (green transects in Fig. 2) and saw no meaningful difference in intrusion distance for u_f=0.5 cm s⁻¹. For the low and high freshwater velocity cases, including turbulent mixing (blue transects in Fig. 2b, c) increases intrusion distance when compared to the laminar test case (green transects in Fig. 2b, c). For the slowest freshwater case, the intrusion is reduced in the absence of turbulent mixing due to the lack of entrainment between the seawater intrusion and the buoyant subglacial discharge. When modeled, this entrainment extends the intrusion by generating a tail of relatively low-density water. This is not to say that turbulent mixing is unimportant in the dynamics of seawater intrusion but rather that intrusion distance is not strongly sensitive to the strength of turbulent mixing (over the range of discharge velocities and C_μ values considered realistic for subglacial and estuarine environments), particularly when compared to other factors such as the geometry of the subglacial environment and freshwater discharge flux. In the flow regimes simulated here, the turbulent viscosity does not get large enough to greatly affect the flow dynamics as shown in Fig. D1. It may be that, for much higher discharge velocities (𝒪 (m s⁻¹)) encountered at times of high subglacial discharge, turbulent mixing plays a more important role than the cases considered here. However, under those conditions, the freshwater flux is likely to be supercritical and prohibit any intrusion development as predicted by theory in Wilson et al. (2020) and Robel et al. (2022).

Figure 2 Depth-averaged density transects along the subglacial environment. The transect termini represent the intrusion distance or the point at which density is less than or equal to the density of freshwater at the pressure-dependent freezing point. Blue lines represent simulations without melting. Dashed lines represent higher turbulent mixing (C_μ=0.18), dotted lines are low turbulent mixing (C_μ=0.045), and solid lines are medium turbulent mixing (C_μ = 0.09). Dashed red transects are for simulations without melting and have free-slip kinematic boundary conditions on the horizontal ice face. The solid red transects are simulations with melting and medium turbulent mixing. The orange and magenta transects in panel (a) are scenarios with different geometries, a taller subglacial environment (H=7.5 cm) and a retrograde slope (θ=0.5°), respectively. The green transects represent laminar flow cases with no turbulent mixing. Shading depicts the first temporal standard deviation or the maximum/minimum value depending on the smallest absolute difference with the average value. Note the varying x axis across the panels.

Table 3 Simulations with varied geometry.

Changing subglacial geometry has a large effect on intrusion distance as demonstrated by the experiment with a taller subglacial channel (H=7.5 cm) plotted as an orange transect in Fig. 2a. Increasing the channel height by 50 % increased the intrusion distance by nearly a factor of 3. This indicates an even stronger sensitivity to the height of the subglacial opening in these more realistic simulations than predicted by Robel et al. (2022), who find a quadratic dependence on the subglacial conduit height. We also tested the effects of a retrograde bed slope (θ=0.5°) on intrusion distance, which increased the intrusion distance by about a factor of 1.5 relative to the flat cases tested (magenta transect in Fig. 2a). We did not find any evidence of an unbounded increase in the intrusion distance under these retrograde slopes, as predicted by Wilson et al. (2020) and Robel et al. (2022). That being said, our finite domain length limits the intrusion distances achievable in this model configuration.

Figure 3 Time-averaged vertical profiles of thermal forcing (a, e, i), salinity (b, f, j), x component of velocity (c, g, k), and buoyancy frequency (d, h, l) along the seawater intrusion for each freshwater velocity tested in the standard geometry with medium turbulence C_μ=0.09. The fraction of intrusion distance represents evenly spaced vertical slices at relative distances along each intrusion. The top row (a–d) is for the fastest freshwater velocity, u_f = 5 cm s⁻¹. The middle row (e–h) is for the middle freshwater velocity, u_f=0.5 cm s⁻¹. Lastly, the third row (i–l) is for the slowest freshwater velocity, u_f=0.05 cm s⁻¹.

The time-averaged vertical profiles of temperature, salinity, and velocity along the intrusion for simulations without basal melting (Fig. 3) depict a strongly stratified two-layered exchange flow, with variable thermohaline gradients for each freshwater flux tested. The fastest freshwater flux (Fig. 3a–d) has the strongest stratification and fastest exchange flow of all the cases tested. For all cases, the height and slope of the thermocline and halocline decrease upstream into the subglacial space. The strength in stratification and exchange flow decreases with decreased freshwater flux. This trend holds for the retrograde slope case and increased channel height case (Fig. E1). The retrograde slope has a nearly identical vertical structure to the similar flat channel test case, while the simulation with increased channel height sees a faster exchange flow due to an increased flux of freshwater. The system decays to a classic pipe flow further upstream into the subglacial space, out of the regime of intrusion. The degree of vertical mixing increases with decreasing freshwater velocity, which is similar to observed mixing dynamics in estuaries (Montagna et al., 2013). For high freshwater fluxes, a well-developed salt wedge forms and high stratification persists (Fig. 3a–d). For the lowest freshwater flux, the diffusive and advective timescales are of the same order, which leads to homogeneity in the vertical distribution of heat and salt (Fig. 3i–l). In estuaries, strong tides generate shear at the bed, increasing turbulence and vertical mixing (Montagna et al., 2013). Tides may play a similar role in subglacial seawater intrusion regimes; however, that is not simulated here.

Buoyancy frequency, N, is a measure of the degree of stratification in the water column:

$$\begin{array}{}\text{(5)}& N=\sqrt{{\displaystyle \frac{-g}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}y}}}.\end{array}$$

Higher N indicates stronger density gradients and stratification, tracking the presence of seawater intrusion. Buoyant forcing, as represented by N, suppresses mixing in the presence of a flat horizontal ice boundary (as we have here), competing with the velocity shear and turbulence that drive mixing over the length of the intrusion. In regimes of seawater intrusion, buoyancy frequency is large due to the large vertical density gradients and small vertical length scale (Fig. 3d, h, l). Buoyancy frequency is largest at the grounding line and decays to zero upstream of the intrusion. This strong stratification can act to “shield” the ice from melting, generating a near-ice layer of fresh, cold water. However, the horizontal density gradient introduced by the characteristic wedge shape of seawater intrusion will drive vertical baroclinic convective motion to flatten isopycnals. Such baroclinic adjustment may be an important source of interfacial mixing, working in tandem with turbulence and double-diffusive convection to reduce stratification within the subglacial environment. This convective-driven mixing mechanism differs from convective mixing caused by a sloping ice boundary, in which a buoyant plume may form. For the retrograde slope geometry tested here, there is a slight reduction in stratification strength relative to the comparable flat geometry simulation (Fig. E1d) likely attributed to natural convective mixing. Where subglacial openings have complex geometry, we anticipate buoyant-driven convection from sloping ice boundaries to aid in driving mixing on small scales, which will reduce the stratification from the intrusion.

Figure 4 Drag coefficient for each freshwater velocity tested with medium turbulence (C_μ=0.09). The dashed line represents the middle freshwater velocity with a retrograde slope, and the dotted line is for the same u_f but with a taller subglacial environment. The scattered points represent the intrusion distance for each case.

Drag from the ice shears fluid flow in the subglacial environment to have zero velocity at the ice, per the no-slip kinematic boundary condition. This shearing determines boundary layer thickness and therefore influences heat and salt transport towards the ice. Within the model framework here, the wall drag coefficient is not a free parameter to be set but rather one diagnosed from the simulations via the following relationship (Pope, 2000):

$$\begin{array}{}\text{(6)}& C_{\mathrm{d}}={\displaystyle \frac{\mathrm{2}\mathit{\nu }}{u^{\mathrm{2}}}}\left({\displaystyle \frac{\partial u}{\partial y}}\right),\end{array}$$

where $\partial u/\partial y$ is the velocity gradient normal to the wall and is evaluated with a linear fit above the peak in velocity of the upper freshwater layer. The freestream current speed, u, is obtained from the local centerline flow. The kinematic viscosity is represented by ν. The full derivation of calculating the drag coefficient is given in Appendix C. In previous work (Wilson et al., 2020; Robel et al., 2022), C_d is set to values around $\sim {\mathrm{10}}^{-\mathrm{3}}$ derived from observed drag coefficients under sea ice; however, here we diagnose C_d's in the range of 10⁻² to 10¹ (Fig. 4) when a no-slip kinematic boundary condition is applied. The high drag coefficients simulated here upstream of the intrusion regime are in line with the expected values for laminar flows with these Reynolds numbers. However, over the regime of intrusion, the drag coefficient increases by nearly an order of magnitude in all cases tested. The increased drag coefficients over the intrusion are more difficult to estimate and likely arise from enhanced turbulence in the interfacial shear layer between the intrusion and buoyant subglacial outflow. Where intrusion occurs, the near-wall velocity gradient grows (darker profiles in Fig. 3c, g, k), increasing the shear velocity and drag. This velocity gradient is an important driver of heat flux to the ice, as discussed in later sections. We also simulated a free-slip ice boundary condition, meaning no drag at the ice boundary (i.e., C_d=0). This resulted in the largest change in intrusion distance (red dashed transects in Fig. 2) and in some cases allowed for seawater to fill the whole simulation domain. This indicates incredible sensitivity of intrusion distance to the applied kinematic boundary conditions. The high drag coefficients simulated here likely explain the much smaller intrusion distances simulated here compared to prior studies (Wilson et al., 2020; Robel et al., 2022) and are discussed in more detail in Sect. 4.

Figure 5 Simulated melt rates for all freshwater velocities with medium turbulent mixing alongside the retrograde slope simulation and taller subglacial geometry simulation. The dots indicate the respective intrusion distance. Two versions of parameterized melt rates are plotted for the middle freshwater velocity. The dashed line is the melt rate using a constant “best-fit” turbulent transfer coefficient in the S99 parameterization (S99_a). Alternatively, the dotted line uses a variable turbulent transfer coefficient calculated in Eq. (8) (S99_b).

3.2 Dynamics of intrusion-induced melt

Intrusion distance does not vary significantly in simulations where basal melt of ice is included (red lines in Fig. 2), when compared to simulations with no-slip kinematic boundary cases (blue transects in Fig. 2). However, in simulating an interface that allows for an added freshwater source from “melting”, we lose the shearing effect of the boundary. Therefore, we find it appropriate to compare the intrusion distances of the simulations with melting to those without melting that have a free-slip boundary imposed (dashed red transects in Fig. 2). Under this comparison, nearly all intrusion distances simulated with melting are reduced by over 50 % relative to their corresponding without melting free-slip simulations (Table 2).

The highest simulated melt rate occurs in the sloped subglacial environment with u_f=0.5 cm s⁻¹, peaking at approximately 90 m yr⁻¹. The taller subglacial environment under the same u_f follows closely, exhibiting a slightly lower peak melt rate but sustaining elevated melting over an extended region. By contrast, the standard flat geometry with u_f=0.5 cm s⁻¹ experiences a peak melt rate of about 60 m yr⁻¹. In scenarios with retrograde slopes and thicker subglacial environments, the area of ice experiencing significant melt is considerably larger than in the standard flat configuration. These results indicate that the sloped and taller geometries not only expand the region affected by seawater intrusion-induced melting but also increase the peak melt rate by nearly 50 %.

In simulations with the fastest freshwater discharge (u_f=5.0 cm s⁻¹) and standard flat geometry, melting is negligible due to minimal intrusion and strong stratification. Conversely, for the slowest discharge (u_f=0.05 cm s⁻¹), the melt rate is not the highest but extends over the largest region, suggesting a trade-off between peak melt intensity and spatial extent.

The melt distribution aligns with the intrusion wedge, with maximum melting concentrated near the grounding line and tapering to negligible rates of millimeters per year at approximately half the intrusion distance. This secondary source of freshwater discharge from melting at the horizontal ice–water interface does not substantially alter the vertical structure of the flow regime. The stratification strength and exchange flow within the intrusion are primarily governed by the magnitude of freshwater discharge, which outweighs the momentum contribution from meltwater input.

3.3 Parameterization of basal melting in regions of seawater intrusion

Large-scale ocean and ice sheet models typically have grid resolutions much coarser than what is necessary to resolve the ice–ocean boundary layer and directly calculate heat and salt fluxes from the ocean towards the ice. In practice, parameterizations are used in large-scale models to approximate the melt rate based on ocean temperatures and salinities modeled outside the boundary layer (10¹–10⁴ m) from the ice. However, such parameterizations have not previously been tested in the flow regimes relevant to seawater intrusion below grounded ice. Here, we simulate intrusion-induced basal melt by employing the conservation of energy (Eq. 4) at the horizontal ice boundary and compare these model results directly to a traditionally used parameterization scheme. In our simulations, we can accurately capture relevant boundary layer heat and salt transport processes and calculate basal melt rates because our model resolution of millimeters directly resolves the near-ice boundary layer.

The most widely used parameterization for modeling ocean-induced ice melt under ice shelves is the three-equation parameterization of Holland and Jenkins (1999) (referred to as S99), which assumes heat and salt transport occurs via shear-driven mixing. Here, we utilize model output to obtain the boundary temperature, T_b, which is constrained by the boundary salinity as described in the methods section. The water temperature, T_w, is set to the temperature at the center of the subglacial space. The melt rate, as parameterized by S99, is therefore given by

$$\begin{array}{}\text{(7)}& {\dot{m}}_{\mathrm{T}}{L}_{\mathrm{i}}{\mathit{\rho }}_{\mathrm{i}}={\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{w}}{\mathit{\gamma }}_{\mathrm{T}}(T_{\mathrm{w}}-T_{\mathrm{b}})-{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}{\mathit{\kappa }}_{\mathrm{i}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{i}}}{\mathrm{d}y}}.\end{array}$$

Here, L_i is the latent heat of ice, ρ_w is the seawater density, and c_w is the heat capacity of seawater. The product of the vertical temperature gradient of ice $\mathrm{d}{T}_{\mathrm{i}}/\mathrm{d}y$, ice density ρ_i, the heat capacity of ice c_i, and the thermal conductivity of ice κ_i represents heat conduction into the ice. All values are presented in Table A2. We set the ice vertical temperature gradient to zero since we do not model any heat conduction within the ice and instead assume all energy is used to cause melting. The turbulent transfer velocity, γ_T, describes the transfer of heat across the outer portion of the boundary layer and into the viscous sublayer adjacent to the ice. This transfer velocity is further parameterized via

$$\begin{array}{}\text{(8)}& {\mathit{\gamma }}_{\mathrm{T}}={\displaystyle \frac{u_{*}}{\mathrm{2.12}\ln \left(u_{*}h{\mathit{\nu }}^{-\mathrm{1}}\right)+\mathrm{12.5}{\mathit{Pr}}^{\mathrm{2}/\mathrm{3}}-\mathrm{9}}},\end{array}$$

where h is the distance from the ice of the chosen reference for T_w, usually taken to be the thickness of the boundary layer (Holland and Jenkins, 1999). The Prandtl number (Pr) is the ratio of viscous forces to diffusive forces for temperature (Kader and Yaglom, 1972; McPhee et al., 1987). The kinematic viscosity is represented by ν, and the shear velocity, u_*, is defined as

$$\begin{array}{}\text{(9)}& u_{*}=\sqrt{\mathit{\nu }{\displaystyle \frac{\partial u}{\partial y}}}\end{array}$$

(Schlichting and Gersten, 2016). However, the shear velocity is generally solved as a quadratic function of the wall drag coefficient using a form of Eq. (6). This means shear-driven parameterized melt rates are sensitive to the choice of drag coefficient (C_d). In this study, C_d is found to range 10⁻² to 10¹ with variability of over a magnitude in the regime of intrusion, which is larger than the observed range for sea ice of 10⁻³ to 10⁻² (Mcphee, 1980; Randelhoff et al., 2014; Brenner et al., 2021). In using shear velocity to drive the parameterization, this framework assumes the system has enough momentum to dictate boundary layer transport processes. When momentum is low, other processes like diffusion or natural convection may contribute to setting the complex boundary layer structure and vertical thermohaline gradients. For the weak flow scenarios tested here, we anticipate shear-driven turbulence to be relatively unimportant (Rosevear et al., 2022).

In this study, we assessed the S99 parameterization for a simulation with melting with u_f=0.5 cm s⁻¹ using two approaches. In the first approach, we evaluated the turbulent transfer coefficient (γ_T) using the relationship defined in Eq. (8) (dotted line in Fig. 5). With these values of γ_T, the S99 parameterization yielded peak melt rates approximately half those simulated directly in our study, underestimating melt in the downstream region of the intrusion. The γ_T values derived from this approach ranged from 10⁻⁵ m s⁻¹ within the intrusion to a background value of 10⁻⁷ m s⁻¹ upstream of the intrusion. In the second approach, we computed the value of γ_T that minimized the root mean squared error (RMSE) between parameterized and simulated melt rates. The optimal γ_T, corresponding to the best agreement with simulated melt rates, was approximately $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ (dashed line in Fig. 5). This optimized turbulent transfer coefficient underestimated peak melt rates in the downstream region of the intrusion while overestimating melt rates upstream near the intrusion nose. Both the tuned and calculated versions of γ_T produced melting along the entire length of the intrusion. In contrast, the simulated melt rates were concentrated along the first half of the intrusion, with a marked drop-off in melting further downstream before the intrusion terminated. Neither parameterization approach was able to replicate this downstream reduction in melt rates in the simulations.

Values of γ_T from our simulations generally fall on the lower end of the range reported in the literature, reflecting the low velocities and thin geometries tested (i.e., small Reynolds numbers). However, the corresponding thermal Stanton numbers ($\sqrt{C_{\mathrm{d}}}{\mathrm{\Gamma }}_{\mathrm{T}}$) are in closer agreement within the range found in the literature (Holland and Jenkins, 1999; Jenkins et al., 2010; Washam et al., 2023). Thermal Stanton numbers represent the ratio of heat transfer to the thermal capacity of a fluid and indicate the balance of these two processes (Jenkins et al., 2010), whereas transfer velocities represent the efficacy of heat transport via turbulent and molecular processes throughout the boundary layer. Stanton numbers are often used as tuning coefficients when the boundary layer cannot be resolved. The thermal Stanton number is equivalent to ${\mathit{\gamma }}_{\mathrm{T}}/u$, and the back-of-envelope calculations for the middle freshwater velocity where γ_T is $\sim {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ and upper layer velocity within the intrusion is ∼4 cm s⁻¹ give $\sqrt{C_{\mathrm{d}}}{\mathrm{\Gamma }}_{\mathrm{T}}\sim \mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{4}}$.

Three distinct issues make it difficult to apply existing parameterizations in regions of seawater intrusion: (1) stratification, (2) the interaction of two boundary layers, and (3) change in water flow direction near the ice boundary. Most melt parameterizations assume a well-mixed and fully developed boundary layer, with reference temperature and salinity taken beyond or at the boundary layer's edge. However, in the simulations presented here, some degree of stratification exists due to insufficient boundary shear mixing between the upper layer of subglacial discharge and the lower layer of seawater intrusion. Intrusion beyond a grounding line entails fluid flow between two boundaries, which is intrinsically different than the geometries considered for other ocean-induced melting regimes (flow bounded by a singular wall on one side). The upper fresh glacial water will have a boundary layer with the ice, and the lower saline ocean water will generate a boundary layer with the bed. In the subglacial environment, the opposing fresh and saline flows meet and create an interfacial boundary with zero velocity. The sharp transition in opposing flows and strong interfacial shear is associated with two accelerating fluid layers that are increasing drag on their respective boundaries. Here, where subglacial environments have thicknesses of 5 or 7.5 cm, and freshwater velocities ranging from 5 to 0.05 cm s⁻¹, a complete and stable boundary layer never develops, indicating the entire subglacial domain feels the effect of at least one boundary. This means multiple transport processes (turbulent and viscous) operate at the same relative importance, rather than one mechanism dominating the other. These characteristics of seawater intrusion within narrow gaps violate the assumptions inherent in traditionally used parameterizations which rely on well-mixed and fully developed boundary layers (e.g., Holland and Jenkins, 1999).

4 Discussion

In the context of previous work, these simulations confirm that freshwater flux and the geometry of the subglacial environment are both strong controls on seawater intrusion. Our simulated intrusions follow the general trend and scale sensitivity to those identified in previous laboratory experiments (Fig. 6) (gray markers; Wilson et al., 2020), which are within a factor of 10 of the theoretical prediction (dashed line) from Robel et al. (2022). Previous studies estimate possible seawater intrusion distances of kilometers to tens of kilometers (Wilson et al., 2020; Robel et al., 2022), which give many more orders of magnitude difference than any of the intrusions simulated in this study. This difference is likely due to the drag coefficient at the ice wall (C_d), which field studies of ocean flow under sea ice (Mcphee, 1980; Randelhoff et al., 2014; Brenner et al., 2021) and geometric parameterizations (Lu et al., 2011) estimate to be of the order of 10⁻³ to 10⁻². In this study, the drag coefficient of the wall cannot be prescribed but rather is an emergent property arising from the no-slip kinematic boundary condition and momentum dissipation within the model from the κ–ϵ closure scheme. Calculating the drag coefficient using model output gives C_d with values of the order of 10⁻² to 10¹ with variability of over a magnitude within the intrusion regime. The analytical theory of intrusion distance (L) for an unobstructed water sheet from Robel et al. (2022) is

$$\begin{array}{}\text{(10)}& L={\displaystyle \frac{H^{\mathrm{2}}{g}^{\prime }}{\mathrm{4}{C}_{\mathrm{d}}^{\mathrm{2}}{u}^{\mathrm{2}}}},\end{array}$$

where H=0.05 m is the height of the subglacial environment and $g^{\prime }=\mathrm{0.20}$ m s⁻² is the reduced gravity. The drag coefficient, C_d, and the fluid velocity, u, are both set to the average values within the intrusion. The reduced gravity is referenced to the density difference between the prescribed pure freshwater and pure seawater. Assuming an unobstructed sheet with negligible drag at the salt wedge interface gives a predicted intrusion distance on the order of 10⁰ (i.e., comparable to the blue transects in Fig. 2a). The theory thus gets the magnitude of intrusion distance correct in the simulations of this study. In some cases with free-slip kinematic boundary conditions, seawater fills the entire subglacial portion of the model domain. Extending the model domain to capture the full seawater intrusion distance is computationally prohibitive, hence why we have mainly focused on no-slip cases in this study. Indeed, Eq. (10) predicts unbounded intrusions for C_d=0. Since both the drag coefficient and fluid velocity exhibit great variability along the intrusion regime, and the theory is sensitive to changes in either, it could make choosing which values to use difficult. Another explanation for the disagreement is the assumption of negligible interfacial drag. Robel et al. (2022) find that including interfacial drag of the order of 10⁻⁴ reduces predicted intrusion distance by about a factor of 2. A factor to consider in the future development of more realistic theories of seawater intrusion is the potential role of meltwater feedbacks. In addition, drag from the bottom of the channel will have different mechanical properties than the ice base and will require independent consideration.

Figure 6 Comparison between Wilson et al. (2020) experimental data (gray markers) and intrusion characteristics found in this study (red and blue markers). The red markers represent simulations with melting, and the blue markers represent simulations without melting. The black dashed line is the numerical solution to Robel et al. (2022) with γ=2.

There are two limitations to consider within ANSYS Fluent when melting is enabled: (1) the geometry of the subglacial environment does not change, providing no geometric feedback like that considered in Bradley and Hewitt (2024), and (2) changing the ice wall to a velocity inlet means there is no boundary drag acting on the fluid. Both of these limitations will influence overall intrusion distance. However, we can conclude that the addition of buoyant melt feedback does not greatly affect intrusion distance or set vertical mixing dynamics, which is primarily controlled by the upstream freshwater fluxes, similar to estuaries.

Figure 5 demonstrates the disagreements between the melt calculated directly in our high-resolution simulations and the predicted basal melt rates from parameterizations. The S99 parameterization assumes transport of heat, and salt is dictated by transfer velocities. The thermal turbulent transfer velocity calculated directly from our simulations is on the lower end of those previously documented due to the low velocities tested here and strong stratification within the intrusion regime. Turbulent transfer velocities are derived from theory and experimental results of pipe flow that assume the fluid is steady and parallel to the wall. In seawater intrusions flow is spatially varying, with opposing directions of flow parallel to the wall and some flow perpendicular to the wall. The combined effect of low velocities and narrow gaps under grounded ice (we test subglacial domains with 5 and 7.5 cm gap height in this study) means that the viscous sublayer and buffer layer together occupy the entire domain considered, and drag from either boundary affects our entire subglacial domain. The turbulent transfer velocities in the S99 parameterization are intended to be used with ocean properties at the edge or beyond the outer boundary layer – if the domain of interest does not develop the full boundary layer, then this theory of heat and mass transfer is not appropriate. Since the entire subglacial environment is within the viscous sublayer (where viscous effects dominate) or the buffer layer (where viscous effects and turbulence both occur), the S99 parameterization is not appropriate for calculating melt from seawater intrusions. This is expected since S99 was not formulated considering domains of seawater intrusion and strong stratification.

If we anticipated viscous effects to dominate in seawater intrusions under grounded ice, then using the thermal and haline molecular diffusivities as so-called “transport velocities” would be appropriate. However, this would require knowing the width of the diffusive sublayer, which, given computational constraints for coupled ice–ocean models, cannot be resolved. In polar settings where seawater is warmer than subglacial discharge, double-diffusive convection could be the fitting heat transfer mechanism, acting against the stratification that persists from seawater intrusion. However, this requirement limits this mechanism's relative importance in Antarctica where seawater can be colder than the pressure-dependent freezing point for subglacial discharge (Davis et al., 2023). Double-diffusive convection is dominant where shear is weak, traditionally in low-flow, highly stratified environments, where the geometry does not favor strong natural convection (Rosevear et al., 2021). Observations beneath Ross Ice Shelf near the grounding zone hint at the potential for double-diffusive convection to contribute to vertical mixing when the water column is structured by an upper layer of cold, fresh water adjacent to the ice and a mixed homogenized layer at depth (Begeman et al., 2018). Conceptually, we may expect double-diffusive convection to be a good description of the source of mixing and basal melt in these simulations (as hypothesized by Robel et al., 2022), given the right oceanographic settings. Due to model limitations, the results presented here do not include melt from salt diffusion, hindering a direct comparison to a double-diffusive convection parameterization (e.g., Eq. B1 in Rosevear et al., 2022). Future work should prioritize parameterizing this form of ice loss in regions of seawater intrusion using high-resolution models that can resolve the boundary layer structure.

Currently, only Bradley and Hewitt (2024) have considered the problem of intrusion-induced melting, which allows the subglacial gap to grow in height due to melting in the region of seawater intrusion. Since a taller subglacial gap leads to further seawater intrusion (as predicted in Wilson et al., 2020; Robel et al., 2022), this geometric effect constitutes a positive feedback on seawater intrusion. Our study includes feedbacks from basal melting on flow within the subglacial environment but not this geometric feedback. We find that intrusion-induced melt does not significantly change the structure of the intrusion, but it may act to reduce the regime experiencing intrusion-induced melt. Further, we find that the region of ice experiencing intrusion-induced melt is approximately half of the total intrusion length. Here, we represented the ocean thermal forcing to the ice by the ice-adjacent thermal gradient, as opposed to an average of the water column's temperature as done in Bradley and Hewitt (2024). In this study, there is non-negligible melt for ∼ 50 % of the intrusion distance, whereas averaging the water column temperature results in melt across the entire intrusion. This latter approach potentially enhances the strength of the geometric feedback over a wider area of the ice. The widening of the subglacial space will compete with the excess buoyant forcing to influence seawater intrusion development. Furthermore, tidal pumping of intrusions will add another scale of complexity in resolving these feedbacks, since they provide a temporally asymmetric melt rate as demonstrated in Gadi et al. (2023). Further work is needed to consider both feedbacks within a single modeling framework to determine which ultimately dominates and under what conditions, since they operate on different temporal and spatial scales.

Due to computational limitations, our simulations cannot span the full range of discharge velocities and geometric sizes of potential seawater intrusion domains. However, the cases considered here are comparable to several cases where subglacial properties are constrained from observations (Carter et al., 2017; Davis et al., 2023; Washam et al., 2020), indicating the need for a parameterization that operates as a function of boundary layer processes, near-ice (order of centimeters) seawater properties, and buoyant stability of the water column. Scaling the simulations presented here to ice-shelf cavity size domains will limit our ability to resolve boundary layer processes. This scaling issue is not independent of the constraints within large-scale coupled models, where the coarse resolution will not allow for seawater intrusions to form beyond grounding lines. For such coupled models, further experimentation is needed to identify a subgrid basal melt parameterization (i.e., similar to the purely numerical schemes identified in Seroussi et al., 2019) that can be applied to grounded ice as a function of local bed characteristics, freshwater flux, and thermal forcing. Complex-sloping ice wall boundaries likely couple the seawater intrusion dynamics to the proglacial plume, providing an avenue to combine these processes in a single framework for coupled models with resolution 𝒪 (km).

The challenges associated with complex fluid dynamics and computational limitations are not limited to domains of seawater intrusion. These challenges also apply to the hundreds of meters downstream of the grounding line where a small vertical gap exists between the ice and seafloor. The transition from a two-wall bounded flow to a single-wall bounded flow influenced by buoyant, shear, and diffusive instabilities is currently not represented in coupled models. It is likely this transition couples plume dynamics to intrusion dynamics, which may provide a baseline framework to unify the transition. This regime transition, which may be akin to those described in Rosevear et al. (2022), is critical to understanding what sets grounding zone melt rates and therefore what drives retreat.

A sophisticated implementation of seawater intrusion in coarse-grid coupled ice–ocean models is possible. It would require simulating subglacial hydrology, estimating discharge fluxes, and deducing the required scale at which the ocean grid should extend under “grounded” glacier ice. From there, a modified basal melt parameterization could be added to calculate an intrusion melt rate. If a closed-form parameterization based on first-principles theory is not possible in the complex case of seawater intrusion, a modified parameterization based on a comprehensive suite of high-resolution CFD simulations, similar to this study, could be used to constrain an effective parameterization. Alternatively, intrusion-induced melt could be included in a more ad hoc manner by extending the ocean grid to some prescribed distance inland set by a modified form of the theory proposed by Robel et al. (2022). A prior ice sheet model intercomparison (Seroussi et al., 2019) study found that when subgrid melt schemes are applied to partially floating grid points (cells that cross the grounding line), projected sea level contributions from the Antarctic ice sheet can be up to twice as high. Such subgrid melt schemes have a similar effect to intrusion-induced basal melt.

5 Conclusions

Seawater intrusion beneath grounded portions of ice sheets leads to complex interactions between the subglacial hydrologic system and the ocean. The ability for seawater to intrude under grounded ice is dependent on the subglacial discharge velocity and geometry of the subglacial environment. Both of these controls are consistent with studies of estuaries (Krvavica et al., 2016), bolstering the idea that grounding zones are subglacial estuaries (Horgan et al., 2013). Turbulent mixing has a negligible effect on intrusion distance, since turbulent viscosity is minimal for the subcritical freshwater velocities required for an intrusion to develop. Here, we found that intrusion-induced basal melting does not drastically alter the intrusion mixing dynamics, which is set by the freshwater velocity. Excess buoyancy forcing from basal melting may constitute a negative feedback that will compete with the geometric positive feedback identified in another study (Bradley and Hewitt, 2024). Our simulations show that seawater intrusion can produce enhanced melting at, and just upstream of, the grounding line where warm seawater occupies most of the subglacial environment, with melt rates on the order of tens of meters per year for the configurations considered in this study. Current basal melt parameterizations that assume a fully developed boundary layer underpredict the peak in basal melt rates due to seawater intrusion and fail to capture the drop-off in melting due to stratification at the nose of the intrusion. Diffusive melting frameworks may do a better job of estimating intrusion-induced basal melt near the grounding line where melt rates are high. However, for this to occur in reality, the seawater has to be warmer than the subglacial discharge, which we expect to limit its applicability in Antarctica. Theories of heat and salt transfer that consider the complex overlap of transport phenomena occurring in the buffer layer should be further considered to better incorporate basal melt from seawater intrusion in coarse models.

Coupled ice–ocean models do a poor job of reproducing observed patterns of enhanced basal melt near grounding lines due to resolution and the assumptions inherent in melt parameterizations (Adusumilli, 2021; Ciracì et al., 2023). Furthermore, the assumption of a hydraulic barrier between subglacial hydrology and the ocean leads to an incomplete consideration of the fluid dynamics at grounding zones. Improving projections of the ice sheet contributions to sea level rise requires a more complete representation of grounding zones as dynamic estuaries where subglacial hydrology, the ocean, and glacier ice all interact.

Appendix A: Variable descriptions

Table A1 Variable descriptions and values for ANSYS Fluent simulations.

Table A2 Variable descriptions and values for processing simulation output.

Appendix B: Model settings and equations

B1 Reynolds-averaged Navier–Stokes equations

ANSYS Fluent solves the Reynolds-averaged Navier–Stokes (RANS) equations. This formulation of the Navier–Stokes equations decomposes the flow into a mean state and fluctuation about the mean:

$$\begin{array}{}\text{(B1)}& u_{\mathrm{i}}=\stackrel{\mathrm{‾}}{u_{\mathrm{i}}}+u_{\mathrm{i}}^{\prime },\end{array}$$

where $\stackrel{\mathrm{‾}}{u_{\mathrm{i}}}$ is the time-averaged velocity, and $u_{\mathrm{i}}^{\prime }$ is perturbations about the mean velocity. Scalar quantities are also decomposed into their mean and fluctuations. Substituting these decomposed values into the momentum equations yields

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+{\displaystyle \frac{\partial }{\partial x_i}}\left(\mathit{\rho }{u}_{\mathrm{i}}\right)=\mathrm{0}\text{(B3)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }{u}_{\mathrm{i}}\right)& +{\displaystyle \frac{\partial }{\partial x_i}}\left(\mathit{\rho }{u}_{\mathrm{i}}{u}_j\right)=-{\displaystyle \frac{\partial p}{\partial x_i}}\\ & +{\displaystyle \frac{\partial }{\partial x_j}}\left[\mathit{\mu }\left({\displaystyle \frac{\partial u_{\mathrm{i}}}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\delta }}_{ij}{\displaystyle \frac{u_l}{x_l}}\right)\right]\\ & +{\displaystyle \frac{\partial }{\partial x}}(-\mathit{\rho }\stackrel{\mathrm{‾}}{u_{\mathrm{i}}^{\prime }{u}_j^{\prime }})+\mathit{\rho }{g}_{\mathrm{i}}.\end{array}\end{array}$$

Equations (B2) and (B3) are the RANS equations and have the same general form as the Navier–Stokes equations; however, the velocities and solution variables now represent the time-averaged values. In addition, to represent the effect of turbulence, new terms are added. These terms are the Reynolds stresses, $\stackrel{\mathrm{‾}}{u_{\mathrm{i}}^{\prime }{u}_j^{\prime }}$. To close this system of equations, the Reynolds stresses must be solved via a turbulence closure. More details on turbulence modeling and our closure choice are given in Sect. B2.

B2 Turbulence modeling

The turbulent–viscosity hypothesis (also known as the Boussinesq hypothesis) is used in turbulence modeling to solve for the Reynolds stresses. This hypothesis states the deviatoric Reynolds stresses (those deviating from the mean) are proportional to the mean strain rate tensor by a positive scalar. This scalar represents the eddy viscosity (also referred to as turbulent viscosity). This relationship is

$$\begin{array}{}\text{(B4)}& -\mathit{\rho }\stackrel{\mathrm{‾}}{u_{\mathrm{i}}^{\prime }{u}_j^{\prime }}+{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\mathit{\rho }k{\mathit{\delta }}_{ij}={\mathit{\mu }}_{\mathrm{T}}\left({\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i}}{\partial x_j}}+{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_j}}{\partial x_i}}\right).\end{array}$$

The only unknown in the system of equations is the eddy viscosity (μ_T), which can be solved by a variety of different turbulence closure schemes. Here, we employ the two-equation κ–ϵ closure scheme which solves for eddy viscosity by relating it to the square of turbulent kinetic energy and inverse of turbulent dissipation by a positive scalar C_μ and fluid density (Eq. 3). This closure scheme requires two additional equations to solve for turbulent kinetic energy and turbulent dissipation. These equations are

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }k\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(\mathit{\rho }k{u}_{\mathrm{i}}\right)& ={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mathit{\mu }+{\displaystyle \frac{{\mathit{\mu }}_t}{{\mathit{\sigma }}_k}}\right){\displaystyle \frac{\partial k}{\partial x_j}}\right]\\ & +G_k+G_{\mathrm{b}}-\mathit{\rho }\mathit{ϵ}-Y_M+S_k\end{array}\end{array}$$

for turbulent kinetic energy and

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }\mathit{ϵ}\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(\mathit{\rho }\mathit{ϵ}{u}_{\mathrm{i}}\right)& ={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mathit{\mu }+{\displaystyle \frac{{\mathit{\mu }}_t}{{\mathit{\sigma }}_k}}\right){\displaystyle \frac{\partial \mathit{ϵ}}{\partial x_j}}\right]\\ & +C_{\mathrm{1}\mathit{ϵ}}{\displaystyle \frac{\mathit{ϵ}}{k}}(G_k+C_{\mathrm{3}\mathit{ϵ}}{G}_{\mathrm{b}})\\ & -C_{\mathrm{2}\mathit{ϵ}}\mathit{\rho }{\displaystyle \frac{{\mathit{ϵ}}^{\mathrm{2}}}{k}}+S_{\mathit{ϵ}}\end{array}\end{array}$$

for turbulent dissipation. Buoyancy effects are included in the evolution equation for turbulent kinetic energy by a source term, G_b, which can be solved for by

$$\begin{array}{}\text{(B7)}& G_{\mathrm{b}}=g_{\mathrm{i}}{\displaystyle \frac{{\mathit{\mu }}_{\mathrm{T}}}{\mathit{\rho }{\mathit{Pr}}_{\mathrm{T}}}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial x_i}},\end{array}$$

where g_i is the component of gravity in the ith direction, μ_T is the turbulent viscosity, and Pr_T is the turbulent Prandtl number (0.85 default value). The density gradient, $\partial \mathit{\rho }/\partial x_i$, is taken over the ith direction. Buoyancy effects are not included in the dissipation equation due to a higher degree of uncertainty (ANSYS, 2022). Note that, in the equation for turbulent dissipation, turbulent kinetic energy is in the denominator, resulting in issues when κ approaches zero near wall boundaries. To resolve boundary layer dynamics with the κ–ϵ closure, we use the low-Reynolds formulation, which employs damping functions and fixes the singularity that arises with low values of κ. The low-Reynolds κ–ϵ closure we employ here is the Yang–Shih version (Yang and Shih, 1993). In this formulation, the near-wall turbulence timescale is set to the Kolmogorov timescale (T_k), which is proportional to $(\mathit{\nu }/\mathit{ϵ}{)}^{\mathrm{1}/\mathrm{2}}$. In doing so, the equation for eddy viscosity near the wall becomes

$$\begin{array}{}\text{(B8)}& {\mathit{\mu }}_{\mathrm{T}}=C_{\mathit{\mu }}{f}_{\mathit{\mu }}k\left(T_k+{\displaystyle \frac{k}{\mathit{ϵ}}}\right),\end{array}$$

where f_μ is the “damping function” and equal to

$$\begin{array}{}\text{(B9)}& f_{\mathit{\mu }}=[\mathrm{1}-\exp (-a_{\mathrm{1}}{R}_y-a_{\mathrm{3}}{R}_y^{\mathrm{3}}-a_{\mathrm{5}}{R}_y^{\mathrm{5}}){]}^{(\mathrm{1}/\mathrm{2})}\end{array}$$

and $R_y=k^{(\mathrm{1}/\mathrm{2})}y/\mathit{\nu }$. The constants a₁, a₃, and a₅ are constrained from direct numerical simulation experiments for turbulent channel flow.

The final adjustment to the standard κ–ϵ formulation for near-wall flows is to add an additional source of dissipation, which results from inhomogeneity in the mean flow field. This takes the following form:

$$\begin{array}{}\text{(B10)}& E=\mathit{\nu }{\mathit{\nu }}_{\mathrm{T}}{\left({\displaystyle \frac{\partial u_{\mathrm{i}}}{\partial u_j\partial u_z}}\right)}^{\mathrm{2}}.\end{array}$$

This formulation of the low-Reynolds κ–ϵ turbulence closure allows for solving the freestream portion of the flow regime as well as the near-wall region where viscous effects dominate since the added terms tend to zero when turbulence is high.

We employed the low-Reynolds κ–ϵ turbulence closure model to allow for increased resolution near the ice boundaries. We tested the Yang–Shih formulation used here against the Launder–Sharma low-Reynolds formulation and saw little disagreement between the intrusion distances simulated (Table B2). However, we note that intrusion distance can change based on the choice of turbulence closure formulation. This is distinctly different than the amount of turbulence negligibly changing the intrusion distance as we discuss in the results section of the main text. Different closure schemes will resolve the flow structure differently, leading to a nominal change in intrusion distance, but changing the amount of turbulence in a given turbulence closure scheme does not largely affect intrusion distance.

The κ–ϵ turbulence model provides a simple pathway to modulate the amount of turbulent mixing, by changing the C_μ value. The C_μ value controls the amount of relative turbulent kinetic energy to dissipation represented in the eddy viscosity value and thus directly affects the Reynolds stresses via the turbulent–viscosity hypothesis. The typical value assigned to C_μ is 0.09, which has been empirically found for simple wall-bounded flows (Pope, 2000). However, for complex stratified flows or highly energetic jets, a standard constant value for the whole domain may not be appropriate (e.g., Lai and Socolofsky, 2019). Based on Eq. (3), it is clear the value of C_μ is dependent on turbulence dynamics. However, in a modeling framework, we have to prescribe C_μ. Since a “true” value for C_μ has not been found for the flow regime considered here (stratified and dynamically variable), we deemed it an appropriate approach to modulate this value to induce more or less turbulent mixing. In setting C_μ, we found numerical instability when C_μ→1, therefore limiting our upper-end choice to a factor of 2 of the standard value used. For the lower-end case, we wanted to avoid choosing too low of a value, since the solution relaxes to laminar flow when C_μ→0. This constrained our lower-end choice to a factor of 0.5 of the standard value.

B3 Species transport and energy modeling

To simulate the effects salinity and temperature have on fluid density, we employed the species transport model with a user-defined equation based on a linear equation of state for seawater (Eq. 2).

Fluent automatically turns on energy conservation when species transport is enabled and solves for water temperature via

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}\left({\mathit{\rho }}_{\mathrm{w}}E\right)& +\mathrm{\nabla }\left(\mathit{u}\right({\mathit{\rho }}_{\mathrm{w}}E+p\left)\right)\\ & =\mathrm{\nabla }(k_{\mathrm{eff}}\mathrm{\nabla }T-\sum _{j=\mathrm{1}}{h}_j{J}_j)+S_{\mathrm{h}}\end{array}.\end{array}$$

The first two terms on the right-hand side represent energy transfer due to the conduction of heat (∇(k_eff∇T)) and species diffusion ($\mathrm{\nabla }\left({\sum }_{j=\mathrm{1}}{h}_j{J}_j\right)$). Conduction represents heat transfer due to thermal gradients. As salt diffuses in the medium, it also transfers heat due to its unique thermal properties and therefore must also be included. The contribution to heat transport from salt has negligible impacts since it is not a leading order term in Eq. (B11), and therefore tracking species contribution to energy transport does not significantly change the results. S_h represents chemical sources of heat; here, there are none, so this term is zero (ANSYS, 2009). The left-hand side represents the time evolution of energy and energy advection. This energy transport equation also represents the transport of temperature throughout the domain, with the relation T_wc_wρ_wV=E, where T is temperature, c_w is the heat capacity, ρ_w is the seawater density, and V is the volume.

Salt is transported as an active tracer within the fluid, employing an advection–diffusion–reaction equation:

$$\begin{array}{}\text{(B12)}& {\displaystyle \frac{\partial }{\partial t}}\left(\mathit{\rho }S\right)+\mathrm{\nabla }\left({\mathit{\rho }}_{\mathrm{w}}\mathit{u}S\right)=-{\mathit{\rho }}_{\mathrm{w}}{D}_{\mathrm{s}}\mathrm{\nabla }S+R,\end{array}$$

where u is the velocity vector, D_s is the diffusion coefficient for salt (set to $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{9}}$ m² s⁻¹), and R is the production rate from reactions. The production rate for salt in these simulations is zero. The boundary salinity, S_b, is set to the near-ice grid cell (scale of nanometers), which follows the evolution equation given by Eq. (B12). This evolution equation accounts for advective and diffusive transport mechanisms. The boundary salinity is used in the liquidus condition (Eq. 1) to set the thermal boundary condition for the ice interfaces. In the reference simulations, S_b only evolves due to advection and diffusion from the salt wedge entering the subglacial domain and the freshwater layer exiting the subglacial space as a plume. For the simulations with melting, S_b is diluted as fresh, cold water enters the domain normal to the ice boundary. This dilution occurs because Fluent is a finite volume solver, so it balances inflows and outflows across each control volume in the grid while conserving overall mass balance in the domain via a pressure outlet (red arrow Fig. 1). The thermal boundary condition (given by Eq. 1) adjusts based on this dilution. Similarly, the chilling of near-ice waters occurs due to local injection of meltwater at the freezing point.

B4 Domain

The 2-dimensional domain is defined as a long and thin subglacial environment attached to a tall rectangular ocean basin. The length of the subglacial environment is 20 m for most cases. The slowest freshwater velocity simulations and the taller subglacial environment simulations required a longer domain to reach a quasi-steady state, and therefore their subglacial environment length is 40 m. The height of the subglacial domain is 5 cm for all cases except the taller scenario where the height is 7.5 cm. The ocean tank's height is 5 m, and its width is 2 m for all cases. For analysis, we designate the “grounding line” or the point where the subglacial environment meets the ocean basin to be at x=0 m. However, in the model domain, this point exists at x=2 m.

A freshwater inlet is defined at the rightmost boundary of the subglacial environment. A seawater inlet is defined at the leftmost boundary. A pressure outlet is set at the top of the ocean domain, enforcing a no-gradient flux across the boundary. The ice boundary walls are defined at the top of the subglacial environment and the right side of the ocean domain. For both melting and non-melting simulations, these ice walls have a thermal boundary condition dependent on near-wall salinity and pressure (Eq. 1). Neither ice wall boundary allows for salinity diffusion, which would be another mechanism of melt to account for. The non-melting simulations employ a no-slip kinematic condition, which forces the fluid velocity to be zero at the wall. For cases with melting, the top boundary of the subglacial space is turned into a velocity inlet to simulate melt. In designating this boundary as a velocity inlet, a free-slip kinematic condition is required. The downward vertical velocity set at the melting ice boundary follows Eq. (4). We turn off melting for the vertical ice boundary to isolate the intrusion-induced melt from the vertical plume dynamics that would arise from the vertical ice boundary.

B5 Meshing

For the given turbulence closure scheme employed here (Yang–Shih low-Reynolds k−e turbulence model), a nondimensionalized distance, y⁺, must be ≤5. This follows the law-of-the-wall principle, where mesh thickness near the wall must be coordinated with the closure scheme for near-wall viscosity effects to be rendered correctly. This nondimensionalized distance y⁺ is a function of dynamic viscosity, density, distance normal to the wall, and shear velocity.

$$\begin{array}{}\text{(B13)}& y^{+}={\displaystyle \frac{u_{*}y}{\mathit{\nu }}}\end{array}$$

Given a value for y⁺ and known kinematic viscosity, ν, all that is needed is the shear velocity, u_*, in order to find the necessary dimensional y distance of the first grid cell. To find u_* we need the skin friction C_f. For laminar flows, the skin friction is

$$\begin{array}{}\text{(B14)}& C_{\mathrm{f}}={\displaystyle \frac{\mathrm{16}}{\mathit{Re}}},\end{array}$$

where Re is the Reynolds number with the characteristic length being the channel length:

$$\begin{array}{}\text{(B15)}& {\mathit{Re}}_L={\displaystyle \frac{uL}{\mathit{\nu }}}.\end{array}$$

To find u_*,

$$\begin{array}{}\text{(B16)}& u_{*}=\sqrt{{\displaystyle \frac{\mathit{\tau }}{{\mathit{\rho }}_{\mathrm{w}}}}},\end{array}$$

where $\mathit{\tau }=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{C}_{\mathrm{f}}{u}^{\mathrm{2}}$. Rearranging we then get

$$\begin{array}{}\text{(B17)}& u_{*}=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{f}}{u}^{\mathrm{2}}}.\end{array}$$

We can then plug this value for u_* back into Eq. (B13) to calculate the near-wall mesh thickness. We generate a mesh for each freshwater velocity and geometric combination (five total) following the steps above and use inflation layers to increase resolution near the wall boundary but reduce it near the center of the domain. Generating the mesh as such greatly reduces total simulation run time and data storage requirements.

Figure B1 Initial condition for all model simulations.

Table B1 Simulations with varied ocean velocity.

Table B2 Simulations with varied modeling choices for low-Reynolds turbulence closure and equation of state (E.O.S.).

B6 Run time

To keep the simulation time consistent, each simulation runs for 43 200 s (12 h). After this period of time, we evaluate whether the solution has reached a quasi-steady state by evaluating the change in density of the subglacial environment over 100 time steps. If this change is below 10⁻⁴ kg m⁻³, then a quasi-steady state has been achieved. A few runs required longer simulation times to reach such a steady state and were consequently run for another 12 h.

To set the time step we used the eddy turnover time. This is the largest time step we can take without causing non-physical effects and is represented by

$$\begin{array}{}\text{(B18)}& t={\displaystyle \frac{h}{u_{*}}},\end{array}$$

where h is the height of the channel. For each simulation, a time step of 5 s was chosen in accordance with the fastest freshwater flow (in turn it has the largest shear velocity and therefore smallest eddy turnover time). Putting this together, for all runs, a total of 8640 time steps were taken at a time step of 5 s. For every 20th time step (100 s), data were exported along a middle transect, the ice base, and the domain's entire fluid surface.

Each simulation is initialized with a fresh and cold subglacial environment (S=0 ppt and $T_{\mathrm{w}}=-\mathrm{0.75}$ °C) and a warm and salty ocean basin (S=30 ppt and T_w=0.5 °C) as shown in Fig. B1. The simulation then runs for 12 h of simulation time. Each simulation experiences a period of “spinup” time where the intrusion develops, and almost all simulations reach a quasi-steady state. We define a quasi-steady state to occur once the time change in subglacial environment's average density is less than 0.0001. The results presented in the main body of this work are the time-averaged results once a quasi-steady state is reached and for the following 5000 s of simulation time.

B7 Variable ocean forcing

Previous theoretical work on seawater intrusion neglected a background ocean velocity in order to get a simple relationship between freshwater flux and intrusion distance. To test the appropriateness of this simplification, we varied the seawater inlet velocity across three magnitudes (u_o=5, 0.5, and 0.05 cm s⁻¹). In addition, we also tested the smallest ocean velocity the model would converge at: $u_{\mathrm{o}}={\mathrm{10}}^{-\mathrm{6}}$ cm s⁻¹. Between the ocean velocities tested, there is essentially no difference in intrusion distance, mixing characteristics, or vertical structure. This supports the idea that the exchange flow that develops within the intrusion is set by the flux of buoyant freshwater and is independent of ocean currents, similar to estuaries. A summary of these results is given in Tables B1 and 1.

B8 Equation of state

For all the simulations reported in the main body of this work, a linear equation of state from Roquet et al. (2015) (Eq. 2) is used to solve for fluid density. We also tested a higher-order equation of state from Roquet et al. (2015):

$$\begin{array}{}\text{(B19)}& {\mathit{\rho }}_{\mathrm{w}}=\mathrm{1000}-{\displaystyle \frac{C_{\mathrm{b}}}{\mathrm{2}}}(T_{\mathrm{w}}-{\mathrm{\Theta }}_{\mathrm{0}}{)}^{\mathrm{2}}-T_{\mathrm{h}}Z{T}_{\mathrm{w}}+b_{\mathrm{0}}{S}_A.\end{array}$$

Here, C_b=0.011 kg m⁻³ K⁻², $T_{\mathrm{h}}=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{5}}$ kg m⁻⁴ K⁻¹, b₀=0.77 kg m⁻³ (kg g⁻¹), and ${\mathrm{\Theta }}_{\mathrm{0}}=-\mathrm{4.5}$ °C. With the higher-order equation of state, the baseline densities of the freshwater and the seawater were lower than the linear equation of state. However, the distribution of salt and heat remained the same because they are transported by their respective equations (Eqs. B12, B11). Because of this, the intrusion distance did not vary significantly between the two equation-of-state formulations.

Appendix C: Calculating the drag coefficient

The drag coefficient cannot be prescribed for the simulations but is rather diagnosed from the model output. To do so, we use the relationship between the drag coefficient, C_d, and the wall shear stress, τ_w (Pope, 2000):

$$\begin{array}{}\text{(C1)}& C_{\mathrm{d}}={\displaystyle \frac{\mathrm{2}{\mathit{\tau }}_{\mathrm{w}}}{\mathit{\rho }{u}^{\mathrm{2}}}}.\end{array}$$

The wall shear stress is equal to

$$\begin{array}{}\text{(C2)}& {\mathit{\tau }}_{\mathrm{w}}=\mathit{\mu }\left({\displaystyle \frac{\partial u}{\partial y}}\right).\end{array}$$

We can rearrange and solve for C_d as a function of the near-wall velocity gradient, which produces Eq. (6).

Appendix D: Turbulent viscosity results

Turbulent viscosity is helpful to understand the strength of turbulence-induced mixing within a fluid. In the fluid regimes tested here, with low Reynolds numbers, turbulent viscosity is small and comparable to seawater's kinematic viscosity. It is for this reason that modulating the degree of turbulent mixing does not impact the intrusion distance. Turbulent viscosity is greatly reduced over the length of the intrusion, likely due to enhanced stratification suppressing mixing. For the middle and slowest freshwater velocities, the intrusion distance is near the transition point defined by the critical gradient Richardson number, which relates the relative importance of shear to buoyancy.

Figure D1 Turbulent viscosity for each freshwater velocity and an each turbulent mixing scenario. Solid lines are for medium turbulent mixing (C_μ=0.09), dashed lines are for high turbulent mixing (C_μ=0.18), and dash–dot lines are for low turbulent mixing (C_μ=0.045). The red line represents seawater's kinematic viscosity. The dots represent the intrusion distance for the medium turbulent mixing case for the given freshwater velocity. The cross markers represent the point at which the critical gradient Richardson number is crossed.

Appendix E: Vertical profiles for retrograde and taller subglacial channel geometries

Changing the subglacial environment's geometry did not greatly impact the overall structure of the intrusion even if it did change the intrusion distance. The vertical distribution of heat and salt for both the retrograde and taller geometries are qualitatively similar to their comparable base case (Figs. 3e, f and E1a, b, e, f). Consequently, the degree of stratification, as measured by the buoyancy frequency, is similar (Figs. 3h and E1d, h). The structure of the exchange flow is quite similar as well, albeit with a slightly stronger exchange flow for the taller geometry due to enhanced flux of freshwater at the subglacial inlet.

Figure E1 Time-averaged vertical profiles of thermal forcing (a, e), salinity (b, f), x component of velocity (c, g), and buoyancy frequency (d, h) along the seawater intrusion for the retrograde (a–d) and taller (e–h) subglacial environment geometries with u_f=0.5 cm s⁻¹ and medium turbulent mixing.