# Turbulent ice-ocean boundary layers in the well-mixed regime:
insights from direct numerical simulations

###### Abstract

The meltwater mixing line (MML) model provides a theoretical prediction of
near-ice water mass properties that is useful to compare with observations. If
oceanographic measurements reported in a temperature-salinity diagram overlap
with the MML prediction, then it is usually concluded that the local dynamics
are dominated by the turbulent mixing of an ambient water mass with nearby
melting ice. While the MML model is consistent with numerous observations, it
is built on an assumption that is difficult to test with field measurements,
especially near the ice boundary, namely that the effective (turbulent and
molecular) salt and temperature diffusivities are equal. In this paper, this
assumption is tested via direct numerical simulations of a canonical model for
externally-forced ice-ocean boundary layers in a uniform ambient. We focus on
the well-mixed regime by considering an ambient temperature close to freezing
and run the simulations until a statistical steady state is reached. The
results validate the assumption of equal effective diffusivities across most
of the boundary layer. Importantly, the validity of the MML model implies a
linear correlation between the mean salinity and temperature profiles normal
to the interface that can be leveraged to construct a reduced ice-ocean
boundary layer model based on a single scalar variable called thermal driving.
We demonstrate that the bulk dynamics predicted by the reduced thermal driving
model are in good agreement with the bulk dynamics predicted by the full
temperature-salinity model. Then, we show how the results from the thermal
driving model can be used to estimate the interfacial heat and salt fluxes,
and the melt rate.

This work has been submitted to the Journal of Physical Oceanography.
Copyright in this work may be transferred without further notice.

We investigate the turbulent dynamics and thermodynamical properties of water
masses below ice shelves using new data from high-resolution simulations. This
is important because observations of ice-ocean boundary layers are currently
too scarce to construct reliable models of ice-shelf melting as functions of
ocean conditions. Our results demonstrate that the turbulent diffusivities of
salt and temperature are approximately equal in the well-mixed regime. This
implies a linear correlation between the mean temperature and salinity
profiles, consistent with the meltwater mixing line prediction that is often
used to interpret polar observations. We take advantage of this correlation to
propose a reduced model of ice-ocean boundary layers that can predict ice-
shelf melt rates at relatively low computational cost.

##  1 Introduction

Melting and freezing processes transform water masses flowing next to ice
shelves, icebergs and sea ice (McDougall et al., 2014; Dinniman et al., 2016;
Cenedese and Straneo, 2023). Ice melting produces meltwater that is typically
colder and fresher than the ambient, such that the seawater temperature and
salinity usually decrease toward the ice-ocean boundary. In the presence of a
turbulent external flow, a well-mixed bulk emerges adjacent to a molecular
diffusive boundary layer (cf. figure 1(a)). Multiple oceanographic
observations have shown that the temperature $T$ and salinity $S$ are linearly
correlated (on average) in the turbulent bulk of an ice-ocean boundary layer
(IOBL) flow (Kimura et al., 2016; Stevens et al., 2020; Rosevear et al.,
2022a; Davis et al., 2023). In the $T-S$ diagram, this linear correlation
becomes evident when the mean values of instantaneous measurements cluster on
a straight line (figure 1(b)), which is then referred to as the meltwater
mixing line (MML).

Gade (1979) derived an analytical solution for the slope of the MML assuming
steady state, a uniform ambient and a turbulent bulk. The slope, or
correlation coefficient, relates the local temperature $T$ and salinity $S$
(in a time-averaged sense) to the constant temperature and salinity of the
ambient $T_{\infty}$ and $S_{\infty}$ through

$\displaystyle\mathcal{C}=\frac{T-T_{\infty}}{S-S_{\infty}},$ (1)

and depends on parameters, such as latent heat of fusion, and ice properties,
including ice salinity, whose variations in space and time are limited and
most often negligible.

(a)(b)

Figure 1: Illustration of sub-ice water mass properties as functions of (a)
depth and (b) in the $T-S$ diagram. mol and bulk highlight the molecular
boundary layer and turbulent bulk, respectively. The springs in (a)
illustrates the correlation between the mean $T$ and $S$ profiles predicted by
Gade’s MML theory and obtained when water mass transformations are dominated
by the turbulent mixing of a uniform ambient $(T_{\infty},S_{\infty})$ with
meltwater. The solid vertical gray line in (a) represents a conductivity,
temperature and depth (CTD) profiler, while the gray dots in (b) show an
ensemble of plausible measurements. $m$ denotes the melt rate.

In the MML model (1), the correlation is linear, i.e. $\mathcal{C}$ is a
constant (with value e.g. $\sim 2.7$ K/(g/kg) beaneath Pine Island Glacier ice
shelf, see Kimura et al., 2016), for one key reason, which is that it assumes
equal effective (turbulent and molecular) diffusivities for temperature and
salt (Gade, 1979).111We note that the linear correlation also hinges on the
subtle assumption that the ice melts in a vast amount of seawater, which is
most likely appropriate for all oceanic applications, and refer to McDougall
et al. (2014) for extensions of the theory to finite volumes of seawater. The
hypothesis that different scalars with different molecular diffusivities may
have the same effective diffusivities has been tested in laboratory
experiments (Jackson and Rehmann, 2003; Martin and Rehmann, 2006) and
numerical simulations (Smyth et al., 2005; Ma and Peltier, 2022) replicating
the open ocean environment. The aforementioned studies all show that the ratio
of the effective temperature diffusivity to the effective salt diffusivity is
indeed order unity when $Re_{b}\gg 100$, where $Re_{b}=\varepsilon/(\nu
N^{2})$ is the buoyancy Reynolds number ($\varepsilon$ is the viscous
dissipation rate, $\nu$ is the kinematic viscosity and $N$ is the Brunt-
Väisälä, or buoyancy, frequency), but that it increases (potentially above 2)
with decreasing $Re_{b}\leq O(100)$ (Gregg et al., 2018). In fact, at low
$Re_{b}$, buoyancy becomes dynamically important (Falor et al., 2023) and the
turbulent Prandtl and Schmidt numbers depend on $Re_{\tau}$ but also $Pr$ and
$Sc$, hence differ, such that diffential diffusion prevails (Bouffard and
Boegman, 2013).

The MML model is not valid when turbulence is weak, spatially localized and/or
intermittent, as in such cases molecular diffusivities influence the mean
temperature and salinity profiles with depth. For instance, melting in a
quiescient ambient leads to double-diffusive convection, with or without
layering, and properties in a $T-S$ diagram that deviate from the MML
prediction (Kimura et al., 2015; Rosevear et al., 2022a). The MML model is
also insufficient to describe the properties of the meltwater mixture when the
ambient is spatially stratified and/or changing over time, due to e.g. tidal
flows or episodic subglacial discharge (Davis et al., 2023). In such cases,
water mass properties span a two-dimensional subspace rather than a one-
dimensional line in the $T-S$ diagram. Assessing the importance of melting
relative to other phenomena in the water mass transformation process becomes
limited, requiring, for instance, a composite tracer method for the
calculation of the meltwater fraction (Jenkins, 1999).

The linear correlation between $T$ and $S$, when it is valid, provides an
opportunity to simplify the mathematical model governing the dynamics of the
IOBL. The temperature can be inferred from the salinity and vice versa, such
that it becomes unnecessary to solve the advection-diffusion equation for each
scalar field. A natural (or practical) choice for the scalar field to be
solved for in a reduced model of the IOBL is not $T$ nor $S$ but an aggregate
variable called thermal driving, which is defined as (Jenkins, 2011)

$\displaystyle
T_{*}=T-T_{f}(S,P_{b})=T-(\lambda_{1}S+\lambda_{2}+\lambda_{3}P_{b}).$ (2)

In equation (2), the freezing temperature $T_{f}$ is a function of salinity
$S$ and reference pressure $P_{b}$ (chosen at the ice-ocean boundary), and is
here linearized, as is the case in most oceanographic studies (with constants
$\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ listed in table 1). Thermal
driving is the temperature above the local freezing point. Thus, at the ice-
ocean interface it must satisfy $T_{*}=0$ while in the far field
$T_{*,\infty}=T_{\infty}-T_{f}(S_{\infty})$. Since temperature and salinity
can be inferred from thermal driving using equations (1) and (2), a reduced
thermal driving model, substituting the advection-diffusion equations for
temperature and salt with a single advection-diffusion equation for $T_{*}$,
can be used to reconstruct, or diagnose, the full IOBL properties in the
turbulent bulk. The advantage of the reduced thermal driving model over the
full $T-S$ model is significant. The diffusivity of thermal driving
$\kappa_{*}$ can be taken much larger than the diffusivity of salt
$\kappa_{S}$, such that the smallest-scale features of the reduced thermal
driving model, given by the Batchelor scale $\ell_{B,*}$, become much coarser
than the smallest-scale features of the full $T-S$ model, i.e., given by
$\ell_{B,S}$ – the ratio been given by
$\ell_{B,*}/\ell_{B,S}\propto\sqrt{Le_{*}}$ with
$Le_{*}=\kappa_{*}/\kappa_{S}$ (Zhou et al., 2017). Setting the molecular
diffusivity of thermal driving equal to the diffusivity of temperature
$\kappa_{*}=\kappa_{T}$, which is about 176 times larger than the diffusivity
of salt (Middleton et al., 2021), leads to thermal driving structures roughly
$\sqrt{176}\approx 13$ times coarser than salinity structures. A significant
limitation of the thermal driving model is that $T$ and $S$ can be inferred
from $T_{*}$ in the bulk only, i.e. not in the diffusive sublayers, since the
linear correlation (1) only applies where turbulence dominates over molecular
diffusion. A model for the diffusive sublayers of a turbulent subglacial ocean
current thus becomes necessary to extrapolate the temperature and salinity
profiles from the bulk to the ice-ocean interface (Gade, 1979).

Jenkins (2011) first introduced thermal driving as a useful prognostic
variable in a Reynolds-averaged model of the IOBL, i.e. wherein Reynolds
stresses and fluxes are replaced by turbulent viscosity and diffusivities to
focus on the smooth large-scale flow. The goal was to explore the effect of a
stratified ambient. Thus, a second prognostic scalar field, density deficit in
this case (relative to the ambient density), was kept in the model, which as a
result had the same complexity as the original $T-S$ formulation (see also
Hewitt, 2020). A Reynolds-averaged thermal driving model with reduced
complexity, i.e. with only one scalar equation assuming equal turbulent
diffusivity, was subsequently proposed to study the transient evolution and
stability of tilted IOBL under the effect of rotation, buoyancy, turbulent
mixing and geostrophic currents (Jenkins, 2016, 2021).

In this paper we achieve two goals. First, we use direct numerical simulations
(DNS) of a model sub-ice ocean current to validate the hypothesis of equal
effective diffusivities in the IOBL (§2). Second, we demonstrate that DNS
results of a reduced thermal driving model can be used to reconstruct
temperature and salinity fields consistent in the bulk with DNS results of a
full $T-S$ model, and to estimate the heat and salt fluxes at the ice-ocean
interface (§3). In so doing we provide evidence supporting the key hypothesis
of the MML prediction and proof that the reduced thermal driving model is a
worthwhile lightweight emulator of the IOBL dynamics, which we expect can
expand the contribution of high-resolution simulations to studies of ice-ocean
interactions (Vreugdenhil and Taylor, 2019; Mondal et al., 2019; Middleton et
al., 2021; Rosevear et al., 2021; Begeman et al., 2022; Vreugdenhil et al.,
2022; Rosevear et al., 2022b; Patmore et al., 2023). We focus on the well-
mixed regime for practical reasons (also called shear-dominated regime, see
Rosevear et al., 2022b). DNS are computationally expensive and the shear-
dominated regime (weak density stratification) can be investigated in
numerical domains much smaller than those required in the buoyancy-dominated
regime (strong density stratification). At relatively low $Re_{\tau}$, a
shear-dominated flow is turbulent in both small and large channels (in the
streamwise and spanwise directions), while a strongly-stratified flow is
turbulent in large channels only, i.e. requiring a relatively broad numerical
domain to avoid relaminarization, i.e. turbulence decay (García-Villalba and
del Álamo, 2011). The MML prediction and thermal driving model are also
expected to fail in the limit of strong density stratification, such that
verifying their validity in the limit of weak stratification seems a useful
first milestone.

