Shear, Stability, and Mixing within the Ice Shelf–Ocean Boundary Current

: When the inclined base of an ice shelf melts into the ocean, it induces both a statically stable stratiﬁcation and a buoyancy-forced, sheared ﬂow along the interface. Understanding how those competing effects inﬂuence the dynamical stability of the boundary current is the key to quantifying the turbulent transfer of heat from far-ﬁeld ocean to ice. The implications of the close coupling between shear, stability, and mixing are explored with the aid of a one-dimensional numerical model that simulates density and current proﬁles perpendicular to the ice. Diffusivity and viscosity are determined using a mixing length model within the turbulent boundary layer and empirical functions of the gradient Richardson number in the stratiﬁed layer below. Starting from rest, the boundary current is initially strongly stratiﬁed and dynamically stable, slowly thickening as meltwater diffuses away from the interface. Eventually, the current enters a second phase where dynamicalinstabilitygeneratesarelativelywell-mixed,turbulentlayeradjacenttotheice,whilebeneaththecurrentmaximum, strong stratiﬁcation suppresses mixing in the region of reverse shear. Under weak buoyancy forcing the time scale for development of the initial dynamical instability can be months or longer, but background ﬂows, which are always present in reality, provide additional current shear that greatly accelerates the process. A third phase can be reached when the ice shelf base is sufﬁciently steep, with dynamical instability extending beyond the boundary layer into regions of geostrophic ﬂow, generating a marginally stable pycnocline through which the heat ﬂux is a simple function of ice–ocean interfacial slope.


Introduction
Beneath the Antarctic ice shelves, which together cover 1.6 3 10 6 km 2 of the Southern Ocean (Fretwell et al. 2013), the interaction between ice and ocean removes just over half the mass that initially falls as snow over the Antarctic Ice Sheet (Rignot et al. 2013) and creates water masses that in key locations contribute to Antarctic Bottom Water (AABW) formation (Nicholls et al. 2009).Recent variability in that interaction is thought to be responsible for thinning of many of the ice shelves (Paolo et al. 2015) and to have contributed to freshening of AABW precursors (Jacobs and Giulivi 2010).Since the ice shelves restrain outflow from the inland ice sheet, their thinning has been accompanied by acceleration of outlet glaciers and an overall loss of grounded ice to the ocean (Shepherd et al. 2018).The impacts of such changes on ocean circulation and sea level, now and in the future, have motivated efforts to incorporate ice shelf-ocean interactions into ice sheet, ocean, and climate models.
A key process controlling the interaction between ice shelves and the underlying ocean is the exchange of heat and freshwater across the ice-ocean turbulent boundary layer.Our understanding of that process is based almost exclusively on studies of the analogous boundary layer beneath sea ice (McPhee 2008).The fundamental physics governing the production, transport, and dissipation of turbulence within, and the resulting mixing of momentum and scalars across, the boundary layer should not differ.However, there are subtle but important distinctions in the physical processes that generate turbulence.Beneath melting sea ice, turbulence is generated by wind-forced motion of the hydraulically rough ice across the ocean surface.Thus, the source of energy to generate turbulence and overcome the stabilizing buoyancy flux is external to the ice-ocean boundary layer.In contrast, beneath an ice shelf, the currents driven by the buoyancy forcing associated with the melting ice provide a key source of energy for turbulence in the boundary layer, although generally supplemented by tidal forcing (Makinson et al. 2011;Jourdain et al. 2019).
The close coupling between the ice shelf-ocean boundary layer exchanges and the large-scale circulation driven by the resulting density gradients represents one of the main challenges for models of the sub-ice-shelf circulation.Estimates of melting and freezing at the base of the ice shelves (Rignot et al. 2013) provide the most comprehensive datasets against which to evaluate model performance.However, it is difficult to assess whether model biases result from poor performance of the parameterizations of turbulent mixing within the boundary layer or from a poor representation of the delivery of heat to the boundary layer by the larger-scale circulation.A poor representation of the circulation could result from incomplete knowledge of the sub-ice-cavity geometry or of the density forcing generated beyond the cavity but will be compounded by biases in the buoyancy fluxes simulated within the cavity.
Given the paucity of observations within the ice shelf-ocean boundary layer, and the difficulties in significantly expanding that database when each observation requires an access hole to be made through, or the deployment of an autonomous vehicle beneath, hundreds of meters of ice (Stanton et al. 2013;Kimura et al. 2015;Davis and Nicholls 2019), process-oriented models can provide critical insight.Direct numerical simulation of near-ice turbulent mixing (Gayen et al. 2016) and large-eddy simulations of the boundary layer (Vreugdenhil and Taylor 2019) represent potentially rewarding approaches.However, our knowledge of even the basic current structure beneath an ice shelf is almost completely lacking, so there is much that can be learnt from simpler models.This paper describes one such model and the insight that it gives into the problem.It is a development of an earlier study (Jenkins 2016) that explored the fundamental dynamical balance in a one-dimensional model of the ice shelf-ocean boundary current.That study made the unrealistic assumption of constant eddy viscosity and diffusivity and discussed the flow within and beyond the stratified Ekman layer that resulted.Here a turbulence closure scheme is added to that model to investigate the interaction between the buoyancy-driven flows and the turbulent mixing that results, and that ultimately dictates the density structure.
The following section gives an overview of the earlier model of Jenkins (2016) and describes the extensions introduced for the present study.The impact of the added turbulence closure scheme is then discussed, both in the context of a conventional ice-ocean boundary layer beneath a horizontal interface and the inclined ice shelf-ocean boundary current.The sensitivity of the solutions to changes in far-field forcing and interface slope are presented, followed by solutions that result from combinations of interface slope and background pressure gradient.The fundamental structure of the ice shelf-ocean boundary current and the interfacial fluxes that it generates are then discussed, while overall findings and their applications are summarized in the concluding remarks.

