2 Formulation of the problem and solution techniques

Figure 1.

Schematic of the simulation domain. The ice face of length $L$ is in contact with the seawater beneath at initial temperature $T_{w}$ and salinity $S_{w}$ . At the bottom right, the domain has an open boundary condition using a sponge layer. Interface conditions at the ice ( $T_{int}$ ; $S_{int}$ and dissolution velocity $V$ ) are evaluated from heat and salt flux balances at that boundary.

The flow field is solved in a rectangular domain shown in figure 1 with length $L$ parallel to the slope, depth $W$ normal to the ice face and a width $D$ in the spanwise direction (normal to the plane of the schematic). Ice–water interface conditions are applied at one boundary (the ice face) of the computational domain (figure 1). The domain and coordinate system are rotated relative to gravity in order to represent the ice slope. Gravity is always directed downward. The flow field is represented by $\tilde{\boldsymbol{u}}=[u_{\unicode[STIX]{x1D702}},v,u_{\unicode[STIX]{x1D701}}]$ , where the wall-normal ( $\unicode[STIX]{x1D702}$ ), spanwise ( $y$ ) and slope-parallel ( $\unicode[STIX]{x1D701}$ ) directions are $u_{\unicode[STIX]{x1D702}}$ , $v$ and $u_{\unicode[STIX]{x1D701}}$ , respectively. The coordinates and velocities are relative to a reference frame fixed at the planar ice–water interface. This is the most convenient reference frame and we make no assumption about the relative speeds of the glacier advance and ablation. We solve the incompressible continuity, Navier–Stokes, heat and salt equations:

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\tilde{\boldsymbol{u}}=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D702}}}{\unicode[STIX]{x2202}t}+(\tilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})u_{\unicode[STIX]{x1D702}}=-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}u_{\unicode[STIX]{x1D702}}+\frac{\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x1D70C}_{0}}g\cos \unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D701}}}{\unicode[STIX]{x2202}t}+(\tilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})u_{\unicode[STIX]{x1D701}}=-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}u_{\unicode[STIX]{x1D701}}-\frac{\unicode[STIX]{x1D70C}^{\ast }}{\unicode[STIX]{x1D70C}_{0}}g\sin \unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}v}{\unicode[STIX]{x2202}t}+(\tilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})v=-\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\unicode[STIX]{x2202}p^{\ast }}{\unicode[STIX]{x2202}y}+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{2}v, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}T^{\ast }}{\unicode[STIX]{x2202}t}+(\tilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})T^{\ast }=\unicode[STIX]{x1D705}_{T}\unicode[STIX]{x1D6FB}^{2}T^{\ast }, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}S^{\ast }}{\unicode[STIX]{x2202}t}+(\tilde{\boldsymbol{u}}\boldsymbol{\cdot }\unicode[STIX]{x1D735})S^{\ast }=\unicode[STIX]{x1D705}_{S}\unicode[STIX]{x1D6FB}^{2}S^{\ast }. & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D70C}_{0}$ is the reference density for pure water at $0\,^{\circ }$ C and $p^{\ast },T^{\ast },S^{\ast }$ and $\unicode[STIX]{x1D70C}^{\ast }$ denote the deviation from the ambient hydrostatic pressure $(p_{w})$ , temperature $(T_{w})$ , salinity $(S_{w})$ and density $(\unicode[STIX]{x1D70C}_{w})$ . The saline water has kinematic viscosity $\unicode[STIX]{x1D708}$ , thermal diffusivity $\unicode[STIX]{x1D705}_{T}$ and salinity diffusivity $\unicode[STIX]{x1D705}_{S}$ . As the flow involves only a small range of temperatures, the equation of state is closely approximated as linear without significant effects on the solution:

(2.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}^{\ast }=\unicode[STIX]{x1D70C}_{0}(\unicode[STIX]{x1D6FD}S^{\ast }-\unicode[STIX]{x1D6FC}T^{\ast }), & & \displaystyle\end{eqnarray}$$

with coefficient of thermal expansion $\unicode[STIX]{x1D6FC}$ and coefficient of haline contraction $\unicode[STIX]{x1D6FD}$ . Computational domains of a given length but different slopes have different vertical heights. Thus we compare the flow and melt rate for different slopes at two given values of the global Grashof number ( $Gr$ , which is the relative strength of buoyancy to viscous force) and one value of the Stefan number ( $St$ ),

(2.8a,b ) $$\begin{eqnarray}\displaystyle Gr\equiv \frac{g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}SL^{3}}{\unicode[STIX]{x1D708}^{2}},\quad St\equiv \frac{\unicode[STIX]{x1D70C}_{s}L_{f}}{\unicode[STIX]{x1D70C}_{w}c_{w}(T_{w}-T_{int})}, & & \displaystyle\end{eqnarray}$$

where, [ $S_{int},T_{int}$ ] are the interface and [ $S_{w},T_{w}$ ] are the far-field salinity and temperature, $\unicode[STIX]{x0394}S=(S_{w}-S_{int})$ is the salinity anomaly, $c_{w}$ is the specific heat, $L_{f}$ is the latent heat of fusion for ice and $\unicode[STIX]{x1D70C}_{s}$ is the density of the ice. The values of $Gr$ and $St$ are independent of slope because the simulations in § 4 show that the interface temperature and salinity are independent of slope. The Prandtl number ( $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}_{T}$ ) and Schmidt number ( $Sc=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}_{S}$ ) are fixed.

Three relations are applied at the ice–water interface. The freezing point of saline water is closely approximated by a linear function of salinity and pressure:

(2.9) $$\begin{eqnarray}\displaystyle T_{int}=a_{s}S_{int}+bP_{int}\simeq a_{s}S_{int}. & & \displaystyle\end{eqnarray}$$

For the present study the effect on the freezing point of hydrostatic pressure difference within the limited domain size is negligible and the interface temperature is assumed to be solely dependent on interface salinity. The slope of the liquidus line is fixed at a _s = -6 × 10^-2 °C ‰^-1 (Holland & Jenkins Reference Holland and Jenkins1999).

