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 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. 5 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 10 melting worthy of further study.


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 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).

Overview of the LES model
The PArallelized Large eddy simulation Model (PALM) was developed at the Institute of Meteorology and Climatology at 60 Leibniz Universitat 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). 65 The governing equations of PALM are the following: Volume conservation ∂u j ∂x j = 0 Heat conservation The momentum terms on the right hand side of Equation 1 are, in order: advection, Coriolis forcing, imposed geostrophic flow, buoyancy forcing, a correction for divergence in the flow using the pressure perturbation π * (imposing incompressibility), 75 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, 80 g = g z [sin α sin β, sin α cos β, cos α] (5) f = 2Ω[sin φ sin α sin β, cos φ sin α cos β, sin φ cos α] 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. 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 Equation 1 combines the contributions of 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 .
ρ a varies along slope and with depth from the ice interface as a function of the hydrostatic pressure: ρ a (x, y, z) = f EOS θ i (z), S i (z), p(x, y, z) where the superscript i denotes the initial state. The reference density ρ 0 is evaluated in the center of the x-y plane (ρ 0 (z) = 90 ρ a (x mid , y 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 −1 . This simplification has the advantage of avoiding pressure discontinuities across periodic boundary.
The terms on the right-hand side of the heat and salt conservation equations (Equations 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 95 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 Section 2.2.

Turbulence closure
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 100 test case published in Abkar and Moin (2017). Typical sub-grid scale momentum and scalar diffusivities in our sub-ice-shelf simulations of stratified turbulence range from 10 −5 to 10 −4 m 2 s −1 in the upper two-thirds of the domain where damping is not applied (the damping methodology is discussed in Section 2.3).
Sources and sinks of momentum, heat and salt due to interactions between the flow and the ice base are all parameterized as sub-grid fluxes at the boundary. The resolved vertical fluxes at the top layer of the model go to zero as w goes to zero according 105 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 scalars and horizontal velocity components are located at mid-depth of the top cell (z = − 1 2 ∆ z ). We denote the interface location with subscript b and the middle of the first cell below the interface with subscript 1 2 .
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 w u i b , w θ b , w S b at the center of the first grid cell from the boundary without interpolation (i.e., w X b = w X 1 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 largely a constant flux layer, as hypothesized by McPhee (1983).

Simulation set-up
The domain is a 64 m 3 cube with horizontal resolution of 0.5 m and vertical resolution of 0.25 m. The ice-ocean interface is 150 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 an order of magnitude (roughly -2000 m to -200 m). The potentially dynamically relevant differences between conducting these simulations at surface pressure and 800 dbar is that the first derivative of the freezing temperature with respect to salinity is 20% smaller and the first derivative of density with respect 155 to salinity is about 2 times larger.
Boundary conditions for velocity, temperature and salinity are periodic at the side boundaries, Dirichlet at the bottom boundary, and von Neumann at the top boundary. 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 bot-160 tom boundary because scalar and velocity gradients go to 0. The von Neumann boundary conditions at the top of the domain correspond to the sub-grid fluxes of heat, salt and momentum (Equations 8,14,and 15), as resolved fluxes go to zero at a no-penetration boundary. The roughness length (z 0 in Equations 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). A list of parameter choices can be 165 found in Table S1.
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 and a relatively high slope for Antarctic ice shelves of 1 • and vigorous far-field inertial oscillations of 20 cm s −1 generated by a 0.03 Pa m −1 pressure gradient. These conditions were chosen to favor an energetic regime, for reasons discussed in Section 4.1. To examine the relationship between thermal driving and melt rate, as well as turbulent flow characteristics, 170 we vary the far-field thermal driving from 0.15 • C 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 −1 . Initially, both temperature and salinity increase with depth over the upper two-thirds of the domain. The background salinity stratification dominates the 175 density stratification, with an inverse stability ratio R * ρ of 20.
The simulation duration is 50 h, corresponding to ∼4 inertial periods of 13 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 (Figure 1a,d). Unless stated otherwise, the results are presented as averages over the last simulated inertial period and over the domain excluding the sponge layer. 180 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 Equations 10 and 16 from the simulated melt rate and ocean properties at any depth below z 1 2 , here chosen at -2 m. We substitute the friction velocity computed by Equation 8 with one computed using a quadratic drag law 185 from the velocity simulated 2 m below the boundary and the applied drag coefficient, consistent with drag implementations in coarse-resolution ocean models.

Results
Melt rates decline over the course of the simulations (Figure 1), 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 190 decreases vertical heat fluxes by reducing vertical turbulent fluctuations. We do not continue our simulations beyond 4 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 subgrid turbulent fluxes by at least a factor of two, but subgrid turbulent fluxes dominate within several meters of the ice base where the stratification is strongest ( Figure S1). Resolved turbulent fluxes are also comparable 195 to subgrid turbulent fluxes late in the simulation of low slope cases after significant TKE has been lost (≤ 0.1 • C slopes, Figure S1d).
We present depth-profiles of temperature, salinity, and velocity at the end of all simulations (

215
Shear production of turbulence dominates the TKE budget. The evolution of TKE can be described by three source terms and dissipation: Here, we characterize the evolution of resolved TKE and primes designate resolved fluctuations. Figure 4b-d shows the source terms in this budget for the variable slope simulations; the analogous figure for the thermal driving simulations is 220 Figure S2 and shows similar patterns. Shear production (F shear ) ranges from ∼ 10 −9 -10 −8 m 2 s −3 with a maximum at the ice- Buoyancy production of turbulence (F buoy ) is 1 -2 orders of magnitude smaller than shear production, ∼ 10 −10 -10 −9 m 2 s −3 .
For a sloped ice shelf, the buoyancy term can be broken into two components, the horizontal buoyancy fluxes that increase 225 turbulence when these fluxes are oriented upslope and the vertical buoyancy fluxes that decrease turbulence when ice is melting (Equation 19). Since our simulations only produce melting, the vertical component is always negative (dashed lines, Figure   4c). 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 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 230 lines, Figure 4c). The net effect of both horizontal and vertical buoyancy components is the destruction of turbulence throughout the IOBL (solid lines, Figure 4b), 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 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 240 in steady-state, with an average dissipation rate over the course of a simulation of O(10 −11 ) m 2 s −3 as seen in Figure 1a,c.
Dissipation in the IOBL is on the same order as shear production (10 −9 m 2 s −3 ), as all other terms in the TKE budget are small.
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 Figure  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 ( Figure S3). 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 255 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 6 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.
Given the temporal variability in TKE in these simulations, we use a threshold in dissipation rate to characterize the mixing 260 layer depth, as 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 horizontal-and hourly-averaged dissipation rate exceeds 10 −9 m 2 s −3 . This mixing layer depth shows greater temporal variability than the IOBL depth as defined by scalar concentration, and drops below that IOBL depth during periods of enhanced entrainment ( Figure 6).

265
There are also time periods over which the dissipation rate drops below 10 −9 m 2 s −3 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 (Figure 6). At higher slopes (≥ 0.5 • ), 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 (Figure 6b). For the less stratified, low-slope cases, the intermittency 270 in turbulence is more frequent, as the IOBL flows more slowly and generates less TKE production by shear (Figure 6a). We discuss this intermittency and its possible implications further in Section 4.1. an anomalously high thermal exchange coefficient over the last inertial period, which features high TKE shown in Figure 6b.
We discuss these derived thermal exchange coefficients further in Section 4.2.2.
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 (∼ 1/3 size) opposite effect from decreasing interfacial thermal driving with increasing slope (not shown).

295
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 (Figure 7d and Figure 2f). This behavior is discussed in more detail in Section 4.2.1. The differences in derived thermal exchange coefficient with varying slope are on the order of a 40% change as slope increases from 0.01 • to 1 • (Figure 7d). intermittency discussed below), the sub-grid diffusivities of momentum, heat and salt closely resemble one another ( Figure S4).
The vertical salt flux profiles are shown in Figure S5 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 305 could be formulated. The distance from the interface is scaled by the Ekman depth: where K e is the mean eddy viscosity in the turbulent boundary layer assuming the total fluxes follow Fick's law: K e profiles are shown in Figure S6.

310
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 Figures 1c, 3c). 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 ( Figure   3c). Scalar fluxes decline near the boundary for the 1 • slope cases, a feature that we discuss further in Section 4.2.3.

315
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 α) 1/4 , Figure 3f). 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 320 vertical fluxes and the possibility of their parameterization in Section 4.2.3.