Table 1: List of the dimensional physical parameters relevant to this study. All values are standard and similar to the ones used in previous publications (Jenkins et al., 2010; Vreugdenhil and Taylor, 2019; Middleton et al., 2021). They are appropriate for seawater below ice shelves up to several hundreds of meters thick. Parameter name | Symbol | Value
---|---|---
gravitational acceleration | $g$ | 9.81 m/s2
reference density | $\rho_{0}$ | 1000 g/kg
seawater heat capacity | $c_{p}$ | 3974 J/K/kg
latent heat of fusion | $L_{i}$ |  $3.35\times 10^{5}$ J/kg
salt-induced variation of the freezing temperature | $\lambda_{1}$ |  $-5.73\times 10^{-2}$ K/(g/kg)
surface freezing temperature of freshwater | $\lambda_{2}$ |  $8.32\times 10^{-2}$ ∘C
pressure-induced variation of the freezing temperature | $\lambda_{3}$ |  $-7.53\times 10^{-4}$ K/dbar
kinematic viscosity | $\nu$ |  $1.8\times 10^{-6}$ m2/s
thermal expansion coefficient | $\alpha$ |  $3.87\times 10^{-5}$ 1/K
haline contraction coefficient | $\beta$ |  $7.86\times 10^{-4}$ 1/(g/kg)
domain depth | $H$ | 1 m
far-field salinity | $S_{\infty}$ | 35 g/kg
far-field thermal driving | $T_{*,\infty}$ |  $10^{-5}$ K

##  2 Test of the equal effective diffusivities assumption

###  2.1 Full $T-S$ model equations

To test the equal effective diffusivities hypothesis we consider an
externally-forced stratified sub-ice ocean flow in a horizontally periodic
channel (cf. figure 2). We use a Cartesian coordinates system ($x,y,z$) with
unit vector $\bm{e_{z}}$ (bold variables denote vectors) upward opposite to
gravity $\bm{g}=-g\bm{e_{z}}$ to describe the dynamics and set the bottom of
the channel as the origin of the $z$ axis. The ocean flow is forced by a
uniform body force per unit volume $F$ (in kg.m-2.s-2) in the $\bm{e_{x}}$
direction. We use the Boussinesq approximation. The reference density is
$\rho_{0}$ and the density anomaly is a linear function of temperature and
salinity (Middleton et al., 2021).

The dimensional velocity vector, temperature, salinity and pressure are
denoted by $\bm{u}=(u,v,w)$, $T$, $S$ and $P$, respectively. For convenience,
we solve the governing equations in dimensionless form. The dimensionless
coordinates and variables, which we denote by tildes, are related to their
dimensional counterparts through

$\displaystyle(\tilde{x},\tilde{y},\tilde{z})=\frac{(x,y,z)}{H},\quad\tilde{t}=\frac{t}{t_{\nu}},\quad\tilde{\bm{u}}=\frac{\bm{u}}{u_{\nu}},\quad\tilde{T}=\frac{T-T_{f}(S_{\infty},P_{b})}{T_{*,\infty}},\quad\tilde{S}=\frac{S-S_{\infty}}{S_{\infty}},\quad\tilde{P}=\frac{P-P_{st}}{p_{\nu}},$
(3)

where $H$ is the channel thickness, $t_{\nu}=H^{2}/\nu$ is the viscous time
scale, $u_{\nu}=\nu/H$ is the viscous velocity scale,
$T_{f}(S_{\infty},P_{b})$ is the freezing temperature (cf. equation (2)) based
on the far-field salinity ($S_{\infty}$) and constant hydrostatic pressure at
the ice-ocean interface ($P_{b}$), $T_{*,\infty}=T_{\infty}-T_{f}(S_{\infty})$
is the far-field thermal driving, $P_{st}=P_{b}-\rho_{0}g(z-H)$ is the
hydrostatic pressure and $p_{\nu}=\rho_{0}\nu^{2}/H^{2}$ is the dynamic
pressure scale. Note that $\tilde{S}$ is a dimensionless salinity anomaly,
which is always negative (since it is always lower than in the far field).

Substituting the dimensionless coordinates and variables (3) into the IOBL
bulk equations, which consist of the Navier-Stokes equations in the Boussinesq
approximation and advection-diffusion equations for the temperature and
salinity fields (resulting from the conservation of energy and salt), yields
the dimensionless governing equations

$\displaystyle\partial_{\tilde{t}}\tilde{\bm{u}}+\tilde{\bm{u}}\cdot\tilde{\nabla}\tilde{\bm{u}}=-\tilde{\nabla}\tilde{P}+\tilde{\nabla}^{2}\tilde{\bm{u}}+Re_{\tau}^{2}\bm{e_{x}}-Re_{\tau}^{2}Ri_{\tau}(\tilde{S}-R_{\rho}\tilde{T})\bm{e_{z}},$
(4a) $\displaystyle\tilde{\nabla}\cdot\tilde{\bm{u}}=0,$ (4b)
$\displaystyle\partial_{\tilde{t}}\tilde{T}+\tilde{\bm{u}}\cdot\tilde{\nabla}\tilde{T}=Pr^{-1}\tilde{\nabla}^{2}\tilde{T}+\tilde{\mathcal{R}}_{T},$
(4c)
$\displaystyle\partial_{\tilde{t}}\tilde{S}+\tilde{\bm{u}}\cdot\tilde{\nabla}\tilde{S}=Sc^{-1}\nabla^{2}S+\tilde{\mathcal{R}}_{S},$
(4d)

with $Re_{\tau}$, $Ri_{\tau}$, $R_{\rho}$, $Pr$ and $Sc$ the dimensionless
control parameters (see table 2 for the definition of all control parameters
and the range of values considered). The rightmost terms in the temperature
and salinity equations (4c)-(4d) are damping terms of the form
$\tilde{\mathcal{R}}_{X}=Re_{\tau}^{3/2}\left(\frac{\langle\tilde{X}\rangle_{\perp}-\tilde{X}_{\infty}}{2}\right)\left[1-\tanh\left(\frac{\tilde{z}-0.2}{0.025}\right)\right]$,
where $\langle\tilde{X}\rangle_{\perp}$ denotes the horizontal average and
$\tilde{X}=\tilde{T},\tilde{S}$. The damping terms force the mean temperature
and salinity fields to the prescribed far-field values $\tilde{T}_{\infty}=1$
and $\tilde{S}_{\infty}=0$ where $\tilde{z}\lessapprox 0.2$. The relaxation
time scale is $Re_{\tau}^{-3/2}$ (in viscous time units), such that it is
roughly inversely proportional to the mean ($Re_{m}$) and centreline
($Re_{c}$) Reynolds numbers (not shown; note that this is a steeper scaling
than the classical $Re_{m}\sim Re_{\tau}^{8/7}$ scaling because of confinement
effects, cf. Hwang (2013)). Thus the dimensional relaxation time scale is
inversely proportional to the far-field velocity, similar to previous studies
(Vreugdenhil and Taylor, 2019). Note that there is no damping of the velocity
field.

We limit our analysis to slowly melting ice such that the ice-ocean interface
can be considered fixed at leading order (i.e. the simulation time is short
relative to the melting time scale $\sim StPr$). The boundary conditions are
free-slip (appropriate for a half channel flow), fixed temperature and fixed
salinity on the bottom boundary, and no-slip, freezing temperature and melt-
induced dilution at the top boundary, i.e., in dimensionless form

$\displaystyle\tilde{\bm{u}}=(\tilde{u},\tilde{v},\tilde{w})=\bm{0},\quad\tilde{T}=\gamma\tilde{S},\quad\partial_{\tilde{z}}\tilde{S}=LeSt^{-1}\left(1+\tilde{S}\right)\partial_{\tilde{z}}\tilde{T},\quad$
$\displaystyle\text{at}\quad\tilde{z}=1,$ (5a)
$\displaystyle\partial_{\tilde{z}}\tilde{u}=\partial_{\tilde{z}}\tilde{v}=\tilde{w}=0,\quad\tilde{T}=\tilde{T}_{\infty}=1,\quad\tilde{S}=\tilde{S}_{\infty}=0,\quad$
$\displaystyle\text{at}\quad\tilde{z}=0.$ (5b)

The melt-induced dilution equation (third equation in (5a)) arises from the
combination of the energy and salt conservation equations at the interface,
which, neglecting heat and salinity fluxes from the ice, can be written in
dimensionless form as

$\displaystyle\tilde{m}St=-Pr\partial_{\tilde{z}}\tilde{T},\quad\tilde{m}\left(1+\tilde{S}\right)=-Sc\partial_{\tilde{z}}\tilde{S},\quad$
$\displaystyle\text{at}\quad\tilde{z}=1,$ (6)

with $\tilde{m}$ the dimensionless melt rate. Note that we neglect meltwater
advection due to the density difference between liquid water and solid ice
when we write the no-slip condition (5a) at $\tilde{z}=1$. This approximation
is valid in the limit of infinitesimal thermal driving, which we focus on, but
can lead to significant errors when thermal driving is of order 1∘C or more as
it neglects the vertical transport of sensible heat due to the melt-induced
up-welling (Jenkins et al., 2001). This transport term is negligible in our
simulations because it is approximately equal to $\tilde{m}\tilde{T}$, which
is of order $10^{-7}$ or smaller in all cases, hence negligible compared to
the turbulent heat flux that is of order $O(10)$.

(a)(b)(c)

Figure 2: Snapshots obtained at the final time of the simulations for $Pr=1$, $Le=10$ and (a) $Re_{\tau}=200$, (b) $Re_{\tau}=400$ and (c) $Re_{\tau}=800$. We show the normalized vertical velocity $\tilde{w}/Re_{\tau}$ on the upstream vertical slice ($\tilde{x}=0$), the normalized heat flux $\tilde{q}_{T}=[-\partial_{\tilde{z}}\tilde{T}/(1-\langle\tilde{T}\rangle_{\perp})]|_{\tilde{z}=1}$ at the ice-ocean interface (a local version of $\tilde{Q}_{T}$ in (7)), and the salinity on the right side of the channel ($\tilde{y}=0$). The reduction of the domain size with $Re_{\tau}$ is obvious. The three colorbars apply to all snapshots. Table 2: List of the dimensionless control parameters appearing in the governing equations (4)-(5). We provide typical values for the ocean besides those considered in this study. The list of physical parameters used to construct the dimensionless parameters can be found in table 1. The ocean values consider a friction velocity between 0 and 0.18 m/s and thermal driving between 0 and 10 K (other parameters taken from table 1). Note that the dimensionless domain lengths $L_{x}$ and $L_{y}$ are two other control parameters that are changed across simulations (see table 3). The friction Richardson number is based on a salinity gradient as salinity differences typically dominate density differences. Parameter name | Symbol | Definition | Ocean value(s) | Study value(s) | Comments
---|---|---|---|---|---
friction Reynolds | $Re_{\tau}$ | $\sqrt{\frac{FH^{3}}{\rho_{0}\nu^{2}}}$ |  $0$ to $10^{5}$ | 200 to 800 |
friction Richardson | $Ri_{\tau}$ | $\frac{g\rho_{0}\beta S_{\infty}}{F}$ | 0 to $+\infty$ |  $1.3\times 10^{5}$ to $2.1\times 10^{6}$ | Salt Rayleigh number $Ra_{S}=Re_{\tau}^{2}Ri_{\tau}Sc$
density ratio | $R_{\rho}$ | $\frac{\alpha T_{*,\infty}}{\beta S_{\infty}}$ | 0 to $1.4\times 10^{-2}$ | $1.4\times 10^{-8}$ |
Prandtl | $Pr$ | $\frac{\nu}{\kappa_{T}}$ | $\approx 14$ | 1 to 10 |
Schmidt | $Sc$ | $\frac{\nu}{\kappa_{S}}$ | $\approx 2500$ | 1 to 100 | Lewis number $Le=Sc/Pr$
salt-induced freezing slope | $\gamma$ | $\frac{\lambda_{1}S_{\infty}}{T_{*,\infty}}$ |  $-\infty$ to $-0.2$ | $-2.0\times 10^{5}$ |
Stefan | $St$ | $\frac{L_{i}}{c_{w}T_{*,\infty}}$ |  $8.4$ to $+\infty$ | $8.4\times 10^{6}$ |

