High‐Resolution Simulations of the Plume Dynamics in an Idealized 79°N Glacier Cavity Using Adaptive Vertical Coordinates

For better projections of sea level rise, two things are needed: an improved understanding of the contributing processes and their accurate representation in climate models. A major process is basal melting of ice shelves and glacier tongues by the ocean, which reduces ice sheet stability and increases ice discharge into the ocean. We study marine melting of Greenland's largest floating ice tongue, the 79° North Glacier, using a high‐resolution, 2D‐vertical ocean model. While our fjord model is idealized, the results agree with observations of melt rate and overturning strength. Our setup is the first application of adaptive vertical coordinates to an ice cavity. Their stratification‐zooming allows a vertical resolution finer than 1 m in the entrainment layer of the meltwater plume, which is important for the plume development. We find that the plume development is dominated by entrainment only initially. In the stratified upper part of the cavity, the subglacial plume shows continuous detrainment. It reaches neutral buoyancy near 100 m depth, detaches from the ice, and transports meltwater out of the fjord. Melting almost stops there. In a sensitivity study, we show that the detachment depth depends primarily on stratification. Our results contribute to the understanding of ice–ocean interactions in glacier cavities. Furthermore, we suggest that our modeling approach with stratification‐zooming coordinates will improve the representation of these interactions in global ocean models. Finally, our idealized model topography and forcing are close to a real fjord and completely defined analytically, making the setup an interesting reference case for future model developments.

Northeast Greenland Ice Stream (Kappelsberger et al., 2021;Schaffer et al., 2017), holding 1.1 m sea-level equivalent of ice (i.e., its ice could lift global sea levels by 1.1 m if melted entirely; Christmann et al., 2021).Schaffer et al. (2020) estimated that 89% of the meltwater leaving the 79NG fjord comes from subglacial melting caused by the ocean.Ice melting on land or at the surface only accounts for the remaining 11% of 79NG meltwater (and even less at other glaciers, see Rignot & Steffen, 2008), which is discharged into the fjord as subglacial runoff at the grounding line.Subglacial melting thins the glacier tongue, which can reduce the buttressing of the ice sheet, that is, the support of the grounded glacier that is provided by the friction between the ice tongue and the lateral fjord boundaries (Goldberg et al., 2009).With a thinner ice tongue, thus less buttressing, the glacier can flow faster into the ocean, which contributes to sea level rise (Goldberg et al., 2009;Humbert et al., 2022;Shepherd et al., 2004).Furthermore, basal melting can destabilize the ice tongue, which can lead to its breakup (Rignot & Steffen, 2008); in consequence, a lot of ice would be discharged into the ocean (Shepherd et al., 2004).This exemplifies the big role of the ocean in melting the Greenland Ice Sheet (Schaffer et al., 2017) and shows that it is important to understand ice sheet-ocean interactions in glacier fjords like the one at 79°North.
The general idea of ice-ocean interactions under a glacier tongue in Greenland is as follows: Atlantic Intermediate Water (AIW) flows over a sill at the fjord entrance into the glacier cavity as a dense, saline, and warm bottom plume.AIW brings heat into the ice cavity, which is used for melting.The meltwater forms a buoyant plume on the underside of the glacier tongue.This plume causes subglacial melting, transports glacially modified water toward the open ocean, and constitutes the return flow of an overturning circulation within the fjord (Schaffer et al., 2020;Straneo & Cenedese, 2015).
The dense bottom plume and the buoyant subglacial plume are two major processes in a glacier cavity.However, they are difficult to study, because measurements in Greenland's fjords are generally sparse (Straneo & Cenedese, 2015), particularly under floating ice tongues, where the ocean is inaccessible to ships and unobservable by satellites.Ice tethered moorings (Lindeman et al., 2020) give some information about processes under the ice, but only at single positions.So numerical models in combination with measurements are necessary to gain a detailed understanding of ice sheet-ocean interactions.This requires that the model formulations properly incorporate the classical relations for stratified boundary layer flow.Baines (2008) distinguished between two such features: gravity currents and plumes.Gravity currents have relatively gentle slopes; they are characterized by sharp interfaces and a balance between buoyancy force and bed friction.These gravity currents show detrainment and intrude into the ambient water.In contrast to that, plumes exist on steeper topography; the buoyancy force is balanced by strong entrainment of ambient water.We will show that both states, gravity currents and plumes, exist under the 79NG ice tongue at different locations.For the turbulence closure model used here, Arneborg et al. (2007) showed that it well resolves entrainment rates of dense bottom currents in the Baltic Sea.This is due to the fact that the turbulence closure model has been properly calibrated to reproduce a steady-state Richardson number of 0.25 (Burchard & Baumert, 1995;Umlauf & Burchard, 2005) and a mixing efficiency of 0.2 (Burchard & Hetland, 2010;Umlauf, 2009).
A challenge for ocean models is to provide sufficiently high resolution in a glacier fjord to accurately simulate the two plumes.The melt rate computed by the model also depends strongly on the vertical resolution (Gwyther et al., 2020).It has been shown that the subglacial plume and particularly its entrainment layer require a vertical resolution of about 1 m or better to correctly model the plume development and the associated melting (Burchard et al., 2022).This is hard to achieve in most ocean models, because of the stark contrast in vertical scales between the fjord depth of several hundred meters and the plume thickness on the order of 1-10 m.
With the vertical coordinates that are commonly used in ocean models, it is unfeasible to achieve a resolution of 1 m along the whole subglacial plume.At 79NG, the plume starts at the grounding line at 600 m depth, so z-coordinate models (Losch, 2008) would require at least 600 vertical layers to resolve the top 600 m of the water column with a 1 m-resolution-much more than can typically be afforded in global models.Such a resolution is currently only feasible at the fjord scale, as shown in the semi-realistic model by Xu et al. (2013) for a Greenlandic glacier without an ice tongue.With σ-coordinates (Gwyther et al., 2020;Timmermann et al., 2012), a high resolution along the whole ice tongue is possible with less layers by activating a zooming toward the ice-ocean interface, thereby obtaining a finer vertical grid in the boundary layer.However, these terrain-following coordinates have problems when calculating the internal pressure gradient over steep topographic slopes (Burchard & Petersen, 1997;Haney, 1991), which are a typical feature in glacier fjords.
With adaptive vertical coordinates (AVC; Hofmeister et al., 2010), the described problems can be considerably reduced.AVC are terrain-following coordinates that allow with a moderate number of layers a high vertical 10.1029/2023MS003721 3 of 28 resolution in places of interest.By minimizing a cost function, AVC adapt automatically to features like stratification, shear, and interfaces (Burchard & Beckers, 2004).This reduces numerical mixing (Gräwe et al., 2015;Klingbeil et al., 2014) and puts more layers in places where more details need to be resolved, while permitting less vertical resolution in more uniform areas.These coordinates have been used successfully for simulating dense and buoyant plumes in various conditions (e.g., Chegini et al., 2020;Umlauf et al., 2010), but not yet for glacier fjord modeling.
In this paper, we will show that the stratification-zooming feature of AVC is useful for modeling the ocean circulation under ice shelves, because this provides a high vertical resolution of 1 m in the subglacial plume and the bottom plume with feasible computational cost.Furthermore, we will present the new insights into the plume dynamics that were gained by using a model that provides such resolution.
We created an idealized 2D-vertical simulation of the 79NG fjord using AVC together with a melt parametrization (Burchard et al., 2022) that is suitable for high vertical resolution at the ice-ocean interface.To our knowledge, this is the first model to use stratification-zooming coordinates like AVC in a glacier cavity.In addition to testing the performance of AVC under an ice tongue, we use our model to study the sensitivity of the 79NG system to environmental influences.With the 20 scenarios of our sensitivity study, we analyze the effect of the salinity and temperature stratification of the ambient ocean, test the importance of the subglacial discharge, examine the role of the sill, and investigate the influence of roughness or smoothness of the ice tongue.This paper is organized as follows.The following Section 2 describes our model setup, compares it to the real 79NG fjord, explains our modeling choices including AVC, and describes our analysis methods.Section 3 shows the results of our default simulation (Section 3.1), the performance of AVC (Section 3.2), and the results of our sensitivity study (Section 3.3).This is followed in Section 4 by a discussion of the physical processes we observe in all our numerical experiments and what we learn from these findings about ice-ocean interactions in glacier cavities.We also discuss the role of AVC in obtaining the presented results.Some conclusions and an outlook are given in Section 5. Appendix A lists the mathematical expressions used to build our idealized setup, so that our model can serve as a reference test case for future model developments.