The second interface relation expresses the balance between latent heat flux $Q_{m}^{H}$ of melting and the divergence of conductive heat fluxes at the interface,

(2.10) $$\begin{eqnarray}\displaystyle Q_{ice}^{H}-Q_{w}^{H}=Q_{m}^{H}, & & \displaystyle\end{eqnarray}$$

where $Q_{ice}^{H}$ and $Q_{w}^{H}$ are the heat fluxes to the interface in the ice and water, respectively. The conductive transfer of heat into the ice is very small compared to the total heat being used to melt the ice shelf under typical Antarctic conditions (Holland & Jenkins Reference Holland and Jenkins1999; Kerr & McConnochie Reference Kerr and Mcconnochie2015). We therefore neglect the diffusive heat flux into the ice (setting $Q_{ice}^{H}\sim 0$ ). It is also assumed that the diffusion of heat in water at the interface is much faster than the advection of heat by the ablation velocity (i.e. $\unicode[STIX]{x1D705}_{T}\unicode[STIX]{x2202}^{2}T/\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}\gg V\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}\unicode[STIX]{x1D702}|_{\unicode[STIX]{x1D702}=0}$ , as $V\ll \unicode[STIX]{x1D705}_{T}/\unicode[STIX]{x1D6FF}_{T}$ , where $\unicode[STIX]{x1D6FF}_{T}$ is the diffusive thermal boundary layer thickness), so that (2.10) becomes

(2.11) $$\begin{eqnarray}\displaystyle \left.\unicode[STIX]{x1D70C}_{w}c_{w}\unicode[STIX]{x1D705}_{T}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\right|_{\unicode[STIX]{x1D702}=0}=\unicode[STIX]{x1D70C}_{s}VL_{f}, & & \displaystyle\end{eqnarray}$$

where $V$ is the ablation velocity. We adopt the convention that positive ablation velocity indicates melting and retreat of the interface in the negative $\unicode[STIX]{x1D702}$ -direction (at speed $V$ ) relative to the ice mass, or equivalently, translation of the ice mass in the positive $\unicode[STIX]{x1D702}$ -direction relative to the interface reference frame used here.

An analogous equation is used to describe the salt flux balance at the interface due to fresh water release and salt flux divergence:

(2.12) $$\begin{eqnarray}\displaystyle Q_{ice}^{S}-Q_{w}^{S}=Q_{m}^{S}. & & \displaystyle\end{eqnarray}$$

Here $Q_{ice}^{S}$ and $Q_{w}^{S}$ are the diffusive salt fluxes at the interface on the ice and water sides, respectively, and $Q_{m}^{S}$ is the advective fresh water flux associated with melting. The latter can be expressed as $Q_{m}^{S}=\unicode[STIX]{x1D70C}_{s}V(S_{ice}-S_{int})$ . The diffusive salt flux in the ice is neglected and (2.12) becomes

(2.13) $$\begin{eqnarray}\displaystyle \left.\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D705}_{S}\frac{\unicode[STIX]{x2202}S}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\right|_{\unicode[STIX]{x1D702}=0}=\unicode[STIX]{x1D70C}_{s}V(S_{int}-S_{ice}). & & \displaystyle\end{eqnarray}$$

Thus there is a balance between the diffusive flux of salt toward the ice and the advective flux of salt away from the ice face (due to the fresh water release by melting). This is in contrast to the neglect of the heat advection as discussed above. For the solutions found here we will assume that the ice is free of salt ( $S_{ice}=0$ ).

In the DNS we impose wall-normal velocity $u_{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D70C}_{S}V/\unicode[STIX]{x1D70C}_{w}$ (Wells & Worster Reference Wells and Worster2011) at the ice face. Effects of any mean gradient $\text{d}u_{\unicode[STIX]{x1D702}}/\text{d}\unicode[STIX]{x1D701}|_{\unicode[STIX]{x1D702}=0}$ (which occurs in the laminar boundary layer cases, where $V$ is dependent on $\unicode[STIX]{x1D701}$ ) are neglected as justified by Carey & Gebhart (Reference Carey and Gebhart1982). Also neglected are the effects of spatial (in $\unicode[STIX]{x1D701}$ and $y$ ) and temporal variations of $u_{\unicode[STIX]{x1D701}}|_{\unicode[STIX]{x1D702}=0}$ associated with the flow fluctuations. It was verified that the DNS solutions are unchanged even when $u_{\unicode[STIX]{x1D702}}|_{\unicode[STIX]{x1D702}=0}$ is everywhere set to zero.

The ‘open ocean’ side of the computational domain is maintained as an open boundary by relaxing temperature and salinity back to its background temperature $T_{w}$ and salinity $S_{w}$ , respectively, through a ‘sponge’ region (Gayen & Sarkar Reference Gayen and Sarkar2011) at $0.5W\leqslant \unicode[STIX]{x1D702}\leqslant W$ . The along-slope, spanwise velocities and scalar fields are relaxed towards the background state in the sponge region by adding damping functions $-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D702})u_{i}(\unicode[STIX]{x1D702},y,\unicode[STIX]{x1D701},t)$ ( $i=2,3$ ), $-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D702})T^{\ast }(\unicode[STIX]{x1D702},y,\unicode[STIX]{x1D701},t)$ and $-\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D702})S^{\ast }(\unicode[STIX]{x1D702},y,\unicode[STIX]{x1D701},t)$ to the right-hand side of the momentum and scalar equations, respectively, where $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D702})$ is based on the time step $\unicode[STIX]{x0394}t$ and changes from 0 at $\unicode[STIX]{x1D702}=0.5W$ to $1/\unicode[STIX]{x0394}t$   $\text{s}^{-1}$ at $\unicode[STIX]{x1D702}=W$ . At the down-slope and up-slope boundaries of the domain, no-slip conditions are imposed for velocities and no-flux conditions are maintained for the temperature and salinity. Both $T_{int}$ and $S_{int}$ vary with location and time on the ice interface due to turbulent variations in the local heat and solute transport to the interface.

