Ice-shelf ocean boundary layer dynamics from large-eddy simulations

Small-scale turbulent flow below ice shelves is regionally isolated and difficult to measure and simulate. Yet these small-scale processes, which regulate heat and salt transfer between the ocean and ice shelves, can affect sea-level rise by altering the ability of Antarctic ice shelves to “buttress” ice flux to the ocean. In this study, we improve our understanding of turbulence below ice shelves by means of large-eddy simulations at sub-meter resolution, capturing boundary layer mixing at scales intermediate between laboratory experiments or direct numerical simulations and regional or global ocean circulation models. Our simulations feature the development of an ice-shelf ocean boundary layer through dynamic ice melting in a regime with low thermal driving, low ice-shelf basal slope, and strong shear driven by the geostrophic flow. We present a preliminary assessment of existing ice-shelf basal melt parameterizations adopted in single component or coupled ice-sheet and ocean models on the basis of a small parameter study. While the parameterized linear relationship between ice-shelf melt rate and far-field ocean temperature appears to be robust, we point out a little-considered relationship between ice-shelf basal slope and melting worthy of further study.

1 Introduction

The largest source of uncertainty in future sea-level rise is the potential loss of ice from the Antarctic Ice Sheet (IPCC, 2014). The rate of grounded ice loss is highly sensitive to the melting of ice shelves, which drain over 80 % of Antarctica’s grounded ice (Reese et al., 2018; Rignot et al., 2013). The ice-shelf ocean boundary layer (IOBL) controls ice-shelf melting by regulating oceanic heat and salt fluxes to the ice-shelf base. Thus, accurate predictions of ice-shelf melting depend on representing the turbulent dynamics of the IOBL. This representation is also critical for evaluating the sensitivity of the coupled land-ice–ocean system to changes in ocean conditions. Of particular concern is the sensitivity of ice-shelf melting to increasing seawater temperature at depth, a trend observed along a wide swath of the West Antarctic coastline and a potential trigger for West Antarctic Ice Sheet collapse (Purkey et al., 2018; Ruan et al., 2021; Schmidtko et al., 2014; Wåhlin et al., 2021).

One indication that ocean models do not capture turbulent dynamics in the IOBL is that the simulated thermal driving, the difference between ocean temperature and the local freezing point, and consequently the simulated melt rate differ significantly between ocean models and as resolution is varied within a model (Gwyther et al., 2020). Furthermore, ocean models predict ice-shelf melting using parameterizations that neglect the effects of the lateral buoyancy gradient of the IOBL, and often of the vertical buoyancy gradient, on the efficiency of vertical mixing near the boundary. This model deficiency likely biases turbulent fluxes through the IOBL and ocean cavity circulation, which is driven in large part by the buoyant flow of water freshened by ice-shelf melting, an overturning circulation also known as the “ice pump” (Webber et al., 2018). A new parameterization of ice-shelf melting that accounts for the dynamics of IOBL turbulence is needed to achieve a more physically based, accurate coupling of ice sheets and oceans in climate models (Dinniman et al., 2016; Edwards et al., 2021; Gwyther et al., 2020; Naughten et al., 2018). Supercooling of the IOBL and frazil ice accretion to the ice-shelf base are also regionally important processes for cold cavities but are not the focus of this work (Galton-Fenzi et al., 2012; Jordan et al., 2015).

IOBLs present unique conditions in the global ocean, involving a stabilizing buoyancy flux from melting ice and a boundary layer that is positively buoyant against a sloping boundary. There is a rich literature on stably stratified boundary layers (typically under constant stabilizing flux boundary conditions), but the dependence of heat, salt, and momentum fluxes on stratification remains a difficult problem, especially for strongly stratified regimes (Zonta and Soldati, 2018). IOBL turbulence has been explored through laboratory experiments, direct numerical simulations, and large-eddy simulations (Middleton et al., 2021; Mondal et al., 2019; Vreugdenhil and Taylor, 2019; McConnochie and Kerr, 2018; Rosevear et al., 2021). However, this body of work has not yet matured to setting a new standard for ice-shelf melt parameterization in ocean models. Today, the most commonly used ice-shelf melt parameterization is still derived from sea-ice conditions (i.e., in the absence of a slope), with some parameters tuned for ice-shelf settings (Holland and Jenkins, 1999; Jenkins et al., 2010; McPhee et al., 1987). A knowledge gap exists in bridging IOBL dynamics across scales and characterizing the structure of buoyant plumes. Recently, this has been addressed through non-turbulence-resolving 2-d or 2.5-d models (Cheng et al., 2020; Jenkins, 2016, 2021). Notably, Jenkins (2021) found that vertical mixing in IOBL settings was represented poorly by the commonly used K-profile parameterization (KPP) in comparison to a higher-order turbulence closure scheme, suggesting that the application of KPP in ocean models may be inappropriate in ice-shelf cavities.

In this study, we model turbulent heat, salt, and momentum fluxes through the IOBL using large-eddy simulation (LES). Whereas the ocean models typically used to model sub-ice-shelf circulation lack both the resolution and appropriate parameterizations needed to capture the relevant turbulent scales for boundary layer dynamics, LES captures the dominant energy-containing scales of turbulence and represents smaller, unresolved scales with varying degrees of complexity. An effective parameterization of ice-shelf melting is likely to rest on an understanding of how turbulent mixing in the IOBL depends on stratification and shear forcing. We vary far-field ocean temperature and ice-shelf slope between model runs and characterize turbulent fluxes and ice-shelf melt rates. In this study we focus on the high-shear, low-thermal-driving regime. The target of previous LES studies has been on low-shear settings (Middleton et al., 2021; Vreugdenhil and Taylor, 2019). Section 2 describes the LES model, its turbulence closure scheme, and the setup of our simulations. Section 3 presents the results of our simulations, with a focus on comparisons across the sampled range of thermal driving and slope. The discussion is given in two parts: Sect. 4.1 contextualizes the features of our simulated IOBL turbulence and discusses limitations of this study, and Sect. 4.2 focuses on how this study informs parameterizations of ice-shelf melting and IOBL turbulence. We provide some closing thoughts in Sect. 5.

2 Methods

2.1 Overview of the LES model

The Parallelized Large-Eddy Simulation Model (PALM) was developed at the Institute of Meteorology and Climatology at Leibniz University Hannover, Germany (Raasch and Schröter, 2001), for simulation of atmospheric and ocean boundary layers. For this application to sub-ice-shelf settings, we developed a new version based on PALM v5.0 (Maronga et al., 2015) with added features, including a boundary flux scheme for the ice–ocean interface, rotating the gravity and Coriolis vectors for sloped domains, and a different turbulence closure scheme. Here, we provide a brief overview of PALM with a focus on our additions. For more details we refer readers to Maronga et al. (2015).

The governing equations of PALM are the following:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}\text{Momentum conservation}{\displaystyle \frac{\partial u_i}{\partial t}}& =-{\displaystyle \frac{\partial u_i{u}_j}{\partial x_j}}-{\mathit{\epsilon }}_{ijk}{f}_j{u}_k\\ & +{\mathit{\epsilon }}_{i\mathrm{3}j}{f}_{\mathrm{3}}{u}_{g,j}+g_i{\displaystyle \frac{\mathit{\rho }-{\mathit{\rho }}_{\mathrm{a}}}{{\mathit{\rho }}_{\mathrm{0}}}}\\ & -{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial {\mathit{\pi }}^{*}}{\partial x_i}}-{\displaystyle \frac{\partial }{\partial x_j}}\\ & \times \left(\stackrel{\mathrm{‾}}{u_i^{\prime \prime }{u}_j^{\prime \prime }}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\stackrel{\mathrm{‾}}{u_i^{\prime \prime }{u}_i^{\prime \prime }}{\mathit{\delta }}_{ij}\right)\end{array},\text{(2)}& {\displaystyle }\text{Volume conservation}{\displaystyle \frac{\partial u_j}{\partial x_j}}=\mathrm{0},\text{(3)}& {\displaystyle }\text{Heat conservation}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}=-{\displaystyle \frac{\partial u_j\mathit{\theta }}{\partial x_j}}-{\displaystyle \frac{\partial }{\partial x_j}}\left(\stackrel{\mathrm{‾}}{u_j^{\prime \prime }{\mathit{\theta }}^{\prime \prime }}\right)+{\mathrm{\Phi }}_{\mathit{\theta }},\text{(4)}& {\displaystyle }\text{Salt conservation}{\displaystyle \frac{\partial S}{\partial t}}=-{\displaystyle \frac{\partial u_{j}S}{\partial x_j}}-{\displaystyle \frac{\partial }{\partial x_j}}\left(\stackrel{\mathrm{‾}}{u_j^{\prime \prime }{S}^{\prime \prime }}\right)+{\mathrm{\Phi }}_S.\end{array}$$