Model
The starting point is the model of Jenkins (2016) that describes the cooling and dilution of seawater by interaction with an overlying ice shelf and the resulting flow of the modified waters along the ice shelf base through one-dimensional diffusion equations for momentum and thermal driving: (2) The equations have been transformed into a translated and rotated coordinate system (Fig. 1), in which z is the coordinate axis perpendicular to a planar ice-ocean interface that slopes at an angle a to the horizontal, and x, y lie in the plane of the interface with x pointing directly upslope.The momentum equation in (1) uses both Boussinesq and hydrostatic approximations, with the latter applied perpendicular to the sloping interface, so its validity is not limited to small slopes.For convenience, the velocity vector parallel to the interface is expressed as a complex number: u 5 u 1 iy , and temperature and salinity have been combined into a single scalar, the thermal driving, defined as the temperature relative to the freezing point at the interface: The combination of temperature and salinity into a single scalar requires the implicit assumption that the Lewis number is one everywhere.In the rotated coordinate system, the Coriolis parameter is given by Variables and physical constants are defined in Table 1.Time derivatives [first terms in (1) and ( 2)], Coriolis acceleration [second term in (1)], and terms for turbulent diffusion perpendicular to the interface [last terms in (1) and ( 2)] take their conventional forms, but the pressure gradient parallel to the interface has two components.The large-scale slope of the interface is set such that there is no forcing when the density deficit Dr relative to the uniform ambient fluid is zero.Thus, when the ice-ocean interface is at z 5 0, the ice shelf is passively afloat in stationary ambient fluid.Flow is forced either by a nonzero density deficit, which creates an upslope buoyancy force [first term on the right-hand side of (1)], or nonuniform displacement of the interface from is equilibrium position (z 5 0) to z 5 h, which creates a depth-independent pressure gradient parallel to the gradient vector of h [second term of on the righthand side of (1)], which is expressed as a complex number: If the ice-ocean interface were horizontal (sina 5 0), the first of the above terms would vanish, irrespective of the density deficit, and the second would be a conventional sea surface slope forcing.With the exception of the prescribed =h term, all gradients parallel to the interface are assumed to be zero, although, in principle, fixed gradients could be prescribed (Jenkins 2016).The neglect of gradients parallel to the interface implies that the only source of buoyancy is diffusion of the melt signal away from the interface.A more traditional route to a reduced-physics model of buoyancy-driven flow along an ice shelf base would retain either one or both of the spatial dimensions parallel to the iceocean interface that have been dropped from (1) and (2), while depth-integrating over the dimension perpendicular to the interface that has been retained above (MacAyeal 1985;Jenkins 1991;Holland and Feltham 2006;Jenkins 2011).Depthintegration of the last terms on the right-hand sides of ( 1) and (2) reduces them to fluxes at the upper and lower limits of the integration, generally assumed to the boundaries of a distinct turbulent layer.Fluxes at the ice-ocean interface are defined in a manner analogous to that described below (Jenkins 1991), while those at the outer edge of the turbulent layer are parameterized as if the layer could be considered an inclined, entraining plume (Ellison and Turner 1959).The thickness of the layer is calculated from a depth-integrated version of the continuity equation, which in the model presented here reduces to the trivial form ›w/›z 5 0 because of the neglect of horizontal gradients.The key advantage of the depth-integrated layer approach is the inclusion of the horizontal advection terms that are essential for simulating processes such as the accumulation of marine ice at the ice shelf base (Jenkins and Bombosch 1995).The main disadvantage is the loss of information on the vertical structure of the turbulent layer that is critical for representing the processes of heat transfer across that layer from far-field ocean to ice.The model presented here focuses on that vertical structure at the expense of losing the advection terms, making it most appropriate to the early stages of turbulent layer growth when those terms are small (Lane-Serff 1995).A possible strategy for combining the best aspects of both approaches is discussed in the concluding remarks.
The model equations, (1) and ( 2), are coupled via the dependence of the dimensionless density deficit on the thermal driving deficit (Jenkins 2016): where the subscript a indicates the defined far-field properties of the ambient ocean and ''i'' indicates ice properties (Table 1).Equation ( 4) is derived from the combination of a linear equation of state: which can be rewritten as the thermal driving definition (3) rewritten similarly: and an expression that relates the cooling and dilution of waters that interact with a melting ice shelf, analogous to that discussed by Gade (1979): Application of (5) requires the assumption that the ambient fluid has uniform properties and that the Lewis number is one FIG.1. Schematic illustrating the transformed coordinate system (x, y, z) in which the model is formulated, with x upslope and y cross slope in the plane of the ice shelf base, and its relationship with the conventional system (x 0 , y 0 , z 0 ), which has x 0 zonal, y 0 meridional, and z 0 vertical, with positive up and the origin at sea level.Hereinafter, it will be denoted by Dr/DT * .
The ice-ocean interface is assumed to be at the freezing point at all times and a no-slip condition is applied, leading to boundary conditions of zero velocity and thermal driving at the interface, while ambient properties and geostrophic flow are imposed in the far field: u 5 0, T * 5 0 at z 5 0, and The above model could now be supplemented with a turbulence closure scheme of arbitrary complexity.However, the aim of the present study is to explore the qualitative structure of the boundary current that emerges when the viscosity and diffusivity vary in a physically meaningful way with distance from the ice-ocean interface.That aim motivates a relatively simple scheme.The local turbulence closure (LTC) of McPhee (1994) fits that requirement and moreover has been shown to perform as well as much more complex closures when compared with observations from the sea ice-ocean boundary layer (McPhee 1999).
The LTC expresses the eddy viscosity and diffusivity as products of the local friction velocity and a mixing length: where the subscript mol indicates molecular values and the friction velocity is derived from the local shear stress: The mixing length scales with distance from the interface according to Monin-Obukhov similarity theory until it reaches a limiting value: where k is von Kármán's constant, and is the Monin-Obukhov length scale estimated from the interfacial friction velocity and buoyancy flux.The interfacial buoyancy flux is related to the thermal driving flux as before: and the problem of estimating thermal driving and momentum fluxes at the interface will be returned to shortly.Under a stabilizing buoyancy flux the limiting length scale is where R c is a critical flux Richardson number, but as L 0 / ' under nearly neutral conditions, rotation restricts the growth of the mixing length beyond a second limit given by where L * is a similarity constant.McPhee (1994) combines these two limiting length scales into one expression that is onehalf of their harmonic mean: The above expression asymptotically approaches the smaller of the limiting length scales in the case that one is much larger than the other.While the LTC as presented thus far is suitable for most of the boundary layer, it cannot be applied to the region very close to the interface, where the presence of the solid boundary suppresses eddies to such an extent that molecular diffusion becomes the dominant process of heat and mass transfer within an interfacial sublayer (McPhee et al. 2008).Such processes are not captured by the above equations, and the model grid is anyway too coarse to resolve them explicitly, so they must be parameterized using an interface submodel (McPhee 2008).
The parameterization of momentum transfer uses a conventional quadratic drag law to express the magnitude of the interfacial shear stress as with the drag coefficient derived from the ''law of the wall'' for a boundary with a specified roughness length of z * 0 : where z is the distance from the interface at which u is defined (McPhee 2008).In the case of a hydraulically smooth interface, the effective roughness length becomes z * 0 5 n mol e 22 u * 0 , while for a rough interface it is related to the characteristic physical height of the roughness elements by (McPhee 2008) z * 0 5 z r /30.
To apply this parameterization within a fixed grid model, the viscosity across the upper-most grid box adjacent to the iceocean boundary is chosen to ensure that the correct turbulent stress is recovered: where dz is the size of the first grid box.
The interfacial sublayer is particularly influential for scalar transfer within the ice-ocean boundary layer because of seawater's large molecular Prandtl and Schmidt numbers.Large parts of the temperature and salinity changes across the boundary layer therefore occur within the interfacial sublayer.Observations beneath sea ice (McPhee et al. 1999) suggest that transfer of heat across the inner part of the boundary layer is well represented by an equation of the form where G is a Stanton number based on the friction velocity that has been found to be constant over a wide range of conditions (McPhee 2008).A theoretically more justifiable formulation would include separate expressions for the heat and salt fluxes, using thermal and haline Stanton numbers that differ in magnitude, reflecting the role of molecular diffusion in setting the transfer rates within the interfacial sublayer.The differences between the two approaches are discussed in Holland and Jenkins (1999), where it is shown that the nonlinearity associated with the more complex formulation leads to interfacial heat fluxes that differ by no more than 10% from those given by ( 7) for typical oceanic conditions.Here, the simpler expression ( 7) is favored for two reasons: it is the one conventionally used in the analysis of ice-ocean heat flux observations, so that the appropriate Stanton number is better constrained by data; it is in a convenient form for the model described above, where the equations for temperature and salinity have been combined into one for T * .As before, the parameterization is used to set the diffusivity across the upper-most grid box such that the correct heat flux is recovered: The above approach requires an implicit assumption that the first grid point lies beyond the thermal boundary layer, but within the logarithmic layer.It converges numerically as the grid size is reduced but would become physically unrealistic if the first grid point were to lie within the interfacial sublayer.A more general alternative would be to place the upper boundary a set distance from the interface and apply a Neumann boundary condition at that grid point, but the simpler approach here is justifiable given that the standard grid resolution is ;0.5 m.Application of the LTC to the turbulent ice-ocean boundary layer beneath a horizontal interface produces a well-mixed boundary layer that is separated from the far-field ocean by a pycnocline.That motivates application of a second set of mixing length scales in the pycnocline, derived as above, but using the friction velocity and buoyancy flux estimated at the base of the mixed layer, and defining the transition between the mixed layer and pycnocline based on a density or stratification criterion (McPhee 2008).In this study, the pycnocline is defined using a Richardson number criterion, while two approaches to estimating mixing there are adopted.The former follows the LTC mixing length scheme, while the latter defines viscosity and diffusivity as empirical functions of the gradient Richardson number following Pacanowski and Philander (1981): where the gradient Richardson number is given by Ri 5 and values used for the empirical constants are given in Table 1.Parameterizations of this form provide the simplest solution to the problem of simulating mixing due to resolved shear in the stratified ocean interior (Fox-Kemper et al. 2019) and are used to estimate diffusivity in the pycnocline in other turbulence closures (Large et al. 1994).
The two pycnocline schemes, the latter denoted PP, are merged with the LTC boundary layer scheme by defining depth limits based on where the gradient Richardson number first exceeds values of 0.25 and 1.The boundary layer scheme is applied up to the first limit, the pycnocline scheme beyond the second limit, and between the limits a linear combination of the two estimates of viscosity and diffusivity is used, such that For convenience, the scheme that merges LTC in the turbulent boundary layer with PP in the pycnocline is referred to hereinafter as the hybrid turbulence closure (HTC) model.The above equations are solved numerically on a staggered grid with viscosity and diffusivity defined on intermediate points and calculated using velocity and thermal driving data from the previous time step.An explicit time-stepping scheme is used with an adaptive step size, based on grid size and maximum diffusivity, to ensure stability.The domain extends to 15(2n n /jfj) 1/2 , sufficient for the far-field boundary condition to play no role in the final solution.Initial conditions have farfield properties everywhere except at the ice-ocean interface, where the boundary condition is applied.Initial investigations use four options to specify eddy viscosity and diffusivity: constant values; PP; LTC, which includes the interface submodel and modified mixing lengths in the pycnocline; and HTC, which combines LTC with the interface submodel in the boundary layer, and PP in the pycnocline.Subsequent discussion focuses on HTC, as it is found to give the most robust results.