Idealized 2D Model of the 79°N Glacier Fjord
We built an idealized numerical ocean model of the 79° North Glacier (79NG) fjord located in Northeast Greenland, using GETM, the General Estuarine Transport Model (Burchard & Bolding, 2002).This model is suitable for our purpose, because 1. GETM comes with AVC that allow high vertical resolution in areas of interest for low computational cost (Section 2.3); 2. GETM includes state-of-the-art vertical turbulence closure with GOTM (Burchard et al., 1999;Li et al., 2021;Umlauf & Burchard, 2005); 3. GETM has been developed specifically for the coastal ocean and estuaries (Klingbeil et al., 2018).
A glacier fjord is a special type of estuary, in which the subglacial discharge plays the role of a river in a classical estuary (Muilwijk et al., 2022;Straneo & Cenedese, 2015).However, the main source of freshwater is not the subglacial discharge, but the subglacial melting of the floating ice tongue (Schaffer et al., 2020).Since this is the first time that GETM is used for simulating a glacier fjord, we extended the model to include ice tongues and basal melting.The details of this new GETM feature are explained in Section 2.2.
The GETM setup presented here is a two-dimensional (x, z) representation of the 79NG fjord with high resolution along the fjord (x) and in the vertical (z), but integrated in cross-fjord direction (y).The fjord circulation is expected to vary also across the fjord (Lindeman et al., 2020), so a 2D model is a simplification and we discuss its implications in Section 4.However, our 2D approach is a useful starting point, as it has the necessary complexity to learn a lot about the plume dynamics and the overturning circulation in the ice cavity.
We consider the main glacier terminus of 79NG, without the adjacent Dijmphna Sund (Figure 1a).The ice tongue is about 75 km long and 20 km wide; our model has the same width (L y = 20 km) and twice the length (L x = 150 km), to have a sufficiently large "buffer" between the glacier cavity-which is our main interest-and the open ocean boundary.We resolve the domain with 300 water columns in x-direction (Δx = 500 m) and one 10.1029/2023MS003721 4 of 28 grid point in y-direction; the resolution in z-direction with 100 adaptive layers is explained in Section 2.3.At this horizontal resolution, neither nonhydrostatic effects associated with the plumes nor nonhydrostatic internal waves are resolved, so it is appropriate to use GETM in hydrostatic mode (Klingbeil & Burchard, 2013).
To construct the bottom topography of our idealized 79NG model, we look at two data sets (Figure 1).The seismic depth soundings by Mayer et al. (2000) are the most accurate measurements of the bathymetry in the part of the fjord that is covered by the ice tongue.The retreat of the ice tongue in recent decades facilitated more detailed bathymetry measurements near the fjord entrance.In their bathymetric survey, Schaffer et al. (2020) showed that the fjord is separated from the open ocean by a sill that is 325 m deep on its deepest point.Since this sill depth is not representative for the whole width of the fjord (Figure 1a), we use a shallower sill in our idealized 2D model (Figure 1b).It is at 300 m depth in our default setup; in our sensitivity study, we analyze the effect of the sill by varying its depth from 200 to 450 m (Section 3.3.4).
The bathymetry of our idealized model is a smooth concatenation of simple, analytical functions (Figure 1b): A third-order polynomial connects the grounding line (600 m depth) with the deepest point in the trough (900 m) and continues until it reaches a slope of 2.5%.It is then connected linearly to the parabola forming the sill with its maximum (300 m) at 80 km from the grounding line.The parabolic sill goes over into an exponentially decreasing shelf that converges toward a depth of 450 m far offshore.The mathematical details are given in A1 in Appendix A. In our sensitivity experiment without a sill, the linear slope is directly connected with the exponential shelf.Apart from the sill, our model bathymetry only differs markedly from the measured section between the grounding line and the trough (Figure 1b).Despite this difference, we think that a simpler bathymetry with fewer parameters is preferable to a perfect fit to a single transect for an idealized model such as ours.Also, this deep part of the fjord is mostly inactive in our simulations.At the grounding line, which forms the left/western boundary of our model (x = 0), subglacial discharge enters the glacier fjord.This runoff is implemented in our GETM setup like river input.It is added as freshwater at the local freezing point (−0.45°C, which is less than 0°C due to pressure) to the first water column.We take a constant discharge rate of 70 m 3 s −1 (equivalent to 0.07 mSv reported by Schaffer et al., 2020) in our default setup and varied this value in our sensitivity study (Section 3.3.3).The discharge is distributed uniformly over the whole water column, which is about 6.3 m thick at the first grid center.
At the open boundary on the right/eastern end of our model domain (x = L x = 150 km), we prescribe the surface elevation η and the ambient ocean stratification.For the former, we use a constant zero elevation.We also tested  2021), but our experiments showed that the melt rate is relatively unaffected by the tidal forcing.This is because in the deep part of the fjord, where the subglacial plume causes melting, the speed of the tidal currents is less than 0.01 m s −1 in absolute value, thus much smaller than the plume velocity of about 0.2 m s −1 .Therefore, the tide is neglected in the present study.Regarding the open boundary stratification, we use idealized and constant-in-time profiles of temperature and salinity.They are specified by T-and S-values at sea level, at 100 m depth, and at 300 m depth (shown in Figure 2 and listed in A3 in Appendix A), using linear interpolation in between and constant extrapolation below.In our default setup, the resulting profiles are close to CTD measurements by Schaffer et al. (2020), see the comparison in Figure 2. We also perform a sensitivity study with modified stratifications (Sections 3.3.1 and 3.3.2).
The model is initialized at rest with a homogeneous stratification equal to the stratification at the open boundary.We run the model with a timestep for the barotropic mode of Δt = 5 s, in accordance with the CFL stability criterion, demanding  Δ ≤ Δ∕ √ max ≈ 5.3 s (using H max = 900 m as the maximum depth of the fjord, see Figure 1, ignoring the ice cover).We use a split factor of M = 3, so that the baroclinic mode is computed every Δt 3D = MΔt = 15 s.While our default setup can be run with a larger baroclinic timestep, the highmelt scenarios give smoother results with a higher temporal resolution, so we decided to use this split factor for all our runs.For the turbulence closure, our setup uses the k-ϵ model with quasi-equilibrium second-moment closure (Cheng et al., 2002), implemented in GOTM.We activated divergence damping with a diffusion of A n = 50 m 2 s −1 on barotropic transports for a conservative smoothing of the sea surface (Vallis, 1992).After a few simulation months, our model approaches a quasi-steady state, in which melting and circulation are almost time-independent.The results shown in this paper are 24 hr-averages taken at the end of a 6-month simulation and represent the steady state.

Implementation of Glacier Ice in GETM
For this study, we added a new feature to GETM that allows simulations of glacier fjords covered by an ice tongue.Where the ice tongue is present, it adds additional pressure (Section 2.2.1), friction (Section 2.2.2), and melt fluxes (Section 2.2.3) to the sea surface.Our implementation allows the ice to move freely vertically, for example, with long waves, but it is fixed horizontally.Calving is not included in our model.
In this paper, we use the term sea surface to refer to the (moving) upper boundary of the ocean, denoted η = η(x, t) and measured from z = 0 with positive values upwards.Depending on the x-position, the sea surface can be the ice-ocean interface or the atmosphere-ocean interface.Furthermore, we use the term sea level to refer to the level z = 0, which is the initial position of the atmosphere-ocean interface.