The momentum terms on the right hand side of Eq. (1) are, in order, advection, Coriolis forcing, imposed geostrophic flow with a geostrophic velocity u_g, buoyancy forcing, a correction for divergence in the flow using the pressure perturbation π^* (imposing incompressibility), and sub-grid-scale momentum flux. All prognostic variables are considered filtered at the grid scale. This is typically denoted by the overbar, which we omit except for the sub-grid turbulent flux terms to emphasize that we only represent averaged effects with a turbulence closure scheme. Double primes denote sub-grid fluctuations.

We represent a sloping ice base by rotating the gravity (g) and Coriolis (f) vectors while keeping the domain a rectangular prism as in Vreugdenhil and Taylor (2019). Specifically,

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathit{g}& {\displaystyle }=g_z[\sin \mathit{\alpha }\sin \mathit{\beta },\sin \mathit{\alpha }\cos \mathit{\beta },\cos \mathit{\alpha }],\text{(6)}& {\displaystyle }\mathit{f}& {\displaystyle }=\mathrm{2}\mathrm{\Omega }[\sin \mathit{\varphi }\sin \mathit{\alpha }\sin \mathit{\beta },\cos \mathit{\varphi }\sin \mathit{\alpha }\cos \mathit{\beta },\sin \mathit{\varphi }\cos \mathit{\alpha }],\end{array}$$

where g_z is the magnitude of gravitational acceleration in the geopotential z direction, β is the angle the up-slope direction makes with north, and α is the slope angle from horizontal or equivalently the angle between the z axis and g. For the rotated Coriolis vector f, Ω is the rotation rate, and ϕ is the latitude. Unless otherwise specified, quantities are oriented with the simulated grid.

The buoyancy term in Eq. (1) combines the contributions of the along-slope pressure gradient due to the slope of the ice shelf in hydrostatic equilibrium and buoyancy due to changes in density relative to the ambient density of the water column ρ_a. We use the nonlinear equation of state from Jackett et al. (2006) to compute densities, and ρ_a varies along slope and with depth from the ice interface as a function of the hydrostatic pressure:

$$\begin{array}{}\text{(7)}& {\mathit{\rho }}_{\mathrm{a}}(x,y,z)=f_{\mathrm{EOS}}\left({\mathit{\theta }}^i\left(z\right),S^i\left(z\right),p(x,y,z)\right),\end{array}$$

where the superscript i denotes the initial state. The reference density ρ₀ is evaluated in the center of the x–y plane (${\mathit{\rho }}_{\mathrm{0}}\left(z\right)={\mathit{\rho }}_{\mathrm{a}}(x_{\mathrm{mid}},y_{\mathrm{mid}},z)$). We neglect hydrostatic pressure gradients along slope except through ρ_a in the buoyancy term. For the maximum slope simulated in this study, the maximum hydrostatic pressure gradient is on the order of 100 Pa m⁻¹. This simplification has the advantage of avoiding pressure discontinuities across the periodic boundary.

The terms on the right-hand side of the heat and salt conservation equations (Eqs. 3 and 4) are advection by the resolved flow, turbulent transport by the sub-grid-scale fluctuations, and the source and sink terms Φ_θ,S. The source and sink terms are zero in our simulations because we treat heat and salt fluxes due to melting at the boundaries through the sub-grid vertical flux term, as discussed in Sect. 2.2.

We employ PALM's implementation of the Piacsek and Williams (1970) second-order advection scheme for momentum and scalars and a third-order Runge–Kutta time-stepping scheme. We also use the Temperton (1992) fast Fourier transform algorithm for the pressure solver.

The turbulence closure scheme is anisotropic minimum dissipation (AMD; Rozema et al., 2015) employing the extension to scalars introduced by Abkar et al. (2016). We validated our implementation against the stable atmospheric boundary layer test case published in Abkar and Moin (2017) and Stoll and Porté-Agel (2008). We follow the Abkar and Moin (2017) model setup exactly and chose their moderate resolution with 72 × 72 × 72 grid points. Our results are in good agreement with the Abkar and Moin (2017) solution employing another LES model with AMD (boundary layer height differs by <10 %) and other LES models with different turbulence closures (Fig. S1 in our Supplement and their Fig. 1 and references therein). Typical sub-grid-scale momentum and scalar diffusivities in our sub-ice-shelf simulations of stratified turbulence range from 10⁻⁵ to 10⁻⁴ m² s⁻¹ in the upper two-thirds of the domain where damping is not applied (the damping methodology is discussed in Sect. 2.3).

2.2 Boundary conditions for ice-shelf melting

Sources and sinks of momentum, as well as of heat and salt due to interactions between the flow and the ice base, are parameterized as sub-grid fluxes at the top boundary. The resolved vertical fluxes at the top layer of the model go to zero as w goes to zero according to the no-flux boundary condition, and the sub-grid fluxes are determined using the scheme described below in place of the AMD scheme at the top layer. PALM is vertically discretized such that the ice boundary (z=0) is located at the edge of the top cell where the vertical velocity component resides, and scalar and horizontal velocity components are located at mid-depth of the top cell ($z=-\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\Delta }}_z$). We denote the interface location with subscript b and the middle of the first cell below the interface with subscript $\frac{\mathrm{1}}{\mathrm{2}}$.

Sub-grid momentum fluxes are parameterized according to the law of the wall following a linear stability function for stabilizing buoyancy forcing as in Vreugdenhil and Taylor (2019):

$$\begin{array}{}\text{(8)}& {\stackrel{\mathrm{‾}}{w^{\prime \prime }{u}_i^{\prime \prime }}}_b=u_{*}{\displaystyle \frac{\mathit{\kappa }\left[u_i\left(z_{\frac{\mathrm{1}}{\mathrm{2}}}\right)-u_{i,\mathrm{b}}\right]}{\ln \left(z_{\frac{\mathrm{1}}{\mathrm{2}}}/z_{\mathrm{0}}\right)-{\mathrm{\Psi }}_m\left({\mathit{\zeta }}_{\frac{\mathrm{1}}{\mathrm{2}}}\right)}}\end{array}$$

for the horizontal velocity components (i.e., $i=\mathrm{1},\mathrm{2}$). ${\mathit{\zeta }}_k=z_k/L_{\mathrm{O}}$ is the depth from the boundary scaled by the Monin–Obukhov length, u_* is the friction velocity, and z₀ is the roughness length. The horizontal velocity at the boundary u_i,b is always 0. The Monin–Obukhov length L_O, computed following McPhee et al. (1987), is a function of both u_* and the melt rate. The stability function Ψ_m is linear with the scaled depth.

$$\begin{array}{}\text{(9)}& {\mathrm{\Psi }}_m\left(\mathit{\zeta }\right)=\mathrm{1}+{\mathit{\beta }}_m\mathit{\zeta }\end{array}$$

This parameterization assumes that Coriolis forces are not locally dominant in the momentum balance.

Similarly, the sub-grid fluxes of scalars are parameterized by

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{‾}}{w^{\prime \prime }{\mathit{\theta }}^{\prime \prime }}}_b={\mathrm{\Gamma }}_{\mathit{\theta }}{u}_{*}\left({\mathit{\theta }}_{\frac{\mathrm{1}}{\mathrm{2}}}-{\mathit{\theta }}_b\right),{\stackrel{\mathrm{‾}}{w^{\prime \prime }{S}^{\prime \prime }}}_b={\mathrm{\Gamma }}_S{u}_{*}\left(S_{\frac{\mathrm{1}}{\mathrm{2}}}-S_b\right),\end{array}$$

where the exchange coefficients are defined with two terms, one representing turbulent transfer within the turbulent surface layer and one representing molecular transfer within the viscous sublayer following McPhee et al. (1987) (their Eq. 11).

