Turbulence in the Ice Shelf–Ocean Boundary Current and Its Sensitivity to Model Resolution

: The ice shelf – ocean boundary current has an important control on heat delivery to the base of an ice shelf. Climate and regional models that include a representation of ice shelf cavities often use a coarse grid, and results have a strong dependence on resolution near the ice shelf – ocean interface. This study models the ice shelf – ocean boundary current with a nonhydrostatic z -level con ﬁ guration at turbulence-permitting resolution (1 m). The z -level model performs well when compared against state-of-the-art large-eddy simulations, showing its capability in representing the correct physics. We show that theoretical results from a one-dimensional model with parameterized turbulence reproduce the z -level model results to a good degree, indicating possible utility as a turbulence closure. The one-dimensional model evolves to a state of marginal instability, and we use the z -level model to demonstrate how this is represented in three dimensions. In-stabilities emerge that regulate the strength of the pycnocline and coexist with persistent Ekman rolls, which are identi ﬁ ed prior to the ﬂ ow becoming intermittently unstable. When resolution of the z -level model is degraded to understand the gridscale dependencies, the degradation is dominated by the established problem of excessive numerical diffusion. We show that at intermediate resolutions (2 – 4 m), the boundary layer structure can be partially recovered by tuning diffusivities. Last, we compare replacing prescribed melting with interactive melting that is dependent on the local ocean conditions. Interactive melting results in a feedback such that the system evolves more slowly, which is exaggerated at lower resolution.


Introduction
The Antarctic Ice Sheet is a significant reservoir of water, holding the equivalent of 57.9 6 0.9 m in sea level rise (Morlighem et al. 2020).Since the early 1990s, West Antarctica has experienced considerable ice sheet mass loss, which has accelerated over time (Paolo et al. 2015;Shepherd et al. 2018).This ice loss is predominately accounted for by oceandriven melting beneath floating ice shelf extensions of the ice sheet (Rignot et al. 2013).Enhanced ablation at the ice shelf-ocean interface leads to reductions in ice sheet buttressing and subsequent ice sheet retreat (Pritchard et al. 2012;Gudmundsson et al. 2019).If we are to address growing concern around ice loss and global sea level rise, it is critical that we understand the processes behind ocean-driven ice shelf melt and use this insight to constrain estimates of melt in ocean models.
Ice shelf melting is dependent on the transfer of heat and salt across the ice shelf-ocean interface.The transfer of heat and salt is, in part, determined by large-scale external factors such as the time-variable delivery of heat from the open ocean (Dutrieux et al. 2014;Jenkins et al. 2018).However, the near-boundary region is also important and often exhibits a buoyancy-driven boundary current that regulates the transfer of scalars across the interface.The dynamics within the boundary current dictate the turbulent heat flux and therefore control the level of heat and salt that reaches the ice.
Despite recent advances in our understanding of the ice shelf-ocean boundary current (e.g., Jenkins 2016;Mondal et al. 2019;Jenkins 2021), much remains unknown about the dynamics of the current.Among the unknowns are the influence of Ekman rolls, which arise from a shear instability in the Ekman layer and have long been observed in the atmosphere (Brown 1980).Although there is yet to be observational evidence of their existence within the ice shelf-ocean boundary current, Ekman rolls were recently identified in simulations applicable to Larsen C ice shelf (Vreugdenhil et al. 2022).
Denotes content that is immediately available upon publication as open access.
The boundary current is suggested to emerge as a wellmixed turbulent layer, separated from the far field by a marginally stable pycnocline (Jenkins 2021).Within the boundary current there is a frictional sublayer where the flow interacts with the solid boundary.Melt rates and turbulent transfer within this ice shelf-ocean boundary layer are commonly represented by the three-equation model (Holland and Jenkins 1999), adapted from theory of the sea ice-ocean boundary layer.This model uses simultaneous equations to describe the conservation of heat and salt across the iceocean interface, and represents phase changes via the liquidus condition: where l 1 , l 2 , and l 3 are constant coefficients; P b is pressure; T and S are temperature and salinity; L is the latent heat of fusion; c w is the specific heat capacity; a b is the ablation rate; subscript i, w, and b represent ice, ambient mixed layer, and ice-ocean boundary, respectively; is the friction velocity; C d is the drag coefficient; and u h w is the horizontal velocity vector in the mixed layer.The three-equation model estimates interfacial fluxes of heat and salt based on temperature and salinity differences across the ice shelf-ocean boundary layer and assumes that T w and S w are sampled from within a well-mixed region outside of the boundary layer.
In ocean models, grid size limits the ability to resolve boundary processes at the spatial scales important for the delivery of heat and salt to the ice shelf-ocean interface.The rate and distribution of melting in models strongly depends on the grid layout and resolution (Losch 2008;Schodlok et al. 2016;Gwyther et al. 2020).Models are typically configured with the three-equation parameterization at the boundary, designed to account for the boundary layer.The ability of the model to capture processes outside of the boundary layer is reliant on a combination of the resolved dynamics and generic subgrid-scale mixing schemes, with no explicit parameterization for the boundary current.In a z-level model, vertical resolution dictates the length scale over which the three equations act and the degree to which the boundary current is resolved.In many cases, constraints on resolution mean that grid size can approach or exceed the length scale of the entire boundary current.As a result, representation of the boundary current is omitted, and it is instead parameterized by the three equations.Although tuning of coefficients can compensate for the deficiencies of this approach, the implementation of the three-equation model to represent the entire boundary current is not appropriate, because this parameterization is formulated strictly for the boundary layer.Until models consistently resolve the boundary current, uncertainty will remain for estimates of heat transport toward the ice shelf-ocean interface.A resolution-independent parameterization of the full boundary current would mitigate these issues.
To accurately parameterize the ice shelf-ocean boundary current, there is a need to understand both the real physics and the degree to which model simulations capture these physics.The task is then to parameterize the difference between the two.Since model representation of physical processes varies with grid resolution, the difference between model and reality is not necessarily uniform across model resolution.We seek to address this by providing a basis on which to develop a resolutionindependent ice shelf-ocean boundary current parameterization.Our investigation is complementary to a recent study that explores resolution dependence of a one-dimensional model, which combines the three-equation model with a resolutiondependent stress parameterization at the ice shelf-ocean boundary (Burchard et al. 2022).Here we focus more directly on representation of the boundary current within threedimensional z-level models.
We begin with intermodel comparisons of different models against a z-level model at turbulence-permitting scales (1-m resolution).The z-level simulations use the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al. 1997a,b;Adcroft et al. 2022).The first comparison is made against a recently developed one-dimensional model (Jenkins 2016(Jenkins , 2021)), where mixing is fully parameterized.We then compare against a three-dimensional large-eddy simulation (LES; Vreugdenhil et al. 2022).Both comparison models use coordinate systems that are rotated parallel to the ice base, removing any effect of stepped topography, and use state-of-the-art mixing schemes designed for the setting.The mixing schemes of the one-dimensional model are derived from previous efforts on the sea ice-ocean boundary layer (McPhee 1994(McPhee , 1999)), while the LES calculates subgrid-scale mixing using the anisotropic minimum dissipation model (AMD; Rozema et al. 2015).None of the three approaches are regarded as representing the truth, but we use the comparison to identify discrepancies and commonalities that arise from each approach.The LES is assumed to give a best estimate of the dynamics, and we use these simulations to investigate the ability of the z-level model to represent the ice shelf-ocean boundary current.The chosen resolution of the three-dimensional models is fine enough to observe Ekman rolls similar to those identified by Vreugdenhil et al. (2022).
Once the fidelity of the z-level model is assessed at high resolution, it is possible to diagnose the degradation of dynamics as horizontal and vertical resolution is reduced from 1 to 16 m, providing insight into the resolution dependency associated with model approximation of subgrid-scale processes.In comparison, previous evaluations of resolution dependence have been on cavity-scale environments, where vertical resolutions of 10-450 m (Losch 2008) and 2-20 m (Gwyther et al. 2020) are used.Turbulent boundary dynamics are not the focus of these studies, and their horizontal resolutions of 6-11 km (Losch 2008) and 2 km (Gwyther et al. 2020) do not permit turbulence.A dedicated turbulence-permitting examination of resolution dependence is yet to be done.We investigate the resolution dependence with two representations of the ice shelf-ocean boundary condition: The three-equation model and a constantflux boundary condition.This allows us to delineate the effects of including meltwater feedbacks in the system.
Section 2 outlines the methods, detailing the z-level model configuration and novel adjustments made to model code that are necessary for this study.Section 3 introduces the LES and one-dimensional configurations.The models are then compared and contrasted against the z-level model, and the underlying physics of the boundary layer are described.Section 4 presents the results of varying model resolution and discusses the grid-dependent dynamics.Section 5 provides the conclusions and discussion on further avenues of exploration.

