Improved Parameterizations of Vertical Ice‐Ocean Boundary Layers and Melt Rates

Buoyancy fluxes and submarine melt rates at vertical ice‐ocean interfaces are commonly parameterized using theories derived for unbounded free plumes. A Large Eddy Simulation is used to analyze the disparate dynamics of free plumes and wall‐bounded plumes; the distinctions between the two are supported by recent theoretical and experimental results. Modifications to parameterizations consistent with these simulations are tested and compared to results from numerical and laboratory experiments of meltwater plumes. These modifications include 50% weaker entrainment and a distinct plume‐driven friction velocity in the shear boundary layer up to 8 times greater than the externally‐driven friction velocity. Using these updated plume parameter modifications leads to 40 times the ambient melt rate predicted by commonly used parameterizations at vertical glacier faces, which is consistent with observed melt rates at LeConte Glacier, Alaska.


Introduction
Recent warming of deep ocean currents that come into contact with the termini of tidewater glaciers has led to enhanced submarine melt (Cowton et al., 2018;Holland et al., 2008;Straneo & Heimbach, 2013;Wood et al., 2018).This is postulated to have led to a speedup of marine-terminating glaciers outflowing at the margins of the Greenland Ice Sheet and Antarctic Ice Sheet (van den Broeke et al., 2016).
At vertical or near-vertical glacier faces, submarine melt is primarily driven by a combination of three dynamical processes: subglacial discharge plumes, ambient melt plumes, and horizontal circulation (Jackson et al., 2020;Straneo & Cenedese, 2015).The first two types of melt are driven by buoyant plume convection of different strengths with vertical velocities reaching 2-3 m/s for subglacial discharge plumes and up to 10 cm/s for melt plumes (Jackson et al., 2020).The horizontal near-glacier velocity is driven by vorticity generation from water mass transformation and has a magnitude of up to tens of cm/s, and varies significantly between fjords and seasons with limited direct observations (Jackson et al., 2020;Sutherland et al., 2014;Straneo & Cenedese, 2015;Zhao et al., 2021Zhao et al., , 2022)).Although subglacial discharge plumes have the potential to drive the fastest melt rates locally, they are often observed to occupy a small fraction of the glacial face, while the other two melt processes occur across the entire glacial face (Cowton et al., 2015;Slater et al., 2018).
Recent studies have discussed the relative importance of these three melt processes (Jackson et al., 2020;Slater et al., 2018) and most studies using parameterizations with standard parameter values predict ambient melt rates outside of subglacial discharge plumes to be low (much less than a meter/day, often cms per day; Carroll et al., 2016;Fried et al., 2015;Zhao et al., 2021).However, there is a mismatch between recent observations and these parameterizations; measured ambient melt rates at LeConte Glacier, Alaska are an order of magnitude greater (1-10 m/day) across the entire submarine terminus, even in parts of the glacier face far from discharge plumes (Jackson et al., 2020;Sutherland et al., 2019).In addition to the discrepancy for ambient melt rates, because discharge plumes often cover a small fraction of the total glacier area, the face-averaged and total observed melt rates are also much higher than those predicted by parameterizations with standard parameter values.
In this study, we extend results from a recently proposed parameterization for vertical glacial ice fronts, which proposed modifying unbounded plume theory using empirical constraints for the efficiency of turbulent heat and salt transfer to match observational data (Schulz et al., 2022).Schulz et al. (2022) also proposed a transfer function that merges the velocity-dependent (shear-dominated) and velocity-independent (buoyancy-dominated) melt regimes, based on best empirical fits to observations.This results in a significantly higher baseline buoyancy-dominated melt rate (in the absence of external flow) than previous literature (e.g., Kerr and McConnochie (2015)) and does not provide a physical explanation why these baseline melt rates should be higher.In this study, we propose a physically-motivated melt parameterization that includes both buoyancydriven and shear-dominated melt regimes that is consistent with existing theories, observations, and laboratory experiments.
In Section 2, we present an updated and integrated overview of free plumes, wall-bounded plumes, and horizontal circulation-driven melt, and how each drive the boundary layer dynamics and melt at a vertical ice face.In Section 3, we present a set of Large Eddy Simulations of a subglacial discharge plume with and without a vertical glacier wall to compare the horizontal and vertical profiles of vertical momentum for unbounded free plumes and wall-bounded plumes.This is compared with existing theories for discharge plumes.In Section 4, we compare the parameterizations with standard parameter values of glacial melt rates at a rapidly melting vertical ice face (LeConte, Alaska) with the updated melt plume theory (Section 2) and discharge plume theory (Section 3).This shows that wall-bounded plume theory (after accounting for buoyancy due to melt from horizontal circulation) is consistent with recent observations, while the commonly used free plume theory underpredicts the melt rate outside of discharge plumes by a factor of 40.Finally, we summarize our key findings and proposed changes to vertical ice-ocean interface parameterizations and conclude.