Pressure Due To Ice and Initial Sea Surface Elevation
Under glacier ice, the pressure at the ice-ocean interface is the atmospheric pressure (constant in our model) plus the contribution from the weight of the ice tongue (Losch, 2008).We can represent this pressure due to floating ice as p i = gρ i h i , where h i is the thickness of the ice column and ρ i its (homogeneous) density (Table 1).Both h i (x) and ρ i are constant-in-time in our implementation and serve as input parameters to the model.
To initialize our model in an equilibrium state, we must prescribe the initial surface elevation η 0 = η(t = 0) such that the ocean with the floating glacier ice is in hydrostatic balance.This is the case if the water displaced by the ice tongue has the same weight as the ice tongue (Archimedes' principle).For an initially horizontally homogeneous stratification with (water) density ρ(z), this condition can be expressed as: North Glacier (79NG) fjord.The shaded area marks the minimum and maximum values tested in our sensitivity study.The CTD profile was taken in 2017 on RV Polarstern (Kanzow et al., 2018) and represents a typical ambient ocean stratification for 79NG (Schaffer et al., 2020, see their Figure 1a for the location of the profile).The freezing point of saline water in (b) corresponds to the shown CTD profile.We used the Python package gsw (TEOS-10;IOC et al., 2010) to convert from the CTD data pressure to depth, practical salinity to Absolute Salinity, and potential temperature to Conservative Temperature, as well as to compute the freezing temperature.
In our setup, we prescribe the lower ice edge η 0 (see below) and determine h i such that Equation 1 is fulfilled, which has the consequence that we have slightly different ice thicknesses h i for different stratifications ρ(z) (difference to the default setup is always less than 20 cm).Note that a corollary of Equation 1 is the handy rule-of-thumb η 0 ≈ −0.9h i , which says that 90% of an ice column is below sea level and 10% is above.
Given the initially horizontally homogeneous (and stable) ocean stratification, we initialize the ice in the chosen equilibrium position by evaluation of the integral in Equation 1, which yields the ice thickness and thus the pressure loading.This pressure loading is maintained throughout the simulation.However, as the simulation runs, the stratification changes due to basal melting, subglacial discharge, ambient water inflow, and mixing, so the equilibrium position of the ice changes as well.Since the ice in our model can move freely vertically with the convergence and divergence of transports, it will adapt to the changing stratification.The setups presented here reach a quasi-steady state, in which the glacier tongue has found a new equilibrium position, which is slightly (on the order of millimeters) different from the initial position.
In our idealized 79NG fjord model, we prescribe a smooth ice-ocean interface between the grounding line at 600 m depth and the calving front at x = 75 km, where the ice-ocean interface is 75 m below sea level.For the idealized ice shape, we choose a hyperbolic tangent with a maximum slope of 2.5% at the grounding line (see A2 in Appendix A for the mathematical details).This fits well with the measured ice slope near the grounding line (see the reference section in Figure 1b).Since subglacial melting is strongest in this area (Schaffer et al., 2020), we believe it is important to reproduce the ice topography well near the grounding line and accept that the idealized shape differs from observations at mid-depths, as we prefer a simple, analytical ice shape over a perfect fit to a single transect.
The calving front, which in reality is an almost vertical wall, is a challenging part of the model domain.If it was modeled as a vertical wall or as a steep slope, the uppermost terrain-following coordinate levels would follow this slope.However, the water near the calving front is strongly stratified (Figure 2a), so individual grid cells would span a large density range.This could lead to numerical mixing and spurious flows (Gwyther et al., 2020).Therefore, we extend the ice-ocean interface with a linear 1%-slope until sea level is reached.We also tested higher slopes at the calving front, but the model results were poorer, because the horizontal flow below the calving front was too much diluted by passing through too many cells.Thus, we use a slope instead of a vertical wall at the calving front.This is a deviation from the real system, but an acceptable one, since our focus lies on processes inside the glacier cavity, which we assume to be not much affected by this difference.

Surface Friction
Where the ocean is covered by glacier ice, there is a no-slip boundary condition at the sea surface (Burchard et al., 2022).This friction at the ice-ocean interface is implemented in GETM according to the law-of-the-wall with a roughness length z 0,ice , similar to bottom friction.In our default scenario, we use the value z 0,ice = 0.01 m.The effects of smoother or rougher ice are tested in our sensitivity analysis (Section 3.3.5).

Parametrization of Subglacial Melting
We implemented the subglacial melt formulation by Burchard et al. (2022).This parametrization, based on the three-equation model (D.M. Holland & Jenkins, 1999), is suitable for high vertical resolutions under the ice.In our free-surface model, meltwater is added like precipitation as a real freshwater flux (Huang, 1993) to the uppermost grid cell of the water column with a melt rate v b (in m s −1 ).There is no salt flux,     = 0 , because the melted glacier ice is assumed to have zero salinity.There is, however, a temperature flux at the ice-ocean interface: (2) Heat capacity of glacial ice

Table 1 Settings and Parameters of Our Model in the Default Scenario
In the squared bracket, the first term corresponds to the energy necessary for heating up the glacial ice from its core temperature T i to the melt layer temperature θ b ; the second term is the latent heat of the phase change from ice to water; the last term appears because water is exchanged between ice and ocean, that is, the ice-ocean interface is a nonmaterial interface in our model (Jenkins et al., 2001).The values of the constants in Equation 2 are given in Table 1.The melt layer is a thin layer at the ice-ocean interface, that is not resolved but parametrized in our model.For a detailed discussion, see Burchard et al. (2022).
The here-described implementation of melting differs from that used by Burchard et al. (2022), because their 1D model has a rigid lid.In a rigid lid model, the water volume cannot increase, so a virtual salt flux through the ice-ocean interface is needed to get the diluting effect of basal melting on salinity, and a virtual temperature flux is needed instead of Equation 2. However, the more realistic approach is adding meltwater explicitly (Huang, 1993;Jenkins et al., 2001), without a salt flux and with only a real temperature flux, as we do it in this study.Even though melting increases the water volume in our model, the ice volume does not decrease.To allow for a decreasing ice volume and a thinning ice tongue, ice dynamics would have to be modeled as well.Instead, we assume that there is a balance between basal melting of the ice tongue and the discharge of glacier ice from land into the ocean.

Adaptive Vertical Coordinates
Our GETM setup uses AVC described by Burchard and Beckers (2004) and Hofmeister et al. (2010).These coordinates are well-suited for representing surface-attached buoyant plumes (Chegini et al., 2020) and dense bottom currents (Hofmeister et al., 2010;Umlauf et al., 2010).AVC are topography-following coordinates, in which the vertical distribution of the model layers changes with time.The temporal change of model layers is implemented by minimizing a cost function depending on the model state, particularly the stratification.The coordinates adapt in a way that there are more layers in parts of the water column with higher stratification.This ensures high vertical resolution in areas of strong vertical density gradients and minimizes numerical mixing (Gräwe et al., 2015;Hofmeister et al., 2010;Klingbeil et al., 2014).
In the 79NG fjord, important density differences exist in two locations: (a) Between the meltwater plume at the ice-ocean interface and the ambient water below, and (b) between the bottom current and the cavity water above (Schaffer et al., 2020).With AVC we can obtain high resolutions in both of these plumes and particularly in their entrainment layers, without a large increase in computational cost (<10% more computation time compared to σ-coordinates).For this, we configured AVC so that they zoom toward stratification and toward the sea surface.
An explicit bottom-zooming is not required, because the stratification-zooming itself provides sufficiently high resolution in the bottom plume (Section 3.2).Activating the zooming toward the sea floor would also result in high resolution in the deep trough of the glacier fjord and on the continental shelf outside the ice cavity, even though these parts are mostly inactive in our simulations.Thus, we do not activate it and opt instead for an even higher resolutions near the ice-ocean interface, which is important for the accurate representation of melting (Burchard et al., 2022).While 50 coordinate levels would be sufficient to achieve a vertical resolution better than 2 m in the plume under the ice, we present in this paper simulations with 100 AVC layers to show the plumes and the circulation in great detail.

Analysis of Plume-Averaged Quantities
To analyze the entrainment of the subglacial plume, we compute its bulk properties, that is, the vertically averaged plume characteristics, in particular the plume thickness D, its buoyancy  b , and its velocity    .We want to diagnose the bulk values following the ideas by Arneborg et al. (2007) in the modified form for plumes under ice shelves (Burchard et al., 2022): where z′ = η − z is the distance from the ice-ocean interface η, b(z) = −g[ρ(z) − ρ 0 ]/ρ 0 is the buoyancy, and ρ 0 is the ambient ocean density.However, the above equations have been derived in a 1D setting with the assumptions that the ambient water below the plume is homogeneous (with density ρ 0 ) and stagnant (u = 0), which is not the case in our 2D model.So an integration to −∞ or to the sea floor at z = −H would not make sense, because it would include several different water masses in the plume analysis.Instead, we choose an integration depth h 0 > 0, consider the water mass at z = η − h 0 as the ambient water, and use the following modified formulas: Dividing Equation 7by Equation 6gives the plume thickness D, dividing Equation 8 by Equation 6gives the plume velocity    , and dividing Equation 6 by D gives the plume buoyancy  b .We take as ρ 0 the density linearly interpolated from cell centers to z = η − h 0 ; a vertical interpolation gives considerably smoother graphs for the bulk values than taking the density of the grid cell containing z = η − h 0 .The factors of b(z) and  b in Equation 8ensure that the integral gives more weight inside the plume than outside, where b(z) is smaller since the local density ρ(z) is closer to that of the ambient water, ρ 0 .We use velocities horizontally interpolated to cell centers (instead of cell interfaces) in Equation 8, so that all bulk values are defined on cell centers.
Following P. R. Holland and Feltham (2006), the bulk values can be used to formulate a conservation equation for the plume volume: where the terms on the right-hand side are the melt rate v b and the entrainment velocity v e .For our 2D system (∂ y = 0) in steady state (∂ t = 0), Equation 9implies which means that the plume thickness increases in x-direction by flow convergence  (− ū) , melting, and entrainment (Jenkins, 1991).Since the melting is computed by our numerical model, we can reformulate (10) to diagnose the entrainment (Burchard et al., 2022): (11) To further analyze the dynamics of the plume, we compute the Froude number which is a non-dimensional number relating the velocity of the plume to the phase speed of long waves at the plume interface (Arneborg et al., 2007;Burchard et al., 2022).In flows that are dominated by a balance between buoyancy and friction, with little acceleration, the approximation holds (Arneborg et al., 2007), where tan α = ∂η/∂x is the slope of the topography and 10.1029/2023MS003721 9 of 28 is the drag coefficient of the subglacial plume (Burchard et al., 2022) with the van Karman constant κ = 0.4.
The choice of the integration depth h 0 requires some considerations.It must be chosen such that (as long as the plume is attached to the ice) z = η − h 0 lies always outside the plume in a weakly stratified region, but not too far away, so that ρ 0 = ρ(z = η − h 0 ) is actually the density of the water surrounding the plume.To find a suitable integration depth, a visual inspection of the model result is helpful.The identified value is a good choice if the computed bulk values are insensitive to small variations of h 0 .In our default scenario, this is the case for h 0 = 10 m.However, the precise choice of h 0 is not critical for the results.
For the analysis of the dense bottom plume, we use an analogous approach, but with integration from the sea floor at z = −H to z = −H + h 0 , and with z′ = H + z being the distance from the sea floor in Equation 7. As integration height h 0 , we take h 0 = 30 m downstream of the sill and h 0 = h 0 (x) = 30 m + H(x) − H(x sill ) upstream of the sill, where x sill = 80 km is the position and H(x sill ) = 300 m the depth of the sill.This way, on the upstream side of the sill, the integration goes from the sea floor to a constant level of z = −270 m, which is the depth that separates the inflowing water mass below from the outflowing water mass above.Like for the subglacial plume, the precise choice of h 0 is not critical.

Analysis of the Overturning Circulation
A key property of a glacier fjord is the strength of its overturning circulation, often reported in milli-Sverdrup (1 mSv = 1,000 m 3 s −1 ).We take as a measure of the overturning strength the maximum (in absolute value) of the (volume) stream function over the sill (x = 80 km).This value is smaller than the overall maximum of the stream function, which is reached in the interior of the cavity, but it allows the comparison of our results with measurements near the calving front (Schaffer et al., 2020).Since the overturning in the cavity is stronger than over the sill, the term exchange flow might be more suitable than overturning strength, but we use the latter for consistency with the literature.The stream function ψ is defined by and the condition that ψ = 0 on the sea floor; L y is the (constant) width of the fjord (Table 1).Numerically, we diagnose ψ by summing the horizontal transports uΔzL y (defined on cell edges) from the sea floor to the sea surface, which follows from Equation 15 and naturally satisfies ψ = 0 at the bottom.Then Equation 16 is automatically fulfilled thanks to the 2D continuity equation, ∂ x u + ∂ z w = 0. Since the model results shown in this paper are in steady state, the contour lines of the stream function ψ are trajectories.

Results
In this section, we present at first the steady state of our default scenario (Sections 3.1 and 3.2), then we perform a sensitivity study with varying physical parameters (Section 3.3).

Circulation and Melting in the Default Scenario
In our default model setup, which is an idealized representation of the present day situation at 79NG as observed by Schaffer et al. (2020), we find an estuarine-like circulation in the glacier cavity (Figures 3a-3d).This circulation is made up of two gravity plumes: strong, turbulent, and focused currents that are driven by density differences.
One is a buoyant plume at the lower ice edge, driving the melting of the ice tongue and transporting glacially modified water out of the fjord into the ambient ocean (blue in Figure 3a).The other plume-a dense bottom current-brings warm and salty AIW from the open ocean over the sill into the glacier cavity (red in Figure 3a).The strength of the overturning circulation is 39 mSv (Figure 3b), consistent with the value of (46 ± 11) mSv obtained from hydrographic measurements (Schaffer et al., 2020).
Subglacial melting creates a layer of cold water just below the ice along the whole glacier tongue (Figure 3c).This meltwater is transported away from the glacier and introduces a layer of cold water into the ambient ocean at depths of around 90-95 m below sea level.Minimum temperatures offshore the calving front are below −1.5°C at 94 m depth.Apart from this layer and its immediate surroundings, the temperature stratification offshore the sill is mostly in equilibrium with the imposed open ocean conditions.As the flow of AIW from the open ocean into the glacier cavity is hindered by the sill, the cavity water becomes colder than the open ocean water by mixing with meltwater (inset of Figure 3c).
Salinity differences are the main drivers of the circulation in the 79NG fjord (Figure 3d).On the one side, the subglacial plume rises along the ice tongue because it is fresher, thus lighter than the water inside the cavity.On the other side, AIW flows down the bottom slope into the glacier cavity because it is saltier, thus denser than the cavity water.Comparing the water at the same depth on both sides of the sill, we see that the cavity water, which is a mixture of AIW with meltwater, is at least 0.1 g kg −1 fresher than AIW (inset of Figure 3d).Offshore the sill, the salinity stratification is almost horizontally homogeneous and in equilibrium with the imposed conditions of the open ocean.
Along the whole ice tongue of 79NG, the basal melt rate is positive, that is, no freezing appears in our simulation (Figure 3e).We find the strongest melting of 58 m yr −1 close to the grounding line and a mostly monotonic decrease of the melt rate afterward.The melt rate reaches practically zero (<0.1 m yr −1 ) at around 42 km from the grounding line.The rest of the ice tongue has an average melt rate of less than 0.01 m yr −1 .The position where the melting stops is the place where the subglacial plume detaches from the ice tongue (see Section 3.1.1).The melt rate averaged over the whole ice tongue is 12.3 m yr −1 (corresponding to 20.3 km 3 yr −1 ) in our model, consistent with the value of (10.4 ± 3.1) m yr −1 , or (17.8 ± 5.2) km 3 yr −1 , estimated by Schaffer et al. ( 2020) based on measurements.Accordingly, also the percentage of subglacial discharge in the total meltwater production at 79NG is similar between our model (9.8%) and observations (11%; Schaffer et al., 2020).This shows that basal melting is by far the dominant freshwater source in the glacier fjord.

The Buoyant Subglacial Plume
The subglacial plume starts at the grounding line (x = 0), where subglacial runoff is discharged into the cavity.Since this discharge is fresher than the water in the fjord, it is positively buoyant and rises along the lower ice edge.We observe in our model that two opposing processes modify the plume water while rising.On its upper side, the plume causes melting of the ice tongue due to the turbulent heat flux, parameterized as a function of the friction velocity, which adds cold and fresh meltwater to the plume.On its lower side, ambient water is entrained upwards into the plume by turbulent mixing, thus making it saltier and warmer.This way, entrainment transports heat toward the ice and amplifies the melting (Burchard et al., 2022;Jenkins, 2011).As the plume rises, it passes through ever lighter surrounding water and reaches a point where its density equals that of the ambient water (Figure 4a).This is between 95 and 100 m below sea level.At this level, the subglacial plume detaches from the ice tongue, propagates horizontally away from the glacier, and transports glacially modified water out of the fjord (Figures 3a-3c).This observation is qualitatively consistent with the plume detachment and cold-water export at mid-depth in the model of an Antarctic ice shelf by Hellmer and Olbers (1989).
Before its detachment, the plume splits up a number of times.The first splitting occurs at 18 km from the grounding line (Figure 5).Until there, the plume was rising through well-mixed water, allowing it to grow and thicken rapidly by entrainment.However, around the depth of the sill (300 m), the ambient water changes from almost unstratified to stably stratified (Figure 5d).The lower part of the plume consisting of denser water that has been advected with the buoyant meltwater overshoots its neutral level.It falls about 70 m down, rises slightly again, and finds its neutral level near z = −290 m, where it propagates away from the ice (Figures 5a and 5b).This creates a buoyancy oscillation visible in the streamlines (Figure 3b).However, the oscillation is strongly damped, because the plume mixes with ambient water during its ascent and descent (Figure 5c), thereby reaching neutral buoyancy quickly (Figure 5d).Similar though smaller splits of the plume can be observed several times until the plume detachment.This creates a vertical velocity profile with a number of velocity peaks between the depth of the calving front and the depth of the sill (Figure 3a).A similar velocity profile has been observed in reality.Velocity measurements at the calving front of 79NG show the main outflow of glacially modified water near 100 m depth, in addition to weaker outflows at greater depths (Schaffer et al., 2020).These deeper outflows may be caused by the splitting of the subglacial plume.
Prior to the splitting of the plume, its thickness increases from D = 3 m at a distance of 5 km from the grounding line to about D = 5 m at x = 18 km (Figure 4a).Over this distance, the plume becomes more buoyant and increases its vertically-averaged velocity    to a maximum of 0.22 m s −1 (Figures 4b and 4c).When the plume splits, its velocity drops and so does its buoyancy  b , because the ambient water below the plume becomes lighter.
After the splitting, the plume thickens more slowly and reaches D = 6 m at x = 40 km, just before its detachment from the ice.When it detaches, the plume buoyancy drops again (Figure 4c), meaning that the plume density is similar to the ambient density, which is the reason for the plume detachment.Note that the buoyancy does not go to zero because the formulas to compute  b (Section 2.4) are only applicable while the plume is within 10 m from the ice edge; afterward, the thin lines in Figures 4b-4d represent the properties of the water just below the ice.
Entrainment at the plume base is only positive until the plume splits for the first time (Figure 6a).The plume thickening afterward is mainly due to flow convergence (Figure 6a) in consequence of the plume slowing down (Figure 4b).It is not due to entrainment, because after the initial phase, the entrainment velocity v e is negative and detrainment appears (Figure 6a).So instead of taking up ambient water, the plume in total loses water to the stratified interior of the cavity (Figure 3).Correspondingly, the vertical velocity under the plume is negative, that is, downward (Figure 6c).The detrained water forms an outflowing layer below the plume (Figure 6b).
Our interpretation of the detrainment is that initially, the weakly stratified water in the deep part of the cavity allows strong turbulence to develop (Figure 6d), leading to high entrainment rates of   = ∕ ū =  ( 2 × 10 −4 ) and rapid plume thickening (Figure 4a), consistent with the initial plume development and entrainment reported by Burchard et al. (2022).When the plume arrives in the more stratified upper part of the cavity, the reduced turbulence is insufficient to sustain the thick plume.Comparing turbulent kinetic energy (TKE) in the entrainment part with the detrainment part, we see that in the latter case, TKE is clearly reduced at the ice-ocean interface, at the plume base, and below the plume (Figure 6d).So the turbulence might be too weak to further entrain ambient water against gravity, and instead the plume detrains water.This manifests in the first plume splitting near x = 18 km and the subsequent smaller splits as described above.
The Froude number, Fr, of the subglacial plume (Equation 12) is very close to its approximation (Equation 13) based on the ice slope and the drag coefficient (Figure 4d).This indicates that the plume is dominated by friction at the ice-ocean interface (Arneborg et al., 2007), which is plausible, as the plume is a rather thin boundary layer.The decreasing Froude number can thus be considered a consequence of the decreasing ice slope.Since the Froude number decreases gradually (Figure 4d), there is no abrupt change in the flow and no hydraulic jump.From classical hydraulic theory, a hydraulic jump might be expected at the position where Fr = 1.However, the situation here seems to be more complicated, presumably because of the detrainment that the plume experiences.

The Dense Bottom Plume
The bottom plume in the 79NG fjord consists of AIW coming from the open ocean.With a density of 1,028.0kg m −3 (Figure 4e), this is the densest water mass in our system, as well as the warmest and saltiest (Figures 3c and 3d).It flows from the sill at x = 80 km down into the cavity, following the bathymetry.As long as the bottom slope increases, the plume accelerates up to a vertically-averaged velocity of    = −0.16m s −1 (Figure 4f).Due to this flow divergence, the plume thins from 17 m over the sill to 10 m thickness 6 km downstream (Figure 4e).The rapid plume thinning is associated with a transition from subcritical flow (Fr < 1) in the plume before it passes the sill to supercritical flow (Fr > 1) as the plume flows down the slope (Figure 4h).Just downstream of the sill, the Froude number becomes equal to one, which means that the sill acts as a hydraulic control for the bottom plume and limits the inflow of AIW into the cavity.This is consistent with hydrographic measurements around the sill at 79NG, which also indicated hydraulic control (Schaffer et al., 2020).While flowing down the bottom slope, the plume entrains ambient cavity water, which has a lower density since it contains meltwater (Figure 4e).In consequence, the plume density and buoyancy (in absolute value) decrease (Figure 4g).Similar to the subglacial plume, the bottom plume transports water below its neutral depth.The water  , c) and turbulent kinetic energy (d) in the 25 m under the ice (default scenario).The colored graphs in (a) represent the processes acting on the plume thickness: flow convergence (orange), subglacial melt rate (green, close to zero), and entrainment velocity at the plume base (purple); summed together, they give the thicker black line, see Equation 10.For the calculation of the graphs in (a), plume thickness D and bulk velocity    were smoothed with a running average of window size ±1 km.
then rises again and adjusts in an oscillating way to its level of neutral buoyancy (Figure 3b), before propagating horizontally away from the bathymetry.This way, the bottom plume fills the cavity with (partially mixed) AIW over a depth range of 450-600 m (Figure 3a).At about 600 m below sea level, the plume has detached completely from the bottom.It cannot propagate further down, because the entrainment of cavity water made the plume lighter than the water in the trough below 600 m depth.The water in the deep trough is dense because it consists of almost pure AIW (from the initialization) with only little meltwater.This is because (a) meltwater enters the cavity only at depths where the ice tongue is present, and (b) the meltwater is not mixed far below the grounding line (600 m) due to the absence of strong motion there.
Outside the cavity, just offshore the sill, even some AIW below the sill level moves upward and flows over the sill (Figure 3b).This overflow is driven by an internal pressure gradient that is vertically homogeneous, since the water on the upstream side of the sill is unstratified.The phenomenon of upward acceleration of dense water against gravity is called aspiration and commonly observed in fjords (Inall & Gillibrand, 2010).