Modeling
This study investigates the initial-value problem of the adjustment of a sub-ice water column to specified initial conditions, with a sloped ice base conducive to a buoyancy-driven current (e.g., Jenkins 2016Jenkins , 2021)).To configure a turbulencepermitting model, simulations must be high resolution, O (1) m.Modeling entire ice shelf cavities at this resolution is not possible, and thus, we chose small domain sizes and focus on the boundary layer at the base of an ice shelf.Incorporating a buoyant current requires periodic boundary conditions in the along-slope direction with an infinite sloping ice base.In both the one-dimensional and LES models, this is achieved by rotating the coordinate system to align with the ice shelf base (Fig. 1a).For the z-level model, we keep the coordinate system aligned with gravity and instead vertically shift the periodic boundary conditions (Fig. 1b).This is dynamically equivalent to rotating the domain but has the advantage that the ice base and mixing are represented exactly as they would be in a cavity-scale model.In this section, we outline the z-level model along with details of the boundary shift.Particulars of the one-dimensional and LES configurations are given in sections 3a and 3c.

a. Shifted lateral boundary conditions
In the unmodified z-level model, communication between cells at the lateral boundaries occurs at a particular z level.This means, with a sloping ice base, as fluid crosses the lateral boundary, it moves vertically closer/further from the iceocean boundary (Fig. 1b).To overcome this problem, we have adjusted the boundary conditions such that lateral communication between cells is set according to the slope.Figure 1b gives an example of this adjustment.With standard boundary conditions, at z 5 k, the cell bounding the ice on the right-hand side would be connected to the ice cell on the left-hand side.The adjusted boundary conditions instead dictate that this cell is connected to the ocean cell at z 5 k 1 1.As a result, the fluid remains at the same distance from the ice shelf-ocean interface as it crosses the lateral boundary and is representative of an infinite planar slope.This is a Boussinesq (incompressible) model using a linear equation of state.Therefore, the density is only dependent upon temperature and salinity, which are continuous across the boundary.The shift of one cell in Fig. 1b is used for illustration purposes, and this offset can be chosen to be any integer number of grid cells.

b. z-level model
The z-level model, MITgcm, solves the nonhydrostatic Boussinesq Navier-Stokes equations: with conservation of volume where Dr 5 r 2 r 0 , r is density, and r 0 is the reference density; Du/Dt 5 u/t 1 u • =u is the material derivative; t is time; u 5 {u, y, w} is the velocity vector where u, y, and w are the velocity components in the up-slope x, across-slope y, and vertical z directions; f is the Coriolis parameter; k is the unit vector pointing upward; p is the pressure; and n is the kinematic viscosity coefficient.Nonhydrostatic components of the Coriolis term are neglected from (3) due to an f-plane approximation.The scalar equations are defined, where k T and k S are the coefficients of diffusivity.Unless stated otherwise, simulations use a linear equation of state, where a and b are thermal expansion and haline contraction coefficients and T 0 and S 0 are the reference temperature and salinity set to the initial conditions.For reference, parameter values used in the model are listed in Table 1.

c. z-level model configuration
Here we outline all configuration choices relevant to the z-level model.A summary of model design is given in Table 2, and Fig. 2a provides a visualization of the domain.Model configurations are referenced using fields in Table 2.For example, MITgcm_Refer refers to the "Refer" class, whereas MITgcm_ Refer_1m refers to the 1-m configuration within that class.All MITgcm configurations make use of the partial-cell method for representation of topography (see Adcroft et al. 1997).There is a sloped ice shelf topography at the upper boundary with slope angle of 0.01 rad.At the lower boundary, the flow is bounded by a wall with the same slope of 0.01 rad.This lower boundary is considered an arbitrary far field rather than a seabed and, as such, a free-slip condition is applied.The lateral boundary conditions are periodic with the shifted adjustment in the x direction described above.At the upper boundary, we represent processes at the ice shelf-ocean interface with a melt condition and a drag parameterization, which are both described in more detail below.
Resolution is varied across configurations with isotropic grid cells of 1, 2, 4, 8, and 16 m in size, and the time step is adjusted linearly with grid resolution.The boundary shift dictates that the ice slope must cover at least one vertical cell over the domain, and so there is a minimum requirement on domain size for a given slope.To address this, we do not use the same domain size across all simulations.We use a 400 m horizontal 3 100 m vertical domain with 1-, 2-, and 4-m resolutions and a 1600 m horizontal 3 400 m vertical domain with 8-and 16-m resolution.Domain-size dependencies have been tested by comparing results for the small and large domain at 4-m resolution.There was little difference in results between the domain sizes.Comparison of MITgcm with the LES is made with the 400 m horizontal 3 100 m vertical domain at 1-m resolution.To reduce computational cost, comparisons with the one-dimensional model are made with a 100-m cubed domain with 1-m resolution.
The majority of simulations use a second-order upwind flux-limited advection scheme (MITgcm option 77).This scheme avoids issues of overshooting solutions but can be overly diffusive at low resolution.In section 4b, we explore whether the results are sensitive to a less-diffusive advection scheme.Our choice is the second-order centered method (MITgcm option 2).These configurations are labeled "AdvOpt" in Table 2.
Viscosity and diffusivity coefficients are chosen as isotropic and constant in space and time.There are two purposes of the background viscosity; one is to keep numerical stability, and the other is to represent unresolved subgrid-scale processes.For our experiments that use the stable second-order upwind advection scheme (option 77), background diffusivity is present to represent unresolved subgrid-scale processes only.Where the second-order centered method (option 2) is used, a small amount of background diffusivity is required for numerical stability.Our choice of diffusivity coefficients is not constrained by this minimum value.We assume that the unresolved turbulence mixes momentum and scalars equally, maintaining a Prandtl number of one, so diffusivity is equal to viscosity unless otherwise stated.For model stability, viscosity is increased linearly with decreases in resolution.In all simulations, the Lewis number (Le) is one, meaning k T 5 k S , so we simply refer to both as k from here on.
In section 4b, we examine sources of diffusivity in the model via the "Diffu" and "AdvOpt" experiments.Unlike the reference case, Refer, these simulations reduce any influence of the background diffusivity by keeping the diffusivity coefficients at a constant value of 1 3 10 24 m 2 s 21 for all choices of resolution.
To aid development of a turbulent field, the horizontal velocities are initialized with small-magnitude random noise centered on zero with limits 61 3 10 23 m s 21 .By default, the model is initialized with uniform temperature and salinity representative of a warm ice shelf cavity environment in the Amundsen Sea, with where T * is the thermal driving.The initial in situ temperature T w is found by rearranging the thermal driving relation to give T w 5 T * 1 T b 5 T * 1 l 1 S b 1 l 2 1 l 3 P b .The pressure P b is set constant throughout, justified by the small variations in pressure over distances examined and necessary to avoid pressure jumps across the shifted boundaries.The pressure P b is calculated separately from the dynamical pressure p, which varies in space and time.With S b 5 34.5, and the pressure dependence removed by setting P 5 300 dbar, T w 5 20:12 1958C: (8) Melting at the ice-ocean interface is prescribed using either the three-equation model or a constant-flux condition.Models described in Table 2 as 1.Fluxes at the ice-ocean interface in the three-equation simulations are computed according to u h w , T w , and S w are sampled at the center of the ocean portion of the cells bounding the ice.This sampling distance is influenced by the partial cells, which can affect the magnitude of melting.A common correction for this influence supplies a weighted mean of these quantities over a given thickness of ocean beneath the ice (Losch 2008).This study seeks to understand the simplest case, and we therefore do not use this correction method.
In all MITgcm simulations, while a free-slip condition is implemented on the lower boundary to represent an arbitrary far field, the upper momentum boundary condition is The drag coefficient C d accounts for the effects of roughness at a given distance from the wall and is commonly assumed constant in space and time.In a constant-stress boundary layer, C d has a logarithmic profile and varies with the distance at which velocities are sampled dz.The choice of constant C d is based on the assumption that samples are far away from the wall, where variations in dz have little impact on C d .Where samples are sufficiently close to the wall, C d is highly dependent on dz and it is no longer appropriate to assume that C d is constant.In this case, a functional form can be used to account for variation in C d with sampling distance.A simple functional relationship for the drag coefficient is the law-ofthe-wall (LoW) scaling (McPhee 2008; Jenkins 2021): where k K is von Kármán's constant and z 0 is the roughness length.
In z-level model simulations, there are multiple sources of variations in dz.Resolution that is nonuniform with depth, changing resolution between configurations, evolving ice in coupled ocean-ice models, and the implementation of partial cells all influence dz.At low resolution, this is generally not a concern due to the weak sensitivity of C d to variation in dz when sampling distance is large.As resolution increases, the assumption of a constant C d has the potential to bias the results.
The choice of C d for each configuration in this study is noted by the wall stress column in Table 2.The default drag coefficient is set to a constant value of C d 5 3 3 10 23 , which is the choice for the majority of the z-level simulations.A selection of configurations deviates from this choice and use the law-of-the-wall scaling instead of a constant value (see LoW in Table 2).This is implemented for consistency with the onedimensional model.The scaling is applied with a roughness length of z 0 5 3.3 3 10 24 m, chosen to match Jenkins (2021).The drag coefficient C d is always consistent across ( 11), ( 12), and ( 13).
The turbulent transfer coefficients G T and G S are associated with the flux of scalar properties across the boundary layer.The transfer is dominated by molecular processes within the viscous sublayer [O(1) mm from the boundary], where turbulence is suppressed and scalar transfer is governed by molecular diffusivity, meaning G T Þ G S .Turbulent transfer coefficients are set constant on the basis that dominant molecular diffusion renders fluxes that have a quasi-linear relationship with u * (Holland and Jenkins 1999;Jenkins et al. 2010).We use constant values of G T 5 1.1 3 10 22 and G S 5 3.1 3 10 24 , following Jenkins et al. (2010).
Figure 2b highlights the transient progression of the density field as meltwater is fluxed into the domain at the upper boundary, and the inset shows the stepped boundaries that exist within a z-level model.In general, the use of partial cells creates a gradual transition in topography and minimizes step size to the extent that most steps are not visible in the inset (Fig. 2b) at this horizontal resolution (see Adcroft et al. 1997).However, the partial-cell formulation requires a minimum cell thickness in order to avoid numerical instability.
Where the layer would be thinner than the minimum thickness, the cell is either collapsed to zero thickness or expanded to the minimum thickness, whichever is closer (Adcroft et al. 1997).Larger topographic steps appear where this occurs (without the use of partial cells, these steps would be significantly larger).The size of the step can be adjusted by changing minimum cell thickness.A lower minimum creates a smoother slope at the cost of requiring a shorter time step to maintain stability.The topographic steps appear to generate waves in the isopycnals that propagate away from the boundary.We examined the sensitivity of boundary current to the prescribed minimum cell thickness at 1-m resolution (not shown).The time evolution of r, u, and y is similar between these simulations, suggesting the waves are of minor influence on the dynamics.