$$\begin{array}{}\text{(11)}& {\mathrm{\Gamma }}_{\mathit{\theta }}=\left({\mathrm{\Gamma }}_{\mathit{\theta },\mathrm{turb}}+{\mathrm{\Gamma }}_{\mathit{\theta },\mathrm{mol}}\right),{\mathrm{\Gamma }}_S=\left({\mathrm{\Gamma }}_{S,\mathrm{turb}}+{\mathrm{\Gamma }}_{S,\mathrm{mol}}\right)\end{array}$$

The turbulent flux component follows a shape function consistent with the parameterization of momentum fluxes in Eq. (8):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}{\mathrm{\Gamma }}_{\mathit{\theta },\mathrm{turb}}& ={\mathit{\kappa }}^{-\mathrm{1}}\left[\ln \left(z_{\frac{\mathrm{1}}{\mathrm{2}}}/z_{\mathrm{0}}\right)-{\mathrm{\Psi }}_{\mathit{\theta }}\left({\mathit{\zeta }}_{\frac{\mathrm{1}}{\mathrm{2}}}\right)\right],\\ & {\mathrm{\Gamma }}_{S,\mathrm{turb}}={\mathit{\kappa }}^{-\mathrm{1}}\left[\ln \left(z_{\frac{\mathrm{1}}{\mathrm{2}}}/z_{\mathrm{0}}\right)-{\mathrm{\Psi }}_S\left({\mathit{\zeta }}_{\frac{\mathrm{1}}{\mathrm{2}}}\right)\right],\end{array}\end{array}$$

with analogous stability functions to momentum, Eq. (9),

$$\begin{array}{}\text{(13)}& {\mathrm{\Psi }}_{\mathit{\theta }}\left(\mathit{\zeta }\right)=\mathrm{1}+{\mathit{\beta }}_{\mathit{\theta }}\mathit{\zeta },{\mathrm{\Psi }}_S\left(\mathit{\zeta }\right)=\mathrm{1}+{\mathit{\beta }}_S\mathit{\zeta }.\end{array}$$

The coefficients of these stability functions are chosen as in Zhou et al. (2017) and are given in Table S1. The molecular flux components are constant in our simulations following McPhee et al. (1987) and are also given in Table S1. If a portion of the ice face experiences freezing, Γ_θ and Γ_S are set to a higher value, indicating destabilizing fluxes (Γ_f, Table S1).

Taking into account the role of meltwater advection, Eq. (10) is replaced by the following “virtual” freshwater flux form in our implementation (Asay-Davis et al., 2016; Jenkins et al., 2001).

$$\begin{array}{}\text{(14)}& {\displaystyle }{\stackrel{\mathrm{‾}}{w^{\prime \prime }{\mathit{\theta }}^{\prime \prime }}}_b=& {\displaystyle }-u_{*}\left[{\mathrm{\Gamma }}_{\mathit{\theta }}-{\mathrm{\Gamma }}_S{\displaystyle \frac{S_b-S_{\frac{\mathrm{1}}{\mathrm{2}}}}{S_b}}\right]\left({\mathit{\theta }}_b-{\mathit{\theta }}_{\frac{\mathrm{1}}{\mathrm{2}}}\right)\text{(15)}& {\displaystyle }{\stackrel{\mathrm{‾}}{w^{\prime \prime }{S}^{\prime \prime }}}_b=& {\displaystyle }-u_{*}\left[{\mathrm{\Gamma }}_S-{\mathrm{\Gamma }}_S{\displaystyle \frac{S_b-S_{\frac{\mathrm{1}}{\mathrm{2}}}}{S_b}}\right]\left(S_b-S_{\frac{\mathrm{1}}{\mathrm{2}}}\right)\end{array}$$

The temperature and salinity at the ice–ocean interface, θ_b and S_b, are unknown. Three equations are used to solve for these quantities and the melt rate m, the so-called three-equation parameterization.

$$\begin{array}{}\text{(16)}& {\displaystyle }\mathit{\rho }{c}_p{\stackrel{\mathrm{‾}}{w^{\prime \prime }{\mathit{\theta }}^{\prime \prime }}}_b& {\displaystyle }=-{\mathit{\rho }}_{w}mL\text{(17)}& {\displaystyle }\mathit{\rho }{\stackrel{\mathrm{‾}}{w^{\prime \prime }{S}^{\prime \prime }}}_b& {\displaystyle }=\mathrm{-}{\mathit{\rho }}_{w}m{S}_b\text{(18)}& {\displaystyle }{\mathit{\theta }}_b& {\displaystyle }={\mathit{\theta }}_{\mathrm{f}}(p,S_b)\end{array}$$

These equations specify that heat and salt are conserved at the ice–ocean interface (Eqs. 16 and 17), and the interface temperature is fixed at the local freezing temperature θ_f (Eq. 18). Equation (16) assumes that the conductive heat flux into ice is negligible. ρ_w denotes the density of freshwater. The freezing point is calculated using the polynomial function from Jackett et al. (2006).

The PALM implementation applies the fluxes ${\stackrel{\mathrm{‾}}{w^{\prime \prime }{u}_i^{\prime \prime }}}_b$, ${\stackrel{\mathrm{‾}}{w^{\prime \prime }{\mathit{\theta }}^{\prime \prime }}}_b$, and ${\stackrel{\mathrm{‾}}{w^{\prime \prime }{S}^{\prime \prime }}}_b$ at the center of the first grid cell from the boundary without interpolation (i.e., ${\stackrel{\mathrm{‾}}{w^{\prime \prime }{X}^{\prime \prime }}}_b={\stackrel{\mathrm{‾}}{w^{\prime \prime }{X}^{\prime \prime }}}_{\frac{\mathrm{1}}{\mathrm{2}}}$). It is noted by the PALM developers that this error was found to be small, but we have not confirmed this for the sub-ice case. This error would be small if the first ∼10 cm were nearly a constant flux layer; McPhee (1981) hypothesized that the sub-ice-surface layer would be nearly but not exactly a constant flux layer.

2.3 Simulation setup

A schematic of the simulation domain is shown in Fig. 1. A list of parameter choices relevant to our simulations can be found in Table S1. The domain is a 64 m³ cube with horizontal resolution of 0.5 m and vertical resolution of 0.25 m. The ice–ocean interface is located on the top boundary of the domain. Large-scale horizontal pressure gradients drive the mean flow, which is geostrophic in the far field where buoyancy does not modify the flow. We choose a pressure of 800 dbar at the top of the domain, an intermediate choice given that the depth of ice-shelf bottoms ranges over an order of magnitude (roughly −2000 to −200 m). The potentially dynamically relevant differences between conducting these simulations at surface pressure and 800 dbar are that the first derivative of the freezing temperature with respect to salinity is 20 % smaller and the first derivative of density with respect to salinity is about 2 times larger.

Figure 1 Schematic of the simulated ocean domain with background pressure gradient $\mathrm{d}p/\mathrm{d}y$. The purple arrow is oriented north, and the green arrow is aligned with gravitational acceleration. The ^* signifies that the bottom boundary condition is Dirichlet, but there is also no flux as a result of damping.

We set von Neumann boundary conditions at the top boundary corresponding to the dynamic sub-grid momentum and scalar fluxes (Eqs. 8, 14, and 15) as resolved fluxes go to zero at a no-penetration boundary. The roughness length (z₀ in Eqs. 8 and 12) is chosen such that the equivalent drag coefficient assuming quadratic drag is 0.003, an intermediate value for sea ice and ice-shelf bottoms, though poorly constrained (Holland and Jenkins, 1999; Holland and Feltham, 2006; MacAyeal, 1984; Nicholls et al., 2006). Boundary conditions are Dirichlet at the bottom boundary, assigned to the far-field temperature, salinity, and geostrophic velocity. The bottom third of the domain is a sponge layer (Rayleigh damping) within which velocity, temperature, and salinity are relaxed toward their assigned values at the bottom of the domain (Klemp and Lilly, 1978; Maronga et al., 2015). The sponge layer results in negligible vertical fluxes of heat, salt, or momentum across the bottom boundary because scalar and velocity gradients go to 0. The flow is periodic along the x and y dimensions.