Performance of the Adaptive Vertical Coordinates (AVC)
AVC is one feature of our model that has not been employed before in simulations of glacier fjords.Our setup uses 100 vertical layers that adapt automatically to the stratification, as explained in Section 2.3.This way, we reach high vertical resolutions in both plumes.
The vertical resolution in the subglacial plume is everywhere close to 1 m and even better in the entrainment layer at the plume base (white lines in Figure 4a).Thus, AVC achieve the necessary resolution to represent the entrainment into and detrainment out of the plume correctly (Burchard et al., 2022).Since the model layers adapt to and follow the plume, its water is advected mostly along the layers and not across.The plume is always resolved by five layers or more while it is attached to the ice, which allows preserving the plume properties well.Models with z-coordinates usually do not achieve this, which causes the plume to spread out.For example, the layer of cold water under the ice is around 50 m thick in the 2D model of Hellmer and Olbers (1989), much more than the 5 m-thin plume in our setup.This shows an important advantage of stratification-zooming coordinates.
When the plume splits (Figure 5) and when it detaches from the ice (Figure 4a), AVC also attempt to follow the flow of the meltwater by partially bending in the horizontal direction, but cannot follow the plume as well as when it is at the ice.In consequence, the plume must pass through layers that are not fully aligned with its flow direction, increasing the numerical diffusion.The calving front presents another challenge for AVC.As terrain-following coordinates, they must connect the lower ice edge to the sea level, a difference of 75 m in depth.However, the flow under the calving front is horizontal and the density is horizontally homogeneous, so there is necessarily a divergence between coordinates and plume.By stretching the calving front over 7.5 km as explained in Section 2.2.1, the vertical position of the ice-ocean interface changes gradually enough, so that the coordinates manage to adapt to the plume to some extent and preserve its properties well (see the inset of Figure 3c).However, a slight dilution of the plume as it passes under the calving front and through several layers can still be seen (Figures 3a-3c).
Similar to the subglacial plume, also the incoming plume of Atlantic water is resolved by several layers with a thickness on the order of 1 m (Figure 4e).
The high resolution in the vicinity of the ice and the bottom comes at the expense of thicker layers in the interior of the glacier cavity.While the vertical layers are less than 10 m thick in most areas, there are up to 15 m-thick layers in the middle of the water column in places where the fjord is deepest.However, we believe that this is a good trade-off, because (a) the thick layers appear in areas where the velocities are small and the water column is only weakly stratified, and (b) we obtain very thin layers in the dynamically relevant parts.