Theory of Vertical Ice-Ocean Interfaces
In this section, we summarize and integrate recent developments in vertical ice-ocean boundary layer parameterizations by first discussing the thermodynamic coupling of the interfacial boundary layer to the corresponding (plume-or external forcing-driven) outer velocity, temperature, and salinity.We then discuss how these outer properties are parameterized for each of the three types of outer boundary layers: subglacial discharge plumes, ambient face-wide melt plumes, and background/external circulation.We refer the reader to recent reviews of glacial plumes and ice-ocean parameterizations for further details and references (Hewitt, 2020;Malyarenko et al., 2020).

Vertical Ice-Ocean Boundary Layers
Current melt rate parameterizations at vertical ice-ocean interfaces can be classified into being relevant in either the buoyancy-driven regime (Kerr & McConnochie, 2015) or the shear-driven regime (Jenkins, 2011;McPhee et al., 2008) based on whether the rate of turbulent heat flux is primarily driven or constrained by buoyancy flux diffusing away from the wall (buoyancy-driven) or the momentum flux diffusing toward the wall (shear-driven); a transition from the first to the second regime occurs if the buoyant updraft has gained significant vertical momentum (Wells & Worster, 2008).In the absence of externally-forced circulation or turbulence, vertical ice-ocean interfaces start off as wholly laminar boundary layers within the first 10-30 cm above their initiation point before they transition to the buoyancy-driven turbulent regime (Josberger & Martin, 1981;Wells & Worster, 2008); however, we do not discuss the laminar regime further due to its limited relevance to the geophysical scale of glaciers.
Within buoyancy-driven boundary layers, the melt rate is velocity-independent and can be approximated as (Kerr & McConnochie, 2015) where c B = 0.25 μm °C 4/3 is an empirically-derived dimensional constant, T is the ambient temperature, T f is the local freezing temperature (which can be calculated using the liquidus condition similarly to Equation S7c in Supporting Information S1) at ambient salinity S. The buoyancy Rayleigh number is defined as Ra = b(z z 0 ) 3 / (νκ), which represents the plume's convective efficiency with respect to diffusion with height from the source (z z 0 ), where buoyancy b = g(ρ ρ a )/ρ 0 (for plume density ρ, an ambient density ρ a , and reference density ρ 0 ), ν = 1.8 × 10 6 m 2 /s is the viscosity and κ = 7.2 × 10 10 m 2 /s the salt diffusivity of seawater.
For shear boundary layers, and Here, c w = 3974 J kg/°C and c i = 2009 J kg/°C are the specific heat capacity of water and ice, respectively, L = 3.35 × 10 5 J/kg is the latent heat of ice, γ T is the turbulent thermal transfer coefficient (with units of velocity), and T b is the boundary layer temperature predicted by solving the 3-equation thermodynamical balance (see Jenkins et al. (2010), Equations 7a-7c in Supporting Information S1).The turbulent transfer coefficient is dependent on the friction velocity and is discussed in the next subsection.
Here, the threshold between buoyancy-driven and shear-driven boundary layers is set by a critical value of the buoyancy Rayleigh number.The critical Rayleigh number Ra c for the transition and its corresponding transition height z c z 0 is the subject of some debate (Grossmann & Lohse, 2000;Wells & Worster, 2008), partly due to the fact that this transition has not been observed in a natural setting, but it is postulated to occur at Ra c = 10 21 (Kerr & McConnochie, 2015).Recent laboratory experiments suggest that this occurs at a vertical velocity of 0.03-0.05m/s for a discharge plume at 3.5 o C above freezing (McConnochie & Kerr, 2017a).A simple way to combine the two regimes in Equations 1 and 2 is to use a melt rate prediction based on the dominant turbulent transfer process at the boundary layer, resulting in m = max{m S ,m B }. (3)