We present two sets of simulations which have a base case in common. The base case has low far-field thermal driving of 0.15 ^∘C, a slope of 1^∘, and vigorous far-field inertial oscillations of 20 cm s⁻¹ generated by a 0.03 Pa m⁻¹ pressure gradient. These conditions were chosen to favor an energetic regime, for reasons discussed in Sect. 4.1. To examine the relationship between thermal driving and melt rate, as well as turbulent flow characteristics, we vary the far-field thermal driving from 0.15 to 0.6 ^∘C in the first set of simulations. The far-field thermal driving, Δθ_∞, is defined as the difference between the far-field temperature (θ_∞) and the freezing temperature based on the far-field salinity (θ_f(S_∞)). In the second set of simulations, we reduce the slope from 1 to 0.01^∘. The ice base always slopes to the east in the positive x dimension of our domain; thus β=90^∘. The far-field salinity for all runs is 35 g kg⁻¹. Initially, both temperature and salinity increase with depth over the upper two-thirds of the domain. The background salinity stratification dominates the density stratification, with an inverse stability ratio $R_{\mathit{\rho }}^{*}$ of 20.

We initiate turbulence over the first 50 min of the simulation with perturbations to the horizontal velocity components on the order of 0.01 m s⁻¹. The simulation duration is 50 h, corresponding to approximately four inertial periods of 13 h. PALM employs adaptive time-stepping; time-steps for the simulations presented here range from 0.5 to 2.75 s after 2 h. By the end of our simulations, the time-mean kinetic energy of the flow has reached steady state. However, the turbulent intensity for all cases continues to decline, with more pronounced turbulence kinetic energy (TKE) loss at lower slopes (Fig. 2a, d). Unless stated otherwise, the results are presented as averages over the last simulated inertial period and over the domain excluding the sponge layer.

Figure 2 Time evolution of (a, d) domain-averaged, resolved turbulence kinetic energy (TKE), (b, e) melt rate, and (c, f) friction velocity for (a–c) thermal driving simulations and (d–f) variable slope simulations. The black curve represents the same simulation in all panels in this and subsequent figures. The analysis window, the last inertial period, is shaded.

We compute an effective thermal exchange coefficient that differs from that employed by the sub-grid scheme to represent the efficiency of heat exchange that may need to be represented to produce accurate melt rates in an ocean model that does not have the vertical resolution or sophisticated turbulence closure that we employ here. This derived thermal exchange coefficient, Γ_T,der, can be computed using Eqs. (10) and (16) from the simulated melt rate and ocean properties at any depth below z_½, here chosen at −2 m. We substitute the friction velocity computed by Eq. (8) with one computed using a quadratic drag law from the velocity simulated 2 m below the boundary and the applied drag coefficient, consistent with drag implementations in coarse-resolution ocean models.

3 Results

3.1 Overview of the mean simulated state

Melt rates decline over the course of the simulations (Fig. 2), preventing the identification of steady-state melt rates for most runs. This decline is in response to the concomitant increase in stratification during the course of the simulations, which decreases vertical heat fluxes by reducing vertical turbulent fluctuations. We do not continue our simulations beyond four inertial periods in the hope of reaching steady-state conditions because the increase in stratification reduces the turbulent length scales and necessitates higher model resolution than we could computationally afford. Our simulations generally show resolved turbulent fluxes exceeding sub-grid turbulent fluxes by at least a factor of 2, but sub-grid turbulent fluxes dominate within several meters of the ice base where the stratification is strongest (Fig. S2). Resolved turbulent fluxes are also comparable to sub-grid turbulent fluxes late in the simulation of low-slope cases after significant TKE has been lost (≤0.1 ^∘C slopes, Fig. S2d). We evaluated the effects of doubling both horizontal and vertical resolutions in a separate simulation with an otherwise identical setup to the base case, which was run for one inertial period after spin-up. While a greater portion of the vertical heat fluxes was resolved as expected, the melt rate averaged over one inertial period only increased by 7 %, and differences in the mean state were sufficiently small to justify the use of our standard resolution (Fig. S3).

We present depth profiles of temperature, salinity, and velocity at the end of all simulations (Fig. 3). The time-averaged far-field velocity shown in Fig. 3c and f removes a periodic signal from inertial oscillations of the geostrophic velocity with a characteristic magnitude of 20 cm s⁻¹. Ekman rotation near the boundary can be seen in all simulations (Fig. 3c, f), but for more strongly sloped runs, buoyancy plays an increasingly strong role in driving the mean flow near the boundary. This effect can be seen most clearly by comparing the up-slope component of the flow within several meters of the ice–ocean interface across the slope-varying simulations (Fig. 3f). These mean buoyancy effects increase flow velocities on the order of a few centimeters per second as slope varies from 0.01 to 1^∘. Flow velocities near the boundary increase on the order of 10 cm s⁻¹ as far-field thermal driving increases from 0.4 to 0.6 ^∘C at 1^∘ slope (Fig. 3c). This velocity increase is attributed to changes in the magnitude of buoyancy forcing near the boundary which is related to differences in melt rates. The vertical momentum flux profiles shown in Fig. 4a and c reveal that flow is accelerated (negative fluxes) over much of the IOBL, with drag dominating only within the first few meters below the boundary.

Figure 3 Depth profiles of simulated properties as they vary with (a–c) thermal driving and (d–f) slope, averaged over the last inertial period. (a, d) The temperature relative to far-field temperature, (b, e) salinity relative to far-field salinity, and (c, f) velocity in the y direction (solid line, positive north, cross-slope) and in the x direction (dashed line, positive east, up-slope).

Figure 4 Vertical profiles of (a, c) momentum flux, (b, e) heat flux, and (c, f) scaled heat flux averaged over the last inertial period. The first row shows temperature cases, and the second row shows slope cases. Momentum flux is expressed in two components: $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$ (dashed) and $\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}$ (solid). Positive flux denotes upward flux (i.e., drag). The horizontal axis limits vary between panels.

All simulations show an evolution from the weakly stratified initial conditions to more strongly stratified conditions, particularly within the first 5 m of the ice–ocean interface (Fig. 3b, e). In none of the simulations do we observe a well-mixed boundary layer with respect to scalars. Rather, the simulations show varying degrees of stratification over the first 20 m from the boundary. Stratification within the boundary layer increases with thermal driving (Fig. 3a, b). Thus, the effect on stratification of the increase in melt-induced buoyancy fluxes with thermal driving dominates over the increase in shear induced by the buoyant flow. Conversely, stratification decreases with increasing slope, indicating that the increase in shear induced by the buoyant flow dominates over the increase in melt-induced buoyancy fluxes with slope.

3.2 Turbulent kinetic energy budget

Shear production of turbulence dominates the TKE budget. The evolution of TKE can be described by three source terms and dissipation.

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}e}{\mathrm{d}t}}& =\underset{F_{\mathrm{shear}}}{\underbrace{-\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}{\displaystyle \frac{\mathrm{d}u}{\mathrm{d}z}}-\stackrel{\mathrm{‾}}{v^{\prime }{w}^{\prime }}{\displaystyle \frac{\mathrm{d}v}{\mathrm{d}z}}}}-\underset{F_{\mathrm{buoy}}}{\underbrace{{\displaystyle \frac{g_{\mathrm{3}}}{\mathit{\rho }}}\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{\prime }{w}^{\prime }}-{\displaystyle \frac{g_{\mathrm{1}}}{\mathit{\rho }}}\stackrel{\mathrm{‾}}{{\mathit{\rho }}^{\prime }{u}^{\prime }}}}\\ & +\underset{F_{\mathrm{trans}}}{\underbrace{{\displaystyle \frac{d}{\mathrm{d}z}}(\stackrel{\mathrm{‾}}{w^{\prime }{p}^{\prime }}+{\stackrel{\mathrm{‾}}{w}}^{\prime }{e}^{\prime })}}-\underset{\mathrm{diss}}{\underbrace{\mathit{\epsilon }}}\end{array}\end{array}$$

Here, we characterize the evolution of resolved TKE, and primes designate resolved fluctuations. Figure 5b–d show the source terms in this budget for the variable slope simulations; the analogous figure for the thermal driving simulations is Fig. S4 which shows similar patterns. Shear production (F_shear) ranges from $\sim {\mathrm{10}}^{-\mathrm{9}}$–10⁻⁸ m² s⁻³ with a maximum at the ice–ocean interface and local maxima in the upslope flow within 5 m of the boundary (Fig. 5a). The increase in TKE throughout the boundary layer as slope increases reflects this shear-induced turbulence (Fig. 5a).