Sensitivity Studies
We now explore how the results change compared to the default scenario for modified environmental influences.Key properties of all presented scenarios are summarized in Table 2.

Influence of the Ambient Ocean Salinity
The subglacial plume detaches from the ice tongue and transports meltwater out of the fjord toward the open ocean at a depth of around 95 m below sea level in our default scenario.This sensitivity study shows that the depth depends strongly on the salinity stratification of the ambient ocean, which is imposed at the open boundary of the model.When the salinity of the upper water column is increased, the plume propagates further along the ice tongue and detaches higher up.With lower salinities above the sill, the plume does not propagate as far up and detaches earlier.
This relation is exemplified by the two sensitivity experiments shown in Figure 7 in comparison with the default case.For the high salinity scenario, we increased the surface salinity from 29 to 33.5 g kg −1 , so that we obtain a linear salinity stratification in the upper 300 m of the water column (Figure 7b).With this stratification, the plume detaches at around 50 m below sea level (Figure 7a).In the low salinity case, we kept the surface value at 29 g kg −1 but decreased the salinity at 100 m depth from 34 to 33.5 g kg −1 (Figure 7f).Then most of the plume detaches between 125 and 150 m of depth (Figure 7e).These experiments also show that the plume detachment is not caused by the abruptly changing stratification that is in the default scenario at a similar depth as the detachment (Figures 7c and 7d).
In fact, it is the salinity of the open ocean that determines the depth where the plume detaches.The salinity at the detachment level is (33.7 ± 0.1) g kg −1 in all three scenarios.We also tested a stratification with a minimum salinity of 34 g kg −1 (not shown), in which case the plume never detaches from the ice tongue but reaches sea level.The reason that the detachment depth depends strongly on salinity is that at this level, the plume density equals that of the ambient ocean, which is set primarily by salinity in the 79NG fjord.
For the deeper half of the ice tongue, the plume developments and melt rates are basically identical between our sensitivity experiments, but they differ in the upper 300 m.At the plume detachments, the subglacial melt rates drop to almost zero, which shows again that the subglacial plume is responsible for the bulk of basal melting.In the scenario with the plume detachment at great depths, a small second plume develops above the main detachment, causing some more melting with melt rates up to 0.7 m yr −1 before detaching near 100 m depth (Figure 7e).Only in the scenario with a late plume detachment, we observe melt rates above 0.2 m yr −1 along the whole ice tongue up to the calving front.However, note that the plume development as it propagates up the calving front in this scenario (Figure 7a) is not entirely realistic, because the calving front is sloping in our model and not vertical (Section 2.2.1).

Influence of the Ambient Ocean Temperature
We investigate the influence of the imposed temperature stratification at the open ocean boundary by varying the temperatures of Polar Water (PW) and AIW individually as well as together.In our model, PW occupies the upper 100 m of the water column and has in the default scenario a linear temperature profile with −1.5°C at sea level and −1.0°C at 100 m depth (Figure 2b).AIW fills the water column below 300 m depth and has a vertically homogeneous temperature of 1.5°C by default.In between 100 and 300 m, we apply a linear temperature gradient.In our sensitivity study, we increase the temperatures of AIW and/or PW by 0.5 or 1.0 K.We also decrease AIW temperatures by 0.5 and 1.0 K.Note that we cannot make PW colder, because the surface temperature is just above freezing in our default scenario (Figure 2b).
We observe that the AIW temperature has a clearly larger impact on the glacier cavity than variations of PW temperature.With increasing AIW temperature, the subglacial melt rate increases along the whole ice tongue (Table 2) and the point at which the plume detaches moves upward.For AIW temperatures of 0.5°C, the plume detaches below 130 m, for 2.5°C above 90 m depth.This can be explained by the increased temperature forcing, which causes more melting and thereby a lighter plume that rises faster and further.Interestingly, in the deep part of the cavity, the thickness of the subglacial plume is not much altered by temperature differences, although this is the part where AIW is present.Note.Melt rate is the subglacial melt rate averaged over the whole ice tongue.
Overturning is the strength of the circulation measured above the sill.Runoff is the percentage of subglacial discharge in the total meltwater outflow (discharge plus melting) of the fjord.Observation cites the values reported by Schaffer et al. (2020).AIW stands for (the temperature of) Atlantic Intermediate Water, PW for Polar Water.

Table 2 Summary of the Presented Simulations
Our findings are qualitatively consistent with modeling studies of the circulation under Antarctic ice shelves.Hellmer and Olbers (1989) reported a plume detachment at greater depth and a reduced overturning circulation when the inflowing bottom water has a lower temperature, and the opposite effect for a higher temperature.Even though they also modified the inflowing salinity in addition to the temperature, they claimed that the observed effects are actually due to the temperature variation, which is confirmed by our results.Grosfeld and Gerdes (1998) observed that increased temperatures of the water flowing into the cavity led to strongly increased melting, which reduces the salinity of the outflow.This fits with our observations of a lighter plume that detaches later from the ice tongue, at a depth where the salinity is lower.
The parametrization by Slater and Straneo (2022) captures the temperature dependence of the melt rate well, but only in the vicinity of the grounding line.Let us first consider the 15 km of the ice tongue directly after the grounding line, which is the part where the plume rises through a water mass that is similar to AIW (Figures 3c  and 3d).In this area, the average melt rate computed in our simulations is best described by the function (8 ± 1) (Δθ) 1.24±0.09, where Δθ is the temperature forcing, that is, the difference between AIW temperature and the freezing point at the grounding line.The values after ± are 95%-confidence intervals, so our fit is consistent with the (Δθ) 1.19 -proportionality used by Slater and Straneo (2022), though with a larger constant of proportionality.However, if we average over the full length of the ice tongue, the melt rate can be parameterized as (1.3 ± 0.2) (Δθ) 1.69±0.09 .The exponent is significantly larger, but smaller than in the (Δθ) 2 -law found by P. R. Holland than in our default scenario (c, d) as well as lower salinity (e, f).For higher salinities above the sill, the subglacial plume propagates further along the ice tongue and detaches higher up.The salinity at the level of plume detachment is always around 33.7 g kg −1 .When the plume detaches early (e), a weaker secondary plume develops above.et al. (2008).This shows that a close-to linear relation between melting and thermal forcing is only applicable near the grounding line (Slater & Straneo, 2022) and should not be applied to the whole ice tongue.A linear relation between melt rate and Δθ was suggested by Jenkins (2011) and Xu et al. (2012), which fits with our modeled melt rates up to about 4 km from the grounding line.
The effects associated with increased PW temperatures are much smaller.Cavity circulation and both plumes look practically the same as in the default scenario.The only (small) difference we observe is in the detachment point of the subglacial plume.It moves about 2 m down for a PW temperature increase of 0.5 K and about 3 m (compared to default) for a 1.0 K-increase.This makes sense because the upper part of the water column is lighter for warmer PW, so the plume reaches its neutral buoyancy earlier.Since subglacial melting almost stops when the plume detaches, the overall melt rate is slightly lower for higher PW temperatures (Table 2).However, note that our model does not simulate calving, which can be intensified in warmer water.
When we increase the temperatures of both AIW and PW together, thus making the whole water column warmer, we observe a combination of the effects described above.The results look similar to those with only increased AIW temperatures, but the subglacial plume detaches at a slightly deeper level.