Model evaluation and dynamical stability a. The 1D model
We first compare MITgcm against the one-dimensional models presented in Jenkins (2016Jenkins ( , 2021)).A first comparison is made against the laminar solution from Jenkins (2016) with constant coefficients of viscosity and diffusivity.A second comparison is made with the turbulent solution from Jenkins (2021).Coefficients of diffusivity and viscosity in this latter model are based on a functional relationship that incorporates the local turbulence closure (LTC) of McPhee (1994) and the gradient Richardson number mixing scheme (PP) of Pacanowski and Philander (1981).This combined mixing scheme is termed the hybrid turbulence closure (HTC) and is developed specifically for the ice-ocean boundary current.Aside from an adjusted taper for the HTC, which we will describe below, the onedimensional results are identical to Jenkins (2016Jenkins ( , 2021)).
The one-dimensional model is cast in terms of a single scalar, the thermal driving T * .The equation of state is written as a function of the thermal driving: where a * 5 2:5 3 10 24 8C 21 .The initial conditions are set such that thermal driving is initially 28C everywhere.Merging of the LTC scheme and PP to form the HTC is based on a Richardson number criterion.We define the single-scalar Richardson number as where u is the ice-slope angle.Given bounds a and b, a , b, LTC is applied where Ri * # a, PP is applied where Ri * .b, and the scheme linearly tapers LTC and PP between a and b.The second model comparison uses the one-dimensional model with the taper range used in Jenkins (2021) of 0.25-1.A third comparison is made where this range is adjusted to 0.125-0.25,which will be referred to as the adjusted taper.
The adjusted taper was chosen based on insight from previous literature by Miles (1961) and Maslowe and Thompson (1971), and no other choices were tested.Taper choice will be discussed in more detail in section 3b.
In the laminar one-dimensional model, the thermal driving is governed by a Dirichlet boundary condition at the iceocean interface that is held at the initial value of T * 5 08C.In the turbulent one-dimensional model, a Neumann condition is applied that takes the form of the three-equation model, framed in terms of the thermal driving, where G T * 5 0:006.Stress at the boundary is modeled differently in each one-dimensional case: the turbulent case parameterizes stress as ( 13) with the law-of-the-wall scaling ( 14) applied; the laminar case uses the no-slip condition.Two MITgcm cases are configured to match the laminar and turbulent one-dimensional configurations, 1DCST and 1DSGS (see Table 2).Scalars are modeled in terms of thermal driving by defining the equation of state according to (15).Although we match model configurations as closely as possible, we do not use a Dirichlet boundary condition within MITgcm.Initial conditions of T * 5 28C impose a large temperature gradient at the iceocean interface, leaving MITgcm susceptible to stability issues at initialization when using a Dirichlet boundary condition.To account for this, when comparing laminar solutions, MITgcm applies a Neumann condition that restores to the boundary value, where c 5 0.006 m s 21 .Since we define a positive heat flux as upward, this acts to restore the thermal driving to zero.In essence, this performs the same function as ( 11) and ( 12).
However, c is chosen to both minimize numerical instability and maximize the flux to the extent that it is initially equivalent to a melt rate of ;5000 m day 21 .This enormous value dictates that the flux at the interface dominates over any flux from below, limiting the temperature difference between the ice and the ocean.The reduction in temperature difference at the boundary maintains T * 5 08C in the top grid cell, aligning the boundary conditions of the two models despite the different formulations.When comparing against the turbulent onedimensional model, we apply the boundary condition given by ( 17).Both comparative MITgcm simulations use the drag law with the law of the wall, creating further configuration differences in the laminar comparison.These differences are summarized in Table 2.It will be shown in what follows that the configuration differences between MITgcm and the onedimensional laminar case have little influence over the results.