Understanding IOBL turbulence and the limitations of an 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,325 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 Section 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.

330
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 Smagorinky 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)   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 345 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 350 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 −1 and thermal driving is relatively low such that we expected to 355 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 360 simulations the destruction of TKE by the stabilizing buoyancy flux is two 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 Our simulations are certainly missing some sources of TKE present in ice-shelf cavities which could modify IOBL structure 370 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, as well as possible resonance with the cavity geometry (Gwyther et al., 2020;Mueller et al., 2012;375 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).

Insights into melt rate sensitivity to ocean conditions
The decline in IOBL turbulence discussed in Section 4.1 results in declining melt rates. Thus, we could not evaluate the 380 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 385 belief is that the simulated sensitivity of melt rate to ocean temperature remains consistently linear through all inertial periods simulated ( Figure S7). On the other hand, the differently sloped simulations continue to diverge at the end of the simulation as the IOBL accelerates (Figure 1e,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.

390
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 4/3 is also consistent with our data (R 2 = 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 395 et al., 2018). It's worth noting that our simulations do not capture the large scale increase in overturning circulation that accompanies a distant increase in thermal driving (i.e., changes to the water masses entering the ice-shelf cavity). The relationship between melt rate and distant thermal driving is a function of both the relationship between the local thermal driving and melt rate and the relationship between distant thermal driving and friction velocity (Holland et al., 2008). The former relationship is addressed by this study; the latter is only partially captured by the increase in IOBL velocity due to an increase in buoyancy 400 while the far-field velocity remains constant. We note that the often cited quadratic relationship between thermal driving and melting pertains to this distant thermal driving (Holland et al., 2008); studies generally find this exponent 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 405 on the local thermal driving, based on the temperature in the first model layer (Equation 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 α) n ) combines two effects: the effect of ice-shelf slope 410 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 Section 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 415 profile. Thus, we cannot fully assess the 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 = 3/2, and Little et al. (2009), who found n = 0.94 for the range of slopes considered here. Our finding that n = 1 may be applicable to the shear-dominated regime, in contrast to the sensitivity of n = 2/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 420 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 425 IOBL turbulence. This is not evident in our simulations, which do not show significant differences in TKE budgets between the two runs ( Figure 4). 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.

430
The thermal exchange coefficient as computed using Equation 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 m 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 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 Figure 7c,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  length (Equation 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), and 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 445 flux parameterizations is needed. The second contributing factor is declining simulated TKE, which reduces thermal exchange coefficients over the course of these simulations (Figure 7c).
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 not evident in our simulations. The decline in exchange coefficient with thermal driving agrees with Vreugdenhil and Taylor (2019), though 450 their sensitivity is greater than that seen in this study (their Figure 8). As shown in Figure 7c, this relationship holds during all inertial periods, with an exception during an interval of strong turbulence (see Section 3). Due to the limitations of our LES modeling discussed in Section 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 455 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 (Figure 7d). 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 460 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 ( Figure S4).
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 465 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 Equation 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 (Γ T ∝ u 0.5 * , w θ ∝ u 1.5 * ; Ramudu et al., 2018). More studies are required to determine this scaling.

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 475 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 Figure S6 with his Figure 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 (Figure 3).

480
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 10s 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 485 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 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.

490
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 m 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 495 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 (Figure 3c,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.

500
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 20 https://doi.org/10.5194/tc-2021-242 Preprint. Discussion started: 1 September 2021 c Author(s) 2021. CC BY 4.0 License. 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.

505
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 510 further complicates the parameterization of momentum fluxes at the ice-ocean interface. A depth-dependent shape function for momentum fluxes which may prescribe negative (downward) momentum fluxes at the top grid cell below ice shelves for typical ocean model resolutions (see Figure 3a,d at, e.g., 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.

Conclusions
In this study we presented a small parameter study of IOBL turbulence as captured by an 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 10s m 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 520 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 525 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 530 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.