###  2.2 Numerical experiment

The governing equations (4)-(5) of the full $T-S$ model are solved numerically
using the open-source pseudo-spectral code Dedalus (Burns et al., 2020). We
use Fourier mode decompositions in the horizontal directions and Chebyshev
mode decomposition in the vertical direction. We use a third-order multi-stage
explicit-implicit Runge-Kutta scheme for time integration and we set the
Courant-Friedrich-Lewy parameter to 0.5. We run 11 simulations with varying
$Re_{\tau}$, $Pr$ and $Le$ in order to test the equal effective diffusivities
hypothesis underlying the MML theory. Owing to computational limitations, we
cannot run simulations with $1\ll Pr\ll Le$, which is the parameter range
relevant to the natural ocean. In particular, the Schmidt number (or inverse
of the normalized salt diffusivity) does not exceed 100, hence is always less
than the ocean value $Sc\approx 2464$ (Middleton et al., 2021). However, the
MML theory does not depend on any particular value of $Sc$ (or $Pr$ and $Le$)
and we will show that the equal effective diffusivity hypothesis is in fact
valid for all $(Re_{\tau},Pr,Le)$ tested, including cases with $Pr=10$ and
$Le=10$ that have the same hierarchy of scales as the natural ocean (i.e.
large momentum structures, medium-size thermal structures and small salinity
structures). The resolution of the simulations is discussed in the
Supplementary Material. In brief, the simulations adequately resolve the
diffusive sublayers, thanks to the clustering of Chebyshev points near
boundaries, as well as the Kolmogorov length scale. The Batchelor length scale
is not always well resolved, especially at high $Le$. However, running higher-
resolution simulations has shown that the flow statistics are not sensitive to
the relatively coarse resolution of the simulations in the bulk.

We focus on the weakly-stratified limit of the IOBL. We ensure that the
simulations remain turbulent at all times (no relaminarization) by consideing
a very small far-field thermal driving $T_{*,\infty}=10^{-5}$ K (and a far-
field salinity $S_{\infty}=35$ g/kg). All simulations are first run without
the scalar fields until a statistical steady state is reached. Then the
governing equations for the scalar fields are solved for and coupled to the
hydrodynamics. The simulations are run for a relatively long time (up to 100
friction time units) to guarantee that a statistical steady-state is reached
and enough data is accumulated to compute reasonably well converged
statistics. Table 3 lists all simulation cases and parameters.

For computational expediency, we run the simulations in a minimal channel flow
unit (MCU), i.e., considering short and narrow channels (Jiménez and Moin,
1991; Flores and Jiménez, 2010). We hypothesize that the validity of the equal
turbulent diffusivities hypothesis in an MCU implies the validity in an
unconfined channel. This hypothesis is supported by the fact that the small-
scale turbulence remains healthy in an MCU, especially near the walls, despite
well-known horizontal confinement effects in the outer layer (Flores and
Jiménez, 2010; Hwang, 2013; Lozano-Durán and Jiménez, 2014; De Giovanetti et
al., 2016). It is also supported by the fact that one test case run in both an
MCU and a larger domain (two times larger in both horizontal directions)
resulted in equal turbulent diffusivities irrespective of the domain size. The
difference between the drag coefficient $C_{D}$ (defined in the next section)
obtained for a MCU and a domain twice as large when $Re_{\tau}=400$, $Pr=1$
and $Le=10$, which we report in the last column of the eighth and last lines
of table 3, highlights the impact that the (small) horizontal lengths have on
the global dynamics (but not on the equal turbulent diffusivity hypothesis),
in agreement with previous studies (Hwang, 2013; De Giovanetti et al., 2016).

###  2.3 Diagnostics

Our analysis of the simulations focuses on horizontal-plane and volume
averages of the variables and on vertical momentum, temperature and salinity
fluxes. We write $\langle\cdot\rangle_{\perp}$,
$\langle\cdot\rangle_{\tilde{z}}$ and $\langle\cdot\rangle$ the horizontal,
vertical and volume averaging operators, respectively, and
$\overline{\;\cdot\;}$ the time averaging operator (computed at statistical
steady state).

We define the dimensionless horizontally-averaged heat and salt fluxes as

$\displaystyle\tilde{Q}_{T}=\langle\tilde{w}\tilde{T}-Pr^{-1}\partial_{\tilde{z}}\tilde{T}\rangle_{\perp},\quad\tilde{Q}_{S}=\langle\tilde{w}\tilde{S}-Sc^{-1}\partial_{\tilde{z}}\tilde{S}\rangle_{\perp},$
(7)

and note that the convective fluxes in (7) are zero at the no-slip top and
impermeable bottom boundaries. The vertical transport efficiency of momentum,
heat and salinity is evaluated from the temperature Nusselt number, salinity
Nusselt number and drag coefficient, defined as

$\displaystyle
Nu_{T}=\frac{\overline{\tilde{Q}_{T}(\tilde{z}=1)}}{\tilde{Q}_{T}^{diff}},\quad
Nu_{S}=\frac{\overline{\tilde{Q}_{S}(\tilde{z}=1)}}{\tilde{Q}_{S}^{diff}},\quad
C_{D}=2\left(\frac{Re_{\tau}}{\langle\overline{\tilde{u}}\rangle}\right)^{2},$
(8)

where

$\displaystyle\tilde{Q}_{T}^{diff}=-Pr^{-1}\left(\langle\overline{\tilde{T}(\tilde{z}=1)}\rangle_{\perp}-1\right),\quad\tilde{Q}_{S}^{diff}=-Sc^{-1}\langle\overline{\tilde{S}(\tilde{z}=1)}\rangle_{\perp},$
(9)

are the heat and salt fluxes of the equivalent purely-diffusive state
($\tilde{\bm{u}}=\tilde{\bm{0}}$). We recall that $\tilde{T}(\tilde{z}=0)=1$
and $\tilde{S}(\tilde{z}=0)=0$ by definition and that the channel height is
unity.

We test the validity of the meltwater mixing line theory by computing the
effective Prandtl and Lewis numbers (superscripts $e$), defined as

$\displaystyle Pr^{e}=\frac{\tilde{\nu}^{e}}{\tilde{\kappa}_{T}^{e}},\quad
Le^{e}=\frac{\tilde{\kappa}_{T}^{e}}{\tilde{\kappa}_{S}^{e}},$ (10)

with $\tilde{\nu}^{e}$, $\tilde{\kappa}_{T}^{e}$, $\tilde{\kappa}_{S}^{e}$
(resp. $\tilde{\nu}^{t}$, $\tilde{\kappa}_{T}^{t}$, $\tilde{\kappa}_{S}^{t}$)
the effective (resp. turbulent; superscripts $t$) viscosity and temperature
and salt diffusivities, which we write as

$\displaystyle\tilde{\nu}^{e}=\tilde{\nu}^{t}+1=\frac{\langle\overline{\tilde{w}\tilde{u}-\partial_{\tilde{z}}\tilde{u}}\rangle_{\perp}}{-\langle\overline{\partial_{\tilde{z}}\tilde{u}}\rangle_{\perp}},\quad\tilde{\kappa}_{T}^{e}=\tilde{\kappa}_{T}^{t}+Pr^{-1}=\frac{\overline{\tilde{Q}_{T}}}{-\langle\overline{\partial_{\tilde{z}}\tilde{T}}\rangle_{\perp}},\quad\tilde{\kappa}_{S}^{e}=\tilde{\kappa}_{S}^{t}+Sc^{-1}=\frac{\overline{\tilde{Q}_{S}}}{-\langle\overline{\partial_{\tilde{z}}\tilde{S}}\rangle_{\perp}}.$
(11)

The convective terms in (11) involve the full variables; however, they could
be replaced by their fluctuating counterparts since
$\langle\overline{\tilde{w}}\rangle_{\perp}=0$ $\forall\tilde{z}$ at
statistical steady state. Finally, the buoyancy Reynolds number reported in
table 3 is defined as

$\displaystyle
Re_{b}=\langle\frac{\langle\overline{\tilde{\epsilon}}\rangle_{\perp}Sc}{Ra_{S}\langle\overline{-\partial_{\tilde{z}}\tilde{S}+R_{\rho}\partial_{\tilde{z}}\tilde{T}}\rangle_{\perp}}\rangle_{\tilde{z}>0.5},$
(12)

with $\langle\cdot\rangle_{\tilde{z}>0.5}$ denoting the average over the upper
half of the domain and $\langle\overline{\tilde{\epsilon}}\rangle_{\perp}$ the
mean (pseudo) turbulent kinetic energy dissipation rate (cf. Supplementary
Material). Here we find that the buoyancy Reynolds number is always much
larger than O(100) (cf. table 3), such that we might expect equal turbulent
diffusivities in all our simulations.