b. MITgcm versus 1D model
We now compare MITgcm against the one-dimensional configurations described above.This comparison is made in Fig. 3.The columns show the density anomaly, the up-slope x velocities, and the across-slope y velocities.The rows represent the three different one-dimensional model configurations: the laminar case from Jenkins (2016), the turbulent case of Jenkins (2021), and a new configuration of the turbulent case where the HTC scheme in the one-dimensional model uses the adjusted taper.For MITgcm, the profiles are calculated as the mean in the ice-parallel direction, which we term as the ice-plane average.The calculation is made by taking a coordinate transform that re-references the grid in terms of distance from ice base, then averages across the horizontal dimensions of the transformed coordinate system.Each profile represents a snapshot in time, where the number of days and hours elapsed since model initialization are noted by the letters d and h.
The black lines in Figs.3a-c show profiles for MITgcm where constant background diffusivity and viscosity coefficients are chosen to match the one-dimensional model, with n 5 5 3 10 23 m 2 s 21 and k 5 5 3 10 24 m 2 s 21 .These values are high enough to suppress any turbulence in the MITgcm simulations, and a purely laminar solution emerges.In this laminar configuration, there is very good agreement between MITgcm and the one-dimensional simulation, with little difference between the curves for all quantities shown.The small deviations in density near the upper boundary likely stem from restoring to T * 5 0 with a Neumann boundary condition [(18)].Density is initially uniform with a step change at the upper boundary.This diffuses over time and the density deficit generates an across-slope y current via geostrophy.In Fig. 3c, the black dashed lines represent the expected flow for MITgcm based on density gradients if the flow was purely geostrophic.This is calculated in line with Jenkins (2016), where y g 5 2(Dr/r 0 )g sinu/f, with Dr defined by ( 6), and u is the slope angle.The deviations from these geostrophic curves are due to viscous stress.The combination of boundary stress and planetary rotation leads to the production of an Ekman boundary current, and up-slope u velocities develop as a result.The Ekman dynamics result in u being a function of the velocity deficit between y and the dashed geostrophic curves (Jenkins 2016).
For the turbulent solution, we compare the HTC background mixing scheme of the one-dimensional model against a turbulence-permitting configuration of MITgcm.Turbulence is introduced in MITgcm by reducing the diffusivity and viscosity coefficients to the background values of the HTC scheme (Jenkins 2021): n 5 1 3 10 24 m 2 s 21 and k 5 1 3 10 25 m 2 s 21 .These values are chosen to match that of Jenkins (2021) and are small enough to produce a turbulent solution.The transition to turbulence can be identified by the convex/concave geometry of the density anomaly.In the laminar solution (Figs. 3a-c), all curves are concave.In the turbulent solution (Figs. 3d-i), curves begin concave and transition to convex as turbulence begins to enhance mixing.Figures 3d-f show the one-dimensional configuration presented in Jenkins (2021).The profiles are qualitatively similar in comparison to MITgcm, with the main difference being the presence of an enlarged well-mixed layer near the surface in the one-dimensional model.The bold lines in Figs.3d and 3f highlight the region of transition between the two turbulence closures, where 0:25 , Ri * # 1. Below this, the pycnocline in the one-dimensional model is marginally stable with Ri * ' 1 throughout, and the emergence of this condition was used in Jenkins (2021) to provide a final stage of the three-phase description of the flow.The three phases begin with an initial laminar phase where velocities are weak and stratification is strong (d0h12 and d1h13).As the simulation progresses, the flow enters the second phase as a turbulent region appears between the current maximum and the upper boundary (d4h2).During the third phase, a constant gradient pycnocline emerges.For MITgcm, the pycnocline is slightly steeper and 0:25 , Ri * # 1 throughout.Thus, MITgcm exhibits the same dynamical phases, but the third phase emerges with a lower Richardson number, yielding a slightly sharper pycnocline and shallower mixed layer than in the one-dimensional model.
In the one-dimensional model, mixing above the current maximum is generally controlled by the LTC scheme, which is applied using the Richardson number criteria outlined in section 3a.This suggests a possible sensitivity of the HTC to the choice of taper between LTC and PP.Jenkins (2021) tested the effect of changing the size of the taper range, but a shift in the range has not been examined.The default range is 0:25 , Ri * # 1.This choice tends to generate a marginally stable pycnocline in the third phase of dynamics that maintains Ri * 5 1.It could be argued that a more appropriate limit of marginal stability is Ri * 5 0:25 (Miles 1961;Maslowe and Thompson 1971).We show in Figs.3g-i the result of meeting this lower limit by shifting the taper of the HTC to 0:125 , Ri * # 0:25.This adjustment to the HTC scheme leads to a remarkable agreement between the one-dimensional model and MITgcm.The pycnocline in the one-dimensional case is now maintained in a state of marginal stability characterized by Ri * ' 0:25.The upper limit of Ri * for the bold segments in Figs.3g and 3i are chosen to reflect the depth of the pycnocline.While this is imposed by the taper for the one-dimensional model, for MITgcm, it is set by a combination of parameterized diffusive fluxes and resolved mixed-layer instabilities.MITgcm has a slightly stronger stratification above the velocity maximum, and the pycnocline is maintained by a larger range of Ri * , up to 0.6.The presence of Ri * values between 0.25 and 0.6 above the pycnocline in MITgcm can be attributed to the occurrence of resolved instabilities, which do not adhere to a strict limit of Ri * and can occur at a range of Ri * above the neutral stability Richardson number (Kaylor and Faller 1972).

c. LES
We compare MITgcm against two simplified configurations of the LES presented in Vreugdenhil andTaylor (2018, 2019).The LES performs well at representing stratified shear-driven turbulence (Vreugdenhil and Taylor 2018), and we treat this model as our best-available source of information about how turbulence influences ice shelf-ocean boundary layer evolution.The comparisons are made against the 1-m configuration of the MITgcm_Refer class (see Table 2).The term "two scalar" will be used to distinguish this configuration from the single-scalar results presented above.We use the comparison to assess the performance of MITgcm at this fine 1-m spatial resolution.The purpose is also to identify the extent to which MITgcm's partial-cell step topography and constant coefficients of diffusivity and viscosity influence the results.
Similar to the one-dimensional model, the LES uses iceparallel coordinates and does not have stepped topography as is required with z levels.The LES also has the option of a subgrid-scale (SGS) mixing scheme that varies in time and space, the AMD (Rozema et al. 2015).To isolate the effects of step topography, the first LES configuration (LES_Const) is designed such that the coordinate system and therefore the lack of stepped topography is the only major difference with MITgcm.This is achieved by replacing the SGS mixing scheme with constant viscosity and diffusivity coefficients that match those of MITgcm, where k 5 n 5 1 3 10 24 m 2 s 21 .The second LES configuration (LES_SGS) differs from MITgcm in both coordinates and subgrid-scale mixing by employing the AMD.
Aside from the SGS mixing scheme and rotated coordinates, the MITgcm and LES are designed with similar configurations and are aligned with the defaults discussed in section 2c.The grid resolution of both models is 1 m, the domain size is 400 m in the horizontal, and the water-column depth is 100 m.The simulations are initialized with a uniform temperature and salinity associated with a 28C thermal driving and apply a constant heat/salt flux at the upper boundary.Recent iterations of the LES include a Monin-Obukhov near-wall model (Vreugdenhil et al. 2022).To minimize differences between configurations, this near-wall model is not used here, and stress is parameterized at the ice-ocean interface using the drag law [e.g., (13)], with C d 5 3 3 10 23 .