a. Impact of turbulence closure model on stratification and current structure
The evolution of the ice-ocean boundary current beneath a horizontal interface is illustrated in Fig. 2. Solutions for constant viscosity and diffusivity (set equal to n n and 0.1n n , respectively) are analogous to those described in Jenkins (2016).The thermal driving and velocity are uncoupled, so while the former evolves as an error function, the latter adjusts to produce a steady Ekman layer adjacent to the ice shelf base.Within the Ekman layer, friction induces shear in the otherwise depth-independent current, so the application of any turbulence model creates a region of enhanced diffusivity within the Ekman layer.That changes the stratification fundamentally (Figs.2b-d), resulting in a well-mixed layer above a sharp pycnocline.
Using PP only (Fig. 2b), the computed diffusivity rises monotonically toward the interface, where current shear is maximum.There is a positive feedback, in that weakening of the stratification near the ice-ocean boundary further enhances the diffusivity there.The pycnocline is sharpened by the opposite feedback, wherein low diffusivity enhances the stratification, leading to a further reduction in diffusivity.In the velocity structure, the basic form of the Ekman solution is still visible, albeit modified by the depth-varying viscosity.
The LTC (Fig. 2c) includes two important physical controls on mixing that are not considered in PP.The size of the eddies is restricted both by the presence of the solid boundary and by the stabilizing buoyancy flux at the interface.The latter control is there even when mixing is sufficient to destroy the stratification near the boundary.With the inclusion of those processes, the LTC gives a diffusivity that rises rapidly with distance from the boundary, before decaying to background levels beyond the boundary layer where the current shear vanishes and the stratification in the pycnocline restricts the vertical scale of the eddies.The interface submodel yields very low diffusivity next to the ice-ocean boundary, resulting in a sharp jump in properties across the first grid box.That jump, which maintains significant thermal driving within the mixed layer, represents the major difference between the PP and LTC models.As the thermal driving in the mixed layer decays, the stabilizing interfacial buoyancy flux weakens and the diffusivity in the mixed layer grows.The pycnocline is stronger than in the PP result, but its depth is very similar, despite the very different turbulence closure assumptions.
That similarity suggests that, while the near-boundary processes included in LTC are critical to estimating mixed layer properties, the structure of the pycnocline is not sensitive to the choice of turbulence closure.Indeed, the HTC (Fig. 2d) yields results that are very similar to those produced by LTC.The velocity profiles are modified in a similar way, with relative uniformity within the mixed layer and steep gradients through the pycnocline.Because the forcing on the flow is externally imposed and unrelated to the stratification, the form of the current profiles changes little as the mixed layer progressively cools and deepens.
The introduction of an ice-ocean boundary slope (now with zero background pressure gradient, Fig. 3) couples the thermal driving and velocity profiles because the thermal driving deficit now represents a buoyancy forcing on the flow.The resulting current shear enhances the viscosity and diffusivity beyond the frictional boundary layer, so the spread of the boundary current found in the constant diffusivity case is also a feature of the other models.The introduction of stability-dependent mixing in PP (Fig. 3b), LTC (Fig. 3c) and HTC (Fig. 3d) leads to the development of a relatively well-mixed boundary layer above a broad pycnocline.However, the transition between them is more gradual than in the case of a horizontal interface (Fig. 2), and it is no longer possible to define an unambiguous mixed layer depth based on stratification.Instead, the term ''turbulent layer'' is applied to the region where the gradient Richardson number is less than 1, with the transition to the pycnocline then occurring at the bottom of the bold line segments in Fig. 3.
As before, cooling of the turbulent layer leads to enhanced diffusivity near the ice-ocean interface because of the weakening stratification (in the case of PP) and the weakening interfacial buoyancy flux (in the case of LTC and HTC).However, the cooling now also increases the density contrast across the pycnocline, strengthening the flow of the turbulent layer, and further enhancing the diffusivity through the greater frictionally generated current shear.Although the boundary current grows in thickness and the buoyancy-forced cross-slope current increases, the peak speed of the frictionally driven upslope flow remains relatively constant, as discussed by Jenkins (2016) for the constant diffusivity model.
The coupling between the thermal driving deficit and the flow gives rise to three distinct phases in the temporal development of the current at the sloping ice shelf base.Initially the thermal driving deficit is confined to a region very close to the ice.The flow is weak, and the stratification is strong enough that the current is initially stable everywhere.Within the HTC, diffusivity and viscosity are estimated exclusively from the PP model, while within LTC, the pycnocline scaling is used.The shear is maximum near the interface, so the weakening stratification leads to an initial instability there and the resulting introduction of the LTC submodel produces a jump in properties over the first grid box (Figs.3c,d).During the ensuing ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2n n /jfj p , the far-field thermal driving T * a , the geostrophic current at the interface y gi 5 (g cos a/f)=h, and the diffusivity under neutral conditions n n , as defined in Table 1.
second phase of evolution, the flow is sufficiently strong for shear instability to generate turbulence between the current maximum and the ice shelf base.Beyond the current maximum the pycnocline remains stable.Eventually, however, as the flow continues to grow in strength and the stratification weaken, the current enters a third phase.Instability extends beyond the current maximum, while a state of marginal stability is established through the pycnocline, where the shear is now primarily geostrophic.The condition for marginal stability (Ri 5 1) is first met at the lower limit of the thicker lines in Fig. 3 and is maintained through the region of constant vertical gradients (Fig. 3d).Its origin will be discussed later.During the subsequent development of the boundary current, the structure is preserved, with a marginally stable, geostrophic region, through which shear and stratification remain constant, beneath a relatively well-mixed, turbulent layer that grows in thickness as it cools and accelerates.Note that the growth in thickness of the boundary layer is controlled by the increase in viscosity, which leads to a growth in the effective Ekman depth, and hence an expansion of the region where frictionally FIG. 3. As in Fig. 2, but showing profiles obtained using (a) constant viscosity/diffusivity, (b) PP, (c) LTC, and (d) HTC for a flow generated by far-field thermal driving of 28C along an ice-ocean boundary with equilibrium slope sina of 0.01 and interface displacement gradient =h of 0. In the middle row, u (solid) is upslope and y (dashed) is cross-slope.Black dashed lines in (d) (top and middle panels) indicate the thermal driving and geostrophic current profiles for marginal stability (Ri 5 1), based on gradients given in ( 9) and ( 10) plotted at an arbitrary depth chosen to coincide with the pycnocline.Scales used to nondimensionalize the results are as in Fig. 2, except the geostrophic current at the interface, y gi 5 2(g sin a/f)T * a (Dr/DT * ).generated shear enhances mixing.The condition of marginal stability is not well captured by LTC, because perturbations around that state can give a jump in effective mixing length, leading to artificial steps in the pycnocline (Fig. 3c), so the subsequent discussion will focus on HTC.
The above results are qualitatively insensitive to uncertainties in key parameters that define the component models of the HTC (Fig. 4).The characteristic roughness of a melting ice shelf-ocean boundary is largely unknown, yet it controls the magnitude of the interfacial fluxes.Larger roughness lengths lead to higher diffusivities throughout the boundary layer, so while the current grows in thickness more rapidly, the jump in properties across the interfacial layer, and the resulting properties of the turbulent layer, are unaltered (Fig. 4a).Results are more sensitive to changes in the Stanton number, which have a bigger impact on turbulent fluxes at the interface (7) than through the boundary layer (Fig. 4b).As the Stanton number falls the interfacial fluxes decline, so the turbulent layer warms  3d, and scales used to nondimensionalize the results are as in Fig. 3. and decelerates in response.Conversely, as the Stanton number rises the interfacial fluxes grow, so the turbulent layer cools and accelerates in response.Results are insensitive to the parameterization of mixing through the pycnocline.Changing the eddy viscosity under neutral conditions rescales the viscosity and diffusivity profiles obtained from PP, but has almost no impact (Fig. 4c), because the large contrast in diffusivities between turbulent layer and pycnocline is more influential than the absolute values in setting the water column structure.
Changing the arbitrary transition region between the PP and LTC turbulence closures has a similarly small impact (Fig. 4d).A sharper transition leads to a more abrupt decrease in diffusivity at the base of the turbulent layer, but its properties and the strength of the pycnocline are little affected.

