Ice base slope effects on the turbulent ice shelf-ocean boundary current

The majority of Antarctica’s contribution to sea level rise can be attributed to changes in ocean-driven melting at the base of ice shelves. Turbulent ocean currents and melting are strongest where the ice base is steeply sloped, but few studies have systematically examined this effect. Here we use 3-D, turbulence-permitting large-eddy simulations (LES) of an idealised ice shelf-ocean boundary current to examine how variations in ice base slope influence ocean mixing and ice melting. The range of simulated slope angles is appropriate to the grounding zone of small Antarctic ice shelves and to the flanks of relatively wide ice base channels, with far-field ocean conditions representative of warm-water ice shelf cavities. Within this parameter space, we derive formulations for the friction velocity, thermal forcing, and melt rate in terms of total melt-induced buoyancy input and ice base slope. This theory predicts that melt rate varies like the square root of slope, which is consistent with the LES results and differs from a previously proposed linear trend. With the caveat that further simulations with an expanded range of basal slope angles and ocean conditions would be necessary to evaluate the validity of our conclusions across the full Antarctic ice base slope parameter space, the derived scalings provide a potential framework for incorporating slope-dependence into parameterisations of mixing and melting at the base of ice shelves.


Introduction
Mass loss from the Antarctic Ice Sheet contributes significant uncertainty to global sea level predictions (Fox-Kemper et al. 2021).The present-day retreat of the ice sheet is predominantly attributed to an increase in ocean-driven melting of its floating ice shelves (Rignot et al. 2013).
Accurately capturing ice shelf-ocean interactions in climate models has therefore become an essential prerequisite for improving the reliability of sea level forecasting.
While continental shelf processes and bathymetry impact basal melting by influencing the properties of the water masses present in the ice shelf cavity (e.g., Azaneu et al. 2023;Haigh et al. 2023), local melt rates are ultimately governed by the vertical convergence of heat and salt fluxes at the ice-ocean interface, in turn regulated by the exchange of heat, fresh water, and momentum across the turbulent ice shelf-ocean boundary layer (Holland and Jenkins 1999).Constrained by computational cost, regional and circum-Antarctic ocean models do not explicitly resolve small-scale turbulence within the ice-ocean boundary layer.As a result, these models rely on parameterizations to represent the structure of the boundary flow and the fluxes that it generates at the ice-ocean interface.
In this paper we follow the terminology introduced by Jenkins (2016) which refers to the ice shelfocean boundary flow, comprising a boundary layer adjacent to the ice and a pycnocline beneath, as the ice shelf-ocean boundary current (ISOBC).In contrast to the boundary layer in which Ekman dynamics are relevant, the flow in the pycnocline has been shown to be mostly geostrophic (Jenkins 2021;Patmore et al. 2022).However, the most commonly used ice shelf-ocean boundary layer parameterizations (Holland and Jenkins 1999) were originally developed to evaluate melt under flat (2019) noted an increase in melt rate with slope.For small ice base slope angles up to 0.02 • , Jenkins (2021) suggested a near linear relationship between melt rate and slope accounting for Earth's rotation, although with parameterized turbulence.More recently, a study by Begeman et al. (2022) based on turbulence-permitting large-eddy-simulations (with rotation) also found a linear sensitivity of melt rate to slope angle, with a threshold-like behaviour at very shallow slopes.
However, Begeman et al. (2022) only considered four slope angles ranging between 0.01 • and 1 • and the suggested linear trend was inferred by fitting a curve through the data points.While the numerical modeling studies listed above have each highlighted important aspects of slopeinduced dynamics, an expression for the melt rate as a function of ice base slope, derived based on turbulence-permitting simulations accounting for the effects of Earth's rotation and melt-induced stratification has not yet been established.
In order to examine ice base slope effects, we conduct a series of turbulence-permitting large-eddy simulations (LES) where we vary the slope angle and fix other parameters.We use an idealized rectangular model geometry with a uniformly sloped, planar ice base.Periodic lateral boundary conditions are applied to generate a boundary flow that is solely forced by the resultant meltwater buoyancy, thereby neglecting any external drivers of turbulence present under real world ice shelves like for example tidal currents (Davis and Nicholls 2019;Richter et al. 2022).We initialize each simulation with uniform temperature and salinity and no mean velocity.This allows us to isolate processes associated with the buoyancy-driven ISOBC.As a result of this assumption, the water column beneath the ISOBC is not turbulent.It is also worth noting that we do not consider the diffusive-convective regime observed in quiescent ice shelf cavities (e.g., Middleton et al. 2022).
These intentional omissions are justified by our overall aim to consider a controlled system in which we can examine the underlying physics of buoyancy-forced ISOBC dynamics, rather than attempting realistic simulations of any specific ice shelf cavity environment.
It is also important to mention that, due to the periodic boundary conditions in the ice baseparallel dimensions required to generate a buoyancy-driven flow in a domain that is sufficiently small to allow for a turbulence-permitting grid resolution, our simulations do not account for largescale lateral advection of heat and salt that may be present under real ice shelves (Schmidt et al. 2023).Given the absence of heat and salt restoring sources in our model configuration, this results in a boundary flow that continuously freshens and cools, implying that steady state conditions will 5 Accepted for publication in Journal of Physical Oceanography.DOI 10.1175/JPO-D-23-0256.1.not be reached until all the available ocean heat has been consumed (at which point there will be no further melting).We are therefore studying the transient initial value problem, in line with the approach taken in previous ISOBC modeling studies (e.g., Jenkins 2016Jenkins , 2021;;Vreugdenhil et al. 2022;Patmore et al. 2022), as well as in large-eddy simulations of other stratified turbulent flows such as wakes (e.g., Zhou 2022).The time-evolving behaviour is a useful problem to study in itself because in reality, the development of the boundary current would be halted at the stage when a balance between the lateral convergence of advected heat and the vertical divergence of interfacial and pycnocline heat fluxes was attained.This would then result in a steady state with similar characteristics to the transient solution at that stage (Jenkins 2016).
The layout of this paper is as follows: section 2 outlines the LES model set up; section 3 describes the effects of ice base slope on the boundary current structure and dynamics, section 4 formulates an expression for transient melt rate in terms of slope angle; section 5 summarizes the key findings, discusses potential ocean modeling implications, and suggests further research avenues.