d. MITgcm versus LES
Figure 4 shows a comparison of MITgcm against the LES for a transient solution.Figures 4a-c show profiles of iceplane mean quantities at different times provided as snapshots: Dr, up-slope velocity u, and across-slope velocity y.The black lines are for MITgcm, blue lines for the LES applied with constant background viscosity and diffusivity coefficients, and orange lines for the LES with the subgrid-scale mixing scheme.LES profiles are calculated consistently with MITgcm by taking ice-plane averages of the three-dimensional fields.It is possible that statistical variability remains present in these ice-plane-averaged profiles, which has potential to bias the results.The magnitude of statistical variability was tested by averaging the LES profiles over hourly bins.We found negligible difference in the profiles when calculating these additional temporal averages (not shown), demonstrating that statistical variability is sufficiently removed by ice-plane averaging.The constant heat and salt flux at the upper boundary causes a reduction in the density anomaly in the top cell over time.The heat and salt then slowly diffuse over time, and the signal of melt is propagated to depth.The dashed lines in Fig. 4c are the geostrophic curves as described in section 3b and highlight the role of Ekman dynamics.
The results show that MITgcm performs well at 1-m resolution.Nonetheless, subtle differences exist that are associated with the stress and buoyancy fluxes produced by both models.Figures 4d-f show the parameterized stress |kdu h /dz|, resolved turbulent stress |w u h |, and total stress, where u h 5 (u, y).Figures 4g-i show the parameterized buoyancy flux kdb/dz, resolved buoyancy flux w b , and total buoyancy flux, where b 5 gDr/r 0 .Primes are calculated as anomalies from the iceplane mean, and the total in each case is the sum of the parameterized and resolved contributions.
The total stress (Fig. 4f) and buoyancy flux (Fig. 4g) are comparable between models, but there are small differences near the upper boundary, with the LES displaying larger magnitudes for each.When averaged over the top 3 m and measured at the last time step shown (d6h16), the total stress and buoyancy flux is 109% and 25% larger for the LES.This larger buoyancy flux in the top 3 m is associated with a more defined mixed layer with lower density gradients near the boundary.The enhanced near-boundary stress is associated with reduced velocities near the boundary and larger velocity gradients above the velocity maximum.Splitting the total contributions shows that each model has a different proportion of resolved and parameterized components.In general, MITgcm has a larger resolved and smaller parameterized contribution compared to the LES.This difference is most pronounced in the parameterized stress (Fig. 4d) over the top 3 m, with the LES exhibiting 928% more parameterized stress than MITgcm on the final time step.The large percentage differences in the stress and buoyancy fluxes lead to small near-boundary dissimilarity in the density and velocity profiles of the models.
Primarily, there are three differences in the models: the underlying numerical methods, the grid design, and the turbulence closure.For the grid, MITgcm uses geopotential coordinates with partial-cell representation for topography, whereas the two LES models have rotated coordinates with "smooth" topography.For mixing coefficients, MITgcm_ Refer_1m and LES_Const are identical with constant coefficients of viscosity and diffusivity, but LES_SGS uses a subgrid-scale mixing scheme with coefficients that vary both spatially and temporally.Profiles are similar for the LES_Const and LES_ SGS simulations, suggesting that the choice of mixing scheme is not important for the differences seen between MITgcm and the LES.The remaining differences between LES_Const and MITgcm_Refer_1m are the numerical methods and the rotated grid, indicating that discrepancies near the boundary arise from one or both of these components.
The clearly defined mixed layer in the LES suggests the presence of a greater level of scalar mixing than in MITgcm.We hypothesize that either spurious numerical mixing or the topographic steps in MITgcm affect the relative mixing of momentum and scalars, prohibiting the development of a well-defined mixed layer.This leads to a larger density anomaly in MITgcm, which drives faster flow near the boundary and therefore a reduced shear below this and a weaker parameterized stress.Future work is required to validate this hypothesis.

e. Marginal stability and Ekman rolls
Section 3d presents the first 6 days of the MITgcm simulation with constant heat and salt fluxes (MITgcm_Refer_1m).
This configuration is integrated for 16 days in total, and the extended period reveals important details about the evolution of the boundary current in a three-dimensional setting.Here we examine this extended 16-day period of the two-scalar result presented above, which models temperature and salinity independently, rather than the single-scalar thermal driving.
Figure 5a investigates along-ice-averaged turbulent kinetic energy (TKE) as a function of depth and time for the full time series of the MITgcm_Refer_1m configuration.The simulation exhibits the three-phase development outlined in Jenkins (2021).Phases 1 and 2 evolve over the initial 9-day period, where there is a monotonic deepening and strengthening of TKE.As the boundary layer grows, the turbulent region spans a greater depth, and a subtle lower magnitude of TKE appears near the surface, beneath the frictional layer.After 9 days, there is a regime shift whereby enhanced levels of TKE appear at depth in a periodic manner.This dynamical change is representative of a progression into the third phase of the dynamics, which is associated with a series of instabilities.
Within the boundary current, we identify features of alternating vertical velocities suggestive of Ekman rolls.Despite the term rolls, these are wave-like features, since isopycnals do not overturn.Ekman rolls are triggered by a linear shear instability associated with an Ekman layer (Faller 1963;Faller and Kaylor 1966).The Ekman rolls develop at an angle to the geostrophic flow due to Ekman turning (Faller 1963).Under stably stratified conditions, the Richardson number is the key metric for the instability (Kaylor and Faller 1972;Brown 1972).Ekman rolls in the ice shelf-ocean boundary layer were first identified in Vreugdenhil et al. (2022) (Coleman et al. 1990;Dubos et al. 2008).This may suggest the instabilities seen in phase 3 are a signature of Ekman roll destabilization.However, the Ekman rolls show no signs of destabilization, persisting throughout the third phase.Therefore, an alternative mechanism is anticipated for the periodic peaks in TKE during phase 3.
The instabilities in the third phase of the flow represent a three-dimensional analog to the region of marginal stability that develops in the one-dimensional model, regulating the strength of the pycnocline.Depth profiles in Figs.5e-j single out the first event after day 9.In the run-up to the first peak in TKE, the boundary current cools and freshens through mixing of the meltwater, generating a buoyant current, which progresses in a uniform manner.Vertical density gradients in the pycnocline sharpen, enhancing the geostrophic component of the across-slope flow y.In response, the flow becomes unstable, producing a higher rate of mixing at the pycnocline depth, with elevated stress and buoyancy fluxes.Mixing then reduces the density gradients and velocity shear, returning the turbulent stresses and buoyancy fluxes to their previous levels.
Results in Fig. 6 show how stability evolves over each of the three phases of the dynamics in MITgcm_Refer_1m.To reflect the broader range of Ri that is unstable for MITgcm (Figs. 3g-i), for the two scalar results, we define the region of marginal stability as the depths at which Ri 2 (0.125, 0.6], with Ri 5 where N 2 5 2g(dDr/dz)/r 0 is the Brunt-Väisälä frequency.In this configuration, stratification has a dual role.At the same time as stabilizing the flow, it is also responsible for driving the across-slope geostrophic current.Therefore, the region of high TKE (Fig. 5a) is associated with relatively stable conditions, where Ri . .0.6.The MITgcm result shown in Fig. 6 exhibits two regions of marginal stability, bracketed by the contours.In general, the lower region of marginal stability tracks the bottom of the pycnocline, and below this depth, where flow is negligible, Ri is set to zero.The dynamics of the one-dimensional model are characterized by the uppermost region of marginal stability (Jenkins 2021).Similar to the single-scalar MITgcm case (section 3b), the upper region of marginal stability in Fig. 6 suggests the first two phases of the dynamics broadly follow the one-dimensional model.In phase 1, a stable region (Ri .0.6) appears as the pycnocline develops, with corresponding strong stratification (Fig. 6b) and weak velocity shear (Fig. 6c).The second phase emerges after day 3, where a region of marginal stability develops over the top 5 m, associated with turbulent mixing above the current maximum that dramatically decreases the stratification.
The third phase produced by the one-dimensional model is characterized by the region of marginal stability extending beyond the current maximum.In Fig. 6, the region of marginal stability periodically deepens and shallows with each instability and only ever briefly coincides with the current maximum (not shown).Although on average the current maximum remains below the region of marginal stability, both Fig. 6b and the magenta, purple, and black contours in Fig. 5e show that the oscillatory behavior limits stratification in the pycnocline.This modulation of the pycnocline stratification is similar to the one-dimensional case, where the mixing scheme taper imposes a region of marginal stability that moderates the stratification.Despite the differences induced by three-dimensional turbulence with MITgcm in phase 3, the broader dynamical feedbacks of marginal stability parameterized by the one-dimensional model remain present.