b. Sensitivity to forcing
Changing the slope of the ice shelf base (Figs.5a,b) or the far-field temperature of the ocean (Figs.5c,d) alters the buoyancy forcing on the current.Increased buoyancy forcing yields a faster current and enhanced, shear-induced mixing, so the turbulent layer cools, as the greater diffusivity enhances heat loss to the ice, and deepens, as the greater viscosity increases the effective Ekman depth.For an equivalent change in the buoyancy forcing, changes in slope have a slightly greater impact than changes in temperature, because gravity acts more obliquely to the interface as the slope increases.As a result, changes in both buoyancy forcing and static stability, through the cosa terms in the Monin-Obukhov length (6) and gradient Richardson number (8), contribute to changing viscosity and diffusivity, whereas only the former acts for changes in far-field thermal forcing.Nevertheless, the impact is small and only apparent in the early stages (Figs.5a,c, bottom row).Once the boundary current is fully developed the main differences are in the marginally stable pycnocline, where the characteristic gradients are functions of interface slope, but independent of far-field thermal driving (Figs.5b,d), as discussed later.
Since these are all transient solutions, a thicker turbulent layer at any particular time indicates more rapid development of the flow rather than any differences in a long-term steady state, which cannot be attained with the present model.In reality, advection processes, missing in the model, would halt the evolution at the stage where they were sufficient to balance the divergence between the interfacial and pycnocline heat fluxes, giving a steady state that should be qualitatively similar to the transient solution at that time (Jenkins 2016).
Imposing a background pressure gradient in addition to the buoyancy forcing leads to much more rapid deepening of the turbulent layer, which develops over the first inertial period (Figs.6a,c).During that initial phase, the buoyancy forcing on the flow is weak, so the structure of the boundary layer is determined by the magnitude of the background flow and is almost independent of its direction.The current profiles shown in Fig. 6c are practically the same as those in Fig. 6a, except for a 908 cyclonic rotation.However, as the boundary layer cools, the buoyancy forcing on the flow builds and a directional asymmetry develops (Figs.6b,d).
When the background flow enhances the buoyancy-driven flow (Fig. 6b), there is a positive feedback as the developing buoyancy forcing enhances current shear, increasing the diffusivity and promoting deepening of the turbulent layer.In contrast, when the background and buoyancy driven flows are in opposition, the developing buoyancy forcing weakens the shear at the interface, causing the viscosity to drop and the turbulent layer to retreat.The effect is particularly marked when the background and buoyancy driven flows exactly cancel (Fig. 6b).Then the Ekman layer is arrested, the frictional shear at the interface vanishes and the resulting turbulent layer is shallower than that obtained with buoyancy forcing alone.An up-or downslope background flow, driven by a cross-slope pressure gradient, always enhances the diffusivity relative to the result with buoyancy forcing alone (Fig. 6d).The directional asymmetry in the currents is less marked because the background and buoyancy-driven flows interact only via their respective Ekman layers.