Model design a. Simulation set up
The turbulent shear flow beneath a sloping ice base is modeled using idealized large-eddy simulations performed in the computational fluid dynamics solver Diablo (Taylor 2008).A schematic of the cuboid model domain is shown in Fig. 1.The upper boundary ( = 0) represents the base of the melting ice shelf and the bottom boundary an arbitrary far-field.The slope is introduced by tilting the computational domain by an angle  from the horizontal such that the - plane aligns with the ice shelf base, with  up-slope,  across-slope, and  perpendicular to the ice base.
Hereafter, we refer to the slope values in percentage rather than angle.While this range of slopes is appropriate to the near grounding line region of small ice shelves and the flanks of relatively wide basal channels (Stanton et al. 2013;Dutrieux et al. 2014;Alley et al. 2016), mean ice shelf slopes can be as low as 0.01%, which is less than the lower bound tested here.However, for slopes smaller than 3% we found that turbulence was suppressed to an extent that it could not be sufficiently 6 resolved on our grid given available computational resources, which is why these shallower angles are excluded from our analysis.
In order to simulate a buoyancy-driven flow along the sloped ice base, while keeping the domain size sufficiently small to enable the use of turbulence-permitting grid resolution, we apply periodic boundary conditions along the ice plane-parallel dimensions and impenetrable conditions in the -direction.No external forcing is applied, meaning that the simulated flow is driven solely by melt-induced buoyancy forcing.We impose a no-flux, free slip condition on the bottom boundary.
At the top boundary, we implement a near-wall model, which estimates the stress, heat flux, salt flux and melt rate at the ice-ocean interface based on the velocity, temperature and salinity at the top grid point of the computational domain (details in section 2c).As noted in the introduction, the use of periodic lateral boundary conditions coupled with the no-flux bottom boundary condition means that the simulations are an initial value problem.The key limitation to this approach is that our results are quantitatively dependent on the choice of initial conditions.To reach a long-term equilibrated state in the periodic domain, the alternatives would have been either to include a restoring zone at the base of the domain, or to specify vertical heat and salt fluxes at the bottom boundary.However, while removing the dependence on initial conditions, given that the physical processes leading to steady state conditions under ice shelves are not well understood, both of these The temperature and salinity fields are initialized with uniform values of  = −0.12• C and  = 34.5 psu.Setting the hydrostatic background pressure to  = 300 dbar (equivalent to a depth below sea level of approximately 300m), this gives a thermal driving (defined as the difference between the ocean temperature and the freezing temperature   (, ) of 2 • C, representative of conditions in 'warm' Antarctic ice shelf cavities flooded by Circumpolar Deep Water (Dinniman et al. 2016).The velocity field is initialized with small amplitude random noise on the order of 0.01m s −1 to help initiate turbulence.Each simulation is run until an accumulated melt of 3m has been achieved, plus one inertial period (  ) of 12.5 hrs.Due to the variation of melt rate with slope (discussed later), this results in different simulation durations.For the simulations presented here the durations range from 23 inertial periods for the steepest slope to 120 inertial periods for the shallowest slope.
The model domain is 64m in each ice plane-parallel direction and 400m in the ice plane-normal direction.The flow is discretized on a uniform grid with a uniform resolution of Δ  = Δ  = 2 m in the ice plane-parallel dimensions and Δ  = 0.5 m in the wall normal dimension.This resolution is comparable to the highest resolution used in Vreugdenhil et al. (2022).For the simulations presented in this paper, and with the exception of the near-wall layer (between the uppermost grid point and the ice base), this grid resolution is sufficiently high for at least 80 % of the scalar and momentum fluxes to be generated by resolved turbulence (with the remaining unresolved portion parameterized with a subgrid-scale model), ensuring accuracy of the solution obtained using LES.
The chosen grid configuration also meets the isotropy requirement for accurate LES of a stratified wall-bounded flow (vertical-to-horizontal aspect ratio of at least 0.25) identified by Vreugdenhil and Taylor (2018).

b. Governing equations
Diablo solves a low-pass filtered version of the non-hydrostatic Navier-Stokes equations under the Boussinesq approximation, along with conservation equations for mass, heat and salt, and a linear equation of state.The simulations apply the Traditional Approximation of Rotation, which neglects the local horizontal component of the Earth's angular velocity.Under this approximation, Earth's rotation is included with a slope-independent Coriolis parameter (set here to  = −1.4× 10 −4 rad s −1 based on typical Antarctic latitudes).This simplification is justified given that the correction is small for the range of slope angles considered here.Moreover, maintaining  independent of slope enables us to isolate the effect of slope-driven buoyancy force when interpreting the result, which supports our aim to gain insight into fundamental physical processes rather than reproduce flow conditions under a specific ice shelf.In the absence of external forcing, and for a model domain tilted by an angle  from the horizontal, the grid-filtered equations are the following: Δ where u = (, , ) is the 3-D resolved velocity field with the corresponding position vector (, , ), i, k are the unit vectors in the -direction (up-slope) and -direction (ice plane-normal) respectively,  is time,  is pressure,  and  are the resolved scalar fields, and   and   are the reference temperature and salinity set to initial conditions.The density deficit is Δ =  −   with  the resolved density field and   the reference density.The gravitational acceleration is  = 9.81 m s −2 and the coefficients of thermal expansion and haline contraction are   = 3.28 × 10 −5 • C −1 and   = 7.84 × 10 −4 psu −1 , chosen to be consistent with typical values used in ice shelf-ocean interaction studies (e.g., Patmore et al. 2022).We use realistic values of molecular viscosity ( = 1.8 × 10 −6 m 2 s −1 ), and molecular diffusivities for temperature (  = 1.3 × 10 −7 m 2 s −1 , Pr = /  = 13.8) and salt diffusivity (  = 7.2 × 10 −10 m 2 s −1 , Sc = /  = 2500).
Derivatives are calculated with a pseudospectral method in the -and -directions.Derivatives in the -direction are computed using second-order finite differences.The equations are solved using an adaptive time-stepping scheme with a semi-implicit Crank-Nicholson method applied for the viscous and diffusive terms and an explicit third-order Runge-Kutta method for all other terms.Further details on the numerical method implemented in Diablo can be found in Taylor (2008).
(1), ( 3) & (4) are evaluated using the anisotropic minimum dissipation (AMD) turbulence closure scheme (Abkar et al. 2016;Vreugdenhil and Taylor 2018).The AMD model set up employed here is identical to the set up implemented in previous simulations of the ice shelf-ocean boundary current (Vreugdenhil et al. 2022).

c. Boundary conditions for ice shelf melting
Our vertical model resolution of Δ  =0.5 m is too coarse to resolve the very small turbulent motions near the ice base, where the size of the eddies is limited by distance from the wall.We therefore use a near-wall model recently developed by Vreugdenhil and Taylor (2019) and verified for a stratified shear-driven turbulent flow beneath a melting ice base (Vreugdenhil et al. 2022) to estimate the stress, scalar fluxes, and melt rate at the ice-ocean interface from the velocity, temperature, and salinity at the top grid point of the computational domain (located Δ  /2 = 0.25 m beneath the ice base).The near-wall model uses Monin-Obukhov similarity scaling (law-ofthe-wall logarithmic scaling modified to account for melt-induced stratification) coupled with the diffusive melt equations (McPhee et al. 1987;Holland and Jenkins 1999).Originally developed for evaluating melt under sea-ice, the diffusive melt equations describe the balance of heat and salt fluxes at the ice-ocean interface and constrain the temperature at the ice base to be equal to the local depth-and salinity-dependent freezing point.Neglecting heat conduction into the ice and volume input of meltwater (as justified in Holland and Jenkins (1999)) yields the following system of equations: where  0 and  0 correspond to the interface temperature and salinity, and   and   are the heat and salt fluxes.The specific heat capacity of water is   = 3974 Jkg −1 • C −1 , the latent heat of fusion is   = 3.34 × 10 5 Jkg −1 , and  1 = −5.73× 10 −2 • C,  2 = 8.32 × 10 −2 • C, and  3 = −7.53× 10 −4 • C dbar −1 are empirical constants used to express the seawater freezing point as a function of salinity and depth (Hewitt 2020).The near-wall model computes   and   at each time step and in each grid cell, which results in spatially and temporally varying melt rates at the ice base.Full details on the derivation of the near-wall algorithm can be found in Appendix B of Vreugdenhil et al. (2022).
It is also worth mentioning here that most regional and circum-Antarctic ocean models that incorporate ice shelves apply the same diffusive melt equations as implemented in the near-wall model.However these larger scale models generally parameterize the turbulent fluxes of heat and salt via prescribed (and constant) heat and salt transfer coefficients Γ  and Γ  : where  and  correspond to far-field conditions.In practice, these values are sampled either in the first grid cell beneath the ice or averaged over a set distance away from the ice base (e.g., Losch 2008), which may or may not correspond to far-field conditions depending on model resolution.The friction velocity  * gives a measure of the strength of the interfacial shear stress.In ocean models, the friction velocity is commonly modeled as a function of the horizontal velocity magnitude just outside the boundary layer  and a constant drag coefficient   as  * = √   .This quadratic drag law friction velocity parameterization neglects boundary layer stratification effects which have been shown to impact basal melt rates (Vreugdenhil and Taylor 2019).In contrast, the near-wall model implemented here estimates momentum flux from the law-of-the-wall following a linear stability function of stabilizing buoyancy forcing, meaning that the effect of boundary layer stratification on basal melt rates is accounted for.Furthermore, given that the values of   , Γ  and Γ  are poorly constrained beneath Antarctic ice shelves, the main advantage of estimating the shear stress and scalar fluxes at the ice-ocean interface using the near-wall algorithm is that it does not rely on having to make assumptions about the most appropriate values of   , Γ  and Γ  .

a. Mean boundary current structure
The temporal evolution of the simulated mean boundary current structure is illustrated in Fig. 2.
The first three panels present time-evolving ice base-normal profiles of density deficit Δ, up-slope velocity , and across-slope velocity , sampled at different inertial periods   throughout the 7% slope simulation.The profiles have been averaged in space in the slope-parallel directions (this average is referred to as 'ice plane-average' and indicated by angle brackets throughout the rest of the paper) and averaged in time (denoted by an overbar) over 2 inertial periods.The density deficit profiles (Fig. 2a) highlight that the melt-induced release of fresh water at the ice base acts to cool and freshen the near-ice region of the boundary layer.As the meltwater is mixed downward, the thickness of the buoyant layer increases.As the simulation progresses, the stratification at the top of the profile weakens, and the density structure evolves from a nearly uniformly stratified boundary current to a two-layer structure with a relatively well-mixed boundary layer adjacent to the ice base, overlying a broad pycnocline.In line with Jenkins (2021) and Patmore et al. (2022), we find that the mean velocity structure under a sloped ice base can be described by up-slope flow in an Ekman layer adjacent to the ice 12 base, embedded within a thicker across-slope geostrophic current (Figs.2b-d).The dashed lines in Fig. 2c correspond to estimates of the ice plane-averaged geostrophic velocity ⟨  ⟩, calculated by integrating the ice plane-averaged thermal wind balance with respect to .In the tilted coordinate frame, this gives the following expression for the ice plane-averaged geostrophic velocity as a function of ice plane-averaged density deficit and the sine of the slope angle: Varying the slope has two main effects.Firstly, for a given amount of meltwater accumulated in the water column, the buoyancy force driving the current increases with slope due to the sin  factor in the buoyancy term of the momentum equation (second term on the right-hand side of Eq. ( 1)).Secondly, at any particular time in the simulation, the amount of accumulated meltwater (and hence the thickness of the boundary current) increases with slope due to the positive feedback between buoyancy and melt rate (caused by an increase in buoyancy-driven flow, which increases turbulence and hence enhances melting).This more rapid evolution of the boundary current for steeper slopes is illustrated in Fig. 3a, where the base of the boundary current, ℎ, is defined as the distance from the ice base where ⟨⟩ < 0.01 × max(⟨⟩).It clarifies the analysis if we single out the first effect, so following Patmore et al. (2022), we choose to compare simulations with different slopes at the same accumulated melt (which is equivalent to depth-integrated buoyancy) rather than at the same time.This removes the second effect.
The deepening of the ISOBC with time shown in Fig. 3 is associated with a thickening of both the Ekman layer and the TWL (see portions above and beneath the pink markers in Fig. 2) Accepted for publication in Journal of Physical Oceanography.DOI 10.1175/JPO-D-23-0256.1.
due to turbulent mixing of momentum through the TWL.However, the turbulent momentum flux (Reynolds stress) is small compared to the Coriolis and pressure gradient terms in the up-slope () momentum equation, and hence thermal wind balance remains a very good approximation.
Note that the turbulent momentum flux is a leading order term in the across-slope () momentum equation (not shown).Fig. 4 shows density and velocity profiles for 1m accumulated melt in the top row and 3m accumulated melt in the bottom row.The mean boundary current structure described earlier for the 7% slope simulation (Fig. 2) is a generic feature.Comparing the density deficit profiles between the differently sloped simulations (Figs. 4a,b) shows that for the same amount of buoyancy input, increasing the ice slope causes the boundary current depth to increase and the density deficit at the ice base to decrease.This can be explained by a more vigorous turbulent transport from the far-field towards the ice for steeper slopes.Furthermore, as expected based on the slope dependence of the buoyancy force term in the momentum equations, for the same input of buoyancy from melting, the flow is faster under a steeper ice shelf basal slope (Figs.4c-f).
For all the simulated slope values, the up-slope velocity peaks within a few meters of the ice base.
In comparison, the maximum across-slope velocity occurs further from the ice base.Although the depth of maximum across-slope velocity increases with accumulated melt, it is independent of slope when evaluated at the same value of accumulated melt (see black markers in Figs.4e,f).For a given total input of buoyancy, the near-surface stratification increases with slope (Figs.4a,b).This can be explained using the stability parameter ⟨⟩, defined as the ratio between the depth scale of the Ekman layer ⟨  ⟩ and the Obukhov length scale with   von Kármán's constant and  0 the interfacial buoyancy flux (Rosevear et al. 2022).The Obukhov length scale describes the relative strength of shear and stratification produced by melting at the boundary.When the melt-induced buoyancy flux is strong compared to shear production (i.e. when the boundary layer is strongly affected by stabilizing buoyancy), ⟨  ⟩ ≪ ⟨  ⟩, hence ⟨⟩ ≫ 1.As the simulation progresses, ⟨⟩ decreases (not shown).When ⟨⟩ approaches unity, rotation starts to limit the thickness of the boundary layer, leading to the formation of a more clearly defined mixed layer.For a given amount of accumulated melt, ⟨⟩ increases with slope (see Fig. It is also interesting to note that in the TWL, ice plane-normal gradients in across-slope velocity are nearly independent of slope (Figs.4e,f), and ice plane-normal gradients in both the density and the across-slope velocity appear to be steady in time (Figs.2a,c).These observations are examined in more detail in the next section.

b. Coupling between stratification, shear, and mixing
Melting beneath a sloped ice base generates a buoyancy-forced sheared flow parallel to the ice as well as a buoyancy gradient perpendicular to the ice.Whilst the former drives turbulence, the latter tends to suppress it.The relative strength of the shear and stratification can be quantified using the gradient Richardson number.Here, we calculate the ice plane-averaged gradient Richardson number in the tilted reference frame as: The buoyancy frequency squared ⟨ 2 ⟩ = ⟨⟩/, with buoyancy defined as  = −Δ/  , gives a measure of the strength of stratification normal to the ice base, and ⟨⟩ denotes the mean shear of the ice plane-parallel flow.Note that for the range of slopes considered in this study, the difference between ⟨Ri g ⟩ in the tilted and geophysical (i.e.non-tilted) coordinate frames is negligible (0.5% at most).As described in the previous section, our simulations all feature a pycnocline within which the time-averaged flow is approximately in thermal wind balance, implying a coupling between mean shear and stratification via thermal wind balance.Assuming that in the pycnocline, the shear of the mean ice plane-parallel flow is dominated by its across-slope component (i.e.⟨⟩ ≈ |⟨⟩/|), the frictionless thermal wind balance equation in the tilted coordinate frame can then be expressed as: If ⟨⟩ can be approximated by a slope-independent value, as observed in Figs.4e,f expression for the gradient Richardson number in the thermal wind layer: If ⟨⟩ is assumed to be slope-independent, this suggests that the ISOBC beneath a steeper ice base features a lower ⟨Ri g ⟩, hence a more turbulent pycnocline.The decrease of the mean gradient Richardson number with slope in the thermal wind layer, as expressed by Eq. ( 14), is illustrated in the top row of Fig. 5.Here the panels show scatter plots of instantaneous ⟨ 2 ⟩ versus ⟨⟩ 2 sampled in the thermal wind layer over four inertial periods for four differently sloped simulations.In each panel, the solid grey line corresponds to the ⟨ 2 ⟩ − ⟨⟩ 2 combinations required for thermal wind balance from Eq. ( 13) and the dashed black lines indicate where ⟨Ri g ⟩ = 1/4.Any data points above the dashed line are then characterized by ⟨Ri g ⟩ > 1/4, and vice versa.As expected based on Eq. ( 14), under thermal wind balance, ⟨Ri g ⟩ is higher for shallower slopes, as indicated by the position of the cluster of darker markers relative to the dashed black lines.Assuming that stratified shear-driven turbulence requires ⟨Ri g ⟩ ≈ 1/4 (Rohr et al. 1988;Holt and Ferziger 1992) leads us to expect that turbulence cannot be maintained in a steady state for the shallower slope angles.
The top row in Fig. 5 also highlights that the pycnocline is not in thermal wind balance at all times.Moreover, and as noted in Figs.4e,f, the simulations with slopes steeper than 4% feature a well defined cluster of data points centered around a similar ice plane-averaged shear value (as indicated by the cyan marker) with a spread of shear values not exceeding ⟨⟩ 2 ≈ 2 s −2 in all of these cases.In comparison, the 3% (not shown) and 4% simulations feature a larger spread of values, both in terms of ⟨ 2 ⟩ and ⟨⟩ 2 .
The clustering around a nearly slope-independent mean shear value (at least for steeper slopes), as well as the decrease in mean stratification and gradient Richardson number with slope is further emphasized in Fig. 6, presenting statistics of ⟨Ri g ⟩, ⟨⟩, and ⟨⟩, each sampled within the TWL for the entire duration of each of the differently sloped simulations.The ⟨Ri g ⟩ probability distribution curves each feature a distinct peak, indicating that for a given ice base slope, the flow through the pycnocline can be characterized by a most probable ⟨Ri g ⟩.The peak probability ⟨Ri g ⟩ decreases with slope, with a percentage change of around 80% between the 10% slope and 4% slope simulations (see coloured markers in Fig. 6d).The results therefore differ from the previous characterization of the TWL as having a state of marginal instability associated with a single critical gradient Richardson number that is independent of forcing (e.g Jenkins 2021).
Compared to ⟨Ri g ⟩, the ice-plane averaged flux Richardson number, defined as the ratio of the buoyancy flux ⟨⟩ to the turbulent kinetic energy shear production rate ⟨⟩, with primes denoting fluctuations calculated here with respect to ice plane-averaged values, features a much smaller spread in values across the eight differently sloped simulations (|⟨Ri f ⟩| = 0.25 + /−0.030, see Fig. 7a).The weak slope dependency of ⟨Ri f ⟩ is highlighted in the bottom row of  between ⟨Ri g ⟩ and ⟨Ri f ⟩, demonstrating a new result from the LES that the relationship between bulk property gradients and turbulent fluxes changes as a function of slope.The small remaining reduction in the absolute value of ⟨Ri f ⟩ with slope suggests that a smaller fraction of the TKE shear production is used to mix the density profile as the slope angle increases.This is consistent with a decrease in pycnocline stratification with increasing slope.
An alternative way of viewing these results is that the absolute value of the turbulent Prandtl number ⟨Pr t ⟩ =   /  , which can be shown to be equal to ⟨Ri g ⟩/⟨Ri f ⟩ (see Kundu et  (for the 3% slope) to 0.70 (for the 10% slope), which are relatively close to a |⟨Pr t ⟩| ≈ 1 often found in stratified turbulent flows (Caulfield 2020).The slight decrease in simulated |⟨Pr t ⟩| with slope is consistent with the positive correlation between |⟨Pr t ⟩| and ⟨Ri g ⟩ that emerges from commonly used turbulence closure models that parameterize effective viscosity and diffusivity as a function of the gradient Richardson number (e.g., Pacanowski and Philander 1981).However, it is interesting to note that the ⟨Pr t ⟩ values calculated from the model of Pacanowski and Philander (1981) (using the simulated peak probability ⟨Ri g ⟩) are larger than 1 and at least twice as large as the simulated |⟨Pr t ⟩| values (see Fig. 7c).
Comparing the vertical spread of the coloured markers in Figs.6e & 6f shows that the mean shear displays a weaker sensitivity to slope angle variations than the buoyancy frequency (≈ 40% vs. ≈ 80% increase between the 10% slope and 4% slope simulations), as discussed above.
This difference in slope sensitivity can be explained based on the ⟨Ri f ⟩ values from the LES combined with an empirical model for Pr t .Venayagamoorthy and Stretch (2010) give an empirical relationship between Pr t and Ri g based on laboratory experiments and direct numerical simulations  these estimated gradient Richardson number values, ⟨⟩ and ⟨⟩ can then be inferred from thermal wind balance [Eqs.( 13) & ( 14)].The inferred percentage changes in ⟨⟩ and ⟨⟩ agree very well with the simulated relative ranges (compare circle markers and star markers in Fig. 6d).This demonstrates that the weak slope sensitivity of the mean shear found empirically from the LES may be explained by combining the weakly varying ⟨Ri f ⟩ with the known increase in Pr t with Ri g that is found in stratified turbulence (Venayagamoorthy and Stretch 2010).
The results above show that the dependence of the shear on slope angle is considerably weaker than for stratification.This motivates us to approximate the shear in the TWL with a slopeindependent value, which we refer to as the critical shear ⟨  ⟩.For our model configuration, ⟨  ⟩ is found to be approximately equal to 0.0076 s −1 , as indicated by the bold dashed line in Fig. 6b.When combined with thermal wind balance, a slope-independent critical shear implies that the gradient Richardson number must increase with decreasing slope angle.For sufficiently large values of the gradient Richardson number, we might expect that the mean shear and stratification cannot support a sustained turbulent flow and that the flow becomes intermittently turbulent.This is explored in the next section.

c. Turbulence intermittency
The intensity of turbulence can be quantified by the turbulent kinetic energy (TKE): The time history of ice plane-averaged TKE, along with the evolution of ⟨Ri g ⟩ and ⟨⟩ 2 , are shown in Figs.8a-c for the 7% slope case.
Turbulence is strongest near the ice base in the Ekman layer.Within the thermal wind layer (beneath the pink line), ⟨⟩ decreases with distance from the wall and it goes to zero at the base of the boundary current.Near the base of the boundary current there is a region in which mixing appears to be intermittent.As highlighted in Fig. 8d, this deeper layer features periods with no mixing interrupted by bursts of turbulence occurring on an inertial timescale.The ⟨Ri g ⟩ and ⟨⟩ 2 Hovmöller diagrams in Figs.8b,c reveal the same distinction between a uniform and an intermittent mixing region within the thermal wind layer.As expected based on the variation of ⟨⟩ with depth observed in Fig. 8a, ⟨Ri g ⟩ increases with distance from the ice base, which is driven, at least partly, corresponds to the peak probability value of ⟨Ri g ⟩, as diagnosed from Fig. 6a for the 7% simulation.In (f), the solid lines correspond to the simulated ⟨ 2 ⟩ (blue) and ⟨⟩ 2 (black), and the dotted black curve corresponds to ⟨⟩ 2 predicted from thermal wind balance [Eq.( 13)] based on the simulated ⟨ 2 ⟩.
by a decrease in mean shear with depth (Fig. 8c).The other simulations exhibit similar oscillations in the thermal wind layer (not shown here), with the thickness of the intermittent region taking up a larger proportion of the pycnocline as the slope gets shallower.This is consistent with the scatter plots in the top row of Fig. 5, which highlighted a larger deviation from the mean thermal wind balance state for lower slope angles.Although not shown here, it is worth pointing out that the melt rate curves do not feature inertial oscillations.
The intermittency of turbulence in the TWL can be explained using the buoyancy length scale   =  ′ /, which is an estimate of the vertical scale at which the vertical kinetic energy is equal to the potential energy associated with displaced parcels (Smyth and Moum 2000).As expected based on the reduction of TKE with distance from the ice base and the approximately constant stratification through the pycnocline, the buoyancy scale decreases with distance from the ice base and goes to zero at the base of the TWL.When   is sufficiently small the turbulence will be suppressed by viscosity.A similar argument can be made using the buoyancy Reynolds number (also known as the turbulent activity parameter).When the downward momentum flux through a turbulent region of the TWL encounters a non-turbulent region, the momentum flux convergence will lead to a local acceleration of the flow and an increase in the shear.This will decrease Ri g until a turbulent event is triggered.If the transition to turbulence is sufficiently fast, it will excite inertial oscillations, which appear to couple with the intermittency.

Melt rate relationship to slope
Ice sheet models rely on parameterizations to represent the effect of ice shelf-ocean interactions on basal melting.The dependence of melt on ice base slope is not accurately captured in presently used basal melt parameterizations (Burgard et al. 2022).In this section we propose a formulation for the transient melt rate as a function of slope angle, supported by the results of the LES.While the assumptions made to derive this scaling will need to be validated in other settings, the derivation described in this section provides a framework that could be adjusted for incorporating slope effects in basal melt parameterizations.
Substituting the parameterized heat flux expression [Eq.( 9)] into the heat balance at the ice-ocean interface [Eq. ( 6)] gives the following expression for the domain averaged melt rate applicable to the shear-controlled ISOBC regime modeled here: with ⟨ * ⟩ the ice plane-averaged friction velocity and (⟨ 1 ⟩ − ⟨ 0 ⟩) the temperature difference between the first grid point beneath the ice base and the ice base.In what follows we refer to this temperature difference as 'thermal forcing' following the terminology employed in Holland and Jenkins (1999).Also note that for the remainder of this section, the subscripts 0 and 1 are used to indicate the ice base and the first grid point beneath the ice base, respectively.
As demonstrated by the differently colored curves in Figs.9a & 9b, the simulated friction velocity and thermal forcing increase with slope angle, irrespective of the sampling time.Assuming that the heat transfer coefficient Γ  is independent of slope and constant throughout the simulation (as justified in Appendix A), the increase of the friction velocity and thermal forcing with slope translates into an increase in simulated melt rate with slope (Fig. 9c).In order to quantify this slope-induced melt rate increase, we now derive expressions for ⟨ * ⟩ and (⟨ 1 ⟩ − ⟨ 0 ⟩), which are then be substituted into Eq.( 17) to yield an expression for the transient ice plane-averaged melt rate as a function of basal slope angle.

a. Friction velocity
Ocean models typically estimate stress at the ice-ocean interface from the standard drag law, which relates the friction velocity to the flow speed just outside the boundary layer, where frictional effects become negligible (Holland and Jenkins 1999).This implementation of the drag law implies that the pressure gradient driving the flow is constant with depth throughout the boundary layer.While this assumption holds for a flat ice base, it is not always valid in the case of a stratified boundary layer on a slope.In the latter configuration, the geostrophic velocity decreases approximately linearly with distance from the ice (see dashed lines in Fig. 4e) so using the velocity just outside the boundary layer would underestimate the friction velocity.Instead, in a rotating stratified flow near a sloping boundary, the magnitude of the stress at the boundary scales with the geostrophic velocity at the boundary (Jenkins 2016(Jenkins , 2021)).Assuming ⟨ 1 ⟩ ≈ ⟨ 0 ⟩, the friction velocity can then be approximated by with    a constant, which we call here the geostrophic drag coefficient following Jenkins (2021).It is worth clarifying here that in our LES simulations, rather than being calculated from Eq. ( 18), the friction velocity is estimated from the near-wall model (as described in section 2c).However, Eq.( 18) provides a simple shear stress formula that could be implemented in a melt rate parameterization.Computing the geostrophic velocity at the uppermost grid point from Eq. ( 11) and substituting into Eq.( 18) gives the following general expression for the ice plane-averaged friction velocity as a function of the ice-plane averaged buoyancy at the top grid point ⟨ 1 ⟩ = −⟨Δ 1 ⟩/  : In our simulations ⟨ 1 ⟩ evolves in time (Fig. 2a).In principle, any functional form could be used to approximate the shape of the buoyancy profile and relate the interfacial buoyancy to the depthintegrated buoyancy.Assuming that the buoyancy decreases linearly with depth (as seen for steep slopes in Figs.4a,b), the ice plane-averaged buoyancy can be estimated from ⟨⟩ ≈ ⟨ 2 ⟩( + ℎ), leading to the following approximation for ⟨ 1 ⟩: The boundary current depth ℎ depends on the total melt-induced buoyancy accumulated in the water column since the start of the simulation.Using the linear buoyancy profile from Eq. 20 and assuming that ⟨ 2 ⟩ is steady in time (Fig. 2a), the integrated buoyancy in the ISOBC is: Rearranging for ℎ, combining Eq. ( 12) and Eq. ( 13) to express ⟨ 2 ⟩ as a function of ⟨Ri g ⟩ (which implies thermal wind balance), and substituting into Eq.( 20) gives the following transient 25 expression for the buoyancy at the top grid point: Substituting this expression into Eq.( 19) gives the following expression for the transient ice plane-averaged friction velocity as a function of the gradient Richardson number: According to Eq. ( 23), if we assume that the flow in the ISOBC has a slope-independent ⟨Ri g ⟩, then for a given amount of accumulated melt (i.e. a fixed depth-integrated buoyancy), the ice plane-averaged friction velocity is independent of slope angle (illustrated by the horizontal dotted green line in Fig. 9g).However, as described earlier, our simulations suggest instead that the mean shear can be approximated by a critical slope-independent value ⟨  ⟩, so that Eq. ( 14) implies an inverse relationship between ⟨Ri g ⟩ and sin .Substituting Eq. ( 14) into Eq.( 23), we obtain the following slope-dependent estimate for ⟨ * ⟩: suggesting that, for a given amount of accumulated melt, ⟨ * ⟩ ∝ (sin ) 1/2 , which aligns well with the simulated trend (cf.dashed black line and red markers in Fig. 9g).23) with a slope-independent ⟨Ri g ⟩ = 1/4 (dotted green curves), and from (b) Eq. ( 24) with a slope-independent ⟨  ⟩ = 0.0076 s −1 (dashed black curves, collapsed).For the simulated curves, the different line colors correspond to the differently sloped simulations, with the same color-code as in previous figures.The value of

√︃ 𝐶
in Eq. ( 23) and Eq. ( 24) was selected to optimize the match between the simulated and estimated values in panel (b).
Fig. 10 shows the simulated ⟨ * ⟩ for all slopes (colored curves; same as plotted in Fig. 9d), as well as the friction velocity estimates from Eq. ( 23) in the left hand panel (dotted green lines) and from Eq. ( 24) in the right hand panel (dashed black lines).As highlighted in Fig. 10b, Eq. ( 24) works remarkably well for all slopes, justifying the assumption made, namely linear and steady stratification, thermal wind balance, and a slope-independent shear.In contrast to the curves estimated based on a slope-independent ⟨Ri g ⟩, the estimates based on a slope-independent ⟨⟩ collapse onto a single curve.The collapse of the simulations can be interpreted by noting that the term inside the brackets of Eq. ( 24), sin  × ∫ 0 −ℎ ⟨⟩ , corresponds to the buoyancy force generated by the total meltwater accumulated in the water column since the start of the simulation.In what follows we refer to this term as 'integrated buoyancy forcing'.As shown in Eq. ( 1), when the basal slope changes, the fundamental change that occurs in the system is a change in the buoyancy force per unit density deficit.It is therefore not surprising that the integrated buoyancy forcing is the predictor in the ⟨ * ⟩ relationship and that scaling the ⟨ * ⟩ curves by this term collapses the curves.
Empirically, we find that any small deviations from the collapsed friction velocity (see Fig. 9d) occur when the ice plane-averaged stability parameter  approaches unity (as highlighted by the black markers).As introduced in section 3.a, the evolution from ⟨⟩ ≫ 1 at the start of the simulations towards ⟨⟩ < 1 after a certain amount of meltwater has accumulated, indicates the transition from a stratified boundary layer adjacent to the ice base to a relatively homogeneous Ekman layer in which turbulence is not as strongly affected by the melt-induced buoyancy flux.When ⟨⟩ ≫ 1, we can assume that the ice-plane-averaged ice-normal density gradient is approximately constant with depth and hence that the buoyancy profile can be approximated by a linear function of depth.

b. Thermal forcing
Since we have expressed friction velocity as a function of integrated buoyancy, we now seek an analogous expression for thermal forcing.Buoyancy in the ISOBC is dominated by salinity, so this requires a linkage between thermal forcing and salinity anomalies.We therefore introduce the concept of thermal driving, which can be used to relate temperature, salinity, and density anomalies in the ISOBC (e.g., Holland and Jenkins 1999;Jenkins 2016).Thermal driving is defined as the temperature of a water mass relative to its salinity-dependent freezing point at the pressure experienced at the ice base.This enables us to combine the freezing point relation [Eq.( 8)], the equation of state [Eq.( 5)], and assumptions about the relative mixing of heat and salt to derive an expression for thermal forcing as a function of buoyancy.
We first consider processes close to the ice, in order to relate the thermal driving at the top grid point, ⟨ * 1 ⟩, to thermal forcing across the top grid cell, ⟨ 1 ⟩ − ⟨ 0 ⟩.Following Holland and Jenkins (1999) and Eq. ( 8), the difference in thermal driving between the first grid point and the ice base can be written as: Dividing Eq. ( 9) by Eq. ( 10) and substituting for   and   using Eqs.( 6) & (7) gives: Assuming that molecular diffusion dominates the ratio between salt and heat transfer, the transfer coefficients ratio in Eq. ( 26) can be approximated by Γ  /Γ  ≈ (Pr/Sc) 2/3 (Holland and Jenkins 1999).Combining the resulting expression with Eq. ( 25) and setting ⟨ * 0 ⟩ = 0 (which is true by definition in our model given that the ice shelf basal temperature equals the local freezing point; see section 2c), yields: Next, we consider processes in the wider turbulent boundary layer in order to link thermal driving and buoyancy in the top grid point.Assuming a linear equation of state [Eq.( 5)], uniform far-field properties, a turbulent Lewis number close to unity, and negligible dilution, thermal driving and buoyancy at the top grid point can be related via: with  * ∞ denoting far-field thermal driving (see Appendix B for the derivation of this expression).
Substituting Eq. ( 28) into Eq.( 27) we obtain the following estimate of thermal forcing as a function 29 Accepted for publication in Journal of Physical Oceanography.DOI 10.1175/JPO-D-23-0256.1.
of ⟨ 1 ⟩: Substituting the previously derived interface buoyancy approximation [Eq.( 22)] into Eq.( 29), it follows that the transient thermal forcing is given by with and

c. Melt rate
Having derived estimates for the friction velocity Eq. ( 23) and thermal forcing Eq. ( 30), we may now substitute these into Eq.( 17) to obtain an expression for the transient ice plane-averaged melt rate.This yields: which does not depend on any assumed relationship between ⟨Ri g ⟩ and ice base slope angle.
Assuming a slope-independent thermal wind shear ⟨  ⟩, implying that ⟨Ri g ⟩ decreases linearly with slope as per Eq. ( 14), we obtain the following expression for the melt rate: The key conclusions from this section are that the constant stratification assumption seems to provide a good approximation of the mean ISOBC structure, and that the relationship between melt rate and slope is tied to the assumed slope dependence of the gradient Richardson number.
Assuming an inverse relationship between ⟨Ri g ⟩ and sin  as inferred from thermal wind balance with a slope-independent shear [Eq.( 34)] implies that melt rate scales with the square root of sin with () and () representing time-dependent, but slope-independent, values.It is important to note that the values of () and () are controlled by three empirically derived parameters The sensitivity of these values to varying far-field conditions and external forcing has not been tested as part of this study, but would constitute an interesting area of further research.

Summary and discussion
Improving our understanding of the physical processes governing vertical mixing and oceandriven melting underneath ice shelves is a necessary step toward reducing uncertainty of Antarctic ice shelf melt rates predictions, and hence sea level rise projections.The central objective of the present study is to examine the influence of ice base slope on the interplay between ice shelf-ocean boundary current (ISOBC) structure, shear-driven turbulence, and melting.Given the scarcity of observations required to investigate these interdependent processes, the study is instead performed using turbulence-permitting large-eddy simulations.Our analysis considers slopes ranging between 3% and 10%, which is representative of ice base slopes typically found near the grounding line of small ice shelves and in the presence of wide basal channels.A number of simplifying approximations are made, including a planar ice base, a small angle approximation in the Coriolis force, uniform far-field conditions, and the omission of external sources of turbulence (the simulated ISOBC is forced solely by melt-induced buoyancy), and the results should be interpreted in this context.In line with recent numerical modeling studies investigating the ISOBC under a sloped ice base (Jenkins 2021;Patmore et al. 2022), our simulations feature the development of an Ekman boundary layer overlaying a broad pycnocline, within which the mean flow is approximately in thermal wind balance.While the up-slope momentum balance in the pycnocline is dominated by the Coriolis and buoyancy terms, the base of the ISOBC deepens with time due to the vertical transport of momentum, heat, and salt associated with turbulence in the TWL.Previous studies have proposed that the pycnocline has a constant gradient Richardson number, but we show that the mean flow in the pycnocline is instead characterized by a critical mean shear value that is only weakly dependent upon slope.Combined with thermal wind balance conditions, the assumption of a constant critical shear value ⟨  ⟩ yields an inverse relationship between ⟨Ri g ⟩ and sin .This scaling could be applied to account for slope effects in turbulence closure models that rely on a gradient Richardson number criterion to diagnose the base of the boundary layer (e.g., Jenkins 2021;Burchard et al. 2022).However, it is important to emphasize that the value of ⟨  ⟩ is derived empirically rather than analytically, so it may not hold outside the conditions tested here.Investigating the emergence of a slope-independent critical shear value for a wider range of slope angles and far-field ocean conditions is an important avenue for future research.
Accounting for the slope-dependence of ⟨Ri g ⟩ as per Eq. ( 14) and assuming a linear buoyancy profile (which we find to be a valid assumption for ⟨⟩ ≫ 1), we propose a formulation for the transient melt rate in terms of total buoyancy input and ice shelf basal slope angle.The resulting relationship between melt rate and (sin ) 1/2 agrees remarkably well with the simulations, implying that while basal features with steeper local slopes do enhance melt rates, this enhancement weakens as the slope increases.Moreover, the derived square root relationship suggests that the linear scaling found by previous studies (e.g., Jenkins et al. 2018;Jenkins 2021;Begeman et al. 2022) may not always be applicable, particularly in the steepest parts of ice shelves, which dominate their overall melting.However, the melt rate scaling derived here may be specific to our idealized model set up, and it is worth reiterating that we have not considered the shallow slope angles appropriate to the overall slope of an ice shelf (up to ≈ 2%).Further simulations with an expanded range of ocean conditions and basal slope angles, and with the inclusion of external sources of turbulence such as tides, internal waves, and mesoscale eddies, would be necessary to evaluate the validity of our conclusions across the full Antarctic ice base slope parameter space.Moreover, our proposed melt rate estimate is only valid under the assumption of a uniform boundary current stratification.
A natural extension to this study would be to evaluate the impact that the formation of under-ice mixed layers as ⟨⟩ decreases would have on the derived slope dependency.
As a result of periodicity in the ice base-parallel dimensions, lateral heat advection processes are omitted from our model, which prevents a steady state from being achieved in our simulations.This represents a key limitation of our study, as it means that our findings are quantitatively dependent on the selected initial conditions (chosen here to be representative of ocean conditions under warm-cavity Antarctic ice shelves).An important avenue of further research therefore consists in evaluating the processes responsible for arresting the transient evolution of the ice shelf-ocean boundary current, and characterizing the structure achieved in this steady state.
Finally, while the analysis presented in this paper is idealized, it provides hypotheses that could be tested with in-situ data obtained from potential future targeted observational campaigns.For example, current and density profiles through the entire ISOBC combined with local ice base slope measurements could be used to evaluate the emergence of thermal wind balance conditions in the pycnocline.Ideally, these observations would be obtained both in quiescent and tidally forced ice shelf cavities to test the sensitivity to background currents.Furthermore, if basal melt rate data was collected at the measurement sites, the melt rate estimate provided in the present modeling study could be tested for realistic conditions.

Fig. 1 .
Fig. 1.Schematic of the simulated LES domain (not to scale).

7
Accepted for publication in Journal of Physical Oceanography.DOI 10.1175/JPO-D-23-0256.1.model set up options would have introduced a dependence on an arbitrary choice of boundary conditions.

)
Note that ⟨Δ⟩, and hence ⟨  ⟩, vary as a function of .The pink markers in Fig.2indicate the depth above which the simulated across-slope velocity deviates from the estimated geostrophic velocity (with the deviation diagnosed based on ⟨⟩/⟨  ⟩ > 1.002).In what follows, we refer to this depth as the top of the thermal wind layer (TWL) and we use the terms pycnocline and thermal wind layer interchangeably.Fig. 2d further emphasizes the geostrophic nature of the pycnocline and the Ekman layer set up by the ageostrophic component of the flow above the top of the TWL.Although not shown here, the base of the ageostrophic layer (i.e. the top of the TWL) is relatively well approximated by the boundary layer thickness in an unstratified Ekman layer    ≈ 0.4 ⟨  ⟩ (Pope 2000), where   is the Ekman scale defined as  * /|  |.

