The Vertical Structure and Entrainment of Subglacial Melt Water Plumes

Basal melting of marine‐terminating glaciers, through its impact on the forces that control the flow of the glaciers, is one of the major factors determining sea level rise in a world of global warming. Detailed quantitative understanding of dynamic and thermodynamic processes in melt‐water plumes underneath the ice‐ocean interface is essential for calculating the subglacial melt rate. The aim of this study is therefore to develop a numerical model of high spatial and process resolution to consistently reproduce the transports of heat and salt from the ambient water across the plume into the glacial ice. Based on boundary layer relations for momentum and tracers, stationary analytical solutions for the vertical structure of subglacial non‐rotational plumes are derived, including entrainment at the plume base. These solutions are used to develop and test convergent numerical formulations for the momentum and tracer fluxes across the ice‐ocean interface. After implementation of these formulations into a water‐column model coupled to a second‐moment turbulence closure model, simulations of a transient rotational subglacial plume are performed. The simulated entrainment rate of ambient water entering the plume at its base is compared to existing entrainment parameterizations based on bulk properties of the plume. A sensitivity study with variations of interfacial slope, interfacial roughness and ambient water temperature reveals substantial performance differences between these bulk formulations. An existing entrainment parameterization based on the Froude number and the Ekman number proves to have the highest predictive skill. Recalibration to subglacial plumes using a variable drag coefficient further improves its performance.