a. Structure of the fully developed ice shelf-ocean boundary current
Jenkins ( 2016) drew a distinction between the buoyancydriven boundary current and the inner boundary layer where friction plays a role in the force balance.Such a distinction is absent when the ice-ocean interface is horizontal (Fig. 2).In that case, the only deviation from the constant, far-field flow is caused by friction, which is the only source of shear instability to create turbulence.The boundary layer and current are then coincident, composing a relatively well-mixed layer above a sharp pycnocline. Figure 3 shows that when the interface is sloped a relatively well-mixed boundary layer may also develop, given sufficient time.However, it is underlain by a broad pycnocline, which forms the outer part of the boundary current, where the inclined isopycnals drive a sheared, cross-slope, marginally stable, geostrophic current.The conditions necessary for the development of such a flow are discussed in the appendix.
When the pycnocline becomes marginally stable, its nature changes fundamentally from that of the stable pycnocline that forms beneath a horizontal ice-ocean interface.In the latter case (Fig. 2), the current shear is an externally defined parameter, so the pycnocline stability, determined by the gradient Richardson number (8), increases with increasing stratification.There is a positive feedback whereby increasing stratification reduces diffusivity and leads to higher thermal driving gradients, generating the sharp transition between turbulent layer above and strong, stable pycnocline below.However, when the interface is sloped (Fig. 3), the current shear everywhere is influenced by the stratification.For the fully developed boundary current, in its final phase of evolution, the pycnocline lies beyond the boundary layer, so the flow is primarily geostrophic, and stratification and shear are directly related through the thermal wind equation: The gradient Richardson number (8) can then be expressed as Ri ' implying the counterintuitive result that as the stratification increases the stability of the pycnocline decreases, because the increase in the current shear is the dominant effect.There is then a negative feedback, whereby increasing stratification destabilizes the water column and the enhanced mixing weakens the stratification.A balance is established through the pycnocline, in which a constant vertical gradient of thermal driving maintains the gradient Richardson number at the critical value where enhanced mixing starts (defined to be Ri 5 and from the thermal wind equation: The above relationships are plotted in Fig. 3d and explain why the thermal driving and current profiles through the pycnocline are so similar in Figs.4-6.There is only a very weak dependence on the far-field thermal driving, through the density parameter in ( 9), leaving interface slope as the single controlling factor (Fig. 5).The diffusivity and viscosity within the marginally stable pycnocline are constant (Fig. 3d) and given by PP with the gradient Richardson number set to 1: FIG. 6. Results obtained with HTC after (a),(c) 1 and (b),(d) 16 inertial periods applying the standard forcing in Fig. 3d with the addition of (left),(left center) varying upslope and (right center),(right) varying cross-slope interface displacement gradient.Scales used to nondimensionalize all results are as in Fig. 3.Note the differing horizontal axis limits along the bottom two rows.The e-X notation in the legend indicates that the leading number should be multiplied by 10 raised to the negative X. nj Ri51 5 2:4 3 10 24 and Kj Ri51 5 5:0 3 10 25 .Equations ( 9) and (10) thus yield simple, slope-dependent expressions for the heat flux and shear stress through the marginally stable pycnocline.Furthermore, they offer a simple explanation for why the pycnocline on very shallow slopes remains stable.The thermal driving gradient required for instability (9) increases as the slope decreases (Fig. 5b).Where the critical gradient exceeds that generated by the background diffusivity, typified by the stable pycnocline in Fig. 2d, that critical gradient can never be attained and the pycnocline will remain stable.