Figure 5 (a) Simulated turbulence kinetic energy (TKE) and (b–d) TKE source terms for variable slope simulations averaged over the last inertial period. (b) Shear production. (c) Buoyancy production: total (solid lines), vertical component (dashed), and upslope component (dotted). (d) TKE transport. Positive denotes production and negative destruction of TKE. Note that the x-axis scales differ between panels.

Buoyancy production of turbulence (F_buoy) is 1–2 orders of magnitude smaller than shear production: $\sim {\mathrm{10}}^{-\mathrm{10}}$–10⁻⁹ m² s⁻³. For a sloped ice shelf, the buoyancy term can be broken into two components: the horizontal buoyancy fluxes that increase turbulence when these fluxes are oriented upslope and the vertical buoyancy fluxes that decrease turbulence when ice is melting (Eq. 19). Since our simulations only produce melting, the vertical component is always negative (dashed lines, Fig. 5c). Over the coarse of an inertial oscillation, the horizontal component of buoyancy production is generally positive when the mean flow is oriented upslope and negative when the flow is oriented downslope. Interestingly, the time-averaged effect of a slope is destruction of turbulence near the boundary and production of turbulence near the base of the boundary layer (dotted lines, Fig. 5c). The net effect of both horizontal and vertical buoyancy components is the destruction of turbulence throughout the IOBL (solid lines, Fig. 5b), with the exception of the base of the IOBL where the horizontal buoyancy component augments entrainment. While these complexities are intriguing from a dynamical perspective, we do not explore them further here since the buoyancy term is not of leading order in the TKE budget.

The transport term (F_trans) contains two components: advection of TKE due to pressure fluctuations and turbulent advection of TKE. The former is negligible, reaching maximum values on the order of 10⁻¹² m² s⁻³, while the latter increases turbulence at the boundary with oscillations of decreasing amplitude with distance from the boundary (Fig. 5d).

Dissipation (diss), the remaining term in the TKE budget, can be inferred from the remainder of these terms and the rate of change of TKE. Since we evaluate resolved TKE (e), dissipation here represents the transfer of energy to the sub-grid scales; there is no flow of energy from sub-grid to resolved scales in the turbulence closure scheme. In these simulations TKE is not in steady state, with an average dissipation rate over the course of a simulation of 𝒪(10⁻¹¹) m² s⁻³, as seen in Fig. 2a and d. Dissipation in the IOBL is on the same order as shear production (10⁻⁹ m² s⁻³) as all other terms in the TKE budget are small.

3.3 Boundary layer turbulence

To demonstrate the simulated turbulent structures in this regime, we present horizontal and vertical cross-sectional snapshots through the domain for the two slope end-members in Fig. 6. Turbulent structures within the IOBL are consistent with propagating Holmboe shear instabilities under stable stratification (Carpenter et al., 2010). Shear is stronger within the IOBL than at the base of the IOBL due to the concentration of buoyant plume flow near the top of the IOBL. Consequently, the amplitude of these structures increases near the boundary for the more strongly sloped simulations (Fig. 6). The difference in IOBL turbulence with slope is perhaps best seen in the turbulent structures at 1 m from the ice–ocean interface (Fig. 6e, f). The structures become increasingly filamentous (i.e., near-wall streaks, e.g., del Álamo and Jiménez, 2003; Hoyas and Jiménez, 2006) and coherent as slope and thus the velocity of the buoyancy-driven current increase. Since the stratification decreases with increasing slope, the ratio of vertical to horizontal velocity variance also increases; velocity fluctuations are less confined to the slope-parallel direction (Fig. S5).

Figure 6 Instantaneous flow structures observed at 40 h in the (a, c, e) 0.01^∘ slope case and (b, d, f) 1^∘ slope case. (a, b) Temperature in cross section mid-way through the y axis. (c, d) Salinity in cross section mid-way through the y axis. (e, f) Resolved cross-slope velocity fluctuations at 1 m below the ice–ocean interface.

We define the base of the IOBL as the depth where the thermal driving relative to the far-field freezing point is 99 % of the far-field thermal driving. The IOBL depth increases through time, reaching 13–19 m at the end of the simulations. The simulated boundary layer depth increases with the far-field thermal driving (from 13 to 19 m) and with ice-shelf basal slope (from 15 to 19 m), reflecting the increase in flow velocities across those parameter changes which drives entrainment into the IOBL. Figure 7 shows the temporal evolution of IOBL depth for a few simulations, with steady IOBL growth at low slopes and punctuated growth at higher slopes corresponding to trends in TKE.

Figure 7 Horizontally averaged turbulence characteristics for (a) 0.15 ^∘C thermal driving and 0.01^∘ slope case and (b) 0.6 ^∘C thermal driving and 1^∘ slope case. Turbulence is considered intermittent when the dissipation contour of 10⁻⁹ m² s⁻³ (dashed line) reaches the boundary. Significant turbulence kinetic energy (TKE, green-yellow shading) can be present when dissipation is low. Higher TKE at −25 m at later times in (a) is due to turbulent transport, while the higher TKE at later times in (b) is due to shear production. IOBL depth is contoured in black.

Given the temporal variability in TKE in these simulations, we use a threshold in dissipation rate to characterize the mixing layer depth as being distinct from the mixed layer depth. This criterion has been deployed for stratifying ocean boundary layers (Franks, 2015; Sutherland et al., 2014). We consider the mixing layer as the depth interval over which the horizontally and hourly averaged dissipation rate exceeds 10⁻⁹ m² s⁻³. This mixing layer depth shows greater temporal variability than the IOBL depth, as defined by scalar concentration, and drops below the IOBL depth during periods of enhanced entrainment (Fig. 7).

There are also time periods over which the dissipation rate drops below 10⁻⁹ m² s⁻³ throughout the water column, indicating intermittency in turbulence. We note that these intervals of low dissipation do not correspond to complete loss of TKE nor a shoaling of the IOBL depth as we have defined it (Fig. 7). At higher slopes ($\ge \mathrm{0.5}{}^{\circ }$), this intermittency is interrupted as inertial oscillations enhance up-slope IOBL flow, increasing shear-driven mixing. This mixing front begins at the boundary and propagates to the base of the IOBL over a few hours (Fig. 7b). For the less stratified, low-slope cases, the intermittency in turbulence is more frequent as the IOBL flows more slowly and generates less TKE production by shear (Fig. 7a). We discuss this intermittency and its possible implications further in Sect. 4.1.

3.4 Melting and its relation to thermal driving and slope

Average melt rates over the last inertial period range from 0.3 to 2.5 m (Fig. 2c, f), and average Monin–Obukhov length scales are 3–5 m (not shown). The dominant temporal frequency in melt rate is the inertial frequency, and the melt response to those oscillations is highly nonlinear. Maximum melt rates occur when the mean flow is oriented between the up-slope direction and the Coriolis-favored direction, and minimum melt rates are roughly 180^∘ opposed to that. These melt rate fluctuations correspond to fluctuations in TKE which are reflected in the friction velocity shown in Fig. 2b and e. The dominant contribution to turbulence is shear production of TKE, which changes its distribution with depth as the mean flow profile evolves. During high-melt periods, the far-field flow is oriented up-slope, and shear production of TKE is concentrated near the boundary. During low-melt periods, the far-field flow is oriented down-slope, and shear production of TKE is concentrated a few meters away from the boundary. Melt rate fluctuations increase in amplitude as thermal driving increases and as slope increases, which we attribute to the increasing importance of buoyancy forcing in driving a near-boundary plume and thus determining the depth distribution of shear. The two simulations associated with the highest thermal driving values, and the highest near-boundary stratification, experience a dramatic reduction in melt rates during the down-slope flow period. This coincides with reduced friction velocity (i.e., shear stress) at the ice interface (Fig. 2b) and reduced vertical velocity fluctuations (not shown).

Time-averaged melt rates depend fairly linearly on far-field thermal driving (Fig. 8a, R²=0.97). This is mostly attributable to a linear relationship between melt rate and interfacial thermal driving as differences in friction velocity and the thermal exchange coefficient are small and do not have a systematic relationship with far-field thermal driving. On the other hand, the derived, time-averaged thermal exchange coefficient representing the efficiency of heat exchange from −2 m depth to the ice base (Γ_T,der) does have a weakly negative relationship with far-field thermal driving (Fig. 8c). This indicates a decrease in the efficiency of heat exchange with increasing near-interface stratification. The highest thermal driving case shows an anomalously high thermal exchange coefficient over the last inertial period, which features high TKE shown in Fig. 7b. We discuss these derived thermal exchange coefficients further in Sect. 4.2.2.