Table 3: List of key input parameters (physical and numerical) and output variables for all simulation cases. $(n_{x},n_{y},n_{z})$ is the number of Fourier and Chebyshev modes in all three directions (not including the scaling factor 3/2 used for dealiasing), $Le^{e}$ is the effective Lewis number based on the effective diffusivities of temperature and salt in the bulk, i.e., averaged between $\tilde{z}=0.3$ and $\tilde{z}=0.7$ and $Re_{b}$ is the buoyancy Reynolds number given by equation (12). For the reduced $T_{*}$ simulations $Nu_{T}$ is the Nusselt number with $\tilde{T}_{*}$ replacing $\tilde{T}$. The far-field thermal driving and salinity are set to $T_{*,\infty}=10^{-5}$ K and $S_{\infty}=35$ g/kg in all cases. Simulation name Full wide on the last line refers to a Full $T-S$ simulation with double the length in the horizontal directions. Type | $Re_{\tau}$ |  $Pr$ ($Pr_{*}$) | $Le$ | $(L_{x},L_{y},L_{z})$ | $(n_{x},n_{y},n_{z})$ | $Le^{eff}$ | $Re_{b}$ | $Nu_{T}$ | $Nu_{S}$ | $C_{D}$
---|---|---|---|---|---|---|---|---|---|---
Full $T-S$ | 200 | 1 | 1 | (2,1,1) | (64,64,48) | 1.000 | $10^{4}$ | 8.6 | 8.6 | 0.0080
Full $T-S$ | 200 | 1 | 10 | (2,1,1) | (128,128,96) | 0.997 | $10^{4}$ | 8.4 | 27.3 | 0.0078
Full $T-S$ | 200 | 1 | 30 | (2,1,1) | (256,256,192) | 1.002 | $10^{4}$ | 8.1 | 41.5 | 0.0079
Full $T-S$ | 200 | 1 | 100 | (2,1,1) | (256,256,192) | 1.002 | $2\times 10^{4}$ | 8.3 | 63.7 | 0.0079
Full $T-S$ | 200 | 5 | 5 | (2,1,1) | (256,256,192) | 0.991 | $2\times 10^{4}$ | 27.2 | 65.0 | 0.0080
Full $T-S$ | 200 | 10 | 10 | (2,1,1) | (256,256,192) | 0.992 | $4\times 10^{4}$ | 20.1 | 39.9 | 0.0080
Full $T-S$ | 400 | 1 | 1 | (1,0.5,1) | (64,64,96) | 1.000 | $10^{5}$ | 13.9 | 13.9 | 0.0056
Full $T-S$ | 400 | 1 | 10 | (1,0.5,1) | (128,128,192) | 0.996 | $10^{5}$ | 14.3 | 54.1 | 0.0057
Full $T-S$ | 400 | 10 | 10 | (1,0.5,1) | (256,256,384) | 1.000 | $3\times 10^{5}$ | 54.2 | 137.3 | 0.0056
Full $T-S$ | 800 | 1 | 1 | (0.5,0.25,1) | (64,64,128) | 1.000 | $2\times 10^{6}$ | 19.7 | 19.7 | 0.0031
Full $T-S$ | 800 | 1 | 10 | (0.5,0.25,1) | (128,128,256) | 1.011 | $2\times 10^{5}$ | 20.2 | 95.0 | 0.0031
Reduced $T_{*}$ | 200 | 1 | - | (2,1,1) | (64,64,48) | - | - | 8.7 | - | 0.0078
Reduced $T_{*}$ | 200 | 5 | - | (2,1,1) | (128,128,96) | - | - | 20.1 | - | 0.0079
Reduced $T_{*}$ | 200 | 10 | - | (2,1,1) | (128,128,96) | - | - | 28 | - | 0.0080
Reduced $T_{*}$ | 400 | 1 | - | (1,0.5,1) | (64,64,96) | - | - | 14.0 | - | 0.0056
Reduced $T_{*}$ | 400 | 10 | - | (1,0.5,1) | (128,128,192) | - | - | 53 | - | 0.0056
Reduced $T_{*}$ | 800 | 1 | - | (0.5,0.25,1) | (64,64,128) | - | - | 19.5 | - | 0.0032
Full wide | 400 | 1 | 10 | (2,1,1) | (256,256,192) | 0.985 | $10^{5}$ | 16.3 | 55.1 | 0.0062

###  2.4 Results

To illustrate the relatively long time of the simulations, we first show in
figure 3(a) the time series of the normalized and horizontally-averaged heat
flux at the ice-ocean interface
$\tilde{Q}_{T}(\tilde{z}=1)/\tilde{Q}_{T}^{diff}$. We show the results from
five simulations only for clarity. The temporal windows used for time
averaging are highlighted by the horizontal black lines and last at least 25
friction time units. The time-averaged normalized heat flux is then shown as a
function of $\tilde{z}$ in figure 3(b). The heat flux is almost depth
invariant for $\tilde{z}\geq 0.2$, in agreement with the time-averaged heat
equation at statistical steady state (without the relaxation term). For the
most ambitious simulation ($Re_{\tau},Pr,Le$)=(400,10,10) we observe a
$\sim$10% deviation with depth (suggesting that the boundary layer is not
perfectly equilibrated yet), which we expect would decrease with increasing
simulation time and statistical averaging. For $\tilde{z}\leq 0.2$ the
source/sink of heat comes from the relaxation term $\tilde{\mathcal{R}}_{T}$,
which is not included in the definition of ${\tilde{Q}_{T}}$. As expected,
$\overline{\tilde{Q}_{T}}\rightarrow 0$ as $\tilde{z}\rightarrow 0$ because of
the no penetration condition at $\tilde{z}=0$ and the imposed uniform mean
temperature profile in the damping region. Figure 3(c) shows the mean
horizontal velocity as a function of depth. As expected, we observe a boundary
layer close to the (top) ice-ocean interface and a maximum at $\tilde{z}=0$
where the flow is shear free. The dashed lines show the mean horizontal
velocity obtained in reduced $T_{*}$ simulations (described in §3). They
overlap perfectly with the solid lines, which show the full $T-S$ simulations
results. This overlap does not yet prove that the reduced thermal-driving
model can replace the full $T-S$ model as we are running the simulations in
the shear-dominated regime, such that the scalar fields are almost passive and
have little influence on the dynamics. The decrease of the drag coefficient
with $Re_{\tau}$ (reported in table 3) is consistent with the increase of the
normalized velocity in figure 3(c). It is partly due to confinement effects,
the higher $Re_{\tau}$ simulations having smaller domain sizes (see figure 4
in Hwang, 2013).

(a)(b)(c)

Figure 3: (a) Normalized heat flux at the ice-ocean interface as a function of
time (in friction time units) for five cases. The horizontal black lines show
the temporal windows used for time averaging. (b) Normalized time-averaged
heat flux as a function of depth. The vertical dotted lines highlight the
normalized heat flux at the ice-ocean interface ($\tilde{z}=1$) while the
horizontal dashed line highlights the upper limit of the damping region. Note
that the $(200,1,100)$ and $(200,1,10)$ results overlap. (c) Normalized time-
averaged horizontal velocity as a function of depth. The velocity profiles
with the same $Re_{\tau}$ overlap. Results from the reduced $T_{*}$ model are
also shown (dashed lines) but are indistinguishable from the full $T-S$ model
results (solid lines).

We show the vertical profiles of $\tilde{\nu}^{t}$, $\tilde{\kappa}_{T}^{t}$
and $\tilde{\kappa}_{S}^{t}$ in figure 4(a) (masking the damping region
$\tilde{z}\leq 0.2$ for clarity). In all cases, $\tilde{\nu}^{t}\gg 1$,
$\tilde{\kappa}_{T}^{t}\gg 1$ and $\tilde{\kappa}_{S}^{t}\gg 1$ when
$0.2\leq\tilde{z}\leq 0.9$, demonstrating the turbulent nature of the flow and
predominance of convective transports over molecular diffusion in the bulk.
The turbulent viscosity reaches a maximum relatively close to the ice-ocean
interface, where the shear is large, and then decreases with depth, in
agreement with previous studies. The turbulent diffusivities reach similarly
large values close to the top boundary but then decay more slowly with depth.
Consequently, $Pr^{e}$ and $Sc^{e}$ become most often less than one in the
bulk (as observed in figure 4(b) for $Pr^{e}$), in agreement with previous
models and observations of turbulent Prandtl and Schmidt numbers when
$Pr,Sc\geq 1$ (Reynolds, 1975; Alcántara-Ávila and Hoyas, 2021). The key
result of this section is that the effective Lewis number
$Le^{e}=\tilde{\kappa}_{T}^{e}/\tilde{\kappa}_{S}^{e}$ is close to one in most
of the bulk in all cases (see figure 4(c)).

(a)(b)(c)

Figure 4: (a) Turbulent viscosity (thick solid lines) and temperature
diffusivity (thin) as a function of depth, (b) effective Prandtl number as a
function of depth, and (c) effective Lewis number as a function of depth. We
only show five cases for clarity. The vertical dashed lines highlight ratios
equal to 1. Note that the turbulent temperature diffusivity is not shown in
the relaxation region $\tilde{z}\leq 0.2$ (below the horizontal dashed lines)
because the damping terms generate depth-invariant temperature (and salinity)
profiles with unrealistically large turbulent diffusivities. By construction,
$Pr^{e}\rightarrow Pr$ and $Le^{e}\rightarrow Le$ when $\tilde{z}\rightarrow
1$ (see the inset in (c) zooming on the top diffusive sublayer).

To further demonstrate the validity of the MML theory in the weakly-stratified
limit, we show the mean temperature profile
$\langle\overline{\tilde{T}}\rangle_{\perp}$ as a function of the mean
salinity profile $\langle\overline{\tilde{S}}\rangle_{\perp}$ in figure 5 (the
$T-S$ diagram). All simulation results perfectly overlap with the MML
prediction in the bulk.222Note that this is true also near the bottom boundary
(top of the curves) because of the imposed relaxation terms. Without
relaxation terms the uniform far-field is replaced by a boundary layer with
$T-S$ correlations that are different from the MML prediction, unless the
bottom boundary is also an ice-ocean interface and buoyancy effects are
negligible. The correlation between
$\langle\overline{\tilde{T}}\rangle_{\perp}$ and
$\langle\overline{\tilde{S}}\rangle_{\perp}$ is always the same far from the
ice-ocean interface, i.e. the curves all have the same slope, regardless of
the control parameters $Re_{\tau}$, $Pr$ and $Le$, in agreement with the MML
theory (derived in appendix).

Figure 5: Mean temperature shown as a function of mean salinity (both are
depth-dependent variables) for five cases. The thick solid black line shows
the freezing temperature while the dashed black line shows the meltwater
mixing line prediction.

The $\langle\overline{\tilde{T}}\rangle_{\perp}$ profiles as functions of
$\langle\overline{\tilde{S}}\rangle_{\perp}$ differ between cases as they
approach the freezing line. The salinity deficit at the interface, i.e.
$-\overline{\tilde{S}}(\tilde{z}=1)$, increases with $Le$ (compare $Le=10$ and
$Le=100$ for $Pr=1$ and $Re_{\tau}=200$) and $Pr$ (compare $Pr=1$ and $Pr=10$
for $Le=10$ and $Re_{\tau}=400$) but decreases with $Re_{\tau}$ (dark to light
red curves). This parametric dependence is consistent with previous studies,
once we realize that the interfacial salinity is inversely proportional to the
ratio of Nusselt numbers, and thus to the ratio of the salt to heat exchanges
velocities, denoted by $\gamma_{S}$ and $\gamma_{T}$ (m.s-1), which are
parameterized in ISOBL models (Holland and Jenkins, 1999). In fact,
$\gamma_{T}=\frac{Nu_{T}\kappa_{T}}{H}$ and
$\gamma_{S}=\frac{Nu_{S}\kappa_{S}}{H}$ by definition, such that

$\displaystyle\frac{\gamma_{S}}{\gamma_{T}}=\frac{Nu_{S}}{Nu_{T}Le}=\frac{\overline{\tilde{Q}_{T}}(\tilde{z}=1)\left[1-\overline{\tilde{T}}(\tilde{z}=1)\right]}{\overline{\tilde{Q}_{S}}(\tilde{z}=1)\left[-\overline{\tilde{S}}(\tilde{z}=1)\right]}=\frac{St\left[1-\overline{\tilde{T}}(\tilde{z}=1)\right]}{-\left[1+\overline{\tilde{S}}(\tilde{z}=1)\right]\overline{\tilde{S}}(\tilde{z}=1)}\approx\frac{St}{-\overline{\tilde{S}}(\tilde{z}=1)}.$
(13)