The solution is obtained using a mixed spectral/finite difference algorithm (Gayen et al. Reference Gayen, Griffiths and Kerr2016). The wall-normal and slope-parallel spatial derivatives ( $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D701}$ ) are computed with second-order finite difference. The spanwise ( $y$ ) direction is considered periodic and derivatives in this direction are treated with a pseudo-spectral method. Time stepping is accomplished with a mixed implicit/explicit strategy with all terms involving viscous contribution being stepped with the alternating direction implicit (ADI) method. All the other terms are treated with a low storage third-order Runge–Kutta method (Gayen Reference Gayen2012).

The physical dimension of the rectangular domain for one set of simulations is $W=0.4~\text{m}$ , $D=0.05~\text{m}$ and $L=1.8~\text{m}$ . Additional simulations with the across-slope width doubled to $D=0.1~\text{m}$ , and the same $W$ and $L$ , show similar boundary layer properties and melt rates. Thus a width of $D=0.05~\text{m}$ is used for the remainder of the simulations for the sake of computational efficiency. A second set of solutions in a larger domain used $W=4~\text{m}$ , $D=0.05~\text{m}$ and $L=20~\text{m}$ . In the two domain lengths the grids have $256\times 64\times 1150$ and $256\times 64\times 1920$ points in the $\unicode[STIX]{x1D702}$ , $y$ and $\unicode[STIX]{x1D701}$ directions, respectively. In order to resolve the salinity boundary layer and turbulent microscales for salinity, grid stretching is used in the $\unicode[STIX]{x1D702}$ direction. All solutions satisfy rigorous grid resolution and grid convergence criteria that are reported in detail in Gayen, Griffiths & Hughes (Reference Gayen, Griffiths and Hughes2014), Gayen et al. (Reference Gayen, Griffiths and Kerr2016). Variable time stepping with a fixed Courant–Friedrichs–Lewy (CFL) number of $N_{CFL}=0.5$ is used. The time step is calculated to be $\unicode[STIX]{x0394}t=(1/N_{CFL})[\unicode[STIX]{x0394}\unicode[STIX]{x1D702}/u_{\unicode[STIX]{x1D702}},\unicode[STIX]{x0394}y/v,\unicode[STIX]{x0394}\unicode[STIX]{x1D701}/u_{\unicode[STIX]{x1D701}}]_{min}$ and varies significantly, from $O(10^{-3})~\text{s}$ in the turbulent regime to $O(10^{-2})~\text{s}$ in the laminar regime. The fractional step method is used to evaluate dynamic pressure at each time step (Gayen & Sarkar Reference Gayen and Sarkar2011).

The slope angle $\unicode[STIX]{x1D703}$ is varied from $2^{\circ }$ to $90^{\circ }$ . In order to focus on the effect of slope on the melting process, the far-field temperature $T_{w}=2.3^{\circ }\text{C}$ and salinity S _w = 35 ‰ are fixed for all cases. All temperatures are measured with respect to the freezing point of pure water and therefore are quoted in Celsius. We fix $g=10~\text{m}~\text{s}^{-2}$ , $\unicode[STIX]{x1D705}_{T}=1.285\times 10^{-7}~\text{m}^{2}~\text{s}^{-1}$ , $c_{w}=3985~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$ , $\unicode[STIX]{x1D708}=1.8\times 10^{-6}~\text{m}^{2}~\text{s}^{-1}$ , $\unicode[STIX]{x1D6FC}=6\times 10^{-5}~\text{K}^{-1}$ and 𝛽 = 8 × 10^-4 ‰^-1, taken from the physical properties of aqueous NaCl solutions at the far-field temperature and salinity (Washburn Reference Washburn1926; Weast, Astle & Beyer Reference Weast, Astle and Beyer1989). We use $\unicode[STIX]{x1D705}_{S}=7.2\times 10^{-10}~\text{m}^{2}~\text{ s}^{-1}$ at $0^{\circ }\text{C}$  (Josberger & Martin Reference Josberger and Martin1981) with a resultant $Sc=2500$ at grid points adjacent to the interface. In the interior we use $\unicode[STIX]{x1D705}_{S}=3.6\times 10^{-9}~\text{m}^{2}~\text{s}^{-1}$ ( $Sc=500$ , Gayen et al. Reference Gayen, Griffiths and Kerr2016) in order to make the turbulence simulations feasible. The Prandtl number is fixed at $Pr=15$ throughout the domain. The primary solution sets are designed to cover the range of slope angles for both $Gr=7.5\times 10^{11}$ (for a domain length of 1.8 m) and $Gr=10.28\times 10^{14}$ (domain length of 20 m).

The critical Grashof number for transitions from laminar to turbulent flow on a vertical wall varies from $Gr_{c}^{\bot }\sim 10^{9}{-}10^{10}$ (Turner Reference Turner1979; Josberger & Martin Reference Josberger and Martin1981; Holman Reference Holman2010). A first approximation for the critical along-slope distance $L_{c}$ beyond which the boundary layer becomes turbulent, based on the critical vertical height $L_{c}^{\bot }$ , is

(2.14) $$\begin{eqnarray}\displaystyle L_{c}\approx \frac{L_{c}^{\bot }}{\sin \unicode[STIX]{x1D703}}=\frac{1}{\sin \unicode[STIX]{x1D703}}\left(\frac{\unicode[STIX]{x1D708}^{2}Gr_{c}^{\bot }}{g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S}\right)^{1/3}, & & \displaystyle\end{eqnarray}$$

giving a critical Grashof number for the sloping interface as

(2.15) $$\begin{eqnarray}\displaystyle Gr_{c}\approx \frac{Gr_{c}^{\bot }}{(\sin \unicode[STIX]{x1D703})^{3}}. & & \displaystyle\end{eqnarray}$$

Based on this assumption we calculated the transition length for various slope angles bounding the previously established critical $Gr_{c}^{\bot }$ (figure 2) and plotted the laminar to turbulent transition length ( $L_{c}$ ) given by the DNS, where $L_{c}$ is defined as the height where spanwise fluctuations ( $v_{rms}$ , $rms$   $=$  root mean square) reach 10 % or more of the average up-slope flow. The results show that the along-slope length required to achieve turbulence increases rapidly with decreasing slope angles. Hence we require the longer domain of $20$ m in order to simulate turbulent conditions under sloping ice at $\unicode[STIX]{x1D703}\leqslant 20^{\circ }$ .