b. Interfacial fluxes for the fully developed ice shelf-ocean boundary current
The melt rate computed at 30 inertial periods for a wide range of slope and thermal driving conditions is shown in Figs.7a and 7b, along with the dependence of melt rate on thermal driving for the case shown in Fig. 2 with zero slope, but forced by a background pressure gradient.The slope  3 (magenta circles) and for Fig. 2 (green diamonds).In (c), results are shown for varying cross-slope (red crosses) and upslope (cyan squares) pressure gradient with all other parameters fixed at values used for Fig. 3, and for varying crossslope pressure gradient with all other parameters fixed at values used for Fig. 2 (green diamonds).Dashed lines, color coded, were generated using the theoretical relationships in ( 11), ( 12), and (13).
dependence of the melt rate (Fig. 7a) and the temperature dependence for zero slope (Fig. 7b) both show only weak nonlinearity, while the temperature dependence of melting for the sloped interface (Fig. 7b) is clearly nonlinear, as has been found in other studies (Holland et al. 2008).That behavior can be understood by inspection of (7), which expresses the interfacial heat flux as directly proportional to the interfacial friction velocity and the local thermal driving.For the horizontal interface, the friction velocity is approximately constant, while the similarity of the results for varying slope and varying temperature (on a nonzero slope) demonstrate that friction velocity is an approximately linear function of buoyancy forcing (Fig. 7a).That suggests a simple dependence of friction velocity on the geostrophic current within the turbulent layer via an approximately constant geostrophic drag coefficient and turning angle: Furthermore, the near-linearity of the melt rate relationship for zero slope indicates a simple dependence of the local thermal driving on the far-field conditions: Using ( 11) and ( 12), the melt rate can be expressed as where the magnitude of the geostrophic current in the turbulent layer is giving a quadratic dependence of melt rate on thermal driving for the buoyancy forced flow and a linear dependence for the pressure-gradient forced flow.Equations ( 11) and ( 13) are plotted in Figs.7a and 7b using a drag coefficient of 1.3 3 10 23 and a thermal driving coefficient, «, of 0.25 for the buoyancy forced cases and 0.08 for the pressure-gradient forced case.The magnitude of « is clearly time dependent (Figs. 2 and 3), but the fact that a single value can characterize all solutions at a particular time indicates that the underlying processes scale with the magnitude of the buoyancy forcing.The weak dependence of the geostrophic drag coefficient and turning angle (;158) on such factors as flow speed and stability is a result of the relatively smooth interfaces and low melt rates under consideration (McPhee 2012).
Melt rate and interfacial friction velocity are shown in Fig. 7c as functions of background pressure gradient for zero (as in Fig. 2) and standard (as in Figs. 3 and 6) slope cases.Again, the interfacial stress conditions are well represented using the same geostrophic drag coefficient (1.3 3 10 23 ) and turning angle (158).Using the same thermal driving coefficient (0.08) gives a good estimate of the melt rate for the zero slope cases, while a linear combination of buoyancy and pressure-gradient forced coefficients can characterize the melt rate dependence on the cross-slope pressure gradient.However, an analogous combination for the upslope pressure-gradient forced solutions cannot fully capture the asymmetry in the response, particularly the low melt rates when pressure gradient and buoyancy forcing are acting in opposition.It is noteworthy that in Fig. 6b, the configuration with h 0 x 5 1 3 10 25 (black line) drives a melt rate that is only half that of the configuration with h 0 x 5 0 (light brown line), despite the strong far-field current in the former.That is a result of the ''slippery'' boundary layer created when the buoyancy forcing acts to reduce the interfacial shear stress.