Modeling implications
a. Sensitivity to resolution with a constant interfacial heat/salt flux We have shown that MITgcm performs well against other state-of-the-art models at high resolution.To improve model representation of the ice shelf-ocean boundary current at more typical ocean-model resolutions, we must understand how the results change when resolution is decreased.
Our first step investigates resolution dependence with a constant surface heat and salt flux, as used in sections 3c-e.Coefficients of viscosity and diffusivity and time step are increased linearly as the grid is coarsened.An increase in horizontal viscosity coefficients is required for numerical stability, and the diffusivities are kept equal to the viscosity under the assumption that the unresolved turbulence mixes momentum and scalars equally.Vertical and horizontal motion becomes increasingly isotropic with coarsening grid resolution, but horizontal motion will still dominate at the spatial scales modeled here.This means our choice of isotropic viscosity and diffusivity results in larger-than-necessary coefficients in the vertical.The impact of these choices is discussed in section 4b.Constant-flux simulations allow for continuity with results shown in sections 3c-e and are beneficial for isolating the effects of melt feedbacks associated with an interactive melt condition by holding melt rates fixed across resolutions.A comparison with an interactive melt condition will be presented in section 4c.
Figure 7 shows Dr, u, and y plotted at a fixed time for five different spatial resolutions, where the top of each profile is the middepth of the uppermost grid cell.There is a clear resolution dependence in the density and velocity profiles.The 1-, 2-, and 4-m cases are qualitatively similar, but the boundary layer structure is lost when resolution is coarser than 4 m.The low-resolution results appear much more diffuse in comparison to high resolution, with enhanced scalar mixing throughout the water column.Toward higher resolution, the density anomaly is higher near the surface and lower at depth, and there is a general increase in vertical density gradients.Since the across-slope y velocities are primarily balanced by stratification through geostrophy, the story is similar for the velocities.Lower vertical density gradients are associated with lower vertical gradients in y as resolution is reduced.
Stress in the boundary cells behaves differently for each model resolution.The boundary stress is based on the velocity magnitude within the top cell, and it influences the interior of the domain through the presence of turbulence and viscosity.In the low-resolution results, the current maximum occurs at the uppermost grid point, and near-boundary stress is entirely modeled within this top cell.As resolution increases, the near-boundary stress begins to be resolved over multiple grid points, and a current maximum appears away from the boundary cell.
The varying degrees of resolved stress are associated with different profiles of up-slope u velocities.Ekman layer dynamics are the sole driving force for up-slope velocities in a periodic sloping ice shelf-ocean boundary current in the absence of tides, because there is no large-scale pressure gradient in the across-slope direction.Up-slope u and across-slope y Ekman velocities are defined as (Jenkins 2021): These velocities also correspond to the deviations from the geostrophic flow.Stress above the peak maximum in acrossslope y velocities is associated with up-slope Ekman velocities u .Where resolution is high and dy 2 /dz 2 is negative, u is positive.However, dy 2 /dz 2 reduces with decreasing resolution, the Ekman layer is less resolved, and flow is increasingly geostrophic.In this case, stress at the boundary is dominated by parameterization within the boundary cell.Since the stress parameterization does not incorporate the Coriolis force, the up-slope Ekman velocities disappear.This deficiency was addressed by the parameterization of Wilchinsky et al. (2007) and results here suggest support for the use of such developments when modeling the ice shelf-ocean boundary current.

b. Effect of viscosity and diffusivity
We have shown there is a significant degradation of results with decreasing resolution.The solutions become overly diffuse, density gradients are greatly reduced, and this feeds into the velocity profiles.There are three sources of diffusion within the model: (i) parameterized mixing, (ii) resolved turbulent mixing, and (iii) numerical diffusion (truncation errors).The magnitude of parameterized mixing is, in part, governed by the diffusivity coefficient.Resolved turbulent mixing will scale with turbulence and thus increase with resolution.Numerical diffusion is associated with the advection scheme and generates spurious mixing that reduces with resolution (Griffies et al. 2000).In this subsection, we investigate the relative importance of each component.
Results presented in Fig. 7 used diffusivity coefficients that increased with resolution, which is required by numerical stability if we seek to keep diffusivities equal to viscosities.As a result, these simulations do not separate the influence of numerical and parameterized diffusion.We isolate the role of numerical diffusion by lowering the resolution while keeping the diffusivity coefficient constant.Figure 8 compares the same profiles of r, u, and y presented in Fig. 7 against configurations where the 1-m resolution reference diffusivity of 1 3 10 24 m 2 s 21 is applied.The viscosity and time step varies linearly in the same manner as the reference cases.All simulations have a closer agreement with the 1-m case than their higher diffusivity coefficient equivalent, with particular alignment at 2-(Figs.8a-c) and 4-m (Figs.8d-f) resolution.Despite the improvements brought about by holding k constant at high resolution, the 8-(Figs.8g-i) and 16-m (Figs.8j-l) solutions remain overly diffuse, suggesting a significant influence of spurious numerical diffusion.
An additional consideration for the observed resolution sensitivity, which we do not fully investigate in this study, is the vertical component of the viscosity and diffusivity coefficients.The choice of parameterizing diffusivity and viscosity as isotropic results in larger-than-necessary vertical coefficients due to anisotropy between horizontal and vertical motions in stratified flows.Reducing vertical viscosity or diffusivity coefficients could reduce resolution sensitivity of the profiles without compromising the numerical stability.Our results provide some indication of the sensitivity to vertical diffusivity coefficients.There is low sensitivity to an order of magnitude change in diffusivity coefficients for the 16-m case (Figs.8j-l), which suggests vertical diffusivity is not a limiting factor in recovering the vertical structure of the 1-m case.On the other hand, we can only speculate on the role of vertical viscosity coefficients, and it is plausible that vertical viscous damping of horizontal motion influences the resolution sensitivity displayed in Fig. 8.
The role of numerical diffusion is dependent on the choice of advection scheme.In general, this choice is based on a trade-off between the smoothness and the conservation of scalar fields (Griffies et al. 2000).Upwind schemes tend to produce smooth fields but are prone to spurious numerical diffusion (Griffies et al. 2000).Centered schemes can generate near-grid-scale noise (dispersion errors) but minimize numerical diffusion (Griffies et al. 2000).To illustrate the influence of numerical diffusion in our simulations, we compare both types of advection scheme.Figure 9 shows a comparison of Dr for simulations with centered and upwind schemes.In all simulations, the diffusivity is held constant.The default in this paper is the upwind scheme (Fig. 9a), chosen because it is stable, produces a bounded solution, and is commonly used with MITgcm.The results show less sensitivity of the solution to resolution with the linear-centered scheme.The sensitivity of the profiles to advection scheme provides further confirmation of the importance of numerical diffusion.
It is well known that dispersive errors and numerical diffusion both scale with resolution (e.g., Griffies et al. 2000), but their relative influence in comparison to resolved and parameterized mixing, across different resolutions, has not been evaluated for the ice shelf-ocean boundary current.Our investigation identifies that spurious numerical mixing is an important source of resolution sensitivity with commonly used advection schemes, and this sensitivity is greatly reduced at 1-, 2-, and 4-m resolution.We experiment with the use of two advection schemes, but many more choices exist.Furthermore, our experiments are limited to constant coefficients of diffusivity and viscosity, the simplest possible turbulence closure.This was chosen to enable control over the magnitude of coefficients, and thus make it possible to isolate parameterized and numerical diffusion.Schemes such as Smagorinsky (1963) or the AMD (Rozema et al. 2015) would be a more appropriate choice for representation of the boundary current, since they vary coefficients in space and time.These schemes would not eliminate numerical diffusion, but there may be feedbacks between variable coefficients and advection schemes.An investigation into the optimal advection scheme for representing boundary layer dynamics at coarse resolution and how results are influenced by variable coefficients of viscosity and diffusivity could form the basis of a useful follow-up study.