The second to last expression in (13) combines equations (6) and (7) while the
last expression is obtained in the limit
$\overline{\tilde{S}}(\tilde{z}=1)\rightarrow 0$. We show $Nu_{S}/(Nu_{T}Le)$,
i.e. the exchange velocity ratio, as a function of $Le$ in figure 6. We
clearly observe a decrease of $Nu_{S}/(Nu_{T}Le)$ with $Le$ and a (relatively)
more modest increase with $Re_{\tau}$ for $Le=10$ and $Pr=1$, consistent with
the $Le$\- and $Re_{\tau}$-dependence of the interfacial salinity. The $Le$
dependence of $Nu_{S}/(Nu_{T}Le)$ can be modelled with a power law with a -2/3
to -1/2 exponent, in good agreement with the parameterizations for
$\gamma_{S}$ and $\gamma_{T}$ used in Holland and Jenkins (1999) for shear-
driven flows. Starting from equations (11)-(12) in Holland and Jenkins (1999)
we find that $\gamma_{S}/\gamma_{T}$ scales like $Le^{-2/3}$ in the large
$Pr$, large $Sc$ and (relatively) small $Re_{\tau}$ limit. The $Le$ dependence
of $Nu_{S}/(Nu_{T}Le)$ is also consistent with the -2/3 (resp. -1/2) power law
scaling uncovered in recent simulations of multicomponent vertical convection
at relatively low $Le$ (resp. high $Le$) (Howland et al., 2023). The decrease
of $-\overline{\tilde{S}}(\tilde{z}=1)$ with $Re_{\tau}$ and increase with
$Pr$ observed in figure 5 are also in good agreement with the exchange
velocity formula in Holland and Jenkins (1999) since $\gamma_{S}$ (resp.
$\gamma_{T}$) increase with $Re_{\tau}$ and $Pr$ (resp. $Sc$) such that the
ratio $\gamma_{S}/\gamma_{T}$ decreases with $Pr$ (as long as $Le>1$) but
increases with $Re_{\tau}$.

Figure 6: Nusselt number ratio shown as a function of $Le$ and normalized by
$Le$ such that it corresponds to the ratio of salt to temperature transfer
velocities in Holland and Jenkins (1999).

##  3 Test of the reduced thermal driving model

###  3.1 Mathematical derivation

In this section we introduce the governing equations for the reduced thermal
driving model and the formula to reconstruct the temperature and salinity
fields in the bulk using dimensionless variables. We substitute the
temperature $\tilde{T}$ and salinity $\tilde{S}$ with a single scalar field,
called thermal driving, which we define as
$\tilde{T}_{*}=\tilde{T}-\gamma\tilde{S}$, i.e. as the temperature above the
local freezing point. From the MML model, which predicts
$\tilde{T}-1=\tilde{S}St$ (see appendix) and the definition of
$\tilde{T}_{*}$, we can rewrite the buoyancy term
$(\tilde{S}-R_{\rho}\tilde{T})$ as
$\tilde{T}_{*}\left(\frac{1-R_{\rho}St}{St-\gamma}\right)+\frac{\gamma
R_{\rho}-1}{St-\gamma}$ and discard the last term since it is constant
($\alpha$ and $\beta$ are constants, see table 1) and can be balanced by the
hydrostatic pressure. The governing equations for the reduced thermal driving
model become

$\displaystyle\partial_{\tilde{t}}\tilde{\bm{u}}+\tilde{\bm{u}}\cdot\tilde{\nabla}\tilde{\bm{u}}=-\tilde{\nabla}\tilde{P}+\tilde{\nabla}^{2}\tilde{\bm{u}}+Re_{\tau}^{2}\vec{e_{x}}-Re_{\tau}^{2}Ri_{\tau}\tilde{T}_{*}\left(\frac{1-R_{\rho}St}{St-\gamma}\right)\vec{e_{z}},$
(14a) $\displaystyle\tilde{\nabla}\cdot\tilde{\bm{u}}=0,$ (14b)
$\displaystyle\partial_{\tilde{t}}\tilde{T}_{*}+\tilde{\bm{u}}\cdot\tilde{\nabla}\tilde{T}_{*}=\frac{1}{Pr_{*}}\tilde{\nabla}^{2}\tilde{T}_{*}+\tilde{\mathcal{R}}_{T_{*}},$
(14c)

with $Pr_{*}$ the Prandtl number for the scalar field $\tilde{T}_{*}$ and
$\tilde{\mathcal{R}}_{T_{*}}$ a relaxation term similar to
$\tilde{\mathcal{R}}_{T}$ in equation (4c). The kinematic boundary conditions
are the same as for the full $T-S$ model, and the boundary conditions for
thermal driving are simply

$\displaystyle\tilde{T}_{*}(\tilde{z}=1)=0,\quad\tilde{T}_{*}(\tilde{z}=0)=\tilde{T}_{*,\infty}=1.$
(15)

###  3.2 Numerical experiment

We solve equations (14) (together with the boundary conditions) numerically
with a pseudo-spectral solver (again from the Dedalus code) similar to the one
used to solve equations (4)-(5) of the full $T-S$ model (similar spectral
expansions and time stepper). We run simulations of the reduced thermal
driving model for $Re_{\tau}=200,$ 400, 800 and $Pr_{*}=1$, 5, 10 (see the
list of simulation cases in table 3) and compare the results with the output
of the simulations of the full $T-S$ model (with similar control parameters).
The final checkpoint file of each thermal driving simulation was used to
initialize a new $T-S$ simulation (splitting the thermal driving field into
temperature and salinity fields according to equations (16)) with $Le=1$, 10
or 100 – i.e. distinct from the original $T-S$ simulations. This procedure was
designed to test the sensitivity of the $T-S$ simulations to the initial
conditions for temperature and salinity (uniform or from thermal driving at
statistical steady state). However, we did not observe any significant
sensitivity to initial conditions, such that we do not include these split
simulations in table 3.

###  3.3 Diagnostics

We assess the reliability of the reduced $T_{*}$ model by comparing the mean
temperature and salinity profiles, and the mean interfacial fluxes, diagnosed
from $\tilde{T}_{*}$ results, with those obtained in full $T-S$ simulations.
The reconstruction of the temperature and salinity from $\tilde{T}_{*}$
combines the definition of thermal driving and the MML (see details in
appendix) and yields for the diagnosed fields (denoted by superscript diag)

$\displaystyle\tilde{T}^{diag}=\frac{St\langle\overline{\tilde{T}_{*}}\rangle_{\perp}-\gamma}{St-\gamma},\quad\tilde{S}^{diag}=\frac{\langle\overline{\tilde{T}_{*}}\rangle_{\perp}-1}{St-\gamma}.$
(16)

The interfacial heat and salt fluxes are then diagnosed provided that we know
the effective (or turbulent) diffusivity $\tilde{\kappa}_{*}^{e}$
($\tilde{\kappa}_{*}^{t}$) of the $\tilde{T}_{*}$ field, defined as

$\displaystyle\tilde{\kappa}_{*}^{e}=\tilde{\kappa}_{*}^{t}+Pr_{*}^{-1}=\frac{\langle\overline{\tilde{w}\tilde{T}_{*}-Pr_{*}^{-1}\partial_{\tilde{z}}\tilde{T}_{*}}\rangle_{\perp}}{\langle-\overline{\partial_{\tilde{z}}\tilde{T}_{*}}\rangle_{\perp}}.$
(17)

If all scalar effective diffusivities are equal, i.e.
$\tilde{\kappa}_{T}^{e}\approx\tilde{\kappa}_{S}^{e}\approx\tilde{\kappa}_{*}^{e}$
(in the bulk), and the fluxes are depth invariant, we can indeed estimate for
a uniform ambient and at statistical steady state

$\displaystyle\tilde{Q}_{T}^{diag}(\tilde{z}=1)=-\tilde{\kappa}_{T}^{e}\langle\overline{\partial_{\tilde{z}}\tilde{T}^{diag}}\rangle_{bulk}=-\frac{\tilde{\kappa}_{*}^{e}St}{St-\gamma}\langle\overline{\partial_{\tilde{z}}\tilde{T}_{*}}\rangle_{bulk},$
(18a)
$\displaystyle\tilde{Q}_{S}^{diag}(\tilde{z}=1)=-\tilde{\kappa}_{S}^{e}\langle\overline{\partial_{\tilde{z}}\tilde{S}^{diag}}\rangle_{bulk}=-\frac{\tilde{\kappa}_{*}^{e}}{St-\gamma}\langle\overline{\partial_{\tilde{z}}\tilde{T}_{*}}\rangle_{bulk},$
(18b)

where
$\langle\cdot\rangle_{bulk}=\frac{1}{0.7-0.3}\int_{0.3}^{0.7}\langle\cdot\rangle_{\perp}dz$
is the bulk averaging operator. In other words, the heat and salt fluxes can
be diagnosed at the ice-ocean interface simply from $\tilde{T}_{*}$ and
$\tilde{\kappa}_{*}^{e}$ in the bulk, since we assume steady state. Note that
equations (18) are valid in the limit
$\overline{\tilde{S}}(\tilde{z}=1)\approx 0$ only (cf. appendix), such that
the exact dependence of the interfacial salinity on $Re_{\tau}$, $Pr$ and $Le$
(cf. figure 5) cannot be captured by the reduced model.

###  3.4 Results

Figures 7(a)-(b) show the mean temperature and salinity profiles with depth,
respectively, and demonstrate that the reduced $T_{*}$ model can reproduce
some of the key features of the full $T-S$ model. Specifically, the good
overlap of the mean temperature and salinity profiles obtained from $T-S$
simulation results (solid lines) with those diagnosed from $T_{*}$ simulation
results (dashed lines) in the bulk validates the MML theory. Figures 7(c)-(d)
zoom in on the ice-ocean diffusive sublayers and show that the
${\tilde{T}^{diag}}$ and ${\tilde{S}^{diag}}$ curves deviate in the boundary
layer from the true profiles obtained with the full $T-S$ model. This boundary
layer disagreement is expected since the MML theory only applies in the
turbulent bulk and illustrates why the MML prediction cannot be used alone to
estimate the interfacial heat and salt fluxes. In addition, the discrepancy
increases as $Le$ increases since $Le$ effects are missing in the reduced
model (as seen for $Le=100$ from the separation of the solid purple line and
dashed purple line, which overlaps with the brown dashed line, in figures
7(a)-(d)). We note that the thermal driving profiles (shown by the black
dotted lines in figures 7(a),(c)) are close to the temperature profiles
diagnosed from equation (16) in all our cases. This is because
$St\gg|\gamma|$, which should always be the case in the oceanographic context.
This overlap is observed not only in the bulk but also in the boundary layer
when $Pr_{*}=Pr$. However, we note that it is not a good idea to use
$\langle\overline{\partial_{\tilde{z}}\tilde{T}_{*}(\tilde{z}=1)}\rangle_{\perp}$
as a proxy for
$\langle\overline{\partial_{\tilde{z}}\tilde{T}(\tilde{z}=1)}\rangle_{\perp}$.
There is no theoretical basis and a more physically-motivated estimate is
readily available from equation (18).

(a)(b)(c)(d)

Figure 7: (a) Mean horizontal temperature and (b) salinity profiles. The thick
solid and dashed lines show results from full $T-S$ and reduced $T_{*}$ model
simulations (with $\tilde{T}$ and $\tilde{S}$ reconstructed from equations
(16)), respectively, while the dotted black lines in (a) show thermal driving
from the reduced $T_{*}$ model simulations. (c) and (d) are zoomed-in versions
of (a) and (b) with the distance from the ice-ocean interface
$(1-\tilde{z})^{+}=Re_{\tau}(1-\tilde{z})$ (in wall units and increasing from
top to bottom) replacing the depth coordinate $\tilde{z}$ on the $y$ axis.
Note that $(1-\tilde{z})^{+}=10$ corresponds to $\tilde{z}=0.95$ for
$Re_{\tau}=200$.