Summary and concluding remarks
Developing our understanding of ice-ocean heat transfer at the base of ice shelves is an important step toward explaining the observed mass loss from the Antarctic Ice Sheet and projecting its continued evolution under future climate change.This study represents one step in that process.It extends an earlier study (Jenkins 2016) that explored the fundamentals of the current structure adjacent to an inclined, melting ice-ocean interface.The key development is the addition of a simple turbulence closure model to account for the interactions between stratification and current shear and the resulting turbulent diffusivity and viscosity.Although the turbulence closure is simple, it combines elements that are either widely available in ocean models (PP) or have been extensively tested against observations of the turbulent boundary layer beneath sea ice (LTC).While more complex closures might produce quantitatively different results, the underlying qualitative behavior that is the main focus of the paper should be robust.It emerges as a result of the computed viscosity/diffusivity profile that first increases with distance from the interface, reaching a maximum in the boundary layer, before decaying to background levels.Such a profile is a fundamental result that must emerge from any turbulence model, regardless of its complexity.
During the early stages of its evolution, the boundary current is dynamically stable, and the limited changes in viscosity/diffusivity mean that the solutions conform quite closely to those of Jenkins (2016), in which the stratification decays monotonically with depth below the interface (Figs. 2  and 3).However, with the onset of dynamical instability, two inflections appear in the thermal driving profile separating a weakly stratified turbulent layer from much stronger stratification next to the interface and through the pycnocline.The changing stratification influences the buoyancy-forced geostrophic currents generated by the inclined interface, but the ageostrophic currents are controlled by the interfacial shear stress and are largely independent of the stratification (Jenkins 2016).Thus, the picture of a frictional boundary layer embedded within an otherwise geostrophic, cross-slope ice shelf-ocean boundary current described by Jenkins (2016) remains, although the stratification is now stronger across the outer, geostrophic current than within the inner, frictional boundary layer.The interplay between stratification, current shear and diffusivity in the outer, geostrophic part leads to a negative feedback between strengthening stratification and weakening stability that can maintain a state of marginal stability.Such a marginally stable pycnocline cannot develop on very low slopes, where the background diffusivity prevents the stratification reaching the required level (9).
If external forcing is absent, the flow along an inclined, melting interface must develop from an initially motionless state.While the thermal driving deficit is confined close to the interface there can be little flow and diffusivities are consequently low.As the flow builds, shear-induced mixing increases, allowing the thermal driving deficit to spread away from the interface, and that in turn increases the strength of the flow.Three distinct phases of development can be identified: an initially stable flow; a thin turbulent boundary layer underlain by a sharp pycnocline through which friction influences the current profile; a fully developed boundary layer underlain by a marginally stable pycnocline through which the shear is primarily geostrophic.The strength of the buoyancy forcing associated with both the slope of the interface and the far-field thermal driving determines how far the evolution can progress.Estimates of the time scales for transition from the first to the second phase (see the appendix) based on the constant diffusivity model of Jenkins (2016) suggest that under conditions quite commonly encountered beneath Antarctica's larger ice shelves (T * a ; 0.18C; sina ; 10 23 ) the transition would take months to years.However, at such low levels of thermal driving, background currents of ;10 22 m s 21 are sufficient to reduce that time scale to days.Currents associated with tides and even the time-mean circulation are typically an order of magnitude larger, suggesting that a turbulent ice-ocean boundary layer should be a widespread phenomenon.Estimates of the slope required for transition to the third phase (see the appendix) suggest that the fully developed ice shelf-ocean boundary current may be less widespread, but those estimates are conservative, being based on the constant diffusivity model that is applicable only during the first phase.
The boundary current structure described in this paper differs from the plumelike form often assumed (MacAyeal 1985;Jenkins 1991;Holland andFeltham 2006, Jenkins 2011).Although there is a relatively well-mixed core, its depth is set by the thickness of the frictional boundary layer, which is controlled by the flow speed through its role in setting the effective turbulent viscosity.In contrast, the flow speed of a plume sets the rate at which the plume thickness grows through entrainment.Beyond the boundary layer, the outer part of the ice shelf-ocean boundary current is stratified, and if the buoyancy forcing is strong enough to induce marginal stability, the heat flux through the pycnocline becomes a function of interface slope only (9).It is independent of the far-field thermal driving, because, as the thermal driving deficit in the turbulent layer grows, the pycnocline broadens to maintain the constant stratification (Fig. 5d).For a plume, the heat flux from entrainment is often set proportional to interface slope (Jenkins 1991), but it is also directly proportional to the thermal driving deficit within the plume because the pycnocline is treated as a discontinuity.Within the present model, the ice-ocean interfacial heat flux remains a function of far-field thermal driving (Fig. 7) because most of that flux is supplied by cooling and deepening of the turbulent layer.If advection were to be included in the model, a steady state could be obtained once the advective heat flux were sufficient to balance the divergence in the diffusive fluxes between interface and pycnocline.
While the conventional formulation of the depth-integrated plume equations naturally retains the along-slope heat advection, the main drawback is the reliance on an untested entrainment law to quantify the heat flux through the pycnocline.That entrainment law has been shown to hold for a wide range of natural and laboratory flows, including buoyancy-driven currents along inclined interfaces (Ellison and Turner 1959), but it has not been verified for cases where the current is subject to a stabilizing interfacial buoyancy flux.In those cases, turbulent kinetic energy must be expended to maintain wellmixed conditions within the current, and it is not obvious that the application of an identical entrainment law is appropriate.The results of this paper could potentially replace the entrainment law within a revised, depth-integrated layer model.Entrainment fluxes could be replaced with heat and momentum fluxes based on ( 9) and (10), while the layer thickness could be set proportional to the Ekman depth based on an effective viscosity that scales with the interfacial shear stress.Interfacial fluxes could be derived from ( 11) and ( 13), while the more complete heat conservation equation, including along-slope advection, would render (12) unnecessary.Appropriate geostrophic drag parameters are given in Table 2 for a range of surface roughness values.
Thus, while the study described in this paper is entirely idealized, the results potentially offer scope to improve the parameterizations of critical vertical mixing processes in models that simulate the ice shelf-ocean boundary current.It remains to be seen whether those improvements could represent a genuine advance in our ability to model the basal melt rates of ice shelves or merely provide an alternative parameterization of the many unknowns.At present there is a complete absence of observations that can be used to quantify directly the interdependence of the current and density profiles through the ice shelf-ocean boundary current that underpins the above findings.The results, therefore, remain hypothetical, and quantitatively dependent on the choice of turbulence closure.Nevertheless, they provide hypotheses that can be tested either by the application of physically more complete models or by targeted observational campaigns.To consider the development of the ice shelf-ocean boundary current from its initial stable state to its final form with a turbulent boundary layer and marginally stable pycnocline, it is instructive to examine the behavior of the solutions for constant viscosity and diffusivity shown in Fig. 3a.Since the initial condition has zero flow everywhere and a step in thermal driving at the interface, that step diffuses into the interior at a rate determined by the constant background diffusivity defined in PP, and the resulting flow is controlled by the constant background viscosity, even when HTC is used.Deviation from that solution occurs as the flow approaches dynamical instability, causing the viscosity/diffusivity to rise above background levels, but greatly elevated values only occur once Ri , 1 and LTC is introduced into the calculation.So, the constant diffusivity solution effectively captures the early stages of development and can be used to estimate an upper bound for the time taken for the current to become dynamically unstable.The thermal driving profile takes the form T * 5 T * a erf 2z 2 ffiffiffiffiffiffiffi ffi K b t p ! at time t (Fig. A1a), and the resulting buoyancy forcing gives rise to a geostrophic flow: , where is the magnitude of the geostrophic current at the interface (Fig. A1b).In association with the geostrophic current, an ageostrophic component develops over the first inertial periods (Jenkins 2016): , where is the Ekman depth for the background viscosity (Fig. A1c).
The current shear profile can then be expressed as ) , leading to the following expression for the gradient Richardson number (8) as a function of space and time: Ri 5 The basic form of the above expression is governed by the physical constants that determine the relative sizes of the thermal and dynamical boundary layers, while the overall magnitude is set by an inverse length scale, that can be expressed, using (A1) for the geostrophic current at the interface, as where the balance between buoyancy forcing of the shear flow and gravitational stability is encapsulated in the dimensionless factor: b Ri 5 cos a T * a sin 2 a Dr DT *
The evolution of the gradient Richardson number (A2) for the case plotted in Figs.3a and A1 (where b Ri 5 2 3 10 27 ) is illustrated in Fig. A2a.Results are inaccurate in the early stages, when the thermal boundary layer is shallower than the Ekman depth and the Ekman layer is consequently not fully developed.However, the analytical solutions capture the later development well (Fig. A1d).The current is least stable at the interface and the region of very weak stability deepens with time (Fig. A2a).There is a core of slightly greater stability around the current maximum that dissipates over time as the stratification weakens, then a return to lower stability in the region of reverse shear.Importantly, the region of instability extends well beyond the Ekman layer and into the region of primarily geostrophic shear beyond.In contrast, the gradient Richardson number for the case plotted in Figs.2a and A1 [evaluated from (A2) with only the first term in the denominator retained] shows that stability is rapidly established as shear-free stratification develops beyond the boundary layer (Figs.A1 and A2a).
With the introduction of the HTC, the solutions depart substantially from those in Fig. A1.Nevertheless, Fig. A2a provides insight into the early evolution of the boundary current before the onset of instability.The minimum gradient Richardson number occurs at the interface and decreases with time.While that minimum stays greater than one, the current will be stable everywhere and remain in the initial phase of its evolution.The transition to the second phase occurs almost instantaneously in Fig. A2a, but for weaker buoyancy forcing the time needed to reach that transition can be estimated from the long time scale evolution of the interfacial gradient Richardson number, derived from (A2) with only the first term in the denominator retained: The duration of the initial phase, relative to the inertial period t in , can be found by setting Ri 0 to 1 and rearranging: Figure A2b shows the duration of the initial phase as a function of ice shelf basal slope and thermal driving.There is a broad region where the current becomes unstable almost immediately, but the transition to a second broad region where the initial phase of stability is essentially infinite occurs at slopes and temperatures that can be encountered beneath real ice shelves.However, the neglect of advection in the model makes the results more appropriate to regions nearer the ice shelf grounding lines (Lane-Serff 1995), where basal slopes and farfield thermal driving tend to be higher, typically greater than ;10 23 and ;10 21 , respectively.Furthermore, the addition of a background current, which in reality is almost always present, provides an additional source of shear to trigger instability (Fig. A2b).Even when far-field thermal driving rises to a few degrees, background currents of a few centimeters per second are sufficient.In most cases, a rapid transition to the second phase of current evolution, with a turbulent boundary layer, should therefore be expected.Establishment of the third phase of evolution requires that the current be unstable beyond the frictional boundary layer.In Fig. A2a it can be seen that, for the sloped interface, the Richardson number near the outer edge of the boundary layer passes through a minimum about 15-20 inertial periods after the start.If that minimum exceeds one, the pycnocline that forms beyond the boundary layer will remain stable.Inserting appropriate numbers (z 5 pd E ) into the expression for the gradient Richardson number (A2) gives the minimum at the outer edge of the boundary layer as The critical Ri min bl contour is plotted in Fig. A2b.Beneath ice shelves, far-field thermal driving can range from a few tenths of a degree to a few degrees above freezing (Nicholls and Makinson 1998;Jenkins et al. 2010;Begeman et al. 2018), conditions that require minimum slopes from ;10 22 to ;10 23 for the final phase of boundary current evolution described earlier to be reached.Such basal slopes are common on ice shelves, particularly near their grounding lines, although the higher limit for colder waters suggests that the third phase may not be attainable everywhere.A caveat is that the above limits are only a rough guide.Once the boundary layer is unstable, a relatively well-mixed layer forms and the thermal driving profile changes fundamentally from the analytical solution (Fig. 3).Nevertheless, while the pycnocline remains stable, its advance into the ambient fluid will be controlled by the locally lower diffusivity, and its structure will conform more closely to the analytical solution.