3 Scaling analysis

Ablation of the ice face takes place when it is in contact with saline water with a temperature greater than the melting temperature. The addition of melt water to the saline water results in freshening of a layer adjacent to the ice and gives rise to a buoyant boundary layer plume. Working from the momentum and solute transport equations the plume can be divided into multiple layers (figure 3) based on the dominant balances discussed below. The buoyancy force is important only within a solutal boundary layer at the ice face, of thickness $\unicode[STIX]{x1D6FF}_{S}$ , and we refer to this as the ‘inner layer’. Diffusion or turbulent advection of heat in the wall-normal direction is important within a thermal boundary layer of thickness $\unicode[STIX]{x1D6FF}_{T}$ . However, the thermal buoyancy is relatively small and does not enter the momentum balance at leading order. In the inner layer, solutal buoyancy and stress are the dominant terms in the momentum equation. Outside the inner boundary layer there is an outer layer, of thickness $\unicode[STIX]{x1D6FF}_{0}$ , in which buoyancy forces make a negligible contribution.

When the flow is laminar the whole of the inner layer will be dominated by the balance of buoyancy with molecular viscous stress, and by molecular diffusion of heat and solute toward the ice, as previously described in the context of a vertical ice face (Wells & Worster Reference Wells and Worster2011). The laminar outer layer is governed by an inertia–viscous balance. That analysis is adapted here (§ 3.1) to a sloping ice face. We find that this laminar solution may have geophysical relevance for very small ice slopes, as a consequence of a stabilising density stratification produced by the melting.

Figure 2.

Simulated and theoretically predicted (2.14) laminar to turbulent transition length, where the predictions (curves) are based on the critical Grashof number for a vertical wall ( $Gr_{c}^{\bot }=10^{9}{-}10^{10}$ ). The simulated transition length $L_{t}$ is calculated based on turbulent statistics from the simulations for $L=1.8~\text{m}$ and 20 m, respectively.

Figure 3.

Schematic of the theoretical model boundary layer structure, with up-slope velocity and solute concentration profiles in the wall-normal direction, for the case of a turbulent plume (thicknesses not to scale).

When the plume is turbulent, small-scale eddies add significantly to momentum transport in the wall-normal direction within the buoyant inner layer, leaving a much thinner laminar and diffusive sublayer (of thickness $\unicode[STIX]{x1D6FF}_{sub}$ ) on the ice face. Within the sub-layer eddy transport and turbulent kinetic energy production are negligible. The smallness of the sub-layer thickness is confirmed by the DNS solutions reported here, which show it to be $O(10)$ times smaller than the inner layer thickness $\unicode[STIX]{x1D6FF}_{S}$ . In the inner layer the dominant momentum balance is between the solutal buoyancy and Reynolds stress. However, scaling indicates that the molecular diffusion of solute toward the ice remains a dominant term in the salt budget. In the outer layer the turbulent plume is dominated by Reynolds stress and turbulent transport of solute, leaving the solute relatively well mixed, which in turn implies a negligible molecular transport of solute. Molecular diffusion of heat can remain important throughout a region of thickness $\unicode[STIX]{x1D6FF}_{T}$ , which may overlap the outer inertial layer. Ambient fluid is entrained into the turbulent plume from the edge of the outer layer. The layers defined here are consistent with a model previously proposed for thermal convection at a heated vertical boundary (Wells & Worster Reference Wells and Worster2008). In the following sub-sections the thicknesses of the sloping inner solutal layer and thermal boundary layer, and the ablation rate, are analysed for laminar and turbulent cases.

3.1 Laminar boundary layer flows: viscous–buoyancy balance

Wells & Worster (Reference Wells and Worster2011) derived similarity solutions for melting and dissolving with a laminar boundary layer for a vertical ice face and quiescent far field. (The final form of the similarity solution for the case of dissolution of a solute at a vertical wall was also quoted by Husband & Ozsahin (Reference Husband and Ozsahin1967). Here we summarise their solution for the case of dissolving and extend it to sloping ice faces. However, for brevity we use simple scaling arguments, which also assist in contrasting the dynamics of laminar and turbulent cases.) In the inner layer $\unicode[STIX]{x1D6FF}_{S}$ and ${\mathcal{L}}$ are the characteristic length scales for the wall-normal and the along-slope directions, respectively, whereas $\unicode[STIX]{x1D6FF}_{0}$ and ${\mathcal{L}}$ are the characteristic length scales associated with the outer layer. At small Grashof numbers ( $Gr<Gr_{c}$ ) or shallow slopes, the inertia terms inside the diffusive boundary layer are negligible, leaving a balance between viscous drag and buoyancy in the up-slope flow in the momentum equations (2.2)–(2.4). The buoyancy is predominantly supplied by the salinity anomaly ( $\unicode[STIX]{x0394}S$ ) across the diffusive boundary layer:

(3.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}^{2}U_{\unicode[STIX]{x1D701}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}\sim \unicode[STIX]{x1D708}\frac{U_{\unicode[STIX]{x1D701}}}{\unicode[STIX]{x1D6FF}_{S}^{2}}\sim g\sin \unicode[STIX]{x1D703}\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S. & & \displaystyle\end{eqnarray}$$

The scaling of the viscous terms in (3.1) holds for large $Sc$ , as the inflection in the up-slope velocity profile must lie outside the buoyancy-driven inner layer. The wall-normal derivatives in the boundary layer are the dominant contributions (i.e. $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D702}\gg \unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}$ as $\unicode[STIX]{x1D6FF}_{S}\ll {\mathcal{L}}$ ) to both the viscous drag and the diffusion of solute. Outside the solutal diffusive boundary layer the buoyancy force is neglected and hence the inertia in the outer layer balances the residual viscous drag force generated by the up-slope flow (Nilson Reference Nilson1985). This leads to