Figure 8 Melt rate sensitivity to (a) far-field thermal driving and (b) sine of the basal slope. (c) Far-field thermal driving is inversely related to Γ_T,der. Dashed line denotes the value recommended by Jenkins et al. (2010). The largest points correspond to the fourth inertial cycle with progressively smaller points for previous inertial cycles. (d) For variations in slope, the simulated friction velocity (solid points) is linearly related to Γ_T,der. The inferred friction velocities used to compute Γ_T,der are shown with open points. Note the difference in y-axis limits from (c). The 0.01 and 0.1^∘ slope cases are overlapping.

There is also a linear relationship between melt rate and ice-shelf basal slope, with threshold-like behavior at slopes less than 0.01^∘. This linear relationship is due primarily to a linear relationship between friction velocity and melt rate, while there is a smaller ($\sim \mathrm{1}/\mathrm{3}$ size) opposite effect from decreasing interfacial thermal driving with increasing slope (not shown). These differences in friction velocity arise from higher IOBL velocities and turbulence at higher slopes. We observe threshold behavior in melting in the two lowest-slope cases, which can be attributed to similar friction velocities arising from similar IOBL velocities and turbulence (Figs. 8d and 3f). This behavior is discussed in more detail in Sect. 4.2.1. The differences in the derived thermal exchange coefficients with varying slopes are on the order of a 40 % change as slope increases from 0.01 to 1^∘ (Fig. 8d).

3.5 Vertical structure of turbulent fluxes

Vertical heat fluxes are shown in Fig. 4b and e. The vertical heat flux has a maximum at the ice–ocean interface and decreases throughout the IOBL, with small values below the pycnocline. Since the IOBL is fully turbulent (with the exception of some intermittency discussed below), the sub-grid diffusivities of momentum, heat, and salt closely resemble one another (Fig. S6). The vertical salt flux profiles are shown in Fig. S7 and have a very similar shape to the vertical heat flux profiles.

The relatively narrow range of conditions simulated here suggests that a depth-dependent shape function for scalar fluxes could be formulated. The distance from the interface is scaled by the Ekman depth:

$$\begin{array}{}\text{(20)}& d_E=\left(\mathrm{2}{K}_{\mathrm{e}}{)}^{\mathrm{1}/\mathrm{2}}\right|f_{\mathrm{3}}{|}^{-\mathrm{1}/\mathrm{2}},\end{array}$$

where K_e is the mean eddy viscosity in the turbulent boundary layer assuming the total fluxes follow Fick's law:

$$\begin{array}{}\text{(21)}& \stackrel{\mathrm{‾}}{w^{\prime }{\mathit{u}}^{\prime }}=K_{\mathrm{e}}\mathrm{d}\mathit{u}/\mathrm{d}z.\end{array}$$

K_e profiles are shown in Fig. S8.

The linear scaling of melt rate with thermal driving suggests a linear scaling of heat flux profiles with thermal driving. We find that this scaling largely collapses the four thermal driving profiles, with notable deviation from this shape for the run that experiences temporal gaps in shear stress during the analysis period (pink curves in Figs. 2a, 4c). The shape of these profiles can be reasonably approximated by a linear decrease in the scaled heat flux with scaled depth over the boundary layer (Fig. 4c). Scalar fluxes decline near the boundary for the 1^∘ slope cases, a feature that we discuss further in Sect. 4.2.3.

Despite melt rates scaling reasonably well with the basal slope (sin α), the vertical heat flux profile does not, showing a much lower sensitivity to slope ($(\sin \mathit{\alpha }{)}^{\mathrm{1}/\mathrm{4}}$, Fig. 4f). There is strong agreement between the scaled vertical heat flux profiles at different slopes, though the threshold behavior at low slopes noted in melt rates is replicated here. We also find that the Ekman depth is a poor predictor of boundary layer depth. This may be due to the depth-variable shear induced by buoyancy which is not reflected in the depth-mean IOBL eddy viscosity used to compute the Ekman depth. We discuss the scaling of vertical fluxes and the possibility of their parameterization in Sect. 4.2.3.

4 Discussion

4.1 Understanding IOBL turbulence and the limitations of a LES approach

In the simulations presented here, turbulence declines throughout the course of the simulation, becoming intermittent. The relationship between stable stratification, shear, and the persistence of turbulence remains an open question (Zonta and Soldati, 2018). Thus, it is not possible a priori to determine whether the level of turbulence simulated by this LES model is appropriate for the regime space we have sampled. While we did conduct model validation against a stably stratified atmospheric boundary layer test case (see Sect. 2.2), the degree of stratification at the boundary in that case did not approach that simulated in the sub-ice configuration. Fundamentally, we cannot guarantee that our LES is not overly dissipative such that the TKE generated by the resolved dynamics is lost too quickly relative to real-world sub-ice settings.

Excess dissipation could arise either through the sub-grid scheme or the model numerics. We found a rapid loss of turbulence in PALM simulations when a dynamic Smagorinsky turbulence closure was used, which is consistent with previous studies on the limited applicability of the Smagorinsky turbulence closure to strongly stratified flows due to the strong anisotropy of those flows (Flores and Riley, 2011; Jiménez and Cuxart, 2005). This motivated our adoption of the AMD turbulence closure scheme (Abkar et al., 2016). However, we found that the buoyancy term added to this scheme by Abkar and Moin (2017) had an unrealistically high magnitude in the vicinity of the ice base where gradients are large. We removed this term, as it was negligible in the simulations of Vreugdenhil and Taylor (2019) (Catherine Vreugdenhil, personal communication, 2020​​​​​​​), but our unrealistic solution for this term suggests that the AMD scheme may not perform optimally at the resolution employed here. Though the resolution of our simulations is significantly lower than that used by Vreugdenhil and Taylor (2019), it enables us to extend the domain from their 2 m to 64 m to allow for the development of a thick IOBL.

Strongly stratified turbulence has been associated with intermittent turbulence (Nieuwstadt, 2005; Wiel et al., 2012), though there are also numerical experiments that fail to produce intermittency even under strongly stable stratification (Arya, 1975; Komori et al., 1983). This is not the first study to find the emergence of intermittent turbulence in stably stratified, sub-ice settings (Vreugdenhil and Taylor, 2019). Donda et al. (2015) have argued that, in strongly stratified flows, the cessation of turbulence is transient provided there are sufficiently large perturbations. This is consistent with our finding that the temporal variability in shear over inertial oscillations provides sufficient perturbations to reinitiate turbulence. Our simulations approach a gradient Richardson number (Ri_g) of 0.25, which is considered to be the approximate value at which turbulence neither grows nor decays (Holt et al., 1992; Rohr et al., 1988). Thus, fluctuations in TKE are plausible at the simulated levels of stratification and shear. However, when turbulence is intermittent, as it is in this study, the application of LES may be inappropriate due to its inherent horizontal averaging over laminar and turbulent regions (Stoll and Porté-Agel, 2008). Thus, our results should be interpreted with caution.

In this study, we did not attempt to reproduce observed conditions at a particular ice-shelf location due to the difficulties of matching unobserved far-field forcings and the exclusion of tides from our simulations. Nonetheless, it appears as though well-mixed boundary layers are seen for a narrower range of conditions in simulations than in observations. The geostrophic flow chosen in these simulations is quite strong at 20 cm s⁻¹, and thermal driving is relatively low such that we expected to produce a well-mixed boundary layer as is observed in melting regions of the Filchner–Ronne Ice Shelf (Nicholls et al., 2001). However, observations to date are insufficient to fully characterize the regime space for a well-mixed IOBL (Malyarenko et al., 2020). The observational picture is quite nuanced with a range of stratification observed even within one ice shelf (Hattermann et al., 2012).

Shear-driven turbulence within the IOBL plays a central role in determining vertical heat fluxes and thus melt rates. In these simulations the destruction of TKE by the stabilizing buoyancy flux is 2 orders of magnitude smaller than shear production of TKE throughout the IOBL. Our finding that shear production of TKE dominates over the buoyancy term is consistent with Davis and Nicholls (2019) who found that shear production was an order of magnitude greater than buoyancy destruction of TKE in the IOBL below Larsen C Ice Shelf.