FIG. 2 .
FIG. 2. (top) Thermal driving, (middle) velocity, and (bottom) viscosity/diffusivity profiles obtained using (a) constant viscosity/diffusivity, (b) PP, (c) LTC, and (d) HTC for a flow generated by far-field thermal driving of 28C at an ice-ocean boundary with equilibrium slope sina of 0 and interface displacement gradient ›h/›x of 25 3 10 26 .Solutions (color coded) are shown at 1, 3, 8, 16, and 30 inertial periods t in .Thicker lines indicate where the gradient Richardson number falls within the range 0.25-1.In the middle row u is solid and y is dashed, and in the bottom row viscosity is solid and diffusivity is dashed.Scales used to nondimensionalize the results are the Ekman depth d n E 5

FIG. 4 .
FIG. 4. Sensitivity of HTC results in Fig. 3d at 16 inertial periods to the value used for (a) the interfacial roughness length z * 0 , (b) the Stanton number G, (c) the diffusivity for neutral stability in the pycnocline n n , and (d) the range in critical gradient Richardson number used to define the transition between turbulence closures.Brown lines are obtained with the standard parameter values used in Fig. 3d, and scales used to nondimensionalize the results are as in Fig. 3.

FIG. 5 .
FIG. 5. Results obtained with HTC after (a),(c) 1 and (b),(d) 16 inertial periods for (left),(left center) varying ice shelf basal slope sina and (right center),(right) varying far-field thermal driving T a * .Results obtained with the standard forcing used in Fig. 3d are shown in brown.Scales used to nondimensionalize all results are as in Fig. 3.Note the differing horizontal axis limits along the bottom two rows.

FIG. 7
FIG. 7. (top)  Interfacial heat flux (expressed as a melt rate) and (bottom) shear stress (expressed as a friction velocity) computed at 30 inertial periods as a function of (a) interfacial slope, (b) far-field thermal driving, and (c) background pressure gradient.In all panels, the large black symbols indicate standard values used in Figs. 3 (circles) and Fig.2(diamonds).In (a), blue plus signs show results for varying interface slope, with all other parameters fixed at values used for Fig.3.In (b), results are shown for varying thermal driving, with all other parameters fixed at values used for Fig.3(magenta circles) and for Fig.2(green diamonds).In (c), results are shown for varying cross-slope (red crosses) and upslope (cyan squares) pressure gradient with all other parameters fixed at values used for Fig.3, and for varying crossslope pressure gradient with all other parameters fixed at values used for Fig.2(green diamonds).Dashed lines, color coded, were generated using the theoretical relationships in (11), (12), and (13).
FIG. A1.Numerical (symbols plotted for every fifth grid point, color coded by time) and analytical (black lines) solutions for constant viscosity and diffusivity (with Prandtl number 5 10).(top) The results in Fig. 2a; (bottom) the results in Fig. 3a for (a) thermal driving, (b) upslope (solid) and across-slope (dashed) geostrophic currents, (c) upslope (solid) and across-slope (dashed) ageostrophic currents, and (d) total current, where the analytical profile is a simple sum of the geostrophic [in (b)] and ageostrophic [in (c)] components.
FIG. A2.(a) Gradient Richardson number as a function of depth and time, and (b) number of inertial periods spent in the initial, dynamically stable phase of evolution as a function of dynamical forcing and thermal driving for the ice shelf-ocean boundary currents illustrated in Fig. A1, with background viscosity n b and diffusivity K b .In (a), unshaded areas are where the thermal driving remains within 0.5% of the initial condition, the horizontal white line indicates the outer edge of the frictional boundary layer, defined as z=d b E 5 p, and black solid, dashed, and dotted lines highlight contours of Ri 5 1, 0.25, and 0.1825, respectively.In (b), the white asterisks indicate the parameters used for Fig. 3 (top panel) and Fig. 2 (bottom panel), and the white line in the top panel indicates slope and thermal driving combinations required for Ri min bl 5 1.Only in the region above and to the right of that line is the third phase of boundary current evolution attainable.All currents will eventually transition into phase two, but the time scale becomes very long, in excess of 60 inertial periods (;1 month) in the yellow shaded regions in both panels of (b).

TABLE 1 .
Symbols and physical constants.

TABLE 2 .
Geostrophic drag parameters for the fully developed ice shelf-ocean boundary current (based on HTC solutions at 30 inertial periods).