Role of the Subglacial Discharge
The meltwater discharged at the grounding line has an important influence on subglacial melting.In our default scenario, we prescribe a constant subglacial discharge of 70 m 3 s −1 , which is the value reported by a field campaign (Schaffer et al., 2020), and we find a clear, peaked melt rate maximum just after the grounding line.In contrast, if we reduce the discharged water volume in our model by an order of magnitude to 7 m 3 s −1 , we observe a flatter melt distribution after the grounding line with a lower and rather constant melt rate over the first 10 km (Figure 3e).Interestingly, after the splitting of the subglacial plume, the melt distributions look similar for low discharge and normal discharge (Figure 3e).Also, the position of the plume detachment from the ice tongue is not much different.These observations suggest that the subglacial discharge has mostly an impact on the early development of the plume (consistent with Jenkins, 2011), while further away from the grounding line, the plume development is mostly determined by subglacial melting and the ambient ocean stratification.
Due to the decreased subglacial melting in scenarios with lower subglacial discharge, the cavity water is warmer, saltier, and denser.This has the effect that the dense bottom plume does not propagate as far down the slope and detaches earlier from the bottom.Also, both plumes are thinner and slower than in the default scenario.The strength of the overturning circulation is reduced by about one third to 26 mSv for a discharge of 7 m 3 s −1 (Table 2).
We observe the opposite effects when we increase the subglacial discharge: The melt rate increases; the cavity water becomes colder, fresher, and lighter; the plumes are thicker and faster.Doubling the discharge to 140 m 3 s −1 increases the overturning strength by about one fourth (relative to default scenario) to 50 mSv and the average melt rate by about one sixth to 14.3 m yr −1 (Table 2).
Only if we restrict the averaging 〈⋅〉 to the first 1 km of the ice tongue, we find that the melt rate is proportional to Q 0.27±0.04(not shown).This fit would include, in its 95%-confidence interval, the Q 0.31 -proportionality used by Slater and Straneo (2022) and is also close to the Q 1/3 -law suggested by Jenkins (2011) and Xu et al. (2012).But the average contains just two grid cells and is thus not representative of the whole ice tongue.
Considering the melt rate averaged over the full ice length, the best fit to our model data is 〈v b 〉 = (6.4± 0.6) Q 0.16±0.02, suggesting that the melt rate is roughly proportional to Q 1/6 (Figure 8).However, such a law would imply a zero melt rate for zero discharge, which is not plausible as plumes develop also under ice shelves without subglacial discharge.For example, at 79NG, mooring data indicate a year-round outflowing plume, even in the months without subglacial discharge (Lindeman et al., 2020;Schaffer et al., 2020).To allow for plume-induced melting in the absence of subglacial discharge, we fit the function 〈v b 〉 = cQ e + d to seven model runs with different values of Q.A least-squares regression gives the coefficient c = 0.8 ± 0.3, the offset d = 7.2 ± 0.5, and the exponent e = 0.44 ± 0.06 (values after ± are 95%-confidence intervals).The fit clearly describes our model results better than the Q 1/3 -law and also slightly better than the Q 1/6 -law (Figure 8).This means that in the absence of subglacial discharge, the average melt rate is about 7 m yr −1 , and increases approximately with the square root of the subglacial discharge,  √  .
In our model, we cannot reasonably increase the subglacial discharge arbitrarily.For example, with a discharge of 700 m 3 s −1 (ten-times the default), the large amount of meltwater leaving the cavity cannot be transported across the open boundary, because the prescribed conditions at the open boundary correspond to the default scenario, which has lower discharge and melting.This causes a density front near the open boundary, which is physically unstable and prevents the system from reaching a steady state.Nevertheless, the model stays numerically stable, even in such a non-equilibrium situation.

Role of the Sill
Our model allows us to test a hypothesis made by Schaffer et al. (2020) based on their hydrographic measurements.They claim that the bathymetry of the 79NG fjord constrains the heat transport from the open Atlantic Ocean into the glacier cavity.According to Schaffer et al. (2020), the height of the sill at the fjord entrance determines how much warm AIW flows into the fjord, and in turn how much heat is available for subglacial melting.
In our idealized 2D model, we can easily modify the sill height (default: 300 m below sea level) or remove the sill completely and check which impact it has.
We find that the cavity water is clearly colder with a higher sill than with a lower sill or without a sill (Figures 9a  and 9b, see also Figure 3c).The higher the sill, the stronger the temperature contrast between the water in the cavity and the water on the continental shelf.Consequently, the melt rate is larger if the sill is at greater depths and vice versa (Figures 9c and 9d).Interestingly, the melt rate is not larger over the full length of the ice tongue, but mostly in the (20 ± 5) km after the grounding line, where the ice is at great depths.The melting of the thinner part of the ice tongue is not much influenced by the sill, neither is the position of the plume detachment from the ice.When the sill is at 350 m below sea level or deeper, the melt rate is almost independent of the sill depth (Figures 9c and 9d).At this depth, the sill cannot effectively prevent the warm AIW from entering the cavity anymore.
So our simulations show that indeed the sill height constrains the heat transport into the cavity and thereby determines the melt rate of the 79NG ice tongue.This "sill effect" almost ends at a depth of about 350 m, measured from sea level.The Q 1/3 -law proposed by Jenkins (2011) does not give a good fit (dotted, blue).The parametrization proportional to Q 0.31 used by Slater and Straneo (2022) is very similar to the Q 1/3 -law (not shown).A better fit is obtained by using an exponent of 0.16 ± 0.02, which is roughly a Q 1/6 -law (dashed, orange).The best fit is obtained by allowing a non-zero melt rate for zero subglacial discharge, which is close to a  √  law that is shifted upward by about 7 m yr −1 (solid, green).

Roughness Length
In our setup, the smoothness or roughness of the ice tongue on its underside is modeled by a roughness length, z 0,ice (Section 2.2.2).This parameter has the value 0.01 m in our default scenario, but it is poorly known which value is realistic for a given ice shelf (P.R. Holland & Feltham, 2006;Jourdain et al., 2017).To test the sensitivity of the 79NG system on this value, we increased the roughness length by a factor of ten (z 0,ice = 0.1 m, scenario z0p) and decreased it by a factor of ten (z 0,ice = 0.001 m, scenario z0m).We also tested intermediate values to ensure that our observations are actually tendencies as reported below.
Our model results show that the shorter the roughness length, the larger the melt rate and the stronger the overturning circulation (Table 2).Due to the higher melting, the subglacial plume becomes colder (Figures 10c  and 10d), fresher, and more buoyant.It accelerates faster and has a higher velocity under the ice and after its detachment (Figures 10a and 10b).Also the inflowing bottom plume is faster with a shorter ice roughness length (not shown), contributing to the increased overturning strength (Table 2).In the scenario with rougher ice (z0p), most of the plume detaches from the ice tongue already at a depth of 200 m and leaves the fjord at this level, while the outflow at 100 m-depth is much weaker (Figures 10b and 10d).Initially, the plume thickens quickly by entrainment, but detrains strongly from 16 to 24 km behind the grounding line, leading to an almost complete detachment of the plume at around 24 km from the grounding line (Figure 10b).Behind this point, the plume is thinner than in the other two scenarios.In the scenario with smoother ice (z0m), the plume is everywhere thinner than in the default scenario.
The reasons for the observed effects of varying ice roughness are complex.The turbulent heat transfer velocity γ T in the melt formulation depends on the ice roughness length z 0,ice (Burchard et al., 2022).Lower roughness leads to higher heat transfer, which can explain the higher melt rate.In consequence of the increased melting, the subglacial plume becomes more buoyant, which means that the density difference between plume and ambient water is larger.The stronger stratification at the plume interface hinders entrainment, and this lower entrainment in turn leads to a stronger stratification, possibly indicating a positive-feedback loop.This loop eventually breaks as the plume rises along the ice, because the lighter plume accelerates faster and reaches higher velocities, which then increase turbulent mixing and entrainment.The divergence caused by the strong initial acceleration together with the initially weak entrainment explain why the plume under smoother ice is thinner.
Consistent with our findings, also Jenkins (1991) reported higher melting in experiments with a lower drag coefficient, which corresponds to a shorter roughness length.The model used by Jenkins (1991) is a 1D plume Figure 9. Temperature in the glacier cavity in a modified 79° North Glacier (79NG) fjord with a high sill (a) and with no sill (b), as well as subglacial melt rate of the 79NG ice tongue with x-resolution (c) and in spatial average (d) for different sill depths (including no sill).When the sill is higher, that is, with a lower sill depth, less warm water can flow into the cavity, so the melt rate is lower.Note that the continental shelf offshore the cavity is at 450 m below sea level, so a sill depth of 450 m means no sill.
model that-like our model-does not take into account the impact of Earth rotation.Models that include Coriolis show additional effects in response to increased ice roughness.In a realistic setup of the Amundsen Sea, Jourdain et al. (2017) varied the drag and the heat transfer independently of each other.They found that melting increases with the heat exchange coefficient Γ T , with the drag coefficient c d , and also with the Stanton number  St = √ Γ (Jourdain et al., 2017).This does not conflict with our results, because heat transfer and drag are not independent in our model.For smoother ice, the drag coefficient is lower, but the heat exchange is higher, in a way that the Stanton number is generally larger.So, the higher heat transfer compensates for the lower drag and results in stronger melting at smoother ice.However, Earth rotation also deflects the plume away from the direction of the steepest ascent.This effect is stronger, the lower the friction at the ice-ocean interface (P.R. Holland & Feltham, 2006).A plume under smooth ice may be deflected until it is at the side wall of the fjord, where wall drag slows down the plume, which then leads to lower melting (P.R. Holland & Feltham, 2006).This makes the consequences of varying ice roughness more complex (Payne et al., 2007).However, the model by P. R. Holland and Feltham (2006), which was also employed by Payne et al. (2007), cannot represent the plume detachment.The same is true for the water column model of a subglacial plume by Burchard et al. (2022), which they applied to a setting similar to 79NG.They found that in the initial phase of the plume development, the melt rate is higher for smoother ice, while in a later phase, the relation is reversed.The transition from the initial to the later phase occurs after about 1 week (Burchard et al., 2022).Our model shows that the plume detaches from the ice always within 1 week.In a Lagrangian sense, the plume needs less than 4 days to reach the point at which it detaches from the ice, in all analyzed scenarios (Figure 10e).Thus, the plume at 79NG goes only through the initial phase.Models without plume detachment might also simulate the later phase, which does not always occur, as shown by our results for 79NG.This can lead to different conclusions regarding the relation between drag at the ice-ocean interface and subglacial melting.