Figure 8(a) shows the mean interfacial fluxes obtained from the full $T-S$
simulation results (unfilled circles) and diagnosed from the $T_{*}$ model
results (unfilled squares) using equations (18). The overall agreement is
good. For $Le=10$, there is a perfect overlap of the circles and squares. For
higher $Le$, however, differences can be seen between the exact and diagnosed
fluxes. Specifically, the exact fluxes ($\overline{\tilde{Q}_{T}}$ and
$\overline{\tilde{Q}_{S}}$ at $\tilde{z}=1$) decrease with $Le$ while those
diagnosed are, by construction, $Le$ independent (cf. equation (18)). This
difference is highlighted in figure 8(b). The relative difference is typically
10% but jumps to 30% for the high $Le=100$ case. A similar result is obtained
for $\overline{\tilde{Q}_{S}}(\tilde{z}=1)St$ as it is equal to
$\overline{\tilde{Q}_{T}}(\tilde{z}=1)$ when $\overline{\tilde{S}}\ll 1$. We
have verified that the relatively coarse resolution of our high $Le$
simulations is not responsible for the flux discrepancy between the two models
(see Supplementary Material) and found that the $Le^{e}=1$ approximation is
well satisfied in all cases. Thus, we hypothesize that the decrease of the
exact fluxes with $Le$ in the full $T-S$ simulations must be related to the
molecular sublayer dynamics. The molecular sublayers are sensitive to the low
salt diffusivity and may throttle the vertical heat and salt fluxes, which are
quasi linearly correlated at the ice-ocean interface, in the relatively low
$Re_{\tau}$ and high $Le$ limit. The lack of a salt diffusivity bottleneck in
the $T_{*}$ model may put a fundamental limit on the range of $Re_{\tau}$ and
$Le$ parameters for which it can be used to diagnose the interfacial melt
rate. However, the relative error is only 30% for $Re_{\tau}=200$ and
$Le=100$, which is relatively close to the ocean value $Le\approx 170$. The
error may also decrease with higher (more realistic) $Re_{\tau}$ as the
boundary layer bottleneck is thought to vanish when turbulence increases (Iyer
et al., 2020).

(a)(b)

Figure 8: (a) Mean normalized heat flux as a function of the mean normalized
salinity flux at the ice-ocean interface obtained directly from full $T-S$
simulation results (unfilled circles) or diagnosed from $T_{*}$ model results
using equations (18) and a prescribed $Le$ (unfilled squares). The mean heat
and salinity fluxes computed from the gradient of the MML-reconstructed
temperature and salinity profiles (cf. (16)) are also shown (filled triangles)
and unsurprisingly are significantly different from the wall gradients
obtained from the full $T-S$ simulations. The dashed line highlights the
equation $\tilde{Q}_{T}=\tilde{Q}_{S}St$, which is the melt-induced dilution
equation in the limit of small interfacial salt deficit. (b) Difference
between the interfacial heat flux estimated from $T-S$ simulation results and
$T_{*}$ simulation results. The vertical bars show the standard deviation for
the heat flux based on the time series of the full $T-S$ simulation results
(note that the standard deviation measured from the time series of the reduced
$T_{*}$ simulations is typically smaller).

We note that our diagnosis of the interfacial heat and salinity fluxes from
$T_{*}$ model results does not allow to recover the salinity deficit at the
ice-ocean interface. Our derivation of $\tilde{Q}_{T}^{diag}(\tilde{z}=1)$
(cf. equation (18)) makes the approximation
$1+\overline{\tilde{S}}(\tilde{z}=1)\approx 1$, which is valid in the weak
thermal driving limit, such that the melt-induced dilution equation based on
the diagnosed fluxes yields systematically
$\overline{\tilde{S}}(\tilde{z}=1)=0$. The absence of a diagnosis for
$\overline{\tilde{S}}(\tilde{z}=1)$ implies that we cannot diagnose the
Nusselt number for the salinity field. However, we can approximate
$Nu_{T}^{diag}=Pr\tilde{Q}_{T}^{diag}(\tilde{z}=1)\approx Nu_{T_{*}}$
(assuming $Pr=Pr_{*}$) since
$\tilde{T}^{diag}\approx\langle\overline{\tilde{T}_{*}}\rangle_{\perp}$ when
$St\gg|\gamma|$ (reported in table 3). A diffusive sublayer model would be
required to diagnose $\overline{\tilde{S}}(\tilde{z}=1)$ but is beyond the
scope of the present work.

##  4 Conclusions

We have shown that the effective Lewis number $Le^{e}$ is close to unity in
direct numerical simulations of temperature- and salt-stratified sub-ice ocean
currents (§2). We limited our attention to ambient water masses with spatially
uniform and constant-in-time properties, and far-field temperatures close to
freezing (weak far-field thermal driving). Thus, we have provided supporting
evidence for the equal effective diffusivities hypothesis of the MML model in
the well-mixed regime.

Motivated by the observation $Le^{e}\approx 1$, we have shown how the results
from a reduced thermal driving model can be used to diagnose the temperature
and salinity profiles in the turbulent bulk and the heat and salt fluxes at
the ice-ocean interface (§3). The exact and diagnosed water properties are in
good agreement overall, but differ increasingly as the molecular Lewis number
$Le$ increases. The high-$Le$ discrepancy is due to the fact that the full
$T-S$ simulation results are sensitive to $Le$ while the diagnosed variables
are, by construction, $Le$ invariant. The neglect of the interfacial salinity
$\overline{\tilde{S}(\tilde{z}=1)}$ in the derivation of the diagnosed fluxes
(equation (18)) is responsible for the lack of $Le$ effects in the
reconstruction procedure (cf. appendix). Thus, the derivation of a boundary
layer model ($Le$ dependent) predicting the interfacial salinity
$\overline{\tilde{S}(\tilde{z}=1)}$ from the reduced model results would
improve the diagnoses. We note that the $Le$-induced discrepancy is expected
to vanish at high $Re_{\tau}$ (§33.4).

We listed in table 3 key output variables from our simulations. We reiterate
that the global quantities we reported, i.e., $Nu_{T}$, $Nu_{S}$ and $C_{D}$,
are artificially modified by the small horizontal dimensions of our numerical
domains (De Giovanetti et al., 2016). Thus, simulations in much larger domains
will be required to derive scaling laws for $Nu_{T}$, $Nu_{S}$ and $C_{D}$
relevant to the broad (unconfined) ocean. Our aim was not to revisit the
parameterization of ice melting in the IOBL in the shear-dominated regime
(Holland and Jenkins, 1999; Malyarenko et al., 2020), but rather to lay the
groundwork for reduced models that will be cost-effective for this
computationally intensive task.

Future work should examine the minimum far-field thermal driving above which
intermittent relaminarization breaks the equal effective diffusivities
hypothesis. Simulations with large far-field thermal driving will need to
consider large domains as relaminarization is artificially activated in
numerically-confined environments (García-Villalba and del Álamo, 2011).
Previous works suggest that domains with streamwise and spanwise dimensional
lengths $6\pi h$ and $2\pi h$, with $h$ the half channel width, or larger,
will be necessary to avoid confinement effects (Flores and Jiménez, 2010;
Bernardini et al., 2014). The effect of an oscillating current would be also
worth considering, as tides play an important role in several ice-shelf
cavities (Richter et al., 2022). A key question would be whether turbulence
survives long enough when tidal currents are low to ensure approximately equal
turbulent diffusivities at all times. We remark that some of our
approximations may not be valid at high thermal driving. Simulations with
relatively high melt rates might need to model deformations of the ice
boundary explicitly (Couston et al., 2021), or have non-zero velocities normal
to the interface, and consider a nonlinear equation of state to correctly
model density fluctuations in the presencen of large temperature and salinity
differences (Hester et al., 2021).

###### Acknowledgements.

This work was granted access to the HPC resources of IDRIS under the
allocation 2023-A0140114116 made by GENCI. We gratefully acknowledge support
from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon
for the computing resources. We will make publicly available (before the paper
is published) a Github repository containing all analysis scripts (required to
produce the figures) and a Zenodo repository containing both full three-
dimensional checkpoint data (at statistical steady state) and post processed
data. The software Dedalus is an open-source Python package available from the
Github repository https://github.com/DedalusProject/dedalus. Derivation of the
meltwater mixing line model for a stationary ice-ocean interface The MML
prediction can be derived from the advection-diffusion equations (4c)-(4d) for
$\tilde{T}$ and $\tilde{S}$ assuming a spatially uniform and constant-in-time
ambient. First, we integrate equations (4c)-(4d) over
$(\tilde{x},\tilde{y},\tilde{t})$, assume statistical steady state
($\partial_{\tilde{t}}\equiv 0$) and discard the relaxation terms, which
vanish in the bulk, to obtain

$\displaystyle\frac{d}{d\tilde{z}}\langle\overline{\tilde{w}\tilde{T}-Pr^{-1}\partial_{\tilde{z}}\tilde{T}}\rangle_{\perp}=0,$
(19a)
$\displaystyle\frac{d}{d\tilde{z}}\langle\overline{\tilde{w}\tilde{S}-Sc^{-1}\partial_{\tilde{z}}\tilde{S}}\rangle_{\perp}=0.$
(19b)

Then we substitute the definitions for the effective (turbulent) diffusivities
(11) into (19) and integrate from $\tilde{z}=1$ (ice-ocean interface) to
$\tilde{z}<1$ to get

$\displaystyle\tilde{\kappa}_{T}^{e}\langle\overline{\partial_{\tilde{z}}\tilde{T}}\rangle_{\perp}=Pr^{-1}\langle\overline{\partial_{\tilde{z}}\tilde{T}(\tilde{z}=1)}\rangle_{\perp},$
(20a)
$\displaystyle\tilde{\kappa}_{S}^{e}\langle\overline{\partial_{\tilde{z}}\tilde{S}}\rangle_{\perp}=Sc^{-1}\langle\overline{\partial_{\tilde{z}}\tilde{S}(\tilde{z}=1)}\rangle_{\perp}.$
(20b)

Finally, we take the ratio of equations (20a) and (20b), make the assumption
of equal effective diffusivities, and use the melt-induced dilution equation
(5) to arrive at the MML prediction, i.e.

$\displaystyle\frac{\langle\overline{\partial_{\tilde{z}}\tilde{T}}\rangle_{\perp}}{\langle\overline{\partial_{\tilde{z}}\tilde{S}}\rangle_{\perp}}=\frac{St}{1+\langle\overline{\tilde{S}(\tilde{z}=1)}\rangle_{\perp}}.$
(21)

Importantly, the right-hand-side in equation (21) can be simplified as $St$ in
the limit $\langle\overline{\tilde{S}(\tilde{z}=1)}\rangle_{\perp}\ll 1$
(small interfacial salt deficit), which is appropriate for all simulation
cases discussed in this paper. In this limit, integrating (21) from the far-
field $\tilde{z}=0$ upward yields a simple relation between the mean bulk
temperature and salinity, i.e.,

$\displaystyle\overline{\tilde{T}}-1=\overline{\tilde{S}}St.$ (22)

Combining (22) with the definition of thermal driving
$\tilde{T}_{*}=\tilde{T}-\gamma\tilde{S}$ yields expressions (16) for the bulk
temperature and salinity diagnosed from $\tilde{T}_{*}$.

The correlation coefficient $\mathcal{C}$ derived for a fixed ice-ocean
interface (right-hand-side of equation (21)) is slightly different from the
correlation coefficient derived for a moving interface (Gade, 1979). In the
latter case, re-doing the analysis in Gade (1979) with dimensionless variables
yields $\mathcal{C}\approx 1+St$, which is close to the result obtained for a
fixed interface in the limit $St\gg 1$. It is worth noting that our original
MML expression (21) predicts a correlation coefficient that depends on the
interfacial salinity, which itself is a function of the melt rate. This
dependence is not physical. It is absent from the MML prediction derived by
Gade (1979) and is due to the neglect of meltwater advection across the ice-
ocean boundary in our formulation (see a related discussion in section 4 of
Jenkins et al., 2001). However, its effect vanishes in the limit of weak
thermal driving, which is why all our simulations exhibit the same melt-
independent MML slope in figure 5.