The flux Richardson number, Ri_f, the ratio of buoyancy flux to shear-driven TKE production, provides a measure of mixing efficiency. The simulated Ri_f values of 0.05–0.1 are well below the critical Ri_f of ∼0.25, indicating that we are in the regime in which mixing efficiency (Ri_f) decreases with increasing stratification (Ri_g) (Miles, 1961; Howard, 1961; Armenio and Sarkar, 2002; Peltier and Caulfield, 2003). This is consistent with our finding that the derived 2 m depth thermal exchange coefficient decreases for the high thermal driving cases which also achieve stronger IOBL stratification.

Our simulations are certainly missing some sources of TKE present in ice-shelf cavities which could modify the IOBL structure and mixing efficiency. These simulations did not include tides, which provide perturbations to the mean velocity that enhance melt rates and entrainment, especially at Filchner–Ronne (Makinson and Nicholls, 1999; Makinson et al., 2011; Mueller et al., 2018). While internal gravity waves arise in LES, in our LES model they may be of smaller amplitude and play a lesser role in mixing than they would in a real ice-shelf settings due to the absence of large-scale external forcings such as tides and storms, seafloor topography, and possible resonance with the cavity geometry (Gwyther et al., 2020; Mueller et al., 2012; Padman et al., 2018; Robertson, 2013). Enhanced drag at the ice-shelf base could increase shear production of TKE; however, under strong stratification, surface roughness elements may suppress turbulence rather than enhance it (Ohya, 2001).

4.2 Representing the IOBL and projecting ice-shelf melt rates with ocean models

4.2.1 Insights into melt rate sensitivity to ocean conditions

The decline in IOBL turbulence discussed in Sect. 4.1 results in declining melt rates. Thus, we could not evaluate the relationship between melting and far-field conditions at steady state, which would have offered the most direct path to assessing melt rate sensitivity to ocean conditions. As discussed in Jenkins (2016), achieving steady-state solutions in simulations may require prescribing large-scale gradients in temperature. We have not included these large-scale gradients in our simulations, but this may be an avenue for future work. However, we believe that our transient solutions do provide some indications of how melt rates and boundary layer properties depend on ocean temperature and ice-shelf slopes. One justification for this belief is that the simulated sensitivity of melt rate to ocean temperature remains consistently linear through all inertial periods simulated (Fig. S9). On the other hand, the differently sloped simulations continue to diverge at the end of the simulation as the IOBL accelerates (Fig. 2d–f). Nonetheless we are able to find relationships between the vertical heat flux profiles and ocean temperature and slope that suggest predictability of the mean effects of turbulence despite the transience of these simulations.

The linear relationship we find between local, far-field ocean temperature and melt rates is consistent with some previous studies in ice-shelf settings (Holland et al., 2008; Rignot and Jacobs, 2002; Vreugdenhil and Taylor, 2019). A slightly higher exponent of $\mathrm{4}/\mathrm{3}$ is also consistent with our data (R²=0.95), a value reported for the regime in which convective instabilities control melting, while our simulations feature shear instabilities (Kerr and McConnochie, 2015). In contrast, in sea-ice settings this sensitivity of melt rates to local thermal driving was found to be significantly smaller with an exponent of 0.38 (Ramudu et al., 2018). The relationship between average ice-shelf melt rates and the thermal driving for the cavity as a whole can be conceptualized in two components. The first is the relationship between the local thermal driving and melt rate, determined primarily by the local ocean turbulence. The second is the relationship between distant thermal driving (i.e., the water masses entering the ice-shelf cavity) and the strength of sub-ice-shelf circulation (Holland et al., 2008). Only the former is addressed by this study. Our simulations do not capture the large-scale increase in overturning circulation that accompanies an increase in distant thermal driving. Our simulations partially capture the increase in IOBL velocity due to an increase in thermal driving, but the far-field velocity is unaffected. We note that the often cited quadratic relationship between distant thermal driving and melt rate involves both local turbulent processes and large-scale processes (Holland et al., 2008); studies generally find this exponent on distant thermal driving to be between 1.5 and 2 (Favier et al., 2019; Jourdain et al., 2017; Little et al., 2009).

To some extent, this linear relationship between far-field thermal driving and melt rate is embedded in the parameterization of heat fluxes employed at the ice base in these simulations. Specifically, the melt rate is prescribed to have a linear dependence on the local thermal driving, based on the temperature in the first model layer (Eq. 10). However, the level of stratification near the ice base, the buoyant acceleration of the IOBL plume, and the transport of heat from depth to the IOBL all mediate the relationship between far-field thermal driving and melt rates, and yet the dependence of melt rate on far-field thermal driving from these simulations is still fairly linear.

The relationship between melt rate and ice-shelf basal slope (m∝(sin α)ⁿ) combines two effects: the effect of ice-shelf slope on the IOBL's mean velocity profile and the effect of ice-shelf slope on IOBL turbulence. In the three-equation parameterization of ice-shelf melting, the former has the strongest effect on the friction velocity as derived from the mean flow velocity at a given depth, whereas the latter is represented by the scalar exchange coefficients, with higher values indicating more efficient turbulent transport. We return to the implications of our study for exchange coefficients in Sect. 4.2.2. As noted previously in this section, while our simulations capture some of this IOBL acceleration, they do not reach a steady-state mean velocity profile. Thus, we cannot fully assess the effect of ice-shelf slope on the IOBL's mean velocity profile. This effect is addressed by Magorrian and Wells (2016), whose scaling analysis predicted $n=\mathrm{3}/\mathrm{2}$, and Little et al. (2009), who found n=0.94 for the range of slopes considered here. We find that n=1 may be applicable to the shear-dominated regime, in contrast to the sensitivity of $n=\mathrm{2}/\mathrm{3}$ found in the convection-dominated regime at higher slopes than those simulated here (5–90^∘; McConnochie and Kerr, 2018; Mondal et al., 2019). However, we emphasize that further investigation is needed beyond the small number of simulations presented here to validate our results.

The threshold-like behavior in melt rates at very low slopes is not predicted by geostrophic balance between Coriolis and buoyancy forcing, which dictates a linear relationship between sin α and IOBL velocity (Jenkins, 2016), nor the scaling analysis of Magorrian and Wells (2016). This threshold behavior could be produced if any additional buoyant acceleration produced by the small increase in slope from 0.01 to 0.1^∘ increases both shear and dissipation, resulting in a negative feedback on IOBL turbulence. This is not evident in our simulations, which do not show significant differences in TKE budgets between the two runs (Fig. 5). From a dynamical perspective, this may be an interesting target for more highly resolved LES in the near-boundary region. However, this threshold is located at low enough slopes that it is likely not of significance to melt rate parameterization for coarse-resolution ocean models, so we do not devote further attention to it here.

4.2.2 Toward non-constant exchange coefficients in melt parameterization

The thermal exchange coefficient as computed using Eq. (11) at 0.125 m below the ice–ocean interface (the uppermost grid cell) differs from that which would be implemented by coarse-resolution models or derived from oceanographic observations, both of which only know ocean properties meters to tens of meters below the ice–ocean interface. To demonstrate the implications of this study for modeling endeavors, we computed the thermal exchange coefficient that yields the simulated melt rate using ocean properties (temperature and velocity) 2 m below the ice–ocean interface. This depth was chosen to capture the faster portions of IOBL flow, representing a best-case scenario for high-resolution ocean models, though this depth choice is somewhat arbitrary.

These derived thermal exchange coefficients are shown in Fig. 8c and d. The derived thermal exchange coefficients are all less than the value of 0.011 derived from observations (Jenkins et al., 2010). There are two main factors that contribute to this result. The first is the choice of the parameterization of fluxes at the ice–ocean interface (i.e., over sub-meter scales). We chose a stability-dependent parameterization in which the buoyancy forcing from melting enters through the Monin–Obukhov length (Eq. 9). We believe there is a stronger case for the flux parameterization we implement than for a constant thermal exchange coefficient in light of the success of Monin–Obukhov similarity theory (Monin and Obukhov, 1954; McPhee, 2008), as well as the depth dependence of scalar fluxes in a more highly resolved sub-ice LES (Vreugdenhil and Taylor, 2019). Consequently, Vreugdenhil and Taylor (2019) also found that thermal exchange coefficients were less than the Jenkins et al. (2010) value for all but the lowest thermal driving case. However, we acknowledge that more validation of the sub-grid boundary flux parameterizations is needed. The second contributing factor is declining simulated TKE, which reduces thermal exchange coefficients over the course of these simulations (Fig. 8c).