Discussion
In large-scale ocean models without explicitly resolved glacier cavities, meltwater from fjords is often introduced at the sea surface (e.g., Stolzenberger et al., 2022).Our model results show that this is generally not realistic for fjords with an ice tongue.This matches with a similar observation from a high-resolution model of a fjord with a vertical glacier front (Xu et al., 2013).In our default scenario, the bulk of meltwater leaves the 79NG fjord between 90 and 100 m below sea level (Section 3.1 and Figures 3a-3c).This level depends primarily on the stratification of the ambient ocean, which is mainly set by salinity.Even a relatively small change in the upper ocean salinity can alter the outflow depth of glacially modified water by 50 m (Section 3.3.1 and Figure 7).The temperature stratification also influences the outflow depth, but less dramatically, as our sensitivity study shows (Section 3.3.2).On the other hand, the outflow depth is almost unaffected by the subglacial discharge and by the sill at the fjord entrance, despite their big influence on subglacial melting and overturning circulation in the cavity (Sections 3.3.3 and 3.3.4).If the base of the ice tongue had a higher roughness, the outflow around 95 m depth would be weaker but still at the same depth as for smooth ice (Section 3.3.5 and Figures 10a and 10b).We suspect that the outflow depth of meltwater does not change much with seasons, because the fjord properties that have a strong seasonality are the subglacial runoff (Lindeman et al., 2020;Schaffer et al., 2020) and the ocean surface temperature, which both have little impact on the outflow level.Whether the sub-surface stratification at 79NG, which is important for the outflow depth, shows seasonal variability, is still unknown, but the existing mooring data shows no clear signature of a seasonal cycle (Lindeman et al., 2020;Schaffer et al., 2020, and own analysis of their data sets).Longer time series of measurements at 79NG are necessary to answer this question.
Our analysis of entrainment rates (Figure 6) reveals that the dynamics of the subglacial plume at 79NG fall into two different regimes.The first one is analog to the so-called plume regime (Baines, 2008).This is the case over the first 17 km from the grounding line, where the ice slope is high and the ambient stratification is low.The plume shows strong entrainment, partly overshoots its density horizon, falls down again, and intrudes over a wide range of depths (Figure 5).However, at 17 km from the grounding line, the entrainment rate switches sign (Figure 6a).The boundary layer leaves the plume regime and enters that of a gravity current.The gravity current, which exists under a more gently sloping ice, is characterized by detrainment, that is, it loses water to the stratified interior of the cavity (Baines, 2008).When the gravity current finds its neutral density level after around 40 km from the grounding line, it detaches from the ice tongue without overshooting (Figure 3).Both behaviors of the turbulent boundary layer at the ice-ocean interface fit the descriptions by Baines (2008) for dense downslope flow, except that for buoyant upslope flow, everything is upside-down.
Current models of subglacial plumes often employ an assumption of continuous entrainment into the plume (Hewitt, 2020;Lazeroms et al., 2018), a process that has so far not been well constrained by measurements (Anselin et al., 2023).Our results put the validity of this assumption into question.In fact, the subglacial plume in our idealized 79NG fjord model shows entrainment only for about half of its way along the ice tongue, but detrainment afterward.Detrainment is generally not included in current models of meltwater plumes.We thus echo the statement by Hewitt (2020) that these models might not capture all important dynamics and should be revised.
The depth at which meltwater leaves the glacier fjord is not only relevant for the export of glacially modified water but also for the development of the ice tongue.Our simulations show that most subglacial melting occurs while the subglacial plume is at the ice-ocean interface.When the plume detaches, the melt rate drops to almost zero.This happens roughly at the same level as the meltwater outflow.Thus, oceanographic measurements of the depth of glacially modified water near a glacier fjord can be used to infer which part of the glacier tongue is likely to show high basal melt rates.This information can be helpful for a decision of where to install measurement stations on a floating ice tongue to monitor ice thickness changes.
At the depth where the subglacial plume propagates away from the ice tongue, the vertical coordinate levels in our model accumulate.This ensures that the water properties of the plume are preserved over long distances with little spurious mixing.It is achieved automatically by the stratification-zooming of AVC.No a priori knowledge of the position of plume detachment is needed, which is an important difference to non-adaptive coordinates that can achieve high vertical resolutions in pre-defined regions.Moreover, AVC change the vertical layer distribution with time, for example, in simulations with tides or other time-varying forcings that alter the stratification.
With z-coordinates, which are often used to model the ocean under an ice tongue or an ice shelf (e.g., Hellmer & Olbers, 1989;Losch, 2008), it would be difficult to obtain equally detailed simulations of the cavity circulation and in particular of the subglacial plume.Due to their step-wise manner of resolving the ice-ocean interface, z-coordinates are usually too diffusive to preserve the plume over longer distances.Without a well-preserved plume, an analysis of the entrainment rate as shown in Figure 6 would not be feasible.An insufficient representation of the plume development has also implications on the accuracy of the computation of basal melt rates (Burchard et al., 2022).Furthermore, a good simulation of meltwater export from the fjord into the open ocean demands a good preservation of the plume properties with minimal spurious mixing.This can be provided by AVC while the plume is under the ice.Further development of the adaptive coordinates should try to improve also the representation of the outflow after it detached from the ice and as it passes under the calving front.
While AVC (Hofmeister et al., 2010) have a number of characteristics, the main feature used in our setup is their capability to zoom toward stratification.This enables high resolutions in the entrainment layers of both plumes and allows the coordinates to follow the outflow to a reasonable extent, so that glacially modified water can be transported far offshore.This stratification-zooming could be combined with other modeling approaches like vertical Lagrangian remapping or the Arbitrary Lagrangian-Eulerian (ALE) method.In these methods, Lagrangian motion of the model grid is followed by a regrid step, in which the coordinate surfaces are moved back to prescribed target positions; the physical fields are then mapped onto this new grid in a remap step (Griffies et al., 2020).The target coordinate layout could be prescribed based on the ocean stratification in the current model state.Such an approach would combine the advantages of ALE with the advantages of stratification-zooming shown in this paper.
As for terrain-following coordinates in general, the calving front presents a challenge for AVC.Our setup uses a gentle slope instead of an almost vertical wall at the ice front to make sure that the plume is well preserved as it leaves the cavity.This part of the ice tongue could possibly be simulated more realistically by a modification of the cost function that determines the zooming of AVC.Instead of zooming to stratification and the sea surface, it might be advantageous to zoom only to stratification and the ice-ocean interface but not to the atmosphere-ocean interface.This way, more layers could be available at the calving front to allow a high calving front slope as well as a good preservation of plume properties.Since AVC (Hofmeister et al., 2010) have not been developed with glacier tongues in mind, and this paper presents their first application to an ice cavity, such a possibility has not yet been implemented.It should however be kept in mind that processes at the calving front are strongly nonhydrostatic in nature and therefore cannot be sufficiently reproduced with classical ocean models anyway.
While our idealized 79NG fjord model shows qualitatively realistic dynamics and processes under the glacier tongue, its quantitative results should be taken with a grain of salt, as exemplified by our sensitivity study on the sill depth (Section 3.3.4).We observe that the melt rate of the ice tongue (Figure 9) and the strength of the overturning circulation in the cavity (Table 2) are very sensitive to the depth of the sill at the fjord entrance, which is 300 m in our default setup.However, no single value can be entirely realistic, because in the real system, the sill is not at the same depth over the whole fjord width (Figure 1a).The depth of the sill, which is the shallowest point that inflowing water must cross, depends on the path from the open ocean into the cavity.It can be as deep as 325 m below sea level but also shallower (see Figure 1 and Schaffer et al., 2020).Since this cross-fjord variability cannot be reproduced in 2D, the quantitative results of a 2D model can only be approximations.
Another effect that is neglected in the 2D approach is Earth rotation.The internal Rossby radius in the 79NG fjord was estimated to be less than 5 km (Lindeman et al., 2020), so at least four-times smaller than the fjord width (Figure 1a).This suggests that the plumes are deflected to the right by the Coriolis effect.We expect the inflowing plume to follow the northern boundary of the fjord, while the outflowing plume will be rather along the southern wall.Indeed, satellite measurements show higher subglacial melt rates along the southern boundary (Wilson et al., 2017), which can be caused by a more intense subglacial plume in the South.Thus, circulation and melting in the 79NG fjord seem to vary in the transverse direction.However, regarding the sill-controlled inflow of AIW, the situation could be different.The sill is located in a narrow strait of circa 2 km width, so the inflowing plume is thinner than the internal Rossby radius (Schaffer et al., 2020).It is thus not a priori clear, whether rotation plays a dominant role for the sill-controlled inflow.To answer this question, an extension of our setup to a 3D model is necessary.