Shear-Driven Turbulent Transfer
The turbulent transfer coefficient in a shear-dominated regime is commonly expressed in terms of horizontal along-glacier (v) and vertical near-glacier (w) ocean velocities as with a drag coefficient C d ranging from 0.001 to 0.0097.See Figure 1a for a list of references for the quoted drag coefficient range and a schematic of the different boundary layers and corresponding velocities.In the absence of boundary layer observations, a commonly used value of C d = 0.0025 is used in the ice plume literature along with turbulent heat and salinity transfer constants of Γ T = 0.022 and Γ S = 6.2 × 10 4 (Jenkins et al., 2010).However, this does not distinguish between the frictional boundary layer thickness (via C d ) in the horizontal and vertical and more importantly, it does not distinguish between the external velocity-field driven shear boundary layers and the plume-driven boundary layers.For melt plumes and discharge plumes, it is also unclear how v and w should be defined as both far-field velocities for plumes are zero.
The total shear stress at a shear boundary layer is the sum of both the viscous and turbulent shear stresses (5) In most externally-forced wall-bounded shear flows (in either the atmospheric boundary layer or horizontal iceocean boundary layers; Jenkins (1991); Kaimal and Finnigan (1994); Pope (2000)), the turbulent Reynolds stress dominates the momentum dissipation contribution.However, recent laboratory and numerical experiments suggest that plume-driven buoyancy forcing at an ice-ocean interface behaves differently than the external far field-forced velocity field.This is because buoyancy (from melting) is generated directly at the interface itself in melt plumes or close to the wall in the case of subglacial discharge plumes (Gayen et al., 2016;Parker et al., 2020Parker et al., , 2021)).Therefore, it is important to distinguish shear stresses associated with the external velocity field from those of the internal plume-driven shear stresses.For both melt plumes and discharge plumes, more of the shear stress contribution is viscous in Equation 5and thus, more of the kinetic energy is dissipated before becoming turbulent.
However, this has been demonstrated to lead to a melt rate that scales strongly with the friction velocity of the shear boundary layer (Gayen et al., 2016;McConnochie & Kerr 2017b;Parker et al., 2020Parker et al., , 2021)).In order to separate the individual contributions of plume-driven shear and externally-driven shear, we express the turbulent thermal and salinity transfer coefficient as where the plume-driven friction velocity is defined such that w 2 * = ν∂ x w p ⃒ ⃒ x=0 at a vertical wall x = 0 (for x in the horizontal and normal direction from the glacier), for a time-averaged plume vertical velocity w p (which assumes all of the viscous stress is converted to turbulent stress).Equation 6 is a simple way of combining the horizontal and vertical components of friction velocity via the velocity magnitude, which is commonly used when there is a 2D external velocity field (Jenkins, 2011), but such an expression has not been validated for combining plumedriven friction velocity and externally-driven friction velocity (see McConnochie and Kerr (2017b) for further discussion).
In previous studies, the plume-driven shear boundary layer is often expressed using an equivalent skin friction coefficient C p d ≡ w 2 * / W 2 p with empirically derived estimates that are significantly higher than its analogouslydefined externally-forced counterpart  Parker et al. (2020Parker et al. ( , 2021))).The characteristic plume velocity W p used in this parameterization is defined as the horizontally-integrated mass flux divided by momentum flux (see Supporting Information S1 Text S1.2 for further discussion).

1D Ice-Ocean Plume Theory
Buoyant plume convection in the absence of a wall is an extensively studied subject (Morton et al., 1956;Turner, 1979).A 1D theory for the vertical (along-plume) variation in characteristic vertical velocity W (z) = W p (z), buoyancy B(z), and plume width D(z) can be solved using the Boussinesq conservation laws of mass, momentum, and buoyancy, and empirically-derived entrainment assumption (see Text S1 in Supporting Information S1 for further details).Here, W(z) and B(z) are defined as the mean of the vertical velocity w p and buoyancy b = g(ρ ρ a )/ρ 0 over the plume width for a plume density ρ, an ambient density ρ a , and reference density ρ 0 .This theory relies on the empirically-supported assumption originating from Morton et al. (1956), that the local time-mean entrainment at each depth is proportional to the characteristic vertical velocity U(z) = αW(z) for a horizontal inflow velocity magnitude U(z) and constant entrainment coefficient α.
For plumes at an ice-ocean interface, this has been modified to include an accounting of the heat budget and drag at the wall, which are necessary for calculating the melt rate in Equation 1 or 2. Therefore, we have the following commonly-used system of equations for 1D line plume evolution originating from Jenkins (1991).
where m is the melt rate, T a and S a are the ambient temperature and salinity, respectively, S i is the ice salinity and each of the parenthetical terms are integrated numerically in z to generate the 1D plume solution (see Text S1 in Supporting Information S1 and Hewitt (2020) for a derivation).Note that m and the effective thermal gradient T ef (i.e., the boundary layer temperature and salinity) must be calculated using 3-equation thermodynamics (see Text S1 in Supporting Information S1).This 1D plume theory has been used extensively in the parameterizations of vertical ice faces (Cowton et al., 2015;Jackson et al., 2017Jackson et al., , 2020;;Jenkins, 2011;Slater et al., 2018;Sutherland et al., 2019;Zhao et al., 2021) due to its simplicity.
Note that this system of equations is valid for both discharge plumes (for an appropriate initial mass and momentum) and melt plumes (for approximately zero initial mass and momentum).These equations can be modified slightly to describe a point source (by replacing plume width D with plume radius R and deriving the appropriate conservation laws for a radial symmetry).However, line discharge plumes (a finite width buoyancy source instead of a point source) have been shown to better reproduce existing near glacier melt fraction observations, but this may be more attributed to the discharge source width instead of the dynamics (Jackson et al., 2017).