## References

  * Alcántara-Ávila and Hoyas (2021) Alcántara-Ávila, F., and S. Hoyas, 2021: Direct numerical simulation of thermal channel flow for medium–high Prandtl numbers up to Re$\tau$=2000. International Journal of Heat and Mass Transfer, 176, 10.1016/j.ijheatmasstransfer.2021.121412.
  * Begeman et al. (2022) Begeman, C. B., X. Asay-Davis, and L. Van Roekel, 2022: Ice-shelf ocean boundary layer dynamics from large-eddy simulations. The Cryosphere, 16 (1), 277–295, 10.5194/tc-16-277-2022.
  * Bernardini et al. (2014) Bernardini, M., S. Pirozzoli, and P. Orlandi, 2014: Velocity statistics in turbulent channel flow up to Re$\tau$ =4000. Journal of Fluid Mechanics, 742, 171–191, 10.1017/jfm.2013.674.
  * Bouffard and Boegman (2013) Bouffard, D., and L. Boegman, 2013: A diapycnal diffusivity model for stratified environmental flows. Dynamics of Atmospheres and Oceans, 61-62, 14–34, 10.1016/j.dynatmoce.2013.02.002, URL http://dx.doi.org/10.1016/j.dynatmoce.2013.02.002.
  * Burns et al. (2020) Burns, K. J., G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, 2020: Dedalus: A flexible framework for numerical simulations with spectral methods. Phys. Rev. Research, 2 (2), 23 068, 10.1103/PhysRevResearch.2.023068, URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.023068.
  * Cenedese and Straneo (2023) Cenedese, C., and F. Straneo, 2023: Icebergs Melting. Annual Review of Fluid Mechanics, 55, 377–402, 10.1146/annurev-fluid-032522-100734.
  * Couston et al. (2021) Couston, L.-A., E. Hester, B. Favier, J. R. Taylor, P. R. Holland, and A. Jenkins, 2021: Topography generation by melting and freezing in a turbulent shear flow. Journal of Fluid Mechanics, 911 (A44), 1–37, 10.1017/jfm.2020.1064, 2004.09879.
  * Davis et al. (2023) Davis, P. E., and Coauthors, 2023: Suppressed basal melting in the eastern Thwaites Glacier grounding zone. Nature, 614 (7948), 479–485, 10.1038/s41586-022-05586-0.
  * De Giovanetti et al. (2016) De Giovanetti, M., Y. Hwang, and H. Choi, 2016: Skin-friction generation by attached eddies in turbulent channel flow. Journal of Fluid Mechanics, 808, 511–538, 10.1017/jfm.2016.665.
  * Dinniman et al. (2016) Dinniman, M. S., X. S. Asay-Davis, B. K. Galton-Fenzi, P. R. Holland, A. Jenkins, and R. Timmermann, 2016: Modeling ice shelf/ocean interaction in Antarctica: A review. Oceanography, 29 (4), 144–153, 10.5670/oceanog.2016.106.
  * Falor et al. (2023) Falor, D., B. Gayen, D. Sengupta, and G. N. Ivey, 2023: Evaporation induced convection enhances mixing in the upper ocean. Frontiers in Marine Science, 10 (May), 1–9, 10.3389/fmars.2023.1176226.
  * Flores and Jiménez (2010) Flores, O., and J. Jiménez, 2010: Hierarchy of minimal flow units in the logarithmic layer. Physics of Fluids, 22 (7), 1–4, 10.1063/1.3464157.
  * Gade (1979) Gade, H. G., 1979: Melting of Ice in Sea Water: A Primitive Model with Application to the Antarctic Ice Shelf and Icebergs. Journal of Physical Oceanography, 9 (1), 189–198, URL https://doi.org/10.1175/1520-0485(1979)009–“%˝3C0189:MOIISW–“%˝3E2.0.CO;2, arXiv:1011.1669v3.
  * García-Villalba and del Álamo (2011) García-Villalba, M., and J. C. del Álamo, 2011: Turbulence modification by stable stratification in channel flow. Physics of Fluids, 23 (4), 10.1063/1.3560359.
  * Gregg et al. (2018) Gregg, M. C., E. A. D’Asaro, J. J. Riley, and E. Kunze, 2018: Mixing efficiency in the ocean. Annual Review of Marine Science, 10, 443–473, 10.1146/annurev-marine-121916-063643.
  * Hester et al. (2021) Hester, E. W., C. D. McConnochie, C. Cenedese, L.-A. Couston, and G. Vasil, 2021: Aspect ratio affects iceberg melting. Physical Review Fluids, 6 (2), 23 802, 10.1103/PhysRevFluids.6.023802, URL https://doi.org/10.1103/PhysRevFluids.6.023802, 2009.10281.
  * Hewitt (2020) Hewitt, I. J., 2020: Subglacial Plumes. Annual Review of Fluid Mechanics, 52, 145–169, 10.1146/annurev-fluid-010719-060252.
  * Holland and Jenkins (1999) Holland, D. M., and A. Jenkins, 1999: Modeling thermodynamic ice-ocean interactions at the base of an ice shelf. Journal of Physical Oceanography, 29, 1787–1800, 10.1175/1520-0485(1999)029¡1787:mtioia¿2.0.co;2, URL Modelingthermodynamicice-oceaninteractionsatthebaseofaniceshelf.
  * Howland et al. (2023) Howland, C. J., R. Verzicco, and D. Lohse, 2023: Double-diffusive transport in multicomponent vertical convection. Physical Review Fluids, 8 (1), 1–20, 10.1103/PhysRevFluids.8.013501.
  * Hwang (2013) Hwang, Y., 2013: Near-wall turbulent fluctuations in the absence of wide outer motions. Journal of Fluid Mechanics, 723, 264–288, 10.1017/jfm.2013.133.
  * Iyer et al. (2020) Iyer, K. P., J. D. Scheel, J. Schumacher, and K. R. Sreenivasan, 2020: Classical 1/3 scaling of convection holds up to Ra = 1015. Proceedings of the National Academy of Sciences of the United States of America, 117 (14), 7594–7598, 10.1073/pnas.1922794117.
  * Jackson and Rehmann (2003) Jackson, P. R., and C. R. Rehmann, 2003: Laboratory Measurements of Differential Diffusion in a Diffusively Stable, Turbulent Flow. Journal of Physical Oceanography, 33 (8), 1592–1603, 10.1175/1520-0485(2003)033¡1592:lmoddi¿2.0.co;2.
  * Jenkins (1999) Jenkins, A., 1999: The impact of melting ice on ocean waters. Journal of Physical Oceanography, 29 (9), 2370–2381, 10.1175/1520-0485(1999)029¡2370:TIOMIO¿2.0.CO;2.
  * Jenkins (2011) Jenkins, A., 2011: Convection-driven melting near the grounding lines of ice shelves and tidewater glaciers. Journal of Physical Oceanography, 41 (12), 2279–2294, 10.1175/JPO-D-11-03.1.
  * Jenkins (2016) Jenkins, A., 2016: A simple model of the ice shelf-ocean boundary layer and current. Journal of Physical Oceanography, 46 (6), 1785–1803, 10.1175/JPO-D-15-0194.1.
  * Jenkins (2021) Jenkins, A., 2021: Shear, Stability and Mixing within the Ice-Shelf-Ocean Boundary Current. Journal of Physical Oceanography, 2129–2148, 10.1175/jpo-d-20-0096.1.
  * Jenkins et al. (2001) Jenkins, A., H. H. Hellmer, and D. M. Holland, 2001: The role of meltwater advection in the formulation of conservative boundary conditions at an Ice-Ocean interface. Journal of Physical Oceanography, 31 (1), 285–296, 10.1175/1520-0485(2001)031¡0285:TROMAI¿2.0.CO;2.
  * Jenkins et al. (2010) Jenkins, A., K. W. Nicholls, and H. F. Corr, 2010: Observation and parameterization of ablation at the base of Ronne Ice Ahelf, Antarctica. Journal of Physical Oceanography, 40 (10), 2298–2312, 10.1175/2010JPO4317.1.
  * Jiménez and Moin (1991) Jiménez, J., and P. Moin, 1991: The minimal flow unit in near-wall turbulence. Journal of Fluid Mechanics, 225, 213–240, 10.1017/S0022112091002033.
  * Kimura et al. (2016) Kimura, S., A. Jenkins, P. Dutrieux, A. Forryan, A. C. N. Garabato, and Y. Firing, 2016: Ocean mixing beneath Pine Island Glacier ice shelf, West Antarctica. Journal of Geophysical Research: Oceans, 121 (5), 3010–3028, 10.1002/2016JC012149.Received.
  * Kimura et al. (2015) Kimura, S., K. W. Nicholls, and E. Venables, 2015: Estimation of ice shelf melt rate in the presence of a thermohaline staircase. Journal of Physical Oceanography, 45 (1), 133–148, 10.1175/JPO-D-14-0106.1.
  * Lozano-Durán and Jiménez (2014) Lozano-Durán, A., and J. Jiménez, 2014: Effect of the computational domain on direct simulations of turbulent channels up to Re$\tau$ = 4200. Physics of Fluids, 26 (1), 10.1063/1.4862918.
  * Ma and Peltier (2022) Ma, Y., and W. R. Peltier, 2022: Diapycnal diffusivities in Kelvin-Helmholtz engendered turbulent mixing: The diffusive-convection regime in the Arctic Ocean. Journal of Fluid Mechanics, 946, 1–35, 10.1017/jfm.2022.590.
  * Malyarenko et al. (2020) Malyarenko, A., A. J. Wells, P. J. Langhorne, N. J. Robinson, M. J. M. Williams, and K. W. Nicholls, 2020: A synthesis of thermodynamic ablation at ice–ocean interfaces from theory, observations and models. Ocean Modelling, 154, 101 692, https://doi.org/10.1016/j.ocemod.2020.101692, URL https://doi.org/10.1016/j.ocemod.2020.101692.
  * Martin and Rehmann (2006) Martin, J. E., and C. R. Rehmann, 2006: Layering in a flow with diffusively stable temperature and salinity stratification. Journal of Physical Oceanography, 36 (7), 1457–1470, 10.1175/JPO2920.1.
  * McDougall et al. (2014) McDougall, T. J., P. M. Barker, R. Feistel, and B. K. Galton-Fenzi, 2014: Melting of ice and sea ice into seawater and frazil ice formation. Journal of Physical Oceanography, 44 (7), 1751–1775, 10.1175/JPO-D-13-0253.1.
  * Middleton et al. (2021) Middleton, L., C. A. Vreugdenhil, P. R. Holland, and J. R. Taylor, 2021: Numerical Simulations of Melt-Driven Double-Diffusive Fluxes in a Turbulent Boundary Layer beneath an Ice Shelf. Journal of Physical Oceanography, 51 (2), 403–418, 10.1175/jpo-d-20-0114.1.
  * Mondal et al. (2019) Mondal, M., B. Gayen, R. W. Griffiths, and R. C. Kerr, 2019: Ablation of sloping ice faces into polar seawater. Journal of Fluid Mechanics, 863, 545–571, 10.1017/jfm.2018.970.
  * Patmore et al. (2023) Patmore, R. D., P. R. Holland, C. A. Vreugdenhil, and A. Jenkins, 2023: Turbulence in the ice shelf-ocean boundary current and its sensitivity to model resolution. Journal of Physical Oceanography, 53, 613 – 633.
  * Pope (2000) Pope, S. B., 2000: Turbulent Flows. Cambridge University Press.
  * Reynolds (1975) Reynolds, A. J., 1975: The prediction of turbulent Prandtl and Schmidt numbers. International Journal of Heat and Mass Transfer, 18 (9), 1055–1069, 10.1016/0017-9310(75)90223-9.
  * Richter et al. (2022) Richter, O., D. E. Gwyther, M. A. King, and B. K. Galton-Fenzi, 2022: The impact of tides on Antarctic ice shelf melting. Cryosphere, 16 (4), 1409–1429, 10.5194/tc-16-1409-2022.
  * Rosevear et al. (2022a) Rosevear, M., B. Galton-Fenzi, and C. Stevens, 2022a: Evaluation of basal melting parameterisations using in situ ocean and melting observations from the Amery Ice Shelf, East Antarctica. Ocean Science, 18, 1109–1130, URL https://doi.org/10.5194/os-2021-111.
  * Rosevear et al. (2021) Rosevear, M. G., B. Gayen, and B. K. Galton-Fenzi, 2021: The role of double-diffusive convection in basal melting of Antarctic ice shelves. Proceedings of the National Academy of Sciences, 118 (6), 10.1073/pnas.2007541118, URL https://www.pnas.org/content/118/6/e2007541118.
  * Rosevear et al. (2022b) Rosevear, M. G., B. Gayen, and B. K. Galton-Fenzi, 2022b: Regimes and transitions in the basal melting of Antarctic ice shelves. Journal of Physical Oceanography, 52, 2589 – 2608, 10.1175/jpo-d-21-0317.1.
  * Smyth et al. (2005) Smyth, W. D., J. D. Nash, and J. N. Moum, 2005: Differential diffusion in breaking Kelvin-Helmholtz billows. Journal of Physical Oceanography, 35 (6), 1004–1022, 10.1175/JPO2739.1.
  * Stevens et al. (2020) Stevens, C., C. Hulbe, M. Brewer, C. Stewart, N. Robinson, C. Ohneiser, and S. Jendersie, 2020: Ocean mixing and heat transport processes observed under the Ross Ice Shelf control its basal melting. Proceedings of the National Academy of Sciences of the United States of America, 117 (29), 16 799–16 804, 10.1073/pnas.1910760117.
  * Vreugdenhil and Taylor (2019) Vreugdenhil, C. A., and J. R. Taylor, 2019: Stratification effects in the turbulent boundary layer beneath a melting ice shelf: Insights from resolved large-eddy simulations. Journal of Physical Oceanography, 49 (7), 1905–1925, 10.1175/JPO-D-18-0252.1.
  * Vreugdenhil et al. (2022) Vreugdenhil, C. A., J. R. Taylor, P. E. D. Davis, K. W. Nicholls, P. R. Holland, and A. Jenkins, 2022: The ocean boundary layer beneath Larsen C Ice Shelf: insights from large-eddy simulations with a near-wall model. Journal of Physical Oceanography, 10.1175/jpo-d-21-0166.1.
  * Zhou et al. (2017) Zhou, Q., J. R. Taylor, and C. P. Caulfield, 2017: Self-similar mixing in stratified plane Couette flow for varying Prandtl number. Journal of Fluid Mechanics, 820, 86–120, 10.1017/jfm.2017.200.