c. Interactive melting and time scales
So far, we have used constant heat and salt fluxes to represent a fixed melt rate at the upper boundary.In reality, melt rates are variable and dependent on ocean conditions.An interactive melt rate allows for a more realistic feedback between ice and ocean.Accurate interactive melt is the central requirement of modeling an ice shelf cavity.Here we investigate the role of interactive melting (three-equation model) at different resolutions.
Figures 10a-c show the difference between using interactive melting and a constant heat and salt flux.With a constant flux, lower resolution triggers a deeper diffusion of meltwater, reflecting the higher levels of mixing discussed in section 4b.For interactive melting, the relationship with resolution is different, and density anomalies reduce throughout the water column as resolution is decreased.At 1-m resolution, the profiles are similar for each choice of boundary condition, but as the resolution decreases, an increasing disparity emerges.For a given time, interactive melting exhibits reduced density anomalies in comparison to the constant-flux cases, with the largest difference at coarse resolution.This demonstrates that interactive melting leads to lower overall meltwater content, thus indicating reduced melt rates.
Figure 10g shows the melt-rate evolution over time for each case.The constant-flux case is designated with a constant melt rate of 10 m yr 21 .The interactive melt cases follow a similar evolution to each other but on differing time scales.The melt rates for lower-resolution results are much slower to build.The vertical dashed line signifies the time at which the results are sampled in Figs.10a-10c.Melt rates are lower for lower resolution in the run up to this point, which creates a lower accumulated melting (Fig. 10h).
The interactive melt is a function of temperature, salinity, and velocities at the ice base.In the absence of tides, the boundary current beneath ice shelves features a strong feedback whereby melting creates buoyant meltwater, which drives a stronger current and hence more melting (Holland et al. 2008;Jenkins 2016Jenkins , 2021)).Enhanced mixing at lower resolution combined with this feedback in the dynamics causes a delay for lower-resolution cases to generate equivalent melt rates to the higher-resolution counterparts.Lower-resolution cases have lower near-ice density anomalies for a given quantity of melt (Fig. 7), which leads to reduced velocity magnitudes.Since melt rates are dependent on velocities, the acceleration of melting is dampened.
Our investigation into the constant-flux configurations omitted any thermodynamic feedbacks at the ice-ocean interface, since there was a prescribed amount of meltwater in the system at any one time.This is not the case when presenting the interactive-melt model.In Figs.10d-f, rather than plotting results at a given time, we plot results after a given accumulated melting and separate the effect of varying meltwater content.The times chosen are marked by the dots on Fig. 10h.Sampling according to accumulated melt (Figs.10d-f) vastly reduces the differences between profiles for the interactive-melt (dashed) and constant-flux (solid) cases at each resolution in comparison to sampling at a specified time (Figs.10a-c).The similarity between profiles at each resolution in Figs.10d-f shows that disparity between the methods is almost entirely down to the meltwater content, and there is little evidence of additional dynamical feedbacks due to melt rates evolving on different time scales.Nonetheless, interactive melting feeds back onto the density profiles at lower resolution in a way that, when viewed at a fixed time, leads to a larger sensitivity of the vertical density structure to changes in resolution.Getting the time evolution of the boundary current and consequently melt rates correct at lower resolution is imperative for transient solutions.

Discussion, conclusions, and modeling suggestions a. Discretization effects on melting
In most of this study, we prescribe a constant flux of heat and salt at the upper boundary to represent melting, but in section 4c, we investigate the role of resolution on the threeequation melting model.There are various ways to discretize the three-equation model, and this choice has a strong influence on melt rates (Gwyther et al. 2020).Our simulations simply sample tracers and distribute meltwater in the top cell.Other methods involve setting these positions to fixed distances from the ice base or relative to the number of vertical model cells.Our results show that melt rates are dependent on resolution.Results from Gwyther et al. (2020) indicate that this resolution dependence may reduce if sampling and meltwater distribution positions are held constant with changes in resolution.On the other hand, if these positions are set relative to the number of grid cells, resolution dependence would be greater.
The z-level models tend to be biased toward enhanced diapycnal mixing, leading to thicker boundary layers and higher melt rates at coarse resolution (Losch 2008;Gwyther et al. 2020).As vertical resolution increases, cavity-scale simulations can produce strong stratification that suppresses this mixing near the boundary (Gwyther et al. 2020).Distributing meltwater away from the boundary has a similar effect to imposing higher mixing rates near the boundary and can increase melt rates at high resolution (Gwyther et al. 2020).This is a common practice in z-level models, usually implemented over the top two cells.When used over multiple cells, it reduces the differences in melt rates between resolutions toward that of the lower resolution (Gwyther et al. 2020).Although this may suggest a negative outcome, simulations of Gwyther et al. (2020) are laminar, and it is possible that the presence of turbulence would similarly somewhat limit the stratification.We present turbulence-permitting simulations with a 1-m isotropic grid.Stratification increases with resolution in our results, but this does not appear to limit the turbulent transfer of heat toward the boundary to the same extent as Gwyther et al. (2020).A 1 m isotropic grid is not achievable for general ocean models, and turbulence-permitting simulations are not realistic on the cavity scale.Conceptually, although distributing meltwater over a number of cells can counteract the consequences of stratification, it may be an overly simplistic fix.Sophisticated parameterization of the boundary current provides a more targeted method for improving the fidelity of models.

b. Numerical diffusion in previous studies
We show that numerical diffusion can be an important influence in coarse-resolution z-level simulations of the ice shelf-ocean boundary current.Numerical diffusion is known to scale with resolution (e.g., Griffies et al. 2000), but its role is not accounted for in parameterized mixing of the ice shelfocean boundary current.Our results suggest that previous studies of ice shelf cavities with z-level models may have experienced excess diapycnal mixing due to numerical diffusion.This has significant implications in a setting where the pycnocline is a primary barrier to heat delivery and thus melting.Despite taking longer to build, melt rates peak at a higher level with coarser resolution in our transient results (Fig. 10g).Spurious homogenization implies enhanced melt rates due to elevated entrainment of heat from below.Numerical diffusion could partially explain the mixing bias observed in previous studies, which is greatly reduced with enhanced resolution (Losch 2008;Gwyther et al. 2020).

c. Comparison with 1D resolution dependence
Our work is complementary to a recent study that addresses resolution dependence in a one-dimensional model of the ice shelf-ocean boundary current (Burchard et al. 2022).Combining the three-equation model with a turbulence closure for the boundary current, they examine the effects of decreasing resolution (Burchard et al. 2022).Their simulations showed little sensitivity to resolution in regions of weak stratification, but, similar to our results, 2-m resolution was required to represent the full boundary current accurately.Burchard et al. (2022) attribute reduced dependence on resolution to the application of the law of the wall within the drag parameterization.Tests within MITgcm are in contrast with this conclusion, showing negligible changes when comparing the law of the wall with a bulk formula for momentum (not shown).It is worth noting, however, the conclusions of Burchard et al. (2022) are drawn under different experimental conditions.Burchard et al. (2022) use pressure-gradientforced simulations with no rotation and a different choice of roughness length.As discussed in section 2c, for example, a smaller roughness length leads to reduced changes in drag with resolution.This is something to bear in mind with future investigations.