Distinction Between Unbounded Free Plumes and Wall-Bounded Plumes
Free plumes do not have an associated skin friction or melt terms in Equations 7a-7d, but otherwise have the same system of 1D line plume equations.For free plumes, the commonly used velocity profile w p (x) has been determined experimentally and is well-characterized by a Gaussian curve (Ramaprian & Chandrasekhara 1989;Paillat & Kaminski 2014; suggestive of a random walk of water mass parcels) where x = x/ (z z 0 ) is the z-scaled x coordinate, such that the plume profiles w( x)/ W max (z) collapse to a single characteristic profile by similarity with height from the source (z z 0 ) and the maximum plume velocity relative to the characteristic vertical velocity is W max /W ≈ 1.35.
The momentum flux equation in Equation 7b includes a skin friction term at the wall, C p d .This term owes its commonly used value of 0.0025 in ice-ocean applications largely to observations at weakly melting nearly horizontal ice-ocean interfaces (Jenkins et al., 2010;McPhee et al., 1987), which have boundary layers that can be well-approximated by Monin-Obukhov theory (Vreugdenhil & Taylor, 2019); this value of skin friction is also commonly used in most other passive surfaces.
The skin friction coefficient C p d can be diagnosed experimentally or numerically via the balance of momentum terms in Equation 7b.In a wall-bounded plume with a shear layer, the across-plume gradients determined experimentally and numerically differ greatly from the Gaussian-shaped vertical momentum profiles of free plumes and have different distributions of horizontal turbulent fluxes of momentum and buoyancy (Parker et al., 2020(Parker et al., , 2021;;Sangras et al., 2000).
These disparities are owed to two major differences between a plume convecting along a wall versus one without a wall: the impermeability condition (u at the wall is zero) and the no-slip condition (w at the wall is zero).Experiments have shown that the impermeability condition leads to reduced eddy meandering and weaker mixing of buoyant fluid away from the wall.This produces higher near-wall vertical momentum, which together with the no-slip condition contribute to significantly higher shear stresses (Parker et al., 2020(Parker et al., , 2021)).For discharge plumes, the shear stress diagnosed from laboratory and numerical experiments exerts a drag on momentum equivalent to 15% of the buoyancy force for discharge plumes (Ezhova et al., 2018;Parker et al., 2020) and 65% of the buoyancy force for melt plumes in small (1 m tall), unstratified domain heights (Gayen et al., 2016;Parker et al., 2021).These experimentally-derived estimates imply a drag coefficient of C p d = 0.015 for discharge plumes, and C p d = 0.15 for melt plumes along with significantly lower entrainment: α = 0.075 for discharge plumes, and α = 0.068 for melt plumes.See Figure 1b for a list of the drag coefficients, entrainment, and a corresponding list of references.
A significant body of experimental and theoretical work supports an across-plume vertical velocity profile for a freely convecting, heated wall, which was first approximated in Eckert and Jackson (1950) as This velocity profile also approximately matches recent experiments of ice-ocean boundary layers (Parker et al., 2021).However, it is unknown how across-plume vertical velocity profiles behave for much larger Ra and rise heights.

Horizontal Circulation-Driven Melt
The near-glacier horizontal velocity has a very different profile v(x) compared to the across-plume vertical velocity profile w(x) (see Figure 1a), where the commonly used parameterization uses the far-field background velocity v(∞) (at 10-100 m away from the boundary) and a drag coefficient of C ext d = 0.0025 (valid for a meterscale boundary layer in which the law of the wall assumption holds; see Davis and Nicholls (2019)).By comparison, plume-driven shear boundary layers are much thinner (centimeters or less) and they have proportionally larger friction velocities compared to the outer velocity The far-field velocity v ∞ may either be observed directly (estimated to be 20 cm/s near the face of LeConte (Jackson et al., 2020)) or parameterized using the theory from Zhao et al. (2021Zhao et al. ( , 2022Zhao et al. ( , 2023)).

Large Eddy Simulations of Subglacial Discharge Plumes
Although recent experimental and numerical studies (discussed in Section 2.3, Gayen et al. (2016); Parker et al. (2020Parker et al. ( , 2021))) demonstrated a larger wall shear stress and weaker entrainment in wall-bounded discharge plumes compared to free plumes, those experiments were not able to diagnose the relative importance of the two wall effects: impermeability and the no-slip boundary condition.In the following experiments, we test the importance of these two effects separately.We examine the horizontal profiles and vertical acceleration of vertical velocity to reconcile the differences (particularly how it is treated in the glacial context) between the theory of unbounded free plumes and wall-bounded plumes.

Model Setup
To examine the difference between a wall-bounded plume and a free plume, we conduct a series of Large Eddy Simulations (LES) of an idealized near-glacier fjord domain using the Massachusetts Institute of Technology general circulation model (MITgcm; Marshall et al. (1997)).The vertical and horizontal resolution are both 1 m and we use a 3D Smagorinsky viscosity parameterization with a coefficient of 0.03 (Smagorinsky, 1963).The model is forced on the open-ocean side by an idealized temperature/salinity (Figure 2b) based on August 2016 observations at LeConte Glacier, Alaska (Jackson et al., 2020).For the cases with a vertical wall, we parameterized melting at the ice face using a shear boundary layer assumption and 3-equation thermodynamics (Equations 2 and 4 and Equations S7a-S7c in Supporting Information S1).We use an idealized bathtub domain with smooth sidewalls (see Figure 2a for the bathymetric variation in y) for a 200 m deep, 1 km wide (in y) by 2 km long (in x) fjord section with a subglacial discharge of 150 m 3 /s (based on August 2016 observations; see Jackson et al. (2020) and Sutherland et al. (2019)).To help initiate turbulent motions near the source, the mass source is distributed at 10 evenly-spaced outlets over a 100 m extent in y (450 ≤ y ≤ 550 m) at z = 200 m over 0 < x < 3 m, which provides a small degree of along-y asymmetry in vertical momentum without any meaningful influence on the y-integrated vertical momentum.The simulations were run for 2 hr (after reaching steady state at approximately 15 min).The length of this fjord section minimizes the interaction between the plume-generated turbulence and the open ocean boundary at x = 2,000 m.See Text S2 in Supporting Information S1 for additional details of the numerical model.

Figures 2a-2d
show the near-face (x = 3 m) vertical velocity, w = 0.2 m/s velocity surface, and meridionallyaveraged vertical velocity for a reference case of a wall-bounded discharge plume at an ice-ocean boundary without drag.This shows the development of 3D convective turbulence, which develops due to xz-plane vorticity sourced at y = 450, 550 m on the margins of the line plume and throughout the plume due to vorticity aligned in the yz-plane; these regions correspond to high horizontal shear in vertical velocity leading to turbulent shear production.Although a finer grid would better resolve the smaller scales of turbulence at depths especially near the plume source, the turbulence in the upper half of the plume is sufficiently well-developed (evidenced by the constant entrainment coefficient with depth above z = 140 m) to calculate across-plume profiles of vertical velocity.Figure 2d also shows the mean meridionally-averaged density field, whose positive buoyancy anomaly within the plume is the source of the buoyancy flux driving upward acceleration.This also shows the gradual decrease in density along the plume due primarily to the entrainment of ambient water.
To compare the characteristics of free plumes and wall-bounded plumes, we examine three test cases.The first case is a free plume, the second case is a wall plume without drag at the wall, and the third case is a wall plume with a drag coefficient of C d = 0.015 to emulate the no-slip condition (which is parameterized and not resolved).
Note that this value of drag coefficient is determined for discharge plumes (which is distinctly lower than the effective drag felt by a melt plume) and is obtained experimentally (see Section 2b and Figure 1b).
Figure 3a shows the horizontal variation of vertical velocity w(x) and Figure 3b shows the vertical variation of characteristic vertical velocity W(z) for each of the three test cases.These demonstrate that the dynamics are consistent with their respective theories; free plumes w p (x) are well-approximated by the Gaussian profiles in Equation 8 and wall-plume profiles are consistent with the turbulent wall-plume theory in Equation 9.These functions were fit using α as a free parameter (which determines the characteristic width scale D of the plume for a given z z 0 ).For the free plume case, this implies α = 0.14; for the wall-plume case without drag, α = 0.083; for the wall-plume case with drag, α = 0.079.The wall bounded plume simulation is consistent with theory except close to the wall, which has a thin (<1 m) boundary layer that is not resolved.Therefore, the closest 3 gridpoints (3 m) were excluded from the estimate of entrainment.These diagnosed values are consistent with those from previous studies (Figure 1b).In addition to lower entrainment coefficients, these profiles of vertical velocity are also similar to the free and wall plume profiles, respectively, from Kimura et al. (2014).
Figure 3b shows the vertical variation of the characteristic vertical velocity W(z) along with the corresponding theoretical solutions from 1D buoyant plume theory (see Equation 7).This comparison demonstrates that the bulk mean vertical momentum W is consistently well-predicted by plume theory in the MITgcm simulations.In particular, 1D plume theory captures the 17% increase in W for the wall-bounded drag-free plume (due to reducing the entrainment coefficient from 0.14 to 0.083) and the 5% decrease in W when drag is added (primarily due to adjustments to the momentum balance and little influence on the buoyancy balance).Note that the characteristic width of the wall-bounded plume is much narrower (in Figure 3a), partially due to weaker entrainment for the wall-bounded plumes.A notable caveat here is that near the source, the plume is not fully resolved at the 1 m horizontal resolution used in the MITgcm simulations and likely contributes to the small mismatch in vertical momentum there.Near the depth of neutral buoyancy in the MITgcm simulations, the plumes start to decelerate below the pycnocline, which is partially due to the emergence of an inner versus outer plume region (likely caused by additional sources of mixing/instabilities not captured by the theory).The inner plume feels a weaker buoyancy anomaly relative to the outer plume compared to the far field stratification.
In summary, these results show that the presence of a vertical wall strongly alters plume dynamics.However, we demonstrate that parameterizations with standard parameter values (based on theory for unbounded/free plumes), can be adapted to produce the observed variability if a lower entrainment coefficient is used.This lower entrainment is consistent with results from previous wall-plume simulations (Ezhova et al., 2018;Kimura et al., 2014).In addition, the critical difference between the entrainment of free and wall-bounded plumes emerges primarily due to the impermeability condition and to a lesser extent, the no-slip condition, although the latter is implicit in the shear boundary layer parameterization (just not as important for the vertical momentum balance).One caveat of these experiments is that the near-wall horizontal resolution in these LES does not allow the near-glacier plume-driven boundary layer to be fully resolved, which appears to be much more important in melt plumes where the wall shear stress decreases the total plume momentum by approx.65% in both laboratory experiments and Direct Numerical Simulations (DNS) (Gayen et al., 2016;Parker et al., 2021).The resolution and computational cost required to resolve the laminar boundary layer (in e.g., Gayen et al. (2016)) would not currently be affordable at vertical scales larger than a few meters.Therefore, we cannot simultaneously simulate the large-scale discharge plume dynamics while resolving the dynamics of the melt plumes and viscous/diffusive boundary layers.Instead, the drag coefficient and entrainment rates from small-scale laboratory experiments and DNS (Gayen et al., 2016;Parker et al., 2021) are used to supplement the melt parameterization for melt plume-driven boundary layers in the following section.

Application of Theory to Observations at LeConte Glacier, Alaska
In this section, we apply the synthesized melt and plume theories discussed in Section 2 to observations at a rapidly melting vertical ice face.This is motivated by recent estimates of glacial melt rate using repeat multibeam measurements at LeConte Glacier, Alaska (Sutherland et al., 2019), which observed much larger melt rate estimates than those predicted by prior applications of melt plume theory (Jackson et al., 2020).
Figure 4 shows a comparison between the melt rate and plume velocity at Leconte Glacier using the temperature and salinity profiles from the August 2016 field campaign (panel c, Sutherland et al. (2019)).This demonstrates a significant difference in melt rate distribution at LeConte glacier calculated using the traditional free plume melt parameterization (panel (a)) versus the updated wall plume melt parameterization (panel (b)).Note that the free plume melt paramterizations uses the different representation of the turbulent transfer coefficient in Equation 4versus the updated form in Equation 6.For the updated wall-bounded melt plume parameterization (panel (b)), we also use a smaller entrainment coefficient α = 0.068 (based on Parker et al. (2021)) and much larger C p d = 0.15 consistent with wall-bounded melt plumes, commonly-used turbulent heat and salinity transfer constants of Γ T = 0.022 and Γ S = 6.2 × 10 4 (based on Jenkins et al. ( 2010)), and the horizontal velocity using a value of 0.2 m/s from Jackson et al. (2020) (which influences the melt rate directly in Equation 6, but also the 1D plume theory in Equation 7).Including these differences results in a maximum melt rate that increases from <0.1 to 2.2 m/day and much larger separation distance between meltwater intrusions (the darker colors in these panels show where plumes intrude/new plumes form; see Jackson et al. (2020) for a discussion of these intrusions).The melt plumes were initiated at 10 m horizontal intervals along the grounding line and stacked vertically (a new plume was initiated following a previous plume reaching zero momentum/forming a meltwater intrusion).The horizontally-adjacent plumes are also assumed to be non-interacting.See Figure S1 in Supporting Information S1 for additional panels that reflect the melt rates for the free plume theory with horizontal circulation and discharge plume scenarios.
The characteristic vertical velocity is shown in panel (d) for 4 different cases: a free plume (with low drag and high entrainment), a free plume with horizontal circulation-driven melt, a wall plume without horizontal circulation, and a wall plume with horizontal circulation.This shows that the including horizontal velocity-driven melt for a free plume and using wall-bounded plume (drag and entrainment) coefficients have similar effects; they increase the vertical velocities from less than 0.01 m/s to 0.04 m/s and also increase the intrusion separation distance by a factor of 5.If a horizontal velocity-driven melt is included for a wall-bounded plume, their combined effects compound for weakly stratified depths (e.g., z = 170 to 90 m), while they do not differ from their component effects for strong stratified depths (e.g., z = 80 to 40 m) since the increased inertia is still not adequate to overcome the background stratification.Note that the wall-bounded plumes and horizontal circulation cases produce reasonable intrusion separation distances comparable to the observations from Jackson et al. (2020) while the free plume coefficients does not.
In panel (e), the meridionally-averaged melt rates from panel (b) are compared with a discharge plume added (in black; note the entrainment coefficient α = 0.079 based on our LES results, and drag coefficient C p d = 0.015 that is used here that is distinct from the melt plume values) and the melt rate estimates from the repeat multibeam survey.The melt rate estimate for buoyancy-controlled boundary layers is shown for comparison (approx.0.3 m/ day using Equation 1).For melt plumes, modifying drag and entrainment leads to an overall effect of amplifying melt rates at all depths by a factor of 40, which can be attributed to 8x due to increased C p d compounded with 5x from the horizontal circulation-driven melt.The horizontal circulation component contributes directly to the melt and buoyancy input (via Equation 6), which feeds back on the melt plume's vertical velocity (via Equation 7).In the discharge plume case, the buoyancy flux increase is relatively minor (<2% since the buoyancy flux from melting is very small compared to the buoyancy flux from the discharge plume), so there is no feedback between the additional melting and vertical momentum of the plume.These higher melt rates leads to proportional higher buoyancy fluxes and overturning circulation within the fjord.
In summary, the total and local melt rates, and meltwater intrusion depths are consistent between the observations and the updated wall plume theory when horizontal circulation is included, while the melt is significantly underpredicted by using free plume theory alone.The melt rates are somewhat underpredicted by using wallplume theory without horizontal circulation or free plume theory with horizontal circulation.

Summary and Conclusions
In this study, we show that wall-bounded plumes are very different dynamically from free plumes.We propose an updated shear-driven parameterization that uses physically-reasonable values for drag coefficients, melt rates, and entrainment for wall bounded plumes and vertical ice-ocean interfaces.We then test the impact and validate (using Large Eddy Simulations and observations) these updated parameterizations.
These differences are summarized as follows: (a) The plume-driven drag coefficient C p d ) is distinct from the externally forced drag coefficient C ext d ) .Unlike an unstratified flow over a flat plate, C p d is not a drag coefficient in the classical sense as it does not depend on the roughness of the surface: in melt and plume theory, it is used as a means of quantifying the buoyancy-driven turbulence and momentum budget.As such, it is necessary to use a drag coefficient that is relevant to the dynamics in question.Based on recent numerical and laboratory experiments, estimates of the plume-driven drag coefficients have been proposed for discharge plumes C p d = 0.015) and melt plumes C p d = 0.15) .These differences reflect the different types of boundary layers (i.e., v(x), w(x), and W(x) in Figure 1).(b) When wall plumes are parameterized, the entrainment coefficient α should use a much smaller value: α = 0.079 for discharge plumes (based on our LES results) and α = 0.068 for melt plumes (Parker et al., 2021).(c) Horizontal boundary layers v(x) and their melt contribution should still use shear boundary layer parameterizations with drag coefficients of C ext d = 2.5 × 10 3 .However, it is important to include the effect of this melt within the ambient melt plumes as their dynamics are sensitive to horizontal melt rates.
Currently, buoyancy fluxes and glacial melt rates at vertical ice-ocean interfaces are commonly parameterized using theories for unbounded free plumes and assume a universal drag coefficient.However, both Direct Numerical Simulations (Gayen et al., 2016) and laboratory experiments (Parker et al., 2020(Parker et al., , 2021) ) suggest that wallbounded plumes leads to different plume entrainment and vertical velocity profiles (with differences between subglacial discharge and melt plumes) due to the presence of a shear boundary layer.In addition, a recently datasupported parameterization of the turbulent transfer function that merges the velocity-dependent and independent (buoyancy-dominated) melt regimes (Schulz et al., 2022) found a significantly higher baseline buoyancydominated melt rate than previous literature (e.g., Kerr and McConnochie (2015)).Our study reconciles these inconsistencies using a physically-motivated melt parameterization that makes minimal assumptions, includes both buoyancy-driven and shear-dominated melt regimes, and is broadly consistent with existing observations, laboratory experiments, and field data.
We compare the predictions of free plume and wall-bounded plume theories to a discharge plume-resolving LES (MITgcm).We show that these LES results are consistent with previous theories for the along-plume and acrossplume profiles of vertical momentum.Finally, we demonstrate that using the wall-bounded plume modifications leads to a 40x factor increase in melt rate prediction for LeConte Glacier, which is necessary for consistency with existing observations.There are a number of caveats related to the melt rate parameterization we use in this study.We have chosen to examine the dynamics in terms of the limiting cases of discharge and melt plumes, and horizontal flows; our choices and suggestions of various relationships are a preliminary step to merging numerical, laboratory, and observational progress of near-vertical ice-ocean boundary layer dynamics.These include the assumption that the turbulent transfer coefficients (Γ T and Γ S ) are constant between discharge plumes, melt plumes, and horizontal flows (and that the differences in melt rate can be explained by different drag coefficients).This is likely more complicated because the momentum dissipation (constrained by the drag coefficient) may not always be proportional to heat and freshwater flux; this may occur at near-horizontal interfacial slopes where stratification near the viscous layer plays a more important role in the momentum budget (Vreugdenhil & Taylor, 2019).Further, drag coefficients are likely to transition smoothly from discharge plumes (0.015) to melt plumes (0.15) for very weak discharge or discharge plumes with longer lifetimes, that is, beneath ice shelves.The lower viscous shear for discharge plumes is primarily due to the initial buoyancy flux, which is distributed farther away from the interface than melt plumes; over time, the buoyancy will likely evolve to be closer to the ice boundary.The temporal/initial mass flux transition between these types of plumes is likely constrained by the critical Ra number (Equations 1 and 2) and requires further study.
Future work may test these parameterizations for consistency with other direct observations and numerical simulations near vertical ice faces and icebergs.In addition, varying discharge rates, external temperature and velocity, and interface slopes would be necessary to validate this parameterization, which is only constrained here by one warm glacier.Additional modeling studies at both the LES and DNS resolution are needed to understand melt plumes, especially for transitions from buoyancy-dominated to shear-dominated boundary layers and in the presence of both plumes and external velocity forcing.In addition, we may extend these ideas to sloping and geometrically-complex ice-ocean interfaces including ice-shelf cavity geometries, which may also include transition regions from near-vertical interfaces to near-horizontal interfaces.Finally, direct observations of the entrainment rate, melt rate, and the boundary layer profiles of both discharge and melt plumes are necessary to improve our understanding of ice-ocean boundaries.