Fig. 3 .
Fig. 3. Evolution of the boundary current depth for the differently sloped simulations (color-coded), as a function of time (a) and as a function of accumulated melt (b).The red and yellow markers indicate the point at which an accumulated melt of 1m and 3m has been reached, respectively.The black markers indicate when the stability parameter ⟨⟩ is equal to 1.

Fig. 4 .
Fig. 4. Ice plane-averaged profiles of density deficit (a,b), up-slope velocity (c,d), across-slope velocity (e,f) averaged over 2 inertial periods at an accumulated melt of 1m (top row) and 3m (bottom row).The dashed lines in (e,f) correspond to the estimated ice plane-averaged geostrophic velocity, calculated based on Eq. (11).Panels (g,h) show the difference between ⟨⟩ (solid lines in (e,f)) and ⟨  ⟩ (dashed lines in (e,f)).Pink markers indicate the top of the thermal wind layer (TWL), black markers in panels (e-h) indicate the depth of the maximum ice plane-averaged across-slope velocity.Note the horizontal axis limits varies between panels (c-d) and (e-f).
3b), explaining why steeper slopes feature a more stratified boundary layer in the profiles in Fig. 4. 15 Accepted for publication in Journal of Physical Oceanography.DOI 10.1175/JPO-D-23-0256.1.
, then this explains why the density gradient required for maintaining thermal wind balance decreases approximately linearly with slope, as observed in Figs.4a,b.Furthermore, combining Eqs.(12) & (13) yields an 16 Accepted for publication in Journal of Physical Oceanography.DOI 10.1175/JPO-D-23-0256.1.