d. Conclusions
The performance of models in simulating the ice shelfocean boundary current has a significant impact on our ability to constrain estimates of ice shelf melt rates and potential ice sheet mass loss.To improve model representation, there is a need to understand the dynamics of this setting and explore model capability.In this paper, we have investigated the dynamics of the ice shelf-ocean boundary current at turbulencepermitting scales.We then considered how the physics of the ice shelf-ocean boundary current is captured within the framework of a z-level model (MITgcm).
Similar to results in Vreugdenhil et al. (2022), we identify the formation of Ekman rolls in the MITgcm simulations.We also show that MITgcm exhibits the three-phase development of the boundary current outlined in Jenkins (2021).The MITgcm results show how modulation of the pycnocline in the one-dimensional model emerges as a series of instabilities when modeled in three dimensions.The Ekman rolls do not fully destabilize at any point and persist throughout this final phase.
It is not yet understood what influence Ekman rolls have on the third phase of development of the boundary current.We suggest that a dedicated modeling study is required to fully understand the role of Ekman rolls in the ice shelf-ocean boundary current.Moreover, they have so far only been observed within model simulations.Identifying the presence of these features in field observations would be highly beneficial to our understanding, advancing assessment of their relative importance and ubiquity.
In comparing MITgcm with a state-of-the-art LES and onedimensional model at turbulence-permitting scales, we find remarkable agreement between the three models.The LES is proven to perform well at representing wall-bounded stratified flows at turbulent scales (Vreugdenhil and Taylor 2018) and is used as a reference case.MITgcm has geopotential coordinates with partial cells for topography, whereas the LES uses a rotated grid with planar topography.Additionally, the LES has the superior AMD for scalar mixing (Rozema et al. 2015), which is not available for MITgcm.Agreement with the LES, both with and without the AMD, suggests that fundamental differences in the grid and mixing schemes do not have a prohibitively negative impact on the MITgcm results.However, the results show the grid and/or numerical methods may impact MITgcm's ability to reproduce equivalent levels of turbulent mixing near the boundary.For the one-dimensional model, we identify that configuration choice of the HTC scheme is an important factor in the agreement, and by shifting the taper region of the HTC, we show a considerable improvement in the consensus with MITgcm.As MITgcm also matches well with the LES, the implication is that the adjustment represents an improvement of the one-dimensional model.
Given the match in representation of the boundary current between the one-dimensional model, with an adjusted taper, and MITgcm at turbulence-permitting resolution, the HTC scheme of the one-dimensional model could form the basis of an effective parameterization of the boundary current.In its current form, this mixing scheme may face the same issues as MITgcm when reductions in vertical resolution are considered.A future investigation of the one-dimensional model undertaken in a similar manner to Burchard et al. (2022) would provide useful insight on this.
Our turbulence-permitting, z-level model has isotropic grid resolution and coefficients of diffusivity and viscosity.We degrade the model resolution uniformly in the horizontal and vertical.At 1-, 2-, and 4-m resolution, the results are qualitatively similar.However, in general, the boundary current becomes more diffuse with coarser resolution and the Ekman layer disappears.The excess diffusion can be partially remedied by tuning diffusivity coefficients, but this is not sufficient at coarse resolution, where the feature is modeled over a small number of grid cells, and numerical diffusion becomes a limiting factor.Thus, heat and momentum transfer have to be parameterized.Existing boundary layer parameterizations do not provide a solution, and interactive melting exacerbates these issues by introducing a time-scale dependence in the evolution of melt rates between resolutions.The time-scale dependence has implications on the transient response to changing conditions and may impact sensitivities such as those found in Holland (2017).We show this dependence can be mostly removed by resampling according to accumulated melt.
Cavity-scale models do not resolve the full ice shelf-ocean boundary current.Parameterization within MITgcm (threeequation model) is targeted at the ice shelf-ocean boundary layer, which typically spans the first meter of the boundary region.This is not necessarily the best approach for coarserresolution models that fail to represent the boundary current that extends further from the boundary.A z-level model presents the challenge that increasing vertical resolution at the ice base for a sloped ice shelf requires higher resolution for almost the entire domain.To overcome this, we suggest the development of a scale-aware boundary current parameterization that is independent of resolution.Our results provide a basis on which to make progress by highlighting some key weaknesses that exist with current strategies.
FIG. 1. Schematic representation of (a) the rotated domain for the 1D and LES models and (b) the shifted lateral boundary conditions undertaken for the z-level model.Shading represents the ice topography (blue), the solid lower boundary (gray), and the communication cells of the periodic boundary (orange).Dashed lines represent the grid layout.The ocean boundary is marked by the solid black line.In (b), black arrows show the alignment of existing boundary communication in the z-level model, and the red arrow shows an example of shifting this communication by one cell.The number of grid cells used in the model simulations far exceeds what is shown here.

FIG. 2 .
FIG. 2. Example of the z-level reference case (MITgcm_Refer_1m) with 1-m resolution.(a) A 3D view of the domain with vertical velocities plotted on each vertical face.(b) The time evolution of isopycnals over a vertical slice of the top 30 m of the domain.Gray shading represents the bottom/ice topography.The inset shows an enlargement of the ice shelf-ocean interface to highlight the stepped partial cell representation.
"constant flux" use a constant heat and salt flux that is equivalent to a b 5 10 m yr 21 .This value is chosen to approximately match the ablation seen for the MITgcm cases that use the three-equation model.The constant heat and salt fluxes are calculated as where L i 5 3.34 3 10 5 J kg 21 , c w 5 3974 J 8C 21 kg 21 , S b 5 30, r 0 5 1030 kg m 23 , and r i 5 916 kg m 23 .These parameter choices for the constant-flux condition follow Jenkins et al. (2010), whereas configurations in

FIG. 4 .
FIG. 4. Comparison of MITgcm (black lines; MITgcm_Refer_1m) against two LES configurations, where SGS mixing is represented by either constant background coefficients (cyan lines; LES_Const) matching MITgcm, where k 5 n 5 1 3 10 24 m 2 s 21 , or an SGS model (orange lines; LES_SGS).The MITgcm configuration is different from that shown in Fig. 3 (see Table 2 for differences).Profiles are iceplane-averaged snapshots for times shown in (a) over the top 30 m of the domain.(a)-(c) Density anomaly Dr, up-slope velocity u, and across-slope velocity y.Dashed black lines in (c) are estimates of geostrophic velocity y g .(d)-(f) Parameterized, resolved turbulent, and total (parameterized 1 resolved turbulent) stress for MITgcm and LES 1 SGS model configuration.(g)-(i) As in (d)-(f), except for the buoyancy fluxes.The stress is taken as the vector norm of the up-slope and across-slope directions.Profile times are differentiated by solid, dash-dotted, and dotted lines in (d)-(i).
FIG. 5. Evidence of Ekman rolls in the 1-m MITgcm simulation with constant heat and salt fluxes at the upper boundary (MITgcm_ Refer_1m).(a) A Hovm öller diagram of the ice-plane-averaged TKE.(b)-(d) Ice-base-parallel slices of vertical velocity at the depth and times shown by the markers in (a).Ice-plane-averaged profiles at times shown by the vertical lines in (a): (e) density anomaly Dr, (f) upslope velocity u, and (g) across-slope velocity y ; (h) total stress (resolved 1 parameterized) in the up-slope and (i) across-slope directions; (j) total buoyancy flux (resolved 1 parameterized).

FIG. 7 .
FIG. 7. Variation of resolution for MITgcm with a constant surface-flux boundary condition (MITgcm_Refer).Profiles are as in previous figures, except color represents resolution of each profile.Profiles are sampled at time day 8 hour 0. Viscosity and diffusivity coefficients and time step are scaled linearly with resolution.

FIG. 8 .
FIG. 8. Comparison of using constant diffusivity coefficients vs one that varies with resolution (MITgcm_Refer and MITgcm_Diffu).Columns show Dr, u, and y.Profiles from Fig. 7 with resolution-dependent diffusivity coefficients are shown as solid lines.Dashed lines show simulations where the diffusivity coefficient uses the 1-m value of k 5 1 3 10 24 m 2 s 21 for all resolutions.

FIG. 10 .
FIG. 10.Comparison of an interactive melting and the constant-flux boundary condition (MITgcm_3EQ and MITgcm_Refer).(a)-(f) Ice-plane-averaged profiles of Dr, u, and y at different resolutions with and without the three-equation model applied.(g), (h) Melt rates and accumulated melt for each simulation over time.In (a)-(c), all profiles are sampled at day 8 hour 0, as before.In (d)-(f), profiles are sampled according to a chosen accumulated melt highlighted by the dots in (h).

TABLE 2 .
List of model configurations.The "Advec."column displays the choice of advection scheme according to the MITgcm numbering system.The "Scalars" column refers to the number of scalars in the equation of state.SGS stands for subgrid-scale and refers to the use of diffusivity and viscosity coefficients that vary in both time and space."Dims."shows the number of spatial dimensions for each configuration.