•
A modified wall-bounded plume parameterization motivated by recent numerical/lab work is proposed as an alternative to free plume theory • Subglacial discharge plume simulations at a vertical ice face are consistent with entrainment/plume dynamics from wall-bounded plume theory • Melt rate estimates using updated parameters are consistent with observations at LeConte Glacier (40 times greater than standard estimates) Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.(a) An overview schematic of wall-bounded discharge and melt plume and the ice-ocean boundary layer.The horizontal and vertical velocity profiles for the discharge plume and melt plume (v(x), W(x), and w(x), respectively) are illustrative and not to scale.(b) A table of reference ranges of drag coefficient C d and plume entrainment coefficient α corresponding to the three types of ice-ocean boundary layers compared to the commonly-used free plume parameterization and their references.

Figure 2 .
Figure 2. (a) Vertical velocity w (m/s) at x = 3 m away from the ice.(b) Open-ocean boundary condition profiles of conservative temperature θ a and salinity S a .(c) The vertical velocity surface w = 0.2 m/s.(d) Meridionally-averaged vertical velocity contours from 0.0 to 2.3 m/s (orange is 0.1 m/s spacing, black is 1.0 m/s spacing) plotted on density (in color).

Figure 3 .
Figure 3. (a) The horizontal variation of vertical velocity w(x) for three test cases: free plume (black), wall-bounded plume without wall drag (red), and wall-bounded plume with wall drag (blue) with theoretical predictions based on Equations 8 and 9 for the free plume, and wall-bounded plume with drag case.(b) The vertical variation of characteristic vertical velocity W(z) for each of the three cases in (a) with plume-theory solutions using Equations 7a-7d.

Figure 4 .
Figure 4. Melt rates at the LeConte glacier face calculated using (a) free plume parameters (Jackson et al., 2020), and (b) wall plume parameters with an additional horizontal circulation melt contribution driven by a uniform horizontal velocity of v = 0.2 m/s.Note that the color ranges between panels (a) and (b) differ by a factor of 40.(c) Temperature and salinity profiles from Sutherland et al. (2019) used in the calculations.(d) The plume velocity as a function of depth (assuming a starting depth of 200 m, as in the location dotted red line in panel b) for free plume parameters (green), a free plume with horizontal circulation (blue), wall plume parameters (red), and wall plume parameters with horizontal velocity (black).The observed mid-depth intrusion separation of approx.20 m in Jackson et al. (2020) is shown for comparison.(e) The meridionally-averaged melt rate for various theories and approximations is shown and compared to the repeat multibeam survey-based estimates from Sutherland et al. (2019).In addition to the cases considered in (d), an additional line discharge plume (with total discharge rate of 220 m 3 /s imposed between y = 250 and 350 m) and a buoyancy-controlled boundary layer melt estimate (Kerr & McConnochie, 2015) are shown for comparison.