Fig. 5 .
Fig. 5. Ice plane-averaged buoyancy frequency squared and shear squared (top row) and buoyancy flux and shear production (bottom row) in the thermal wind layer sampled over 4 inertial periods centered at an accumulated melt of 1m.The colormap indicates the frequency of events (light = low frequency; dark = high frequency).In the top row the solid grey lines correspond to the thermal wind balance relationship calculated as per Eq.(13) using the relevant slope angle, and the dashed black lines correspond to ⟨Ri g ⟩ = 1/4.The cyan markers indicate where ⟨⟩ = 0.0076 s −1 .In the bottom row the dashed black lines correspond to |⟨Ri f ⟩| = 1/4.

Fig. 6 .
Fig. 6.Statistical distribution of ice plane-averaged ⟨Ri g ⟩, ⟨⟩, and ⟨⟩ sampled at each -direction grid point within the thermal wind layer over the entire duration of the differently sloped simulations (a-c).In panels (d-f), the colored markers indicate the most probable values inferred from the probability curves in panels (a-c), with the colors corresponding to the slope values as labeled in (a).The black star markers correspond to estimated values (expressed as percentage change compared to the value of the sin  = 0.1, i.e. slope=10%, case), with estimates obtained based on the Pr t formulation from Venayagamoorthy and Stretch (2010) combined with thermal wind balance conditions (more detail given in the text).Note that the absolute values on the left -axis only apply to the simulated values (grey markers), while the percentage change values on the right -axis apply to both quantities.