(3.2) $$\begin{eqnarray}\displaystyle \frac{U_{\unicode[STIX]{x1D701}}^{2}}{{\mathcal{L}}}\sim \unicode[STIX]{x1D708}\frac{U_{\unicode[STIX]{x1D701}}}{\unicode[STIX]{x1D6FF}_{0}^{2}}. & & \displaystyle\end{eqnarray}$$

The convective salt transport by mean flow in (2.6) balances solute diffusion through the inner boundary layer:

(3.3) $$\begin{eqnarray}\displaystyle U_{\unicode[STIX]{x1D701}}\frac{\unicode[STIX]{x0394}S}{{\mathcal{L}}}\sim U_{\unicode[STIX]{x1D702}}\frac{\unicode[STIX]{x0394}S}{\unicode[STIX]{x1D6FF}_{S}}\sim \unicode[STIX]{x1D705}_{S}\frac{\unicode[STIX]{x0394}S}{\unicode[STIX]{x1D6FF}_{S}^{2}}. & & \displaystyle\end{eqnarray}$$

Here total velocity $u_{i}=U_{i}+u_{i}^{\prime }$ is decomposed into the mean boundary layer flow $U_{i}$ and fluctuating velocity component $u_{i}^{\prime }$ with the latter assumed to play a negligible role for laminar flow. We also assume that the salinity anomaly in the along-slope direction is equivalent to the salinity difference across the boundary layer. Combining (3.1) and (3.3) leads to the scaling of the thickness of the solutal boundary layer as

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{S}\sim \unicode[STIX]{x1D6FF}_{S0}(\sin \unicode[STIX]{x1D703})^{-(1/4)}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{S0}$ is the salinity boundary layer thickness for the case of a vertical ice face

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{S0}\sim \left(\frac{\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}_{S}{\mathcal{L}}}{g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S}\right)^{1/4}, & & \displaystyle\end{eqnarray}$$

or equivalently,

(3.6) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6FF}_{S0}}{{\mathcal{L}}}\sim Gr_{{\mathcal{L}}}^{-1/4}Sc^{-1/4}, & & \displaystyle\end{eqnarray}$$

with the local Grashof number $Gr_{{\mathcal{L}}}=g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S{\mathcal{L}}^{3}/\unicode[STIX]{x1D708}^{2}$ (which is based on the up-slope distance ${\mathcal{L}}$ ). From (2.13) the ablation rate can be estimated as $V\sim (\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x1D705}_{S}\unicode[STIX]{x0394}S)/(\unicode[STIX]{x1D70C}_{S}S_{int}\unicode[STIX]{x1D6FF}_{S})$ , which gives

(3.7) $$\begin{eqnarray}\displaystyle V\sim V_{0}(\sin \unicode[STIX]{x1D703})^{1/4}, & & \displaystyle\end{eqnarray}$$

where $V_{0}$ is the ablation rate for a vertical ice face given by

(3.8) $$\begin{eqnarray}\displaystyle V_{0}\sim \frac{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x0394}S}{\unicode[STIX]{x1D70C}_{S}S_{int}}\left(\frac{g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S\unicode[STIX]{x1D705}_{S}^{3}}{\unicode[STIX]{x1D708}}\right)^{1/4}{\mathcal{L}}^{-1/4}, & & \displaystyle\end{eqnarray}$$

or

(3.9) $$\begin{eqnarray}\displaystyle \frac{V_{0}{\mathcal{L}}}{\unicode[STIX]{x1D705}_{S}}\sim \left(\frac{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x0394}S}{\unicode[STIX]{x1D70C}_{S}S_{int}}\right)Gr_{{\mathcal{L}}}^{1/4}. & & \displaystyle\end{eqnarray}$$

This scaling predicts that the ablation rate decreases in the up-slope direction as $\unicode[STIX]{x1D701}^{-1/4}$ , consistent with dissolution by laminar flow on a vertical ice face ( $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/2$ ) found by Wells & Worster (Reference Wells and Worster2011). It also predicts that ablation rate decreases with slope angle as $(\sin \unicode[STIX]{x1D703})^{1/4}$ . Considering the diffusive heat transport in the thermal boundary layer, the ablation rate from (2.11) must satisfy

(3.10) $$\begin{eqnarray}\displaystyle V\sim -\frac{\unicode[STIX]{x1D70C}_{w}c_{w}\unicode[STIX]{x1D705}_{T}\unicode[STIX]{x0394}T}{\unicode[STIX]{x1D70C}_{S}L_{f}\unicode[STIX]{x1D6FF}_{T}}. & & \displaystyle\end{eqnarray}$$

Equating (3.10) and (3.7), the thermal boundary layer thickness $\unicode[STIX]{x1D6FF}_{T}$ can then be expressed as

(3.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{T}\sim \unicode[STIX]{x1D6FF}_{T0}(\sin \unicode[STIX]{x1D703})^{-(1/4)}, & & \displaystyle\end{eqnarray}$$

where the thermal boundary layer thickness $\unicode[STIX]{x1D6FF}_{T0}$ for the case of a vertical ice face is

(3.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{T0}\sim \frac{c_{w}\unicode[STIX]{x1D705}_{T}S_{int}}{\unicode[STIX]{x1D705}_{S}L_{f}}\frac{\unicode[STIX]{x0394}T}{\unicode[STIX]{x0394}S}\left(\frac{\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}_{S}{\mathcal{L}}}{g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S}\right)^{1/4}, & & \displaystyle\end{eqnarray}$$

or in dimensionless form

(3.13) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6FF}_{T0}}{{\mathcal{L}}}\sim \left(\frac{\unicode[STIX]{x1D70C}_{S}}{\unicode[STIX]{x1D70C}_{w}}\frac{S_{int}}{\unicode[STIX]{x0394}S}\right)StLeGr_{{\mathcal{L}}}^{-1/4}Sc^{-1/4}. & & \displaystyle\end{eqnarray}$$