Conclusions and Outlook
We developed a numerical ocean model of a glacier fjord in 2D with high horizontal and vertical resolution.The fjord and its forcing were built to resemble 79NG in an idealized, analytical way (Figures 1 and 2).Quantitative results of our default simulation are a good approximation of reality.In particular, the subglacial melt rate and the strength of the overturning circulation are consistent between our model and measurements at the glacier (Table 2).Thanks to the simplicity of the model, its qualitative results (Figure 3), which we explored further in a sensitivity study, will also hold for other glacier cavities.
Our model shows that the buoyant plume, which develops on the underside of the ice tongue, is responsible for the bulk of subglacial melting.When the plume reaches neutral buoyancy and detaches from the ice, basal melting almost stops.At this level, which is about 95 m below sea level in our present-day (default) scenario, the plume transports meltwater out of the fjord toward the open ocean.The detachment depth is set primarily by the stratification of the ambient ocean, particularly its salinity (Figure 7).In between the detachment depth and the sill depth, there are weaker outflows out of the cavity caused by splitting of the subglacial plume (Figure 5).The plume splits at around 18 km from the grounding line, because the turbulence in the plume is too weak to further entrain ambient water, so detrainment occurs (Figure 6).Furthermore, we confirmed that the depth of the sill at the fjord entrance has a big influence on the melt rate and the overturning strength in the fjord.With a deeper sill, the dense bottom plume brings more warm Atlantic water into the cavity and thus more heat is transported toward the ice tongue (Schaffer et al., 2020), which intensifies subglacial melting.In case of 79NG, this sill effect ends at around 350 m depth (Figure 9).The two plumes that make up the estuarine circulation in the glacier cavity are resolved by our model in great detail (Figures 4 and 6), thanks to the stratification-zooming of AVC (Hofmeister et al., 2010).We showed for the first time that with this modeling approach, a vertical resolution of less than 1 m in the entrainment layer of the buoyant plume under an ice tongue can be achieved (Figure 4), which is important for the correct representa tion of subglacial melting and plume development (Burchard et al., 2022).The computational cost compared to non-adaptive σ-coordinates is increased by less than 10% (Section 2.3), which is much cheaper than increasing the number of vertical layers.Further advantages of AVC are that they minimize the pressure gradient error (Gräwe et al., 2015;Hofmeister et al., 2010) and that they follow the plumes to some extent, which preserves the properties of the outflowing water mass quite well (Figure 3).We believe that the application of AVC in more ocean models will mean an improvement to the way processes under ice tongues and ice shelves are simulated.When stratification-zooming is used together with a melt parametrization that is suitable for high vertical resolutions (Burchard et al., 2022), this can refine projections of ice sheet melting and glacier stability.
Given the successful demonstration of AVC in an idealized 2D glacier cavity, a next step should be to extend this setup into a realistic 3D model of the 79NG fjord.This should include resolving the across-fjord dimension with the same high resolution as the along-fjord direction, using the real geometry and topography of the fjord, as well as forcing the regional ocean model with actual observational or reanalysis data.Such a setup will allow to study effects that have been neglected so far, for example, the Coriolis effect, and will back up our qualitative results with accurate quantitative assessments.

Appendix A: Analytical Description of the Setup
Our setup is built to resemble the 79NG fjord in an idealized way that can be completely described by simple, analytical functions.Here we give the mathematical expressions of these functions for the future use of our setup as a reference test case.

Figure 1 .
Figure 1.(a) Map of the 79° North Glacier fjord and its surroundings showing the bottom elevation from the RTopo-2.0.4 data set (Schaffer et al., 2019, resolution 30′′ = 1/120°) together with the positions of seismic depth soundings by Mayer et al. (2000) (measured in 1997/1998 and published by Mayer et al., 2018).The floating ice-tongue extends from the grounding line in the Southwest to the northern calving front in the Dijmphna Sund and to the main calving front in the East.Atlantic water must pass over a 325 m-deep sill (labeled deepest sill) to flow from the open ocean into the cavity.The reference section is a path from the grounding line toward the open ocean that follows the depth soundings up to the calving front and passes over the deepest sill.(b) Bathymetries and ice topographies along the reference section (from RTopo), along the section by Mayer et al. (2000), and in our idealized 2D fjord model.The position where subglacial discharge enters the cavity is marked with a wedge.Note that the deepest sill is the shallowest point along the reference section.The sill depth in our default setup (b) is 300 m, shown as a thin dashed contour in (a).

Figure 2 .
Figure 2. Stratification used in our model as boundary and initial conditions compared with salinity (a) and temperature (b) measurements near the 79° North Glacier (79NG) fjord.The shaded area marks the minimum and maximum values tested in our sensitivity study.The CTD profile was taken in 2017 on RV Polarstern(Kanzow et al., 2018) and represents a typical ambient ocean stratification for 79NG(Schaffer et al., 2020, see their Figure1afor the location of the profile).The freezing point of saline water in (b) corresponds to the shown CTD profile.We used the Python package gsw(TEOS-10; IOC  et al., 2010)  to convert from the CTD data pressure to depth, practical salinity to Absolute Salinity, and potential temperature to Conservative Temperature, as well as to compute the freezing temperature.

Figure 3 .
Figure 3. Model results in steady state for our default scenario of the 79° North Glacier fjord showing horizontal velocity (a), stream function (b), temperature (c), salinity (d), and melt rate (e).Insets in panels (c) and (d) show vertical profiles of temperature and salinity, respectively, at positions on both sides of the sill marked with vertical lines in the same colors as the graphs; conditions at the open boundary are shown in black for comparison.The thin orange line in panel (e) corresponds to a sensitivity experiment, in which the subglacial discharge is reduced by an order of magnitude compared to the default scenario (Section 3.3.3).

Figure 4 .
Figure 4. Details of the buoyant plume (a-d) and of the dense plume (e-h) in the default scenario.Panels (a) and (e) show density in the 15 m just below the ice tongue and in the 15 m just above the sea floor, respectively.The shown areas are marked in white in the overview panel (o).Note the different starting points of the colorbars.White lines in (a) and (e) are coordinate levels (upper/lower cell edges of the model grid) and emphasize the high vertical resolution of about 1 m obtained by adaptive vertical coordinates in the entrainment layers of both plumes.The red lines represent the thicknesses D of the plumes before their detachments, which are marked by dotted vertical lines.After its detachment, the water of the subglacial plume in (a) flows horizontally in parallel to the yellow z = −100 m isobath.Note that the bulk values in panels (f) and (g) have opposite signs than those in (b) and (c), because the plumes go in opposite directions and are oppositely buoyant.After the plume detachments, bulk values are shown as thin lines, because they do not represent the plumes anymore, but are averages of the uppermost 10 m under the ice (b-d) or the lowermost 30 m above the sea floor (f-h), see Section 2.4 for details.Panels (d) and (h) show the Froude numbers (solid) computed by Equation 12 in addition to their approximations (dashed) computed by Equation 13.

Figure 5 .
Figure 5. Zoom to the first splitting of the subglacial plume (default scenario).The arrows represent the flow direction resulting from the combined effects of the horizontal (a) and vertical (b) velocity components: (1) The plume rises along the ice tongue; (2) the lower part of the plume falls down from about 250 m depth to about 320 m depth, while becoming slightly colder and lighter due to mixing with ambient water; (3) the plume rises to about 290 m depth and (4) flows horizontally away from the ice tongue.

Figure 6 .
Figure6.Development of the subglacial plume thickness D before the detachment from the ice tongue, with areas of entrainment and detrainment highlighted (a); vertical profiles at positions near maximum entrainment (red) and detrainment (blue), showing velocity components (b, c) and turbulent kinetic energy (d) in the 25 m under the ice (default scenario).The colored graphs in (a) represent the processes acting on the plume thickness: flow convergence (orange), subglacial melt rate (green, close to zero), and entrainment velocity at the plume base (purple); summed together, they give the thicker black line, see Equation10.For the calculation of the graphs in (a), plume thickness D and bulk velocity    were smoothed with a running average of window size ±1 km.

Figure 7 .
Figure 7. Experiments on the sensitivity of 79° North Glacier to the open ocean stratification, with higher salinity (a, b)than in our default scenario (c, d) as well as lower salinity (e, f).For higher salinities above the sill, the subglacial plume propagates further along the ice tongue and detaches higher up.The salinity at the level of plume detachment is always around 33.7 g kg −1 .When the plume detaches early (e), a weaker secondary plume develops above.

Figure 8 .
Figure 8. Subglacial melt rate averaged over the whole ice tongue, 〈v b 〉, in the default model run and six sensitivity experiments with varying subglacial discharge, Q, (black circles) compared to three simple parametrizations (colored graphs).The Q 1/3 -law proposed byJenkins (2011) does not give a good fit (dotted, blue).The parametrization proportional to Q 0.31 used bySlater and Straneo (2022) is very similar to the Q 1/3 -law (not shown).A better fit is obtained by using an exponent of 0.16 ± 0.02, which is roughly a Q 1/6 -law (dashed, orange).The best fit is obtained by allowing a non-zero melt rate for zero subglacial discharge, which is close to a  √  law that is shifted upward by about 7 m yr −1 (solid, green).

Figure 10 .
Figure 10.Influence of the roughness length of the ice tongue, z 0,ice , on the circulation (a, b) and the temperature (c, d) below the ice tongue as well as on the Lagrangian time of the subglacial plume (e).The Lagrangian time t, defined as the integral of  d = d∕ ū , is shown up to the (principal) plume detachment, see panels (a), (b), and Figure 3; this position is identified by a clear reduction of under-ice velocity u and plume velocity    .