(
their Eq.(3.6)).Using their relationship with the peak probability values of ⟨Ri f ⟩ from the LES, the definition ⟨Pr t ⟩ = ⟨Ri g ⟩/⟨Ri f ⟩, and an asymptotic value of the flux Richardson number of 0.32 (  ∞ in Eq. (3.6) of Venayagamoorthy and Stretch (2010)), gives a prediction for ⟨Ri g ⟩.Using

Fig. 7 .
Fig. 7. Statistical distribution of ice plane-averaged |⟨Ri f ⟩| and |⟨Pr t ⟩| sampled at each -direction grid point within the thermal wind layer over the entire duration of the differently sloped simulations (a,b), and |⟨Pr t ⟩| as a function of slope angle (c).The colored markers in (c) indicate the most probable values inferred from the probability curves in panel (b), with the colors corresponding to the slope values as labeled in (a).The black star markers correspond to estimated values obtained from the turbulence closure model of Pacanowski and Philander (1981) combined with the peak probability ⟨Ri g ⟩ values from the LES.

Fig. 8 .
Fig. 8. Temporal evolution of ice plane-averaged TKE (top row), gradient Richardson number (middle row), and shear squared (bottom row) for the 7% sloped simulation over the entire simulation (a-c), and over 3 inertial periods at a single depth (d-f) indicated by the cyan dotted lines in the bottom right corner of (a-c).The pink lines in (a-c) correspond to the top of the TWL.The colormaps in (a-c) have been deliberately capped to make the oscillatory features near the base of the boundary current visible.In panel (b) ⟨Ri g ⟩ is set to zero beyond the base of the boundary current ℎ, to make the boundary current more clearly visible.The dashed line in (e)