Here $Le=\unicode[STIX]{x1D705}_{T}/\unicode[STIX]{x1D705}_{S}$ is the Lewis number. The thermal boundary layer thickness shows a similar dependence on the up-slope distance and slope angle as does the solutal boundary layer (3.5), (3.6). Scaling of the outer inertial layer thickness $\unicode[STIX]{x1D6FF}_{0}$ (using (3.2) and (3.5)) becomes

(3.14) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6FF}_{0}}{{\mathcal{L}}}\sim Sc^{1/2}\frac{\unicode[STIX]{x1D6FF}_{S}}{{\mathcal{L}}}\sim Gr_{{\mathcal{L}}}^{-1/4}Sc^{1/4}, & & \displaystyle\end{eqnarray}$$

and this too increases with up-slope distance in the same manner as the inner layer thickness.

3.2 Turbulent boundary layer flows: boundary layer inertia–buoyancy balance

At large $Gr$ the boundary layer becomes unstable to both buoyancy- and shear-driven instabilities, leading to small-scale motions (Josberger & Martin Reference Josberger and Martin1981; Holman Reference Holman2010; Gayen et al. Reference Gayen, Griffiths and Kerr2016). In the outer layer diffusive transport is negligible and there is a balance in (2.6) between advection of solute by the mean flow and turbulent solute transport, leading to

(3.15) $$\begin{eqnarray}\displaystyle \frac{U_{\unicode[STIX]{x1D701}}\unicode[STIX]{x0394}S}{{\mathcal{L}}}\sim \frac{\overline{u^{\prime }S^{\prime }}}{\unicode[STIX]{x1D6FF}_{0}}. & & \displaystyle\end{eqnarray}$$

In the inner layer there is potentially a regime, at intermediate $Gr$ , in which turbulent transport of solute dominates over molecular transport while viscous stress remains important relative to Reynolds stress. However, here we consider a regime at very large $Gr$ , in which Reynolds stresses $\overline{u_{i}^{\prime }u_{j}^{\prime }}$ produced by the small-scale motions dominate over the viscous stress and establish a dominant balance with local buoyancy in the inner layer. Thus we assume

(3.16) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\overline{u_{\unicode[STIX]{x1D702}}^{\prime }u_{\unicode[STIX]{x1D701}}^{\prime }}\sim \frac{u_{\unicode[STIX]{x1D702}}^{\prime 2}}{\unicode[STIX]{x1D6FF}_{S}}\sim g\sin \unicode[STIX]{x1D703}\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S & & \displaystyle\end{eqnarray}$$

along with local isotropy ( $u_{\unicode[STIX]{x1D702}}^{\prime }\sim u_{\unicode[STIX]{x1D701}}^{\prime }\sim v^{\prime }$ ).

For this turbulent case small-scale fluctuations contribute to solute transport in the inner layer and the mean convective transport in the wall-normal direction becomes relatively small, as indicated by the ratio of terms $(\overline{u^{\prime }S^{\prime }}/\unicode[STIX]{x1D6FF}_{S})/(U_{\unicode[STIX]{x1D701}}\unicode[STIX]{x0394}S/{\mathcal{L}})\sim \unicode[STIX]{x1D6FF}_{0}/\unicode[STIX]{x1D6FF}_{S}\gg 1$ (using (3.15)). This is consistent with previous studies of natural thermal convection at a heated vertical boundary (George & Capp Reference George and Capp1979; Tsuji & Nagano Reference Tsuji and Nagano1988; Wells & Worster Reference Wells and Worster2008). The solute transport by fluctuations balances, to leading order, the diffusive transport of salt giving

(3.17) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\overline{u_{\unicode[STIX]{x1D702}}^{\prime }S^{\prime }}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}+\frac{\unicode[STIX]{x2202}\overline{u_{\unicode[STIX]{x1D701}}^{\prime }S^{\prime }}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}\sim \unicode[STIX]{x1D705}_{S}\frac{\unicode[STIX]{x2202}^{2}S}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}^{2}}, & & \displaystyle\end{eqnarray}$$

or in scaled form

(3.18) $$\begin{eqnarray}\displaystyle \frac{\overline{u_{\unicode[STIX]{x1D702}}^{\prime }S^{\prime }}}{\unicode[STIX]{x1D6FF}_{S}}\sim \unicode[STIX]{x1D705}_{S}\frac{\unicode[STIX]{x0394}S}{\unicode[STIX]{x1D6FF}_{S}^{2}}, & & \displaystyle\end{eqnarray}$$

where $S^{\prime }$ is the salinity fluctuation and $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D702}\gg \unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D701}$ . For a vertical ice face the fluctuations in the salinity field $S^{\prime }$ and density field $\unicode[STIX]{x1D70C}^{\prime }$ scale with $\unicode[STIX]{x0394}S$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ and (3.16), (3.18) lead to the turbulent inner layer thickness and ablation rate:

(3.19a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{S0}\sim \left(\frac{\unicode[STIX]{x1D705}_{S}^{2}}{g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S}\right)^{1/3}\quad \text{and}\quad V_{0}\sim \frac{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x0394}S}{\unicode[STIX]{x1D70C}_{s}S_{int}}(g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S\unicode[STIX]{x1D705}_{S})^{1/3}. & & \displaystyle\end{eqnarray}$$

Kerr & McConnochie (Reference Kerr and Mcconnochie2015) reported scaling for the ablation velocity in this vertical case having this same dependence on buoyancy, $V_{0}\sim (g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S)^{1/3}$ , but different molecular dependences, $V_{0}\sim (\unicode[STIX]{x1D705}_{S}^{2}/\unicode[STIX]{x1D708})^{1/3}$ . The present theory neglects the effect of viscosity, which is assumed small compared with the Reynolds stress in the inner layer.