Our results also suggest a modest decline in the thermal exchange coefficient at higher thermal driving. This may lead to a sub-linear relationship between thermal driving and melt rate at higher thermal driving values, though this is not evident in our simulations. The decline in exchange coefficient with thermal driving agrees with Vreugdenhil and Taylor (2019), though their sensitivity is greater than that seen in this study (their Fig. 8). As shown in Fig. 8c, this relationship holds during all inertial periods, with an exception during an interval of strong turbulence (see Sect. 3). Due to the limitations of our LES modeling discussed in Sect. 4.1, we cannot recommend a best fit relationship between thermal exchange coefficient and thermal driving. Given the climatic importance of accurately simulating high thermal driving regimes associated with dynamic ice-shelf thinning, LES coupled with observational validation across thermal driving regimes may be a fruitful avenue for future work.

We find a linear relationship between the derived thermal exchange coefficient and the slope of the ice-shelf base, indicating that mixing is more efficient at higher slopes. This relationship also holds for the derived friction velocity from a quadratic drag law (Fig. 8d). Note that these quadratic-drag friction velocities are greater than the model's parameterization of friction velocity as the former neglects stratification effects, while the latter includes them. This enhanced turbulent efficiency is due to enhanced shear instabilities and is reflected in turbulent diffusivities parameterized by the AMD scheme. Much greater variability in sub-grid turbulent diffusivities is seen over the variation in slope than the variation in thermal driving tested here (Fig. S6).

The changes in the thermal exchange coefficient with slope were relatively small, from 0.0045 to 0.007 between 0.1 and 1^∘ cases. For a coarse-resolution simulation, we anticipate that the failure to capture this slope sensitivity will not be a leading source of error in melt projections. Reproducing an accurate friction velocity is likely of a greater concern. Thus, the burden of melt projection accuracy may fall more heavily on parameterizing or resolving buoyant flow than on improving the slope dependence in the parameterization of scalar fluxes (i.e., improving Eq. 10).

These simulations do not reveal whether a similar, linear relationship between friction velocity and thermal exchange coefficient holds when the background flow is varied rather than the slope. LES of sea-ice melting, for which there is no slope, suggests a sublinear relationship between thermal exchange coefficient and friction velocity (${\mathrm{\Gamma }}_T\propto u_{*}^{\mathrm{0.5}},\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}^{\prime }}\propto u_{*}^{\mathrm{1.5}}$; Ramudu et al., 2018). More studies are required to determine this scaling.

4.2.3 Toward a vertical mixing scheme for the IOBL

There is reason to believe that improving ice-shelf basal melt projections using ocean models will require not only an accurate melt parameterization but also an improved vertical mixing scheme. Jenkins (2021) demonstrated significantly different IOBL characteristics when KPP is the turbulence closure scheme in contrast to a low- or high-order scheme. The eddy viscosity simulated by our LES model shows better agreement with the viscosity solutions from Jenkins (2021) employing the low- and high-order turbulence closure schemes than that employing KPP (compare our Fig. S6 with his Fig. 3). Thus, this study offers additional support for the use of a more sophisticated turbulence closure scheme or a modified KPP scheme in sub-ice settings, such as one based on the depth dependence of vertical turbulent fluxes presented here (Fig. 4).

A complication in melt parameterization of slope effects is the differential ability of ocean models and their vertical grid configurations to capture IOBL flow. We find that this buoyant flow is concentrated in the uppermost 10 m, with peak velocities as close as 3 m from the interface. This flow is unlikely to be resolved by most ocean models, which have typical vertical resolutions near the interface of ∼10 m, though some model configurations are reaching ∼1 m resolutions (Gwyther et al., 2020). If the buoyancy-driven boundary current is unresolved in most coarse-resolution ocean models, then the friction velocity as computed from a quadratic drag law and the velocity at the first grid cell is likely to underestimate the true friction velocity. For instance, a simulation that lacks this buoyancy-induced current may look more like our negligible (0.01^∘) slope simulations than those with a slope and may consequently underestimate melt rates. Thus, more accurate melt rate projections may result from employing a depth-dependent parameterization of vertical scalar and momentum fluxes to provide some sensitivity of melt rates to model resolution.

Our simulated depth dependence for vertical scalar fluxes agrees with the shape of the turbulent scalar flux parameterization of McPhee et al. (1987), which features a linear decrease in fluxes with distance from the boundary within the IOBL. The low thermal-driving simulations with slopes of 1^∘ diverge from this linear behavior within a few meters from the boundary. We attribute this to boundary effects that reduce TKE and eddy viscosity close to the boundary, reminiscent of the constant flux layer predicted by Monin–Obukhov similarity theory (Monin and Obukhov, 1954). It is unclear whether a constant flux layer is absent under some conditions or whether our simulations do not resolve it for those highly stratified cases. Since these differences in eddy viscosity are incorporated into the Ekman depth scaling, the scaled heat flux profiles all show quasi-linear slopes (Fig. 4c, f). This linear depth scaling for IOBL fluxes should be explored under a broader range of conditions and at higher resolution to determine whether it is broadly applicable and whether vertical mixing parameterizations need to improve their representation of the constant flux layer in the IOBL.

A few complexities are immediately apparent in pursuing a depth-dependent shape function for vertical scalar flux parameterization. One is that a functional form for the eddy diffusivity is also needed to compute the Ekman depth. Another is that any implementation of these depth-dependent scalar fluxes necessitates an inequality between the vertical heat and salt fluxes from the top model grid cell and the heat and salt fluxes associated with melting. Thus, the closure of those budgets will likely involve horizontal scalar fluxes as well.

A depth-dependent function for momentum is also necessary in ice-shelf settings given that strong stratification causes momentum fluxes to significantly deviate from the quadratic drag function typically employed in ocean models. Near-interface stratification significantly decreases drag (García-Villalba and del Álamo, 2011; McPhee et al., 2008). Tuning the drag coefficient to fit present-day melt rates neglects the functional dependence of drag on melt rate through stratification (as represented by the Monin–Obukhov length). The presence of an ice-shelf slope in ocean models that do not fully resolve a buoyant plume further complicates the parameterization of momentum fluxes at the ice–ocean interface. For instance, a depth-dependent shape function for momentum fluxes may prescribe negative (downward) momentum fluxes at the top grid cell below ice shelves for typical ocean model resolutions (see Fig. 4a and d at, for example, 5 m depth). In other words, when the buoyant plume is unresolved, the momentum boundary condition may need to accelerate the flow at the boundary (at least in the up-slope direction) to produce more accurate shear-driven mixing in the resolved portion of the domain.

5 Conclusions

In this study we presented a small parameter study of IOBL turbulence as captured by a LES model and chiefly explored the sensitivity of ice-shelf melting to ocean temperature and ice-shelf basal slope. To our knowledge, this is the first study to explore turbulent effects of low slopes on ice melting at scales spanning tens of meters from the ice base, allowing for IOBL development. These vertical scales allow us to examine the depth dependence of vertical turbulent fluxes, a key component for the development and validation of vertical mixing parameterizations in ice-shelf settings. Building on Jenkins (2021), our results suggest an improvement on the standard KPP vertical mixing scheme is needed, but a larger simulation campaign and field validation effort would be required to place such a scheme on strong footing.

The thermal exchange coefficient is a key parameter that determines the sensitivity of ice-shelf melting to changing ocean conditions in ocean models. While a constant value is not theoretically supported, it remains the standard choice due to a lack of consensus regarding the appropriate scaling. We concur with a previous study that the thermal exchange coefficient is dependent on thermal driving such that melt rates are less sensitive to changes in ocean temperature (Vreugdenhil and Taylor, 2019). We also suggest that the thermal exchange coefficient should be enhanced at higher slopes to account for shear instabilities that are unresolved by coarse-resolution ocean models, with the caveat that our simulations did not reach steady state.

Ultimately, we acknowledge that there is currently no universal guidance for parameterizing sub-grid scalar and momentum fluxes near the ice base for ocean models, each of which capture IOBL flow differently depending on resolution, vertical coordinate, numerical representations of the pressure gradient, and likely other aspects of the implementation (Gwyther et al., 2020). Thus, progress toward accurate ice-shelf melt projections will require not only theoretical advances in understanding melt rate sensitivity but also advances in numerical modeling to bring resolution (scale) awareness to melt parameterization.