Figs
Figs. 8d-f show ice plane-averaged , Ri g ,  2 and  2 over 2 inertial periods sampled at  = −250m (chosen to be within the region of intermittent turbulence, see cyan dotted line in the bottom right corner of Figs.8a-c).Comparing Figs.8d and 8e highlights that turbulence intensity and the gradient Richardson number are inversely related; when ⟨⟩ is highest, ⟨Ri g ⟩ reaches its minimum, and vice versa.The flow oscillates between laminar flow when the ice plane-averaged TKE is near zero and a state of maximum turbulence intensity characterized by the same ⟨Ri g ⟩ as in the turbulent region of the pycnocline.Moreover, the curves in Figs.8d-f highlight the oscillation of the flow about a state of thermal wind balance, as described in the previous section.During laminar periods, buoyancy and momentum accumulate, which leads to an increase in ⟨ 2 ⟩ and ⟨⟩.The resulting mixing weakens the stratification and shear, and since ⟨⟩ 2 decreases faster than ⟨ 2 ⟩, increases ⟨Ri g ⟩ until the flow becomes laminar again.

Fig. 9 .
Fig. 9. Friction velocity, thermal forcing across the top grid cell, and melt rate sensitivity to ice base slope.The first two columns show the evolution of the ice plane-averaged properties for the differently sloped simulations as a function of simulated inertial periods (a-c) and integrated buoyancy force (d-f).The different line colors correspond to the differently sloped simulations, with the same color-code as in previous figures.The red markers indicate the point at which an accumulated melt of 1m has been reached, the black markers indicate when the stability parameter ⟨⟩ = 1.Panels(g-i) show the simulated ice plane-averaged properties as a function of ice base slope angle, sampled at 1m accumulated melt (red markers) compared with the estimated quantities based on a slope-independent gradient Richardson number of ⟨Ri g ⟩ = 1/4 (dotted green lines), and based on an empirically derived slope-independent critical shear value (dashed black lines).

Fig. 10 .
Fig. 10.Simulated (solid colored curves) and estimated ice plane-averaged friction velocity, with the estimates calculated from (a) Eq. (23) with a slope-independent ⟨Ri g ⟩ = 1/4 (dotted green curves), and from (b) Eq. (24) empirically derived constants ⟨  ⟩ = 0.0076 s −1 , √︃    = 0.0220, and Γ  = 0.0112.As indicated by the black and green dashed lines in Fig.9i, the estimate with the slope-independent shear assumption reproduces the simulated results better than the estimate based on a constant gradient Richardson number.