When the ice face is inclined, the fresh water flux due to ablation generates stratification, which inhibits the turbulent fluctuations. Hence, salinity fluctuations $S^{\prime }$ cannot directly scale with salinity anomaly $\unicode[STIX]{x0394}S$ across the boundary layer. Another independent equation is required to solve for $S^{\prime }$ . We assume that the mean density gradient in the inner stratified layer scales with $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D6FF}_{S}$ and the frequency of turbulent fluctuations (the eddy turnover rate $1/\unicode[STIX]{x0394}t$ ) scales with the buoyancy frequency $N$ . Hence $\unicode[STIX]{x0394}t\sim 1/N\sim [\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D6FF}_{S}/g\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}]^{1/2}$ , where $N^{2}=-(g/\unicode[STIX]{x1D70C}_{0})\text{d}\unicode[STIX]{x1D70C}/\text{d}z\sim (g/\unicode[STIX]{x1D70C}_{0})\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}/\unicode[STIX]{x1D6FF}_{S}$ . This leads to a simple linearised equation for density fluctuation:

(3.20) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}t}\sim u_{\unicode[STIX]{x1D702}}^{\prime }\frac{\unicode[STIX]{x2202}\overline{\unicode[STIX]{x1D70C}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}}\quad \rightarrow \quad \frac{\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x0394}t}\sim \frac{u_{\unicode[STIX]{x1D702}}^{\prime }\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D6FF}_{S}}. & & \displaystyle\end{eqnarray}$$

These assumptions may not hold for near-vertical ice faces, where the effect of local stratification becomes negligible. Substituting the time scale in (3.20), the density and salinity fluctuations ( $S^{\prime }\sim \unicode[STIX]{x1D70C}^{\prime }/\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D70C}_{0}$ ) become

(3.21a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}^{\prime }\sim {u_{\unicode[STIX]{x1D702}}}^{\prime }\sqrt{\frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D70C}_{0}}{g\unicode[STIX]{x1D6FF}_{S}}},\quad S^{\prime }\sim {u_{\unicode[STIX]{x1D702}}}^{\prime }\sqrt{\frac{\unicode[STIX]{x0394}S}{g\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}_{S}}}. & & \displaystyle\end{eqnarray}$$

From (3.16), (3.18) and (3.21) the salinity fluctuations scale as

(3.22) $$\begin{eqnarray}\displaystyle S^{\prime }\sim \unicode[STIX]{x0394}S(\sin \unicode[STIX]{x1D703})^{1/2} & & \displaystyle\end{eqnarray}$$

and the solutal boundary layer thickness becomes

(3.23) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{S}\sim \unicode[STIX]{x1D6FF}_{S0}(\sin \unicode[STIX]{x1D703})^{-2/3}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{S0}$ is again the inner solutal layer thickness for the case of a vertical ice face, this time for turbulent convection. In dimensionless form:

(3.24) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6FF}_{S0}}{{\mathcal{L}}}\sim Gr_{{\mathcal{L}}}^{-1/3}Sc^{-2/3}. & & \displaystyle\end{eqnarray}$$

Hence $\unicode[STIX]{x1D6FF}_{S0}$ and $\unicode[STIX]{x1D6FF}_{S}$ are independent of up-slope distance. From (2.13) the ablation velocity becomes

(3.25) $$\begin{eqnarray}\displaystyle V\sim V_{0}(\sin \unicode[STIX]{x1D703})^{2/3}, & & \displaystyle\end{eqnarray}$$

where the ablation velocity (3.19) for a vertical ice face is

(3.26) $$\begin{eqnarray}\displaystyle \frac{V_{0}{\mathcal{L}}}{\unicode[STIX]{x1D705}_{S}}\sim \left(\frac{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x0394}S}{\unicode[STIX]{x1D70C}_{S}S_{int}}\right)Gr_{{\mathcal{L}}}^{1/3}Sc^{2/3}. & & \displaystyle\end{eqnarray}$$

Thus ablation rates too are predicted to be independent of distance up the slope. Ablation decreases with decreasing slope angle as $(\sin \unicode[STIX]{x1D703})^{2/3}$ .

Using (2.11) and (3.25) we solve for the thermal boundary layer thickness

(3.27) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{T}\sim \unicode[STIX]{x1D6FF}_{T0}(\sin \unicode[STIX]{x1D703})^{-2/3}, & & \displaystyle\end{eqnarray}$$

where for the vertical case

(3.28) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{T0}\sim \frac{c_{w}\unicode[STIX]{x1D705}_{T}S_{int}}{\unicode[STIX]{x1D705}_{S}L_{f}}\frac{\unicode[STIX]{x0394}T}{\unicode[STIX]{x0394}S}\left(\frac{\unicode[STIX]{x1D705}_{S}^{2}}{g\unicode[STIX]{x1D6FD}\unicode[STIX]{x0394}S}\right)^{1/3}, & & \displaystyle\end{eqnarray}$$

or

(3.29) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x1D6FF}_{T0}}{{\mathcal{L}}}\sim \left(\frac{\unicode[STIX]{x1D70C}_{S}S_{int}}{\unicode[STIX]{x1D70C}_{w}\unicode[STIX]{x0394}S}\right)St^{-1}LeGr_{{\mathcal{L}}}^{-1/3}Sc^{2/3}. & & \displaystyle\end{eqnarray}$$

Like the inner solutal boundary layer, the inner thermal boundary layer has thickness independent of distance up the slope.

The outer layer scaling can be established independently using turbulent entrainment parameterisation. Continuity (2.1) for the mean flow in the outer layer shows

(3.30) $$\begin{eqnarray}\displaystyle \frac{U_{\unicode[STIX]{x1D701}}}{{\mathcal{L}}}\sim \frac{U_{\unicode[STIX]{x1D702}}}{\unicode[STIX]{x1D6FF}_{0}}. & & \displaystyle\end{eqnarray}$$

The mean normal velocity $U_{\unicode[STIX]{x1D702}}$ is equivalent to the entrainment velocity and is assumed to be linearly proportional to the up-slope velocity, $U_{\unicode[STIX]{x1D702}}\sim EU_{\unicode[STIX]{x1D701}}$ , where $E$ is the entrainment coefficient (Morton et al. Reference Morton, Taylor and Turner1956). This leads to

(3.31) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{0}\sim E{\mathcal{L}}. & & \displaystyle\end{eqnarray}$$