Supplementary Material

——————————–

Turbulent ice-ocean boundary layers in the well-mixed regime:
insights from direct numerical simulations
——————————–

Louis-Alexandre Coustona

aENSL, UCBL, CNRS, Laboratoire de physique, F-69342 Lyon, France

——————————–

In this Supplementary Material we provide details about the numerical
resolution of the simulations and show that the key results of the main text
are not sensitive to resolution changes.

We compare the grid spacing of our simulations with the turbulent Kolmogorov
and Batchelor length scales (Pope, 2000), defined as

$\displaystyle\tilde{\eta}_{K}=\frac{1}{\langle\overline{\tilde{\epsilon}}\rangle_{\perp,\tilde{z}>0.5}^{1/4}},\quad\tilde{\eta}_{B}=\frac{\tilde{\eta}_{K}}{Sc^{1/2}},$
(23)

The later is based on $Sc$ rather than $Pr$, such that it provides an estimate
of the smallest (salinity) features in the flow (Falor et al., 2023). Here
$\langle\overline{\tilde{\epsilon}}\rangle_{\perp,\tilde{z}>0.5}$ is the
(pseudo) turbulent kinetic energy dissipation rate averaged over the upper
half of the domain, i.e.

$\displaystyle\langle\overline{\tilde{\epsilon}}\rangle_{\perp,\tilde{z}>0.5}=\langle\overline{\frac{\partial\tilde{\bm{u}}^{\prime}}{\partial\tilde{x}}\cdot\frac{\partial\tilde{\bm{u}}^{\prime}}{\partial\tilde{x}}+\frac{\partial\tilde{\bm{u}}^{\prime}}{\partial\tilde{y}}\cdot\frac{\partial\tilde{\bm{u}}^{\prime}}{\partial\tilde{y}}+\frac{\partial\tilde{\bm{u}}^{\prime}}{\partial\tilde{z}}\cdot\frac{\partial\tilde{\bm{u}}^{\prime}}{\partial\tilde{z}}}\rangle_{\perp,\tilde{z}>0.5},$
(24)

where the prime denotes fluctuations relative to the horizontal mean. The
relative difference between
$\langle\overline{\tilde{\epsilon}}\rangle_{\perp,\tilde{z}>0.5}$ and the full
(averaged) turbulent kinetic energy dissipation rate is much less than
$10^{-2}$, hence negligible. Note that $\tilde{\eta}_{K}$ and
$\tilde{\eta}_{B}$ vary only by about 10% and always less than 30% when
replacing the volume averaging operator in (24) with another operator, such as
the horizontal average evaluated at the depth of maximum turbulent kinetic
energy.

To further asses whether our simulations are adequately resolved we count the
number of grid points lying within the momentum, thermal and haline boundary
layers. We define the lower edge of each boundary layer as the position where
the horizontally-averaged diffusive and convective fluxes become equal.
Mathematically, we obtain the thickness of the momentum, thermal and haline
boundary layers, which we denote by $\tilde{\delta}_{u}$, $\tilde{\delta}_{T}$
and $\tilde{\delta}_{S}$, respectively, by ensuring that they satisfy

$\displaystyle\langle\overline{\tilde{w}\tilde{u}}\rangle_{\perp}=$
$\displaystyle-\langle\overline{\partial_{\tilde{z}}\tilde{u}}\rangle_{\perp},\quad$
$\displaystyle\text{at}\quad\tilde{z}=1-\tilde{\delta}_{u},$ (25a)
$\displaystyle\langle\overline{\tilde{w}\tilde{T}}\rangle_{\perp}=$
$\displaystyle-
Pr^{-1}\langle\overline{\partial_{\tilde{z}}\tilde{T}}\rangle_{\perp},\quad$
$\displaystyle\text{at}\quad\tilde{z}=1-\tilde{\delta}_{T},$ (25b)
$\displaystyle\langle\overline{\tilde{w}\tilde{S}}\rangle_{\perp}=$
$\displaystyle-
Sc^{-1}\langle\overline{\partial_{\tilde{z}}\tilde{S}}\rangle_{\perp},\quad$
$\displaystyle\text{at}\quad\tilde{z}=1-\tilde{\delta}_{S}.$ (25c)

As in the main text we denote by subscript $+$ variables pre-multiplied by
$Re_{\tau}$ (i.e. expressed in terms of wall units).

Figure/Table 9 shows for each Full $T-S$ simulation discussed in the main text
the grid resolution (without the 3/2 dealiasing factor), the friction
Kolmogorov and Batchelor length scales, the grid spacing normalized by the
Batchelor length scale, the boundary layer thicknesses, the number of grid
points $\mathcal{N}(Sc^{-1/2})$ and $\mathcal{N}(\tilde{\delta}_{S}^{+})$
satisfying $1-\tilde{z}^{+}<Sc^{-1/2}$ and
$1-\tilde{z}^{+}<\tilde{\delta}_{S}^{+}$, respectively, and the temporal range
of time-averaged statistics (in friction time units).

Figure 9: Table comparing the grid spacing against the smallest length scales
of the flow for all Full $T-S$ simulations discussed in the main text. Note
that the typical time step $Re_{\tau}d\tilde{t}$ varies between $5\times
10^{-5}$ and $6\times 10^{-4}$ across all simulations.

The Chebyshev collocation points cluster sufficiently close to boundaries for
the first grid point to always lie within the distance $1/(Re_{\tau}Sc^{1/2})$
of the ice-ocean interface and for about 10 or more points to lie within the
haline molecular sublayer (6th and 7th columns). However, the fourth column
clearly shows that our simulations only coarsely resolve the Batchelor scale
in the bulk. Most notably, the simulations with the highest Lewis number have
grid spacings in $\tilde{x}$ and $\tilde{z}$ that are 7.7 times the Batchelor
length scale, exceeding by a factor $\sim 2$ the generally-accepted criterion
for resolution requirement (Falor et al., 2023). We note that coarsely
resolving the Batchelor length scale is not uncommon (Middleton et al., 2021).
In fact little kinetic energy is expected at and below the Kolmogorov length
scale such that numerical errors can be expected to remain at small scales.

(a)(b)

Figure 10: Time series of (a) the mean streamwise flow and (b) the normalized
heat flux for the same set of physical parameters, i.e., $Re_{\tau}=200$,
$Pr=1$, $Le=100$ but five different spectral resolutions (see table 11). The
simulations with higher resolution than the reference simulation (solid
magenta) start from the reference simulation’s last checkpoint at time
$Re_{\tau}\tilde{t}=40$ (vertical dotted line). The horizontal lines show the
mean value over the temporal window used for time averaging. Figure 11: Table
showing simulation results with $Re_{\tau}=200$, $Pr=1$, $Le=100$ for a range
of spectral resolutions. The results of the reference simulation, discussed in
the main text, are on the second line ($n_{x},n_{y},n_{z}$=256,256,192).
$\sigma$ denotes the standard deviation in time.

While our simulations do not visibly exhibit non-physical behaviours despite
the coarse resolution, except numerical ringing occasionally, we investigate
the impact of changing the resolution on a handful of statistics for one of
the simulations with the coarsest (relative) resolution, i.e. $Re_{\tau}=200$,
$Pr=1$, $Le=100$. Figure 10 shows the mean flow rate and normalized heat flux
at the ice-ocean interface over time for five different resolutions (cf.
figure/table 11). The curves for the lower resolution and reference
simulations overlap almost perfectly from $Re_{\tau}\tilde{t}=0$ to
$Re_{\tau}\tilde{t}\approx$5, while those at higher resolution (restarted from
the reference simulation at time $Re_{\tau}\tilde{t}=40$) overlap from
$Re_{\tau}\tilde{t}=40$ to $Re_{\tau}\tilde{t}\approx$45\. At later times the
results diverge but the mean value remains close to or within the standard
deviation (in time) of the reference simulation’s result. The relatively
negligible impact of resolution on first order statistics is generalized in
figure/table 11. The mean value of all key statistics is close to or within
the standard deviation of the reference simulation’s result in all cases. Most
importantly, the turbulent Lewis number is always close to 1.

We conclude that the key results of the main text, including $Le^{e}\approx
1$, are not sensitive to resolution changes. Some of the first order
statistics may vary by up to one standard deviation (typically 10% or less).
However, similar changes may be obtained with longer time averaging windows as
most of our simulations are relatively short. Having first order statistics
converged to a sufficiently high level of accuracy to derive scaling laws is
beyond the scope of the manuscript and reserved for future work.