2 of 34 melting and glacier flow is clear, because melt-driven thinning of the floating ice tongue reduces buttressing of the flow across the grounding line (Goldberg et al., 2009).
For the largest floating ice tongue of the 79°N glacier, only 11% of the freshwater enters the fjord directly as subglacial discharge, and about 89% stems from subglacial melting at the ice-ocean interface (Schaffer et al., 2020). Hence, a cold and relatively fresh buoyant water mass composed of contributions from subglacial discharge, melt water and entrained ambient ocean water is produced that propagates upwards toward the calving front as a turbulent plume (Hewitt, 2020). The correct quantification and prediction of the subglacial melt rate under highly variable environmental conditions has been an aim of polar oceanography for decades. Melt processes in the melt layer are typically parameterized based on a three-equation model for melt-layer temperature, melt-layer salinity and melt rate derived from equilibrium fluxes of freshwater and heat across the melt layer (Hellmer & Olbers, 1989). The challenge is to relate these processes to the properties of the underlying subglacial plume. To this end, similarity relations are typically applied, resulting in logarithmic profiles for momentum and tracers in the near-interfacial region of the melt water plume (Kader, 1981;Kader & Yaglom, 1972;Yaglom & Kader, 1974).
After applying the similarity relations for a vertically integrated model of subglacial plumes (Jenkins, 1991(Jenkins, , 1992, plume models have become powerful tools for understanding melt processes underneath floating ice tongues and ice shelves (Holland & Feltham, 2006;Jenkins, 2011;Jenkins et al., 2010;Payne et al., 2007). Their strength is the computational efficiency allowing high horizontal resolution and the reproduction of the plume thickness, but they rely on accurate parameterizations of entrainment of ambient ocean water into the plume. Many different entrainment parameterizations exist Jungclaus & Backhaus, 1994;Wells et al., 2010). These are generally derived for dense bottom currents and typically depend on non-dimensional bulk parameters such as the bulk Richardson number, the Froude number or the Ekman number. A specifically simple and robust parameterization is based on a constant entrainment rate (ratio of entrainment velocity to plume current speed, Jenkins, 1991). The diversity of entrainment parameterizations shows that there is quite an uncertainty in determining the plume dynamics. In spite of their success, plume models have their specific limitations: They do not predict the ocean temperature and salinity underneath the plume such that the amount of entrained heat and salt is highly uncertain.
To overcome the limitations of plume models and to predict better the effects of larger scale ocean processes, three-dimensional ocean models with explicit ice shelf-ocean interfaces were developed.
For ocean models with geopotential coordinates, a result of the typical step-like approximation of slopes is that the sloping ice-ocean interface is poorly resolved (Losch, 2008). To avoid this issue, terrain-following coordinates are often used where the top layer follows the ice-ocean interface and the lower most layer follows the bottom topography, with non-linear zooming of layers toward surface and bottom (Dinniman et al., 2007). Due to the pressure gradient errors in models with terrain-following coordinates (Haney, 1991), large scale ocean models including ice shelves sometimes apply hybrid coordinates with terrain-following properties inside ice shelves and geopotential coordinates elsewhere (Timmermann et al., 2012). Terrain-following coordinates have the clear advantage of smoothly resolving the ice-ocean interface at high vertical resolution. However, their disadvantage is that the vertical resolution near the ice-ocean interface depends directly on the water depth. Typical top-layer resolutions of terrain-following ice-shelve models around Antarctica vary between 0.5 m near the grounding line and 5 m near the calving front (Gwyther et al., 2020). In models with higher vertical resolution near the ice-ocean interface, the insulating effect of subglacial plumes could be better reproduced than coarse-resolution models (Gwyther et al., 2020). As a consequence, coarse resolution models tend to overestimate melt rates at the ice-ocean interface. Vertically adaptive coordinates with specifically high resolution in the entrainment layer (Gräwe et al., 2015;Hofmeister et al., 2010), which would have the potential to resolve subglacial plumes independently of the water depth, have not yet been used in models with an ice shelf-ocean interface.
We can expect that improved strategies for vertical coordinates and available computer resources will allow very high vertical resolution of subglacial plumes and gravity currents, so that related processes can be simulated more accurately with the prospect of higher predictability of the melt rate. A similar emphasis should be placed on realistic turbulence closure schemes in circulation models underneath ice shelves, because the basal melt rates strongly depend on the parameterization of mixing processes and entrainment (Dansereau et al., 2014). Exploring Water column models, also with second-moment turbulence closures, have been used to study melting (and freezing) under sea ice (Mellor et al., 1986;Omstedt & Svensson, 1984;Steele et al., 1989). Analogous studies of the vertical structure of subglacial plumes are in their infancy, but include models with simple (Jenkins, 2016(Jenkins, , 2021 and two-equation turbulence closures (Cheng et al., 2020). In these models, a well-mixed turbulent boundary layer is separated from the ambient water underneath by a stratified layer at marginal stability, across which quiescent ambient water is entrained. The resulting profiles of velocity and eddy diffusivity are very sensitive to parameters such as the roughness of the ice-ocean interface or the transfer velocities for salt and heat. However, the simplicity of the applied turbulence closures in the former case and the particularity of the application in the latter case render general applicability of the results very uncertain specifically in the region of the entrainment layer. While Direct Numerical Simulation (Rosevear et al., 2021) and Large Eddy Simulation (Vreugdenhil & Taylor, 2019) have been applied to the ice-shelf-ocean boundary, the entrainment layer has not yet been studied in these applications because of the limited spatial scales considered. In the present study we overcome these limitations and develop, present, and apply a more general high-resolution water-column model for subglacial plumes that includes realistic second-moment turbulence closures.
Melt processes under floating ice tongues are very difficult to observe in their harsh and barely accessible polar environments. Therefore, the dynamic analogy between buoyant plumes under shelf ice and dense bottom currents due to overflows across sills have been applied to validate plume models (Jenkins, 2016). The main difference is that in ice shelves the buoyancy is mostly produced locally due to subglacial melt, but the (negative) buoyancy in dense bottom currents is a result of upstream processes. While this certainly has substantial effects on larger time and space scales, the vertical structure of both regimes may be comparable. Exploiting this analogy, most formulations for entrainment in plume models are derived from studies of dense bottom currents. In the present study, we apply previous modeling concepts of simulating rotational dense bottom currents in the Western Baltic Sea Umlauf et al., 2010). The subglacial plume model developed here serves the following purposes: 1. Develop a consistent dynamic coupling between parameterized melt layer processes and turbulent processes within the plume and the entrainment layer 2. Develop consistent and convergent discretization methods for melt fluxes that give robust results also for relatively coarse resolution, and 3. Test existing formulations and calibrate a new parameterization of entrainment that can be applied in vertically integrated plume models This paper is structured as follows: First, the underlying mathematical formulations are given, with the water-column equations (Section 2.1), the boundary conditions (Section 2.2), the melt formulations (Section 2.3), the tracer roughness lengths (Section 2.4), and a stationary analytical model of the vertical plume structure (Section 2.5). Afterward, numerical issues are discussed, with discretization methods for velocity (Section 3.1) and tracers (Section 3.2), the numerical treatment of the free surface (Section 3.3), and with numerical convergence experiments (Section 3.4). The transient model simulations with the turbulence closure model are described, with the model setup (Section 4.1), the models results including the default scenario and sensitivity studies are presented (Section 4.2), and a comparison of the model results to the performance of entrainment parameterizations is made, including calibration of a new formulation (Section 4.3). Finally, the main results of the study are discussed (Section 5) and some conclusions are drawn (Section 6). In the appendix, details of the analytical solution (Appendix A) and an analytical dependence of plume speed and friction velocity on the interfacial roughness length (Appendix B) are given.

Water-Column Model Equations
The hydrodynamic and hydrographic water column equations for a buoyant melt water plume under a planar ice-ocean interface with slope ∂ x z b = tan α x , ∂ y z b = tan α y (with the vertical position of the ice-ocean interface z = z b , where z is the upward directed vertical coordinate with the origin at the undisturbed mean sea level) are BURCHARD ET AL.

10.1029/2021MS002925
4 of 34 based on the Reynolds-averaged Navier-Stokes equations with the Boussinesq assumption and the down-gradient parameterization of vertical turbulent fluxes (Umlauf & Burchard, 2005). We assume stagnant and homogeneous ambient water with velocities u = v = 0, potential temperature θ = θ 0 , salinity S = S 0 and potential density ρ = ρ 0 below the plume. The z-axis is assumed to be pointing upwards exactly opposite to the gravitational acceleration. The plume properties are assumed to be homogeneous along the ice-ocean interface, that is, all gradients along the slope vanish: furthermore, the flow is assumed to be parallel to the planar ice-ocean interface, resulting in the kinematic condition = tan + tan .
the geometry of the one-dimensional water column model is sketched in Figure 1. In a one-dimensional hydrostatic water column model the pressure-gradient driven acceleration in x-direction is calculated as with the surface pressure p(z b ) (from atmospheric pressure plus the additional pressure due to glacial ice). Using Equation 1 we obtain for the ambient stagnant water below the plume with z → −∞ and ρ(−∞) = ρ 0 , we demand that the pressure gradient vanishes, that is, such that we obtain with the buoyancy which is positive inside the subglacial plume and vanishes in the ambient water. The pressure gradient in y-direction is calculated accordingly. Due to the geometry of the flow, the total derivative of any property c equals the partial time derivative: where Equations 1 and 2 have been used, such that the dynamic equations for the velocity components u and v read with the eddy viscosity t and the Coriolis frequency f. The second terms on the left hand side represent the stress divergence with the stress vector similar dynamic equations have been used for simulations of subglacial plumes (Jenkins, 2016(Jenkins, , 2021 as well as for dense bottom currents, where less dense ambient water resides above the plume . In both modeling concepts, the coordinate system is defined such that the z-axis is orthogonal to the slope of the model instead of being aligned with the gravitational forcing. However, for mild slopes, the differences to our approach outlined above are negligibly small. In the present study, the formulation of a vertical z-axis is used in order to be consistent with hydrostatic three-dimensional ocean models, which loose their validity for steep slopes where the vertical acceleration becomes relevant.
The budget equations for potential temperature θ and salinity S are formulated as with the eddy diffusivity ′ . The hydrographic Equation 11 are linked to the hydrodynamic Equation 9 by means of an equation of state for potential density, calculated according to Jackett et al. (2006), with the atmospheric pressure at the sea surface, p 0 . Consequently, ρ 0 = ρ(θ 0 , S 0 , p 0 ). Water column stability at a depth with pressure p z = const, that is, ∂ z p z = 0, is then calculated as with the Brunt-Väisälä frequency N.
Eddy viscosity t and eddy diffusivity ′ are calculated in two ways here. For the analytical calculations presented in Section 2.5, 2.5.3, 3 and Appendix A parabolic profiles for t and ′ are chosen that extend over the entire thickness of the plume, see Section A1 and A2. Such parabolic profiles are often used for well-mixed open channel flow, see the recent discussion by Absi (2021), and allow for analytical treatment of velocity and tracer profiles (Burchard et al., 2013;Lange & Burchard, 2019).
For more realistic simulations that do also allow for predictions of entrainment rates at the base of the plumes, eddy viscosity and eddy diffusivity are determined by means of a two-equation turbulence closure model with an algebraic second-moment closure (Umlauf & Burchard, 2005). This closure is based on an equilibrium assumption for the second moments (turbulent transports of momentum and tracers), that is, the transport terms for the second moments are neglected and only the source and sink terms are retained. The two equations of the closure model represent budgets of the turbulent kinetic energy k and its dissipation rate ɛ. The eddy coefficients are then calculated as where c μ and ′ are quasi-equilibrium (assuming an equilibrium condition for the budget of k only for the second-moment closure) non-dimensional stability functions representing the second-moment closure. The argument of the stability functions is the non-dimensional buoyancy number the weak-equilibrium stability functions that additionally depend on a non-dimensional shear number (Umlauf & Burchard, 2005) are not used since they were found to induce some small-scale oscillations in the entrainment layer at the base of the subglacial plume.
The buoyancy term in the ɛ equation is calibrated in a way that for homogeneous shear layers in equilibrium the gradient Richardson number converges toward the steady-state value of one quarter (Burchard & Baumert, 1995).

of 34
This guarantees the correct representation of entrainment rates at the base of surface mixed layers (Umlauf & Burchard, 2005) or on top of dense bottom currents (Umlauf et al., 2010).
In contrast to the second-moment closure used here, Cheng et al. (2020) applied the approach of a standard-k-ɛ model with constant stability functions for their simulations of super-cooled subglacial plumes. In their model, the buoyancy term is not specifically calibrated for reproduction of realistic entrainment rates.

Boundary Conditions
At the upper boundary at z = z b (ice-ocean interface) a no-slip boundary condition for velocity is fulfilled: In the framework of this water-column model, the upper boundary is treated as a rigid lid, that is, melt and freezing processes do not lead to a change in water depth, other than in free-surface models. The dilution of the surface water due to addition of melt water is parameterized here as a virtual salinity flux, see Jenkins et al. (2001) for a discussion of boundary conditions for material (rigid-lid) and immaterial (free-surface) boundary treatment at the ice-ocean interface. Free-surface boundary conditions for freshwater and heat, where the melt water is added to the water column, are given in Section 2.3. The diffusive ocean-to-ice fluxes (orthogonal to the ice-ocean interface) of potential temperature and salinity, and , are located at the same position as the no-slip condition for momentum: note that = 0 is the ocean-to-ice heat flux at the ice-ocean interface, with the heat capacity of ocean water, c. For simplicity, we apply the ocean-to-ice fluxes in the vertical direction, without prior projection from the orthogonal direction. This approximation is valid for small slopes. For example, for a slope of tan α x = 5 ⋅ 10 −3 (i.e., a slope angle of 0.28°), the error is about 5 ⋅ 10 −3 .
Near the boundary, the spatial variation of all momentum and tracer fluxes can be neglected, such that their exact vertical location within the melt layer is not relevant. This plays a role when constructing logarithmic near-boundary profiles based on these fluxes and Dirichlet boundary conditions that are located at slightly different vertical locations (see Section 2.4).
The boundary conditions for the turbulent quantities at the ice-ocean interface are best explained by means of near-boundary profiles as functions of the distance from the interface, z′ = z b − z: where 0 is the equilibrium stability function for unstratified conditions, κ is the van Karman constant, and z 0 is the hydrodynamic roughness of the ice-ocean interface (see Section 2.4 for details). From the turbulence boundary profiles Equation 19 two sets of boundary conditions for z′ = 0 have been derived, Dirichlet conditions and Neumann conditions, of which Burchard and Petersen (1999) could show that the latter are much more accurate. Note that with Equation 14, the near-interface profile of the eddy viscosity is linear:

Melt Rate
To derive formulations for the melt rate and the heat fluxes at the ice-ocean interface, a very thin melt layer at freezing temperature is assumed. The fluxes of potential temperature and salinity across the ice-ocean interface strongly depend on the respective molecular diffusivities, T = /Pr T and S = /Pr S , where is the molecular viscosity, Pr T = 13.8 is the Prandtl number for temperature and Pr S = 2,432 is the Schmidt number for salinity.
For the derivation of the melt rate, v b , that is, the rate at which water is added to the ocean by means of subglacial melting, we largely follow the paper by Holland and Jenkins (1999) who compare various formulations. We adopt the well-known three-equation model that is based on flux equilibria of heat and salt across the melt layer and a linear equation for the freezing temperature. With this, the upward heat flux at the ice-ocean interface is composed of the diffusive heat flux into the ice and the latent heat flux needed to melt the ice: where L i is the latent heat of fusion. Note that ρ i v i is the mass of ice per unit time and unit area that is melted, such that v i is the velocity at which the ice-ocean interface is retreating. The mass of the ocean water that is gained due to melting must be equal to the mass of ice that is melted such that where v b is the melt rate, that is, the increase of sea surface height due to melting per unit time. For the flux into the ice, various formulations are available. We adopt the approach based on an advection-diffusion equation of temperature in the glacial ice, with the vertical advection velocity v i . Based on that, the heat flux into the ice due to diffusion can be formulated as (Holland & Jenkins, 1999 with the heat capacity of ice, c i , and the ice-core temperature, T I . Note that we use here the potential freezing point temperature θ b instead of the in-situ freezing point temperature T b , to allow for an easy comparison with the ocean potential temperature. Combining Equations 22-24, we obtain for the upward flux of temperature at the ice-ocean interface with the heat capacity of sea water, c. Using the potential melt layer temperature instead of the in-situ temperature and comparing it to the ice-core temperature here does not pose a problem, due to the large difference between melt layer and ice core temperatures and the typically large uncertainty in the latter.
Since the total salt flux into the glacial ice must be zero, the diffusive salt flux into the melt layer must be assumed to be opposite to the advective salt flux: we use a linear equation for the freezing temperature, assuming that the melt layer temperature is at the freezing point: with the empirical parameters λ 1 , λ 2 and λ 3 . Note that we use slightly different empirical values than Holland and Jenkins (1999), to apply the potential temperature of the freezing point instead of its in-situ temperature. The BURCHARD ET AL. Information about the plume properties in terms of velocity, temperature and salinity are required to close the meltrate computations. This will be provided by either an analytical solution for the vertical structure of the plume (Section 2.5) or from a numerical model which uses the analytical model to consistently provide the plume information (Section 3.3).
Following Jenkins et al. (2001), the free-surface tracer boundary fluxes can be formulated as

Roughness Lengths for Potential Temperature and Salinity
Similarly to the classical logarithmic law of the wall for velocity profiles, logarithmic profiles are constructed for temperature and salinity in order to derive numerically consistent boundary conditions. These profiles are highly simplified and do not resolve but parameterize the effects of wall roughness and the viscous sublayer by means of surface roughness lengths. While the boundary condition for the velocity profiles in the melt layer (see Section 2.3) is a no-slip condition (u = v = 0), boundary values for temperature and salinity in the melt layer are given by θ = θ b and S = S b . It should however be noted that the locations for the boundary values for velocity, temperature and salinity are slightly different. This is due to the substantially different values for kinematic viscosity and the laminar diffusivities for θ and S. These boundary values are formally located at positions slightly above the interface z = z b : where 0 ≪ 0 and 0 ≪ 0 are formally defined as tiny roughness lengths specific for temperature and salinity fluxes. The formulations for these roughness parameters given below have been taken from Kader and Yaglom (1972), Yaglom and Kader (1974) and Kader (1981). For a hydrodynamically rough interface, the roughness length for momentum is given as where k s is the characteristic height of the roughness elements and B′ = 8.5 is an empirical parameter, such that z 0 ≈ k s /30. For a hydrodynamically smooth interface where B = 5.5 is an empirical parameter, such that 0 ≈ 0.11 ∕ * .
For both, hydrodynamically rough and smooth interfaces, the roughness scale with respect to the flux of c (where c represents any tracer such as T or S) is where the value of β c is calculated differently for rough and smooth interfaces and Pr = ∕ ′ is the turbulent Prandtl number. For a hydrodynamically rough interface, with the non-dimensional roughness scale + 0 = 0 * ∕ . For a hydrodynamically smooth interface Equation 32 holds with The dependence of the tracer roughness length on the Prandtl number is shown in Figure 2. Note that it is only the logarithms of the roughness lengths that are evaluated and not their direct values (which are partially too tiny to be computed). Furthermore Equations 32 and 34 are valid for + 0 < 0.1 and Equations 31 and 33 are valid for + 0 > 3.33 . Here we concentrate on the rough wall conditions and therefore set 0 = rough 0 and = rough .

Analytical Plume Model
The analytical plume model that is derived here serves two purposes: It is used to construct consistent discrete formulations for the boundary conditions for θ and S (see Section 3) and it is used to perform an analytical parameter space study for subglacial plumes. To allow for an analytical solution, Earth rotation is neglected (f = 0 in Equation 9), such that only one velocity component needs to be taken into consideration. Furthermore, the interfacial slope is assumed to be positive, such that * = * > 0 . At the plume base, turbulent entrainment of ambient water is assumed by prescribed values of friction velocity * , turbulent temperature flux and turbulent salinity flux . A further simplification is that for the pressure gradient force, the buoyancy of the plume is assumed to be constant: b = const. The profiles are formulated as a function of prescribed depth-mean values of velocity, potential temperature and salinity, , ̄ and ̄ . The profiles are calculated over the entire plume thickness D, assuming parabolic eddy viscosity and eddy diffusivity. For simplicity the distance from the ice-ocean interface z′ = z b − z is used as vertical reference.

Velocity Profile
The derivation of the analytical stationary velocity profile under a sloping ice-ocean interface is shown in Appendix A1: where * will be calculated by means of solving the quadratic Equation A9, and A is a non-dimensional integration constant defined in Equation A10. The first term in Equation 35 is the classical logarithmic law of the wall written in a form where its vertical average is . The second term represents the effect of the entrainment of ambient water at the plume base. It has a vertical average of zero and diverges to ±∞ for z′ → D, depending on the sign of * . A similar solution had been proposed by Lange and Burchard (2019) for the effect of surface wind stress in estuarine exchange flow.

Tracer Profiles
With Equation A20, neglecting a tiny exponential expression in the term representing the effect of entrainment, we can formulate the analytical profiles for potential temperature (c = θ) and salinity (c = S) as follows: and with the potential temperature θ b and the salinity S b of the melt layer. As for the velocity profile Equation 35, also the profiles of potential temperature and salinity diverge toward ±∞ for z′ → D, but also here the vertical averages are finite.
It should be noted that the boundary values for potential temperature and salinity do not converge to θ b and S b for z′ → 0. This is also the case for the classical logarithmic laws with zero entrainment fluxes and . This inconsistency results from the strong gradients of θ and S in the melt layer due to their small Schmidt numbers. However, since the fluxes of θ and S across the melt layer are applied as boundary conditions and since they are assumed to be constant within the boundary layer, this inconsistency is acceptable. This is clearly seen in the analytical profiles shown in Figures 3b and 3c with the integration constants A T and A S calculated according to Equation A22. By combining Equations 36 and 37 with 38 and 39, a formulation of the profiles is obtained that depends on prescribed values of the depth-averaged potential temperature ̄ and salinity ̄ . The melt layer freezing temperature and melt layer salinity θ b and S b can now be determined, using the melt layer formulation given in Section 2.3. Combining Equations 25 and 26 for the interface fluxes of potential temperature and salinity, the linear freezing temperature formulation Equation 27 with the vertically averaged Equations 38 and 39 for potential temperature and salinity of the plume gives five equations for the five unknowns , , v b , T b and S b . These equations are combined in a way that a quadratic equation for S b results, of which the positive solution is the physically correct one.

Analytical Examples
Two sets of plume profiles are calculated, without entrainment at the plume base (experiment N) and with entrainment at the plume base (experiment E). Results for u, θ and S are shown in Figure 3, using the empirical parameters given in Table 1. For both experiments, bulk values of the plume thickness D = 20 m, the depth-mean velocity = 0.2 m s −1 , the depth-mean temperature ̄= −1.75 • C, the depth-mean salinity ̄= 33.1 g kg −1 , the interfacial depth z b = −300 m, and an interfacial roughness length of z 0 = 10 −2 m are prescribed. These values are similar to those at the end of the transient default experiment that will be presented in Section 4.2.1. Results are shown in linear and logarithmic scale.
For the experiment without entrainment, the profiles of velocity, temperature and salinity are exactly logarithmic (Figures 3d-3f). Due to the small Schmidt numbers, the slopes of the temperature and salinity profiles are very small. In the case of salinity, the difference across the full depth of the plume is about only 6 ⋅ 10 −3 g kg −1 , such that in the non-logarithmic presentation (Figure 3c), vertical salinity gradients can not be detected by visual inspection. A striking feature of the analytical temperature and salinity profiles is the substantial difference between the melt layer values, θ b and S b , and the boundary values, θ(z b ) and S (z b ), see Figures 3b and 3c. This is a direct consequence of Equation 29, and does not pose a practical problem, because it is not the boundary values that are applied, but the fluxes of temperature and salinity which are assumed to be constant on the scale of the melt layer thickness. There is a conceptual issue, because the construction of the free-surface boundary conditions assumes that melt layer values and boundary values are identical (see Jenkins et al., 2001, and Equation 28 of this study). Since these free-surface boundary conditions are not used here, it is beyond the scope of the present study to resolve this inconsistency.
The analytical solution including fluxes of momentum, temperature and salinity across the base of the plume (exp. E, Figures 3g-3l) allows to mimic entrainment of ambient water. The entrainment fluxes are estimated as follows: with the entrainment velocity = 0.036̄sin , see Jenkins (1991) and Section 4.3, and the ambient values are chosen to be the same as in the default scenario of the transient simulations presented in Section 4.2.1: tan α = 0.005, u 0 = 0, θ 0 = 1°C and S 0 = 34.5 g kg −1 . With the above parameters, we obtain v e = 3.6 ⋅ 10 −5 m s −1 , * = −7.6 ⋅ 10 −6 m s −1 , = 9.9 ⋅ 10 −5 K m s −1 and = 5.0 ⋅ 10 −5 g kg −1 m s −1 . As a result, the entrainment has only a minimal influence on the velocity and salinity profiles (Figures 3g and 3j and 3i,3l). However, the temperature profile (Figures 3h-3k) shows slightly larger vertical gradients due to the entrainment of relatively warm water. Due to the stationary character of the analytical solution with fixed values of , ̄ and ̄ , there is no predictive skill. The melt rate (v m = 5.22 m y −1 ) and the ocean-to-ice heat flux ( = 63 W m −2 ) depend only weakly on the entrainment. But due to entrainment of warm and salty ambient water, the instantaneous trends in average plume temperature and salinity are changed from ̄= −6.3 ⋅ 10 −2 K day −1 to ̄= +0.37 K day −1 and from ̄= −2.2 ⋅ 10 −2 g kg −1 day −1 to ̄= +0.19 g kg −1 day −1 .

Numerical Implementation
In numerical models, the analytical logarithmic profiles for u, θ and S derived in Section 2.5 can be used to calculate fluxes of momentum, temperature and salinity at the ice-ocean interface in a consistent and convergent way. To do so, we assume that within the first grid layer underneath the ice-ocean interface the logarithmic laws Equations 35, 36 and 37 are valid, neglecting the fluxes across the plume base , such that the formulations reduce to and where we have extended the analytical solution to the full velocity vector (u, v). Note, that the analytical profiles Equations 41-43 are expected to hold in some vicinity of the interface only, depending on the forcing of the plume.

Momentum Fluxes
Let in a numerical model ( with the numerical drag coefficient This is the numerically consistent discretization for the momentum flux at the ice-ocean interface as defined in Equation 17 layer is applied in many ocean models for the bottom boundary layer (Klingbeil et al., 2018, their Section 7.6). If, however, a constant value for c d is chosen that does not depend on resolution, then a refined resolution near the boundary will not result in a reduction of velocity and consequently in an underestimation of shear and friction velocity (see Section 3.4).

Tracer Fluxes
The calculations of the temperature and salinity fluxes are carried out in a similar way than the momentum fluxes. When for the discrete profiles of temperature and salinity, ( = max with ′ max = 1 2 ℎ max and the upper layer thickness max are known, the following relations can be derived from the logarithmic tracer laws Equations 42 and 43: and max = + Pr * ( ln note that Equations 47 and 48 can be reformulated as with the exchange velocities with the coefficients β T and β S from Equations 33 or 34. With Equations 25, 26, 27, 49 and 50 we have now five equations to solve for the five unknowns T b , S b , v b , and . Note that we solve this system of equations in order to calculate the melt rate v b which is applied to add fresh and cold water to the ocean surface due to subglacial melting, see Section 3.3 for the implementation.

Free-Surface Versus Rigid-Lid Models
For a free-surface model, max and the tracer concentration max should be discretized as follows (assuming zero volume and tracer fluxes across the interface at the bottom of the layer): and BURCHARD ET AL.

10.1029/2021MS002925
15 of 34 which can be combined into where the superscript indicates the number of the time step. For a rigid-lid model as applied in the present study, the numerical scheme Equation 55 will apply, but with a constant layer thickness max : where the respective free-surface fluxes for temperature and salinity are calculated according to Equation 28.

Numerical Convergence Experiment
To test the sensitivity of the numerical methods developed in Section 3.1 -3.3 in terms of spatial resolution and the numerical method for calculating ocean-to-ice fluxes, the momentum Equation 9 for non-rotational flow (f = 0 and v = 0) and the temperature and salinity Equation 11 were discretized. The entire water column was accelerated by a barotropic (constant in space) pressure gradient in a way that the prescribed depth-averaged plume velocity of resulted. Additionally to the melt fluxes, temperature and salinity were forced with a depth-independent source/sink term compensating for the freshening and cooling, see Equation A18. Empirical parameters were chosen identically to Experiment N in Section 2.5.3, such that an analytical solution is available for quantifying the accuracy of the numerical scheme for different vertical resolutions and treatments of the melt flux parameterization.
An Euler-forward central-difference discretization of these diffusion equations was applied with a sufficiently small time step. The numerical scheme was executed until a stationary numerical solution was approximated at high accuracy. Three different vertical discretizations were chosen, ranging from very coarse to very fine (see details in the caption of Table 2).
For each vertical resolution, three different numerical treatments of the melt fluxes were chosen: Note. The very high resolution for k max = 200 layers was obtained by logarithmic zooming toward the ice-ocean interface. still valid because of the high vertical homogeneity of the tracer profiles with max ≈̄ and max ≈̄ and γ T and γ S being largely independent of the layer thickness. In contrast, method 3 diverges from the analytical reference solution (see Figures 4d-4f), since the bulk formulation for the momentum flux is decreasingly representative for increasing resolution with c d not increasing toward ∞.
Since the accurate calculation of the melt rate v b is one of the main goals of refined numerical simulations, it is given in Table 2 for all numerical experiments (three resolutions and three melt flux formulations). For the consistent melt flux formulation 1, melt rates converge toward the analytical value of 5.2236 m y −1 . In the range of the resolutions tested here, method 2 shows an almost indistinguishable behavior, although it should formally not converge toward the analytical solution. Also for the coarse resolutions both methods 1 and 2 give relatively accurate values with only about 3-4 ‰ error for the very coarse resolution. In contrast to this, the bulk method 3 for calculating momentum fluxes diverges substantially, and gives an error of about 1.1 m y −1 (22%) for a top-layer resolution of max = 0.4 m, and an error of 1.6 m y −1 (30%) for the very high resolution of max = 0.015 m. Therefore, it is highly recommended to use methods 1 or 2 for the momentum and tracer fluxes, with a preference for the fully consistent method 1.

Transient Numerical Experiments
The General Ocean Turbulence Model (GOTM, Burchard and Bolding (2001); Umlauf and Burchard (2005); Li et al. (2021), see also www.gotm.net), a one-dimensional water-column model coupled to a library of turbulence closure models, was modified to reproduce the vertical structure of subglacial plumes. For the surface, that is, the ice-ocean interface, no-slip conditions for velocity Equation 16 and fluxes of salinity and heat due to melt processes at the ice-ocean interface Equation 18 were added. Since GOTM treats surface freshwater fluxes like a rigidlid model (i.e., considering changes in concentrations instead of changes in volume), the fluxes of salinity and heat have been implemented according to Equation 56. A further change to GOTM needed to reproduce subglacial plumes was the application of pressure gradients due to vertical buoyancy gradients under a sloping ice-ocean interface, as described in Equation 6. As turbulence closure model, the k-ɛ model with the quasi-equilibrium second-moment closure by Cheng et al. (2002) was used.

Model Setup
The simulations analyzed here start from rest (zero velocity), with an initial plume thickness of D = 5 m. The ice-ocean interface is located at a depth of z b = −300 m and the water column is z b + H = 150 m deep, such that the bottom at z = −H = −450 m is sufficiently far below the ice to allow for an undisturbed plume deepening for all sensitivity studies. The geometry of this one-dimensional plume model is sketched in Figure 1 (not to scale). This set of initial conditions allows for a subglacial plume that is quickly adjusted dynamically to the local conditions. The vertical discretization uses k max = 500 layers with zooming toward the ice-ocean interface such that the resolution is gradually increasing from h 1 = 0.5 m at the bottom to max = 0.09 m at the surface. In a sensitivity study about the effects of coarser vertical resolution, a total of only 50 or 25 layers will be used (Section 4.2.3).
The ambient water is at rest and has a high ocean salinity of S 0 = 34.5 g kg −1 and potential temperatures of θ 0 = 1°C ± 1°C (depending on the scenario). The initial plume salinity and potential temperature are S = 32 g kg −1 and θ = −1°C, such that its potential density ρ is lower than the potential density ρ 0 of the plume water. With this, the initial pressure gradient drives a subglacial plume rising upwards along the slope of the ice-ocean interface. The latitude of the water column location is 79°N, such that the Coriolis parameter has a value of f = 1.43 ⋅ 10 −4 s −1 . The ice-ocean interface is sloping toward the north, while the slope toward the east vanishes (α x = 0). During the simulation time of 14 days, the plume velocity is expected to point toward the north-east (u, v > 0) direction as a consequence of the force balance between northward pressure gradient force, Coriolis force and frictional effects. The plume is subject to cooling and freshening due to melt fluxes at the ice-ocean interface and to warming and salinification due to entrainment of warmer and saltier ambient water. This simulation can be thought of as a plume underneath an infinite plain, where all plume properties are homogeneous along the interfacial slope , see Figure 1.
We analyze one default simulation in detail and carry out six sensitivity simulations with variations in northward slope α y , interfacial roughness z 0 and ambient temperature θ 0 . The parameters for the sensitivity study are given in Table 3.

Analysis of Bulk Values
According to Arneborg et al. (2007), the bulk properties of the plume can be diagnosed from individual profiles as follows:̄= with the plume thickness D, the depth-averaged buoyancy ̄ , and the depth-averaged plume velocity vector (̄̄) . The vertically averaged plume speed is defined as = ( 2 + 2 ) 1∕2 . Depth-averaged salinity and potential temperature are defined accordingly: where Fr > 1 marks supercritical flow and Fr < 1 marks subcritical flow, the Ekman number (ratio of frictional to rotational effects) =̄ (60) and the bulk Richardson number (ratio of bulk stratification to bulk shear) which we need to define for the entrainment formulation by Jungclaus and Backhaus (1994). For the Froude number, the overall slope angle is calculated as = arctan ( tan 2 + tan 2 ) 1∕2 .

Default Scenario
In the default scenario, the plume thickness increases from its initial value of D = 5 m to about D = 20 m within 14 days ( Figure 5). In the initial phase, the plume is accelerated northwards along the slope due to the pressure gradient force, reaching up to a depth-averaged value of = 0.25 m s −1 within 1 hour (Figure 6a). In Figure 5. Simulated cross-slope and along-slope velocity, salinity and potential temperature profiles for the default scenario def during the 14-day simulation period.
BURCHARD ET AL.

10.1029/2021MS002925
19 of 34 this initial phase, the flow is supercritical (Fr > 1) for a short time (Figure 6b). Afterward, due to Earth rotation, the plume velocity veers toward the cross-slope direction (Figures 5a and 5b). This effect is strongest at the plume base where frictional and rotational effects combine in a complex way (Figure 7d), see Figure 11 of Umlauf et al. (2010) for details. During the further development of the plume, frictional effects are reduced due to increased plume thickness and decreased velocity (expressed as strongly decreasing Ekman number, see Figure 6b), such that the depth-averaged velocity vector (̄̄) is further veering toward the downslope direction (Figure 6a), with the downslope velocity peak remaining at the plume base. After the initial adjustment phase, the plume is close to a dynamic balance as indicated by the relatively close agreement between plume velocity diagnosed from the simulation result and the analytical equilibrium velocity diagnosed from the dynamic steadystate assumption of the vertically integrated momentum equations (see Appendix B).
After the initialization, the flow becomes subcritical, with the Froude number slowly decreasing to Fr = 0.6 at the end of the simulation (Figure 6b). The square root of the vertically integrated and thus the shallow-water speed (̄) 1∕2 is slowly increasing during simulation due to the buoyancy fluxes at the ice-ocean interface (Figure 6a).
For no such buoyancy fluxes theory predicts a constant shallow water speed . Potential temperature and salinity are quickly approaching relatively constant values, suggesting a balance between melt fluxes and entrainment fluxes (Figures 5c and 5d).
Since the vertical structure of the plume is almost self-similar after the initial adjustment , the vertical profiles of plumes properties shown in Figure 7 at the end of the simulation time are largely representative for the plume in dynamical balance. When rotating the velocity vector profiles into the direction of the depth mean flow vector, it becomes evident that the velocity peak at the plume base is most pronounced in the cross-flow velocity component (Figure 7a). In the frictionally dominated part of the plume, the cross-flow velocity is negative due to Ekman dynamics (Umlauf & Arneborg, 2009;Umlauf et al., 2010). Potential temperature and salinity are well-mixed in the bulk of the plume, with a gradual increase in stratification at the plume base (Figures 7b and 7c). The entire plume is stably stratified (Figure 7d). Most of the stratification seems to be due to entrainment of denser water from below, but the increase of N 2 toward the ice-ocean interface indicates that some stratification is also induced by the stabilizing ocean-to-ice fluxes. Due to the strong shear at the ice-ocean interface (Figure 7e), the gradient Richardson number decreases continuously in upward direction. It has a maximum of about Ri = 0.75 directly above the entrainment layer. It is characteristic of two-equation turbulence closure models that they allow active mixing at such high stability conditions due to vertical turbulent transport of TKE (Umlauf, 2009). In the entrainment layer itself, Ri attains the value of the steady-state gradient Richardson number of Ri st = 0.25, which is a result of the calibration procedure of the two-equation turbulence closure model Figure 6. Time series of (a) depth-averaged velocity , and = √ 2 + 2 diagnosed from General Ocean Turbulence Model, using Equations 57 and 63, and equilibrium velocity eq , eq and eq = √ ( eq ) 2 + ( eq ) 2 predicted by the steady-state theory presented in Appendix B together with the shallow water wave speed √̄ , as well as the phase velocity and (b) Froude number Fr and Ekman number K for the def scenario.
BURCHARD ET AL.

10.1029/2021MS002925
20 of 34 (Burchard & Baumert, 1995;Umlauf & Burchard, 2005). The rotated stress vector (̂,̂) and the eddy viscosity t are compared to the analytical formulations from Equations A3 and A4 in Figures 7g and 7h, using the simulated surface stress and a mixed-layer depth diagnosed from the first zero-crossing of the simulated shear stress. For the shear stress the agreement is very good, but the parabolic analytical eddy viscosity overestimates the simulated profile by about one third, because it does not take into account effects of stratification that are present at the base of the plume. Still near the ice-ocean interface, the agreement between simulated and analytical eddy viscosity is very good, and both profiles converge to κ(z′ + z 0 ) near the interface. For the simulated profile, this is a consequence of the formulation of the Schmidt number in the ɛ-equation of the turbulence closure model (Burchard & Baumert, 1995;Umlauf & Burchard, 2005) for which the resulting ɛ-profile is shown in Figure 7i.

Sensitivity to Forcing Parameters
In the default scenario described in Section 4.2.1, the entrainment velocity v e is about 1-2 orders of magnitude larger than the melt velocity v m , with peak values of v e reaching almost 2 km y −1 , while the melt velocity has maximum values of about 20 m y −1 (magenta lines in Figure 8). With that, the assumption of a rigid lid does not significantly influence the plume thickness. Highest entrainment and melt velocities are reached in the early adjustment phase of the plume due to maximum Froude and Ekman numbers (Figure 6) after which a steady decrease is observed.
Increasing or decreasing the slope of the ice-ocean interface by a factor of five in the scenarios ayp and aym, has significant effects on the development of the plume. The steeper slope more than triples the final plume thickness to more than D = 75 m, with entrainment velocities of up to v e = 20 km y −1 and melt velocities of up to  (Figures 8a-8g). In contrast, the scenario aym with the strongly decreased slope (Table 3), leads to on a very weak increase in plume thickness during the 14-day simulation, with entrainment velocities of about v e = 20 m y −1 and melt velocities decreasing from initially v m = 10 m y −1 to v m = 1 m y −1 at the end of the simulation.
Also increasing (z0p) or decreasing (z0m) the roughness of the ice-ocean interface by an order of magnitude with respect to the default scenario has a measurable effect on the development of the plume (Figures 8b and 8e, 8h). As shown by the equilibrium theory (Appendix B, Figure B1), an increased roughness should lead to a decreased velocity but to an increased friction velocity. Since conditions for the equilibrium theory are met for the later stages of the plume development (Figure 6a), plume velocity and friction velocity do show this behavior here (not shown). With that, higher interfacial roughness leads to more turbulence inside the plume, and thus a higher entrainment velocity and a higher melt rate (Figures 8e-8h).
In contrast, changes in ambient potential temperature (scenarios tap and tam) have no significant impact on plume thickness and entrainment velocity (Figures 8c-8f). As expected, the melt rate increases by about 1.7 m y −1 for an increase of 1°C in ambient temperature (Figure 8i).

Sensitivity to Vertical Resolution
To study the effect of coarser vertical resolution in three-dimensional numerical models, the water column of 150 m height is discretized with k max = 50 layers (10 layers over the upper 10 m) and k max = 25 layers (5 layers over the upper 10 m) instead of k max = 500 layers (38 layers over the upper 10 m). Such coarser resolutions are in the order of what three-dimensional models of ice-cavities using surface-following coordinates can typically afford. Given the fact that the initial plume thickness of D = 5 m is resolved by only a few discrete values, it is already quite coarse. The resulting vertical profiles (Figure 9) for k max = 500, k max = 50 and k max = 25 layers show that the coarse resolution simulations reproduce the high-resolution profiles with sufficient accuracy. Velocity,  (Figures 9a-9c) are reproduced very accurately, due to the largely resolution-independent ocean-to-ice flux parameterization (Section 3.1 and 3.2). Buoyancy frequency squared N 2 and shear frequency squared M 2 and thus the gradient Richardson number (Figures 9d-9f) are well reproduced at the ice-ocean interface and in the plume interior, but at the plume base, the sharp peaks of N 2 and M 2 are not properly resolved. For the coarse resolution with k max = 25 layers, the gradient Richardson number does not yield a value of Ri st = 0.25 inside the entrainment layer.
Resulting plume thickness, entrainment velocity and melt rate are shown for all three vertical resolutions in Figure 10. It can be seen that for k max = 50 layers the development of the plume thickness is still accurately reproduced with a maximum error of about 1 m (Figure 10a). For k max = 25 layers, due to the resolution of the initial plume thickness with only a few model layers, there is an overall underestimation of the plume thickness by about 3 m, which is however not significantly increasing during the simulation. After some deviations in the initial phase, also entrainment velocity and melt velocity are well-reproduced by both coarse resolution simulations (Figures 10b and 10c).

Comparison to Entrainment Parameterizations
Entrainment is the process of turbulent transport of relatively stagnant ambient water into the turbulent plume layer through its base with the consequence of an increase in plume thickness and density. Despite its complex hydrodynamics in a region of sharp vertical gradients, various parameterizations for the entrainment process have been successfully developed. Moreover, computationally efficient vertically integrated models of subglacial plumes have become a common tool in investigating subglacial melt processes in ice cavities (Hewitt, 2020;Holland & Feltham, 2006;Jenkins, 1991). Here, we first introduce these parameterizations and then compare their performance to the results of the vertically resolved model for the default scenario and the six sensitivity scenarios.
The entrainment velocity v e is defined as the rate of plume thickening due to entrainment of ambient water across the plume base, diagnosed from the plume thickness budget, with E 0 = 0.036. 2. The entrainment model by Jungclaus and Backhaus (1994) calculates the entrainment rate as the present simulations with this modified version showed a substantial underestimation of entrainment rates in comparison to other formulations and to the vertically resolved plume model, such that we retained the original value that Jungclaus and Backhaus (1994) adopted from Kochergin (1987). 3. For subcritical flow, Wells et al. (2010) proposed an entrainment rate that depends on the Froude number Fr, with the bulk mixing coefficient Γ = 0.2 , and 0.1 < D/l h < 0.3 with the characteristic horizontal turbulent length scale l h . 4. The entrainment rate proposed by Arneborg et al. (2007) with a = 0.084, b = 2.65, c = 0.6 and a constant drag coefficient C d = 0.0025 depends on the Ekman number K, additionally. 5. We finally recalibrate the approach by Arneborg et al. (2007) which has been optimized for dense bottom currents in the Baltic Sea using a constant drag coefficient C d =̂Fr̂̂, and formulating the mean-square error between diagnosed and predicted E/c d as with ′ = ln ( ) , where E i , ( ) , Fr i and K i are diagnosed values for all experiments and all time steps, a leastsquare method with ∕ ′ = ∕ = ∕ = 0 gives an optimal parameter set (̂̂̂) with a minimum error. The resulting values ( = 0.052 , ̂= 2.56 and = 0.29 ) are similar to the original values by Arneborg et al. (2007), but with a significantly smaller influence of the Ekman number which includes the variable roughness parameter c d . This model will be denoted as Present study in Figure 11.
As already shown in Figure 8, the entrainment velocity can vary over several orders of magnitude across the different scenarios. In Figure 11, we compare the non-dimensional entrainment rates diagnosed from the vertically resolved plume simulations with the values predicted by the five plume parameterizations given above. As input for the entrainment parameterizations we use bulk values for , ̄ and D diagnosed from the GOTM simulations, using Equation 57. The results for the entrainment rate differ substantially between the parameterizations. For some scenarios, the constant entrainment rate by Jenkins (1991) reproduces correctly the order of magnitude of the diagnosed entrainment, but obviously not its temporal evolution also during later balanced states of the 25 of 34 plume. The scenarios including variations of the slope angle, ayp and aym, show that the concept of formulating the entrainment rate as function of the slope angle in this simple parameterization roughly reproduces the correct order of magnitude of the entrainment process. The bulk Richardson number dependent parameterization by Jungclaus and Backhaus (1994) does largely follow the decreasing trends of the entrainment rate, but is for most scenarios generally significantly overestimating or underestimating the magnitude, a performance that is also seen for the parameterization by Wells et al. (2010) which is based on the fourth power of the Froude number. The performance of the Froude and Ekman number dependent parameterization by Arneborg et al. (2007) is generally better than those previously discussed, and its accuracy could be strongly improved by the recalibration to the present seven scenarios. Therefore, the good performance of the original Arneborg et al. (2007) calibration and the optimal performance of the newly calibrated formulation is not a surprise. This is also because the Arneborg et al. (2007) calibration used the same turbulence closure model as the present study. But the design of this parameterization depending on two non-dimensional plume parameters, Fr and K, that vary independently (see Figure 6) seems to be the most promising for rotational plumes.

Discussion
The discussion of the results of this study will concentrate on five issues that might be of interest for future modeling of ice shelves: the benefits of the analytical solution (Section 5.1), the implications of the numerical analysis (Section 5.2), the consequences of the vertical plume structure and resulting entrainment of ambient water (Section 5.4) and the remaining uncertainties in modeling of subglacial plumes and melt rates (Section 5.5).

Analytical Solution
The analytical solution for the vertical profiles of velocity, temperature and salinity of subglacial plumes that is presented as Equations 35-37 is based on a number of simplifying assumptions: Neglect of Earth rotation, vertically homogeneous acceleration of the plume, parabolic eddy viscosity and diffusivity, and stationarity. Despite its idealized character, the solution can be used for a number of purposes: 1. It can be used as a simple test bed for melt flux parameterizations that is not affected by numerical uncertainties. Despite the high degree of simplification, realistic values for melt rates and ocean-to-ice heat fluxes were calculated (Section 2.5.3). The analytical solution also shows the high degree of vertical homogeneity of the profiles of temperature and salinity and the significant differences between boundary and melt-layer values of temperature and salinity, due to the high Schmidt numbers ( Figure 3). Also, entrainment fluxes can be quantitatively compared to melt fluxes 2. The analytical solution is a basis to construct consistent and convergent formulations for the discrete melt layer fluxes (Section 3.1 and 3.2) 3. Finally, the analytical solution allows for the analysis of numerical convergence for vertically resolving plume models (Section 3.4)

Numerical Accuracy
The convergence analysis of a numerical plume model toward the analytical solution (Section 3.4) shows that ocean-to-ice fluxes and consequently melt rates can accurately be calculated also with a relatively coarse vertical resolution near the ice-ocean interface. The major requirement to achieve this is the proper discretization of the ocean-to-ice momentum flux. Using a drag coefficient that is independent of the upper-layer thickness leads to highly inaccurate and divergent results, due to the strong velocity gradients near the interface. In contrast, discrete tracer flux formulations that are independent of the upper-layer thickness are numerically inconsistent. Still, they do not lead to a measurable loss of accuracy when the vertical tracer profiles are quasi-homogeneous.
Vertical resolution does however matter for the numerical reproduction of the entrainment process at the base of the plume, where strong vertical gradients are present. A typical vertical extent of the entrainment layer is of the order of 2-4 m (Figure 7). Although second-moment turbulence closure schemes are relatively robust with respect to vertical resolution (Li et al., 2021;Umlauf & Burchard, 2005), vertical grid resolutions should be of the order of 2 m or higher in the region of the entrainment layer. This is demonstrated in Figure 9 where an entrainment-layer resolution of about 1.5 m still properly reproduces the vertical plume structure of a high-resolution 27 of 34 model at the end of a 2-week simulation. Also entrainment velocity and melt rates are sufficiently reproduced by this resolution, but for a resolution of 3 m, results start deteriorating ( Figure 10). For typical plume thicknesses of the order of 10 m, this means that more than 5 numerical layers should be present in the plume region.

Implementation Into Climate Models
The parameterizations and numerical schemes developed and tested here can be implemented into climate models, also when coupling between ocean, cryosphere and atmosphere is included. This is feasible, because the turbulence closure schemes used here are in terms of computational efficiency comparable to schemes typically used in climate models such as KPP (Li et al., 2021). The remaining key question is if a vertical resolution of the order of 2 m can be achieved in the vicinity of the ice-ocean interface. For geopotential coordinates, this might mean an overall vertical resolution of 2 m at all depths where the ice-ocean interface is present, that is, typically several 100 m below the undisturbed mean sea level. For surface-following coordinates (where vertical resolution decreases with water depth) such resolutions are still quite a challenge, but using non-linear S-coordinates (Burchard et al., 1997;Shchepetkin & McWilliams, 2005) should make this feasible. The use of surface-following coordinates in global models including ice shelves has already been demonstrated (Timmermann et al., 2012). A more flexible solution could be the concept of vertically adaptive coordinates (Gräwe et al., 2015;Hofmeister et al., 2010), since they allow concentrating the resolution at sharp density interfaces such as in the entrainment layer. This principle has been used by Umlauf et al. (2010) for the simulation of channelized dense bottom currents and could also be applied to subglacial plume simulations.

Vertical Plume Structure and Entrainment
The vertical structure of subglacial plumes is believed to resemble that of dense bottom currents, turned upside down (Jenkins, 2016). That explains why often model tools are applied to subglacial plumes that have been developed for dense bottom currents and oceanic overflows. In both cases, the determination of entrainment rates is important. For oceanic overflows, the entrainment rates determine their potential for ventilating the deep ocean, for subglacial plumes they determine the transport of relatively warm and salty ambient water toward the ice-ocean interface and therefore play an important role in setting the melt rate. The fundamental dynamical difference between dense bottom currents and subglacial plumes is the additional (stabilizing) interfacial buoyancy flux due to melt processes at the ice-ocean interface of subglacial plumes that ultimately drives the flow.
Entrainment velocities analyzed from the present model results are within the range of common entrainment parameterizations typically used in vertically integrated plume models of subglacial plumes ( Figure 11). Those parameterizations however lead to very different estimates of entrainment rates. To provide a robust and reliable entrainment parameterization, the formulation by Arneborg et al. (2007) has been re-calibrated for a depth-dependent drag coefficient Equation 69 and shows an agreeable accuracy over a large range of plume parameters.

Remaining Uncertainties
The vertical structure of dense bottom currents has been well-observed in the ocean, for example, the Faroe Bank Channel overflow (Fer et al., 2010) and overflows into the Baltic Sea . Such high-resolution observations of subglacial plumes are not available due to the thickness of the glacial ice cover and the large surface area of the floating ice tongues. The few available observations (e.g., Washam et al., 2020), do not provide sufficient resolution and coverage of the entire plume thickness. It can therefore only be assumed by analogy how subglacial plumes are vertically structured. In agreement with observations and models of dense overflows in the Western Baltic Sea Umlauf et al., 2007Umlauf et al., , 2010 with similar characteristics in terms of similarity parameters such as the Froude and Ekman numbers, the interior of subglacial plumes can be assumed to be well-mixed and the plumes can be assumed to be bounded by sharp density interfaces in the entrainment layer. However, the uncertainties about the vertical structure of subglacial plumes can only be reduced by very challenging in-situ field observations under glacial ice tongues.
Moreover, the underlying assumption for one-dimensional water column models of a plume that is laterally homogeneous along the ice-ocean interface is highly idealized and unrealistic. The ice-ocean interface may be smooth on the small scale, but it is plausible to assume non-negligible sub-grid scale roughness that exerts a 28 of 34 form drag on the flow. This drag needs to be parameterized in coarse ocean models via an effectively roughness length k s (Section 2.4). Our model experiments show that the roughness length has a large influence on plume thickness, entrainment velocity and melt rate ( Figure 11). Since the determination of this effective roughness is highly uncertain, sensitivity studies with respect to this parameter are recommended.
The biggest uncertainty however remains for the detailed structure of the topography of the ice-ocean interface and the effects on the plume dynamics. Subglacial plumes are often thought of as wide layers of buoyant water propagating across ice tongues with relatively plain subglacial topographies (Holland & Feltham, 2006). It is however known that plumes occur as highly channelized flows, similar to overflows in the ocean (Fer et al., 2010;Umlauf et al., 2007), but probably occurring on even smaller scales (Rignot & Steffen, 2008;Washam et al., 2020). To investigate the impact of these topographically highly diverse plume dynamics on the net basal melt rate, detailed studies of more dimensional flows need to be carried out such as cross-sections models of channelized subglacial plumes, similar to the two-dimensional model of dense bottom currents presented by Umlauf et al. (2010).

Conclusions
A numerical one-dimensional water column model of subglacial plumes has been presented here that should help to constrain ocean models of ice shelves. To our knowledge this is the first high-resolution one-dimensional model that couples the physics of the melt layer to second-moment turbulence closures inside the plume and across the entrainment layer. This modeling strategy allows for quantitative predictions of entrainment processes of ambient water into the plume, such that it can serve as a benchmark for models with simpler physics such as bulk entrainment models.
Using an analytical solution of the plume, accurate and convergent numerical expressions for fluxes across the ice-ocean interface are formulated. Specifically, they do also reproduce these fluxes accurately for relatively coarse near-interface resolution. The probably most critical finding is that the vertical model resolution in the region of the entrainment layer should ideally be order of 2-3 m, which provides a challenge to existing ocean models for ice shelves. Future efforts should be directed at developing flexible numerical schemes that allow this locally high resolution also in large scale ocean models.

A1. Velocity Profile
Under the conditions discussed in Section 2.1, and ignoring Earth rotation, the momentum budget Equation 9 within the plume has the following form: where * is the friction velocity in the entrainment layer at the base of the plume. Assuming stationarity of the velocity profile and combining Equations A1 and A2, a linear stress profile is resulting: the parabolic profile of eddy viscosity is given as see Burchard and Hetland (2010), where for small z′ Equation 20 is retained. Combining Equations A3 and A4 gives BURCHARD ET AL.
the plume-averaged velocitȳ= combining Equations A7 and A9 finally gives the velocity profile fulfilling a prescribed depth-averaged velocity , see Equation 35.

A2. Tracer Profile
Let c be a tracer obeying the following one-dimensional budget equation: with the upward turbulent tracer flux, where the location of the boundary value for c at the ice-ocean interface is slightly shifted with respect to the no-slip boundary condition for velocity (see Section 2.1, 2.2 and 2.4). As flux boundary conditions, we define an eddy diffusivity profile is constructed by dividing the parabolic eddy viscosity profile Equation A4 by the the turbulent Prandtl number Pr t : 30 of 34 assuming that ∂ t c is independent of z, the tracer flux will have the following linear profile: combining Equations A11 and A16 gives the rate of change of c: with the vertically averaged tracer concentration and the integration constant in Equations A20 -A22, the formulations including β c instead of 0 , using Equation 32, are those that should be used for computations, since 0 is typically so small that calculations of its reciprocal would result in overflows. Combining Equations A20 and A21 leads to a tracer profile based on the depth mean tracer concentration instead of the melt layer tracer concentration c b .
BURCHARD ET AL.
The dependence of the non-dimensional current speed and the non-dimensional friction velocity * = 1∕2 on the non-dimensional interfacial roughness is shown in Figure B1 for three different values of the non-dimensional buoyancy forcing ̃s in . Expectedly, a larger roughness length leads to an decreased plume velocity, however, it leads to an increased friction velocity indicating a more turbulent plume. This can be explained by the fact that the decrease of current speed due to increased roughness 0 is smaller than the corresponding increase in 1∕2 . This non-linear effect is not reproduced by models that use a constant drag coefficient.

Data Availability Statement
The model code, the configuration file for the subglacial plume simulation as well as compile and execution instructions can be found at https://doi.org/10.5281/zenodo.6203838.