Similar scaling was suggested for the outer layer in the case of natural thermal convection at a heated vertical boundary (Wells & Worster Reference Wells and Worster2008), for turbulent wall plumes driven by a uniformly distributed wall buoyancy flux (Cooper & Hunt Reference Cooper and Hunt2010) and also for the melt boundary layer at a vertical ice wall (Kerr & McConnochie Reference Kerr and Mcconnochie2015; Gayen et al. Reference Gayen, Griffiths and Kerr2016). Under an inclined ice face $E$ is likely to be dependent on the slope angle, potentially following the result of Ellison & Turner (Reference Ellison and Turner1959) for a dense plume flowing down a sloping boundary. For our case of melting of a sloping ice face, the result is an outer layer thickness (the overall region of turbulent flow) that grows linearly with distance up the slope.

5 Discussion and conclusions

This paper presents the first DNS of ice dissolution due to convection under a sloping ice face. We have simulated fully turbulent flow at geophysically relevant slopes ( $\unicode[STIX]{x1D703}\leqslant 40^{\circ }$ ) at temperature and salinities relevant to Antarctic conditions. The typical Grashof number based on the vertical height of the ice–seawater interface is of the order of $Gr^{\bot }\sim 10^{17}{-}10^{19}$ , for heights in the range 200–800 m. Although this range is not achievable using DNS with present computational capacity, we have chosen large Grashof numbers ( $Gr^{\bot }\sim 10^{10}{-}10^{11}$ ), which are well above the critical $Gr_{c}^{\bot }\sim 10^{9}$ for the transition to turbulent convection on a vertical heated wall. The solutions confirm that these conditions ensure a steady but turbulent mean flow. The boundary layer flow changes significantly with the slope. For steep angles a narrow up-slope flow of relatively fresh buoyant water develops close to the wall. For small slope angles ( $\unicode[STIX]{x1D703}<20^{\circ }$ ), the buoyant up-slope flow is relatively weak due to the reduction of buoyancy force in the along-slope direction. However, the solutions at large Grashof numbers again show a turbulent flow.

Boundary layer properties and ablation rates are dependent on whether the boundary layer is laminar or turbulent. Both the thermal and solutal boundary layer thicknesses increase with decreasing slope angle, resulting in the reduction of heat and salt transport to the interface and consequent reduction of the ablation rate. In the laminar case the ablation rate decreases with along-slope distance as $\unicode[STIX]{x1D701}^{-1/4}$ , consistent with theoretical scaling. In contrast the turbulent cases have ablation rates that are nearly uniform along the slope, as previously found for vertical ice faces (Kerr & McConnochie Reference Kerr and Mcconnochie2015; Gayen et al. Reference Gayen, Griffiths and Kerr2016). We hypothesise that for domain lengths larger than those achieved in the present DNS, and having a Grashof number greater than the critical value, dissolution rate will be given by the asymptotic dynamics of turbulent boundary layer flow, and the present results can be extrapolated to geophysical scales following the scaling presented here. In that case, for a purely convective boundary layer, the ablation of ice faces into the ocean will follow $V\simeq (\sin \unicode[STIX]{x1D703})^{2/3}$ . This leads to a simple modification of the previously derived dependence of melt rate on the driving temperature difference (Kerr & McConnochie Reference Kerr and Mcconnochie2015; Gayen et al. Reference Gayen, Griffiths and Kerr2016), now taking account of the interface slope as:

(5.1) $$\begin{eqnarray}\displaystyle V\simeq 8.98\times (\unicode[STIX]{x0394}T_{L})^{4/3}(\sin \unicode[STIX]{x1D703})^{2/3}\;\text{m}~\text{yr}^{-1}. & & \displaystyle\end{eqnarray}$$

An alternative scaling reported by Magorrian & Wells (Reference Magorrian and Wells2016) is based on a buoyant plume model and gives a greater sensitivity of melt rates on $\unicode[STIX]{x1D703}$ for small slopes ( $V\sim (\sin \unicode[STIX]{x1D703})^{3/2}$ ) and an inverse dependence on slope angle ( $V\sim 1/\sin \unicode[STIX]{x1D703}$ ) for near-vertical interfaces. Those trends are not reflected in the present DNS results. The discrepancy may result from a different regime of the convective boundary layer. The model of Jenkins (Reference Jenkins1991) and Magorrian & Wells (Reference Magorrian and Wells2016) assumes a regime in which the thickness of the laminar sub-layer near the ice face is controlled primarily by shear instability (Grossmann & Lohse Reference Grossmann and Lohse2000; Wells & Worster Reference Wells and Worster2008) rather than convective instability as found here. The transition to this shear-dominated regime at steep slopes for saline convection was predicted to occur at $Gr\sim 10^{20}$ , which occurs for vertical ice heights of hundreds of metres. The present study uses heights/domain lengths where the boundary layer on near-vertical interfaces remains controlled by turbulent convection. Hence it remains to be demonstrated whether the transition can occur. Further discussion of transition between these two regimes can be found in McConnochie & Kerr (Reference McConnochie and Kerr2017b ).

The turbulent kinetic energy budget shows the presence of statistically steady turbulence in the simulated flow fields for slope angles as small as, $\unicode[STIX]{x1D703}=5^{\circ }$ . For near-vertical slopes ( $\unicode[STIX]{x1D703}>80^{\circ }$ ) contributions to turbulent kinetic energy from shear production and buoyancy flux are comparable, with a slightly greater contribution from the buoyancy flux. For small slopes the production of turbulent kinetic energy by buoyancy fluxes is significantly smaller than the turbulent shear production. This potentially implies that the shear associated with large-scale ambient geostrophic currents and barotropic tides in the ocean is more likely to contribute to the turbulent transport at the ice face and enhance the melt rate for small slopes.

The present study has focused on the effect of ice face slope on melting that is driven by natural convection. The natural convection on its own can be viewed as a base, or reference, case given that natural convection will always be present irrespective of the magnitude and influence of shear associated with ambient geostrophic currents, internal waves or sub-glacial discharge plumes. The next step will be to include ambient stratification (McConnochie & Kerr Reference McConnochie and Kerr2016b ), subglacial discharge of fresh water (McConnochie & Kerr Reference McConnochie and Kerr2017a ) and ambient shear.