1. Introduction
Early theoretical and numerical models of the ocean-driven (marine) glacier melt focused on the Antarctic Ice Sheet, where the melt is forced by warm ocean waters circulating under the floating ice shelves connected to the Southern Ocean; several theoretical models for the buoyancy (melt-) driven circulation under the ice, as well as submarine melt parameterizations have been proposed (Hellmer and Olbers, Reference Hellmer and Olbers1989; Jenkins, Reference Jenkins1991, Reference Jenkins2016, Reference Jenkins2021, Losch, Reference Losch2008; Lazeroms and others, Reference Lazeroms, Jenkins, Hilmar Gudmundsson and RS2018, Reference Lazeroms, Jenkins, Rienstra and RS2019). The Greenland Ice Sheet (GrIS) is an important element of the global climate system, storing about six meters of sea level equivalent of ice (Fox-Kemper and others, Reference Fox-Kemper and Masson-Delmotte2021). The development of ocean ice-melt parameterizations for climate ice sheet models is thus a priority (Nowicki and Seroussi, Reference Nowicki and Seroussi2018). Although the basic physics and thermodynamics of submarine melt are the same as in Antarctica, typical GrIS conditions present additional complicating factors.
Greenland’s marine terminating glaciers can be either ice shelves (in the north) or calving quasi-vertical ice cliffs (in the south) and connect to the ocean in long and narrow glacial fjords (Straneo and others, Reference Straneo2013; Straneo and Cenedese, Reference Straneo and Cenedese2015). Many of the calving glaciers typically develop a floating shelf-like ‘tongue’ in the winter (e.g. Moyer and others, Reference Moyer, Nienow, Gourmelen, Sole, Slater, Truffer and Fahnestock2019).
The fjords have often complex bathymetry with sills (submerged moraines from previous more extensive ice margins), exerting controls on the circulation in the fjords (Jakobsson and others, Reference Jakobsson2020; Schaffer and others, Reference Schaffer, Kanzow, von Appen WJ, von Albedyll, Arndt and Roberts2020; Nilsson and others, Reference Nilsson, Van Dongen, Jakobsson, O’Regan and Stranne2023). Around Greenland, hydrographic observations from fjords show a deep inflow of a warm and salty layer of water with Atlantic origin (Atlantic Water, AW) overlain by a cold and fresh layer of Polar Surface Water (PSW) (Straneo and Cenedese, Reference Straneo and Cenedese2015). The AW inflow toward an ice shelf (or quasi-vertical ice cliff) provides a thermal energy source for the submarine glacier melt, resulting in an input of fresh, buoyant melt water (MW) that rises in a turbulent plume along the ice–ocean interface. As the plume rises, it entrains ambient water reducing its buoyancy until it reaches the neutral buoyancy level upon which the plume detaches from the ice and forms an outflow layer of glacially modified water (GMW, the mixture of MW, AW, possibly PSW, and subglacial discharge (SGD), if present) toward the open ocean, closing the estuarine circulation (Straneo and Cenedese, Reference Straneo and Cenedese2015; Wiskandt and others, Reference Wiskandt, Koszalka and Nilsson2023 and Figure 1). The marine terminating glaciers of the GrIS generally have larger subglacial discharges of the surface ice melt at the grounding line than in Antarctica; the SGD provides an additional fresh water source and thus stronger buoyancy forcing, faster circulation close to the ice–ocean interface and increased melt rates (Slater and Straneo, Reference Slater and Straneo2022; Wiskandt and others, Reference Wiskandt, Koszalka and Nilsson2023).
Schematic illustration of the Ryder Glacier and fjord, based on Figure 6 in Jakobsson and others (Reference Jakobsson2020), reused in accordance with Creative Commons CC BY license (https://creativecommons.org/licenses/by/4.0/). It shows the bathymetric features, inflow of Atlantic Water (AW), overlain by outflow of glacially modified water (GMW), and essentially stagnant Polar Surface Water (PSW). The colours delineate temperature, ranging from cold PSW (∼ –1.8 ∘C) to warmer AW (∼ 0.3 ∘C). At Ryder Glacier the sill depth varies across the fjord from ∼400 to ∼200 m (light grey shading; for simplicity, a constant sill depth of 400 m is used in the numerical simulations). Below the ice shelf, there is a buoyant plume (grey arrow) driven by basal (submarine) melt and affected by mixing with fjord waters (mixing is indicated by the circular eddy features). On the glacier side of the sill, the AW inflow accelerates and mixes with colder outflowing GMW, causing reduced temperatures at the grounding line. The downstream acceleration and mixing of inflowing AW are reproduced in numerical simulations (see Figures 2 and 4). The insert map in the upper left corner shows Greenland with the location of Ryder Glacier marked in a magenta circle.
Due to inaccessibility of the area, observations from Greenlandic outlet glaciers and fjords, particularly in the northern sector, are too scarce to resolve the topographic features and intricacies of the ocean circulation. High-resolution numerical model simulations of individual fjords and glaciers based on available observations are used to gain understanding of the physical processes underlying submarine melt and develop melt parameterizations. But high-resolution simulations are computationally expensive, so it is necessary to find a trade-off between the level of complexity and information gain. The ocean density stratification (vertical distribution of sea water density dependent on temperature, salinity and pressure) determines the baroclinic Rossby radius, the length scale at which the Earth’s rotation becomes important for the stratified ocean circulation. At latitudes associated with Greenland, the Rossby radius is typically comparable to or larger than the width of the fjords (of the order of 10 km). This, assuming further that the ice geometry and bathymetric variations in cross-fjord direction have negligible influence on the fjord circulation, motivates using two-dimensional (along fjord-along ice front and vertical) numerical model configurations for the sake of numerical efficiency. Two-dimensional (2D) models have been widely used to examine the submarine melt rates in Greenlandic glacial fjords forced by varying ocean temperatures and subglacial discharge rates, with an ultimate goal to develop submarine melt parameterizations for climate models (Sciascia and others, Reference Sciascia, Cenedese, Nicolì, Heimbach and Straneo2014; Cai and others, Reference Cai, Rignot, Menemenlis and Nakayama2017; Reinert and others, Reference Reinert, Lorenz, Klingbeil, Büchmann and Burchard2023; Wiskandt and others, Reference Wiskandt, Koszalka and Nilsson2023; Bao and Moffat, Reference Bao and Moffat2024).
Only a few previous studies have investigated the three-dimensional ocean-circulation and submarine melt patterns in ice cavities below ice shelves in North Greenland: Millgate and others (Reference Millgate, Holland, Jenkins and Johnson2013) and Prakash and others (Reference Prakash, Zhou, Hattermann, Bao, Graversen and Kirchner2022) have investigated Petermann Glacier and Wekerle and others (Reference Wekerle2024) and Kanzow and others (Reference Kanzow2024) studied the 79∘N Glacier and Slater and coauthors (Reference Slater, Straneo, Das, Richards, Wagner and Nienow2018) studied the Sarqardleq Fjord, West Greenland. There are also studies examining how fjord geometry and Earth’s rotation affect the oceanic heat transport to the vertical faces of Greenlandic marine-terminating tide-water glaciers. These studies have mainly focused on the larger-scale approximately hydrostatic fjord circulations and their impacts of submarine melt (Cowton and others, Reference Cowton, Sole, Nienow, Slater, Wilton and Hanna2016; Fraser and others, Reference Fraser, Inall, Magaldi, Haine and Jones2018; Zhao and others, Reference Zhao, Stewart and McWilliams2021, Reference Zhao, Stewart and McWilliams2022). Finally, there are studies investigating the transient dynamics of Greenlandic fjords of varying width subject to forcing from the shelf that do not address the submarine melt rates at all (e.g. Jackson and others, Reference Jackson, Lentz and Straneo2018).
The present model simulations are partly inspired by the conditions in Sherard Osborne Fjord in North Greenland, where Ryder Glacier, Greenland’s third largest remaining ice shelf, drains (Hill and others, Reference Hill, Rachel Carr, Stokes and Hilmar Gudmundsson2018). Notably, Sherard Osborne Fjord has a prominent sill that partly blocks the inflow of warm Atlantic Water to Ryder Glacier (Jakobsson and others, Reference Jakobsson2020). This paper focuses on how the effects of Earth’s rotation affect the submarine melt rate of an ice shelf by changing the structure of the flow in the cavity below the ice shelf.
2. Methods
To show that the two-dimensional model approach, albeit appealing for its efficiency, can be imperfect, we use the Massachusetts Institute of Technology General Circulation Model (MITgcm) configuration of Wiskandt and others (Reference Wiskandt, Koszalka and Nilsson2023) with a rigid lid, initialized with hydrographic profiles collected during a survey of Ryder Glacier and Sherard Osborn Fjord in northern Greenland by Jakobsson and others (Reference Jakobsson2020).
The model solves the non-hydrostatic Boussinesq form of the Navier–Stokes Equations using a finite volume formulation on an Arakawa C-grid with vertical z-levels employing partial cells (Marshall and others, Reference Marshall, Hill, Perelman and Adcroft1997; Adcroft and others, Reference Adcroft, Hill, Campin, Marshall and Heimbach2004). MITgcm 2D configurations have been used to study ice-sheet ocean interactions in Greenlandic fjords (Xu and others, Reference Xu, Rignot, Menemenlis and Koppes2012; Sciascia and others, Reference Sciascia, Straneo, Cenedese and Heimbach2013; Carroll and others, Reference Carroll, Sutherland, Shroyer, Nash, Catania and Stearns2015; Cai and others, Reference Cai, Rignot, Menemenlis and Nakayama2017; Jordan and others, Reference Jordan2018). The relevant processes that we are aiming to represent with our configuration are shown schematically in Figure 1. In particular, our configuration aims to represent the quiescent features of the plume dynamics and the fjord dynamics explicitly, i.e. without the need to employ a plume parameterization. The 3D domain is varying in width (1, 2, 5 and 10 km wide model experiments), 64 km long and 810 m deep, with resolution of 
$dz=3\,\mathrm{m}$ in the vertical and 
$dx=dy=100\,\mathrm{m}$ in the horizontal. It includes an ice shelf with a grounding line at 800 m and at 
$x=0\,\mathrm{km}$ and a 400 m deep sill centred around x=45 km. From the grounding line, the ice base rises linearly to a 200 m deep vertical front at 
$x=25\,\mathrm{km}$, resulting in a slope of s = 0.024.
Subgrid-scale processes are parameterized using a Laplacian eddy diffusion of temperature, salinity and momentum with constant coefficients. To study the effect of rotation in 3D, we use the exact same model parameters in the 2D and the 3D simulations. We use a horizontal viscosity of 
$\nu_h=2.5$ m2 s−1 and 
$\nu_v=10^{-3}$ m2 s−1 in the vertical.
The MITgcm applies the semi-implicit pressure method for non-hydrostatic equations with a rigid-lid, variables co-located in time and with Adams–Bashforth time stepping. The advective operator for momentum is second-order accurate in space. We apply a third-order direct space–time tracer advection scheme with flux limiter due to Sweby.
On solid boundaries (sea floor, side walls and ice base) no-slip conditions are applied, together with a quadratic drag 
$C_d =1.5 \times 10^{-3}$ as in Holland and Jenkins (Reference Holland and Jenkins1999), while free-slip conditions are applied at the ocean free surface. The northern border of the fjord (at x = 64 km) is the only open boundary. The outflow is balanced at the boundary, yielding a net-zero cross-boundary flow. Initial conditions are taken as the average of CTD-profiles taken glacierward of the sill. The boundary conditions are the average of CTD-profiles outside the first fjord’s sill and are implemented in a 2 km wide sponge at the open boundary. Temperature and salinity are restored to the boundary conditions in this restoring zone with a restoring timescale of 1 d at the innermost grid point (x = 62 km) and 1 h at the outermost point (x = 64 km).
For simplicity and because the nonlinear effects are small in the range of S–T values we are considering, a linear equation of state for the density ρ:
\begin{equation}
\rho = \rho_0\bigl[1-\alpha(T - T_0) + \beta(S - S_0)\bigr],
\end{equation}where the values of coefficients are given in SI.
The melt parametrization in MITgcm (SHELFICE package, Losch (Reference Losch2008)) solves the three equation melt parametrization (Holland and Jenkins, Reference Holland and Jenkins1999; Jenkins and others, Reference Jenkins, Hellmer and Holland2001), as in Hellmer and Olbers (Reference Hellmer and Olbers1989), Holland and Jenkins (Reference Holland and Jenkins1999) and Wiskandt and others (Reference Wiskandt, Koszalka and Nilsson2023):
\begin{align}
T_b &= \lambda_1 S_b + \lambda_2 + \lambda_3 P_b
\end{align}
\begin{align}
c_{p,w} \rho_i \gamma_T(T_w-T_b) &= -L_i q - \rho_i c_{p,i} \kappa_i \frac{(T_s-T_b)}{H_i}
\end{align}
\begin{align}
\rho_i \gamma_S (S_w - S_b) &= -S_b q
\end{align}The interface boundary layer temperature (Tb) is the in-situ freezing point temperature obtained from the boundary layer pressure and salinity (Pb and Sb, respectively) using the linear equation of state (Eqn 1) where λj are constants. Equations (3) and (4), which describe heat and salt balances at the interface, respectively, are used to calculate Sb and q, where q is the upward freshwater flux (negative melt rate, in units of freshwater mass per time) and Li is the latent heat of fusion. Upward heat flux implies submarine melting (a downward freshwater flux), hence the minus sign (Losch, Reference Losch2008). Note that the submarine melt rates depend directly on both ocean temperature and velocity at the ice–ocean boundary (the freezing temperature is also dependent on salinity and pressure). The temperature gradient in the ice depends on the melt rate (Equation 26 and 31 in Holland and Jenkins (Reference Holland and Jenkins1999)), assuming constant vertical advection and vertical diffusion of temperature into the ice, to accurately represent the heat budget at the ice–ocean interface (Wiskandt and Jourdain, Reference Wiskandt and Jourdain2024). The fresh MW input is represented as a virtual flux (locally changing temperature and salinity, no volume flux into the ocean domain). There is no external momentum or heat forcing applied; the model experiments are started from rest, and the circulation in the fjord is driven entirely by the virtual flux of the MW. The simulations for all experiments were run for 200 days and averages over the last 40 days were analysed, when the domain-integrated temperature and kinetic energy are in a statistically-steady state. A complete list of the model parameters is given in the Supplementary Information (SI).
The baroclinic Rossby radius is a length scale quantifying the importance of the Earth’s rotation for a stratified flow (the Coriolis force becomes important on horizontal scales larger than that) and sets the size of fronts and eddies in rotating stratified flows. The first mode Rossby radius is estimated as in Nurser and Bacon (Reference Nurser and Bacon2014):
\begin{equation}
R_1 = \frac{1}{\pi f} \int_{-H}^0 \ N \ dz,
\end{equation}where f is the Coriolis parameter, H is the depth of the water column and N is the buoyancy frequency.
3. Results
When comparing the area-average melt rates of 2D simulations to 3D simulations, we find immediately that there is a large discrepancy, even for fjords with a width of 10 km or less (Table 1): The 2D simulations suggest average melt rates of 15.8 m yr−1 and the rotating 3D simulations suggest 3.1–4.8 m yr−1 depending on the fjord width.
Shelf averaged and shelf integrated melt rate from 2D and 3D simulation of varying fjord widths, without Earth’s rotation (
$\mathrm{f}=0$) and including it (
$\mathrm f=1.44\;\times\;10^{-4}$ s −1)
To quantify the relative importance of friction on the side wall of the fjord and rotation, we first compare the results from the 2D non-rotating simulation to the results from a 3D non-rotating simulation (setting the Coriolis parameter f = 0) before comparing with the results from a 3D rotating simulation 
$f=1.44\;\times\;10^{-4}\,\mathrm s^{-1}$, corresponding to the Ryder Glacier latitude of 82o North). We are investigating differences in circulation patterns and strength and the resulting melt rate distribution underneath the floating ice shelf.
3.1. 2D and 3D without rotation
The 2D simulation exhibits an estuarine circulation of a salty warmer AW inflow and a fresher colder GMW outflow, driven by the buoyant plume due to submarine melt at the ice–ocean interface (Figure 2a). This flow pattern is typical for long and narrow fjords (Straneo and Cenedese, Reference Straneo and Cenedese2015). The inflow over the sill crest is hydraulically controlled, which limits the inflow compared to the case without or with a deeper sill (Nilsson and others, Reference Nilsson, Van Dongen, Jakobsson, O’Regan and Stranne2023). Notably, almost half of the GMW recirculates glacierward of the sill. The non-recirculating fraction of the GMW leaves the domain across the sill balanced by the AW inflow, in combination referred to as the exchange flow in the following (Figure 2a). This leads to cooler temperatures inside the sill, compared to outside the sill, which reduces the submarine ice melt (Figure 2a), compared to a no-sill case. The upper-most layer in the fjord (0–250 m) that exhibits a strong temperature, salinity and density gradient (Figure 2a, only temperature shown) is essentially stagnant, as the buoyant plume does not penetrate through the strong halocline in the upper water column.
Across fjord averages of temperature (background colour) and streamfunction (in m2 s−1, contour lines) for a 2D simulation (a) and 3D simulations without rotation of widths of 10 km (b), 5 km (c), 2 km (d) and 1 km (e). Solid (dashed) streamlines indicate clockwise (counterclockwise) circulation. Coloured lines in the top left corner show the melt rate distribution along the ice base (right y-axis). For 3D simulations, the green shows the across fjord average, blue shows the across fjord maximum and grey shows the 2D melt distribution for reference. The small-scale pattern apparent in the melt rates is due to implementation of the melt parameterization in MITgcm with the partial-cell discretization.
Close to the grounding line, the plume is slow and melt rates are low (Figure 2a). The plume accelerates along the ice, leading to an increase in melt with growing distance away from the grounding line. Once the plume detaches at around 20 km, the melt rate decreases rapidly towards zero.
In the non-rotating 3D simulations, the circulation shows a similar structure to the 2D simulation, for all fjord widths (Figure 2b–e). We note that for decreasing fjord width the across fjord averaged overturning strength (volume transport per unit width) glacierward of the sill decreases from a maximum of 3.5 m2 s−1 in the 10 km wide simulation (Figure 2b) to a maximum of 1.5 m2 s−1 for the 1 km wide simulation (Figure 2e). The exchange flow across the sill decreases from 2.0 m2 s−1 at 10 km width to 1.0 m2 s−1 at 1 km width (Figure 2b–e).
As both the Coriolis force and any geometric cross-fjord variations are absent, we attribute this decrease to lateral friction near the boundary becoming more important in narrower fjords. This suggests that for narrower domains, frictional effects near the side walls influence a greater proportion of domain, which acts to decrease the cross-fjord averaged flow speed.
As a result of the decreasing velocities, the melt rate maximum decreases from almost 30 m yr−1 for 10 km wide simulation to around 20 m yr−1 for 1 km wide simulation (Figure 2b–e). The distribution of melt rate along the ice base is similar for all fjord widths, with low melt rates close to the grounding line that increase with distance away from the grounding line. For narrower domains (Figure 2d,e) we observe that the melt rate decreases more gradually before the plume detaches, consistent with a gradual decrease in average plume speed and temperature in that region (Figure 2d,e). Due to the no-slip condition on side boundaries of the domain, velocities and hence the melt rate decrease close to the fjord walls (Figure 3a–d). When the fjord width decreases, the frictionally influenced near-wall layers with weak melt rates occupy an increasing proportion of the domain, causing the area averaged melt to decrease (Figure 5a). For the smallest fjord widths (1 and 2 km), decreases in maximum melt rates reinforce the reductions of area averaged melt (Figure 2b–e).
Melt rate (background colour) and depth integrated streamfunction (in m3 s−1, contour lines) 3D simulations without (a–d) and with (e–h) rotation of widths of 10 km (a,e), 5 km (b,f), 2 km (c,g) and 1 km (d,h). Solid (dashed) streamlines indicate clockwise (counterclockwise) circulation. The small-scale pattern apparent in the melt rates is due to implementation of the melt parameterization in MITgcm with the partial-cell discretization.
The combination of increasing fractional area of the frictional near-wall layers and decreases in maximum melts (see blue melt lines in Figure 2), as the fjord width narrows, causes a decline in ice shelf averaged melt rate, along with an expected decline of integrated melt rate (Table 1 and Figure 5).
3.2. Effects of Earth’s rotation
In the rotating simulations, the flow develops a clear cross-fjord variation, which is evident in depth-averaged circulation glacierward of the sill: the inflow (toward the grounding line) is concentrated at depth on the western side of the fjord and the outflow is concentrated at the ice base on the eastern side (Figure 3 and SI Figure 2). The Rossby number of the flow is about 0.1, suggesting that the flow in the interior is in near geostrophic balance. Consistent with the small Rossby number, the vertical component of the circulation, i.e. the overturning, is weak compared to the non-rotating case (Figure 4a)—the circulation is primarily horizontal (SI Figure 1).
Across fjord averages of temperature (background colour) and streamfunction of the across fjord averaged flow per unit width (in m2 s−1, contour lines) for 3D simulations with rotation of widths of 10 km (a), 5 km (b), 2 km (c) and 1 km (d). Solid (dashed) streamlines indicate clockwise (counterclockwise) circulation; note the different intervals between contour lines. Coloured lines in the top left corner show the melt rate distribution along the ice base (right y-axis): green shows the across fjord average, blue shows the across fjord maximum and grey shows the 2D melt distribution for reference. The small-scale pattern apparent in the melt rates is due to implementation of the melt parameterization in MITgcm with the partial-cell discretization.
As in the 2D simulations, the fjord exchange flow (AW inflow and GMW outflow) is affected by the presence of the sill (Nilsson and others, Reference Nilsson, Van Dongen, Jakobsson, O’Regan and Stranne2023; Reinert and others, Reference Reinert, Lorenz, Klingbeil, Büchmann and Burchard2023), leading to a recirculation of GMW glacierward of the sill. Furthermore, the exchange flow across the sill is influenced by the rotation. The inflow of AW is concentrated on the western side of the sill (flowing with the fjord wall to the right; SI Figure 2). Notably, the highest outflow velocities are also encountered on the western side of the sill, indicating that some of the outflowing GMW crosses the fjord from the east to the west between the ice shelf edge and the sill crest (SI Figure 1–2).
In the open fjord, away from the ice shelf, Equation (5) yields the Rossby radius of deformation 
$R_1 \approx$ 10 km, which is a typical value for Greenlandic fjords. In the ice cavity, however, where the water column depth decreases and the PSW layer is absent, we find that 
$R_1 \approx$ 2 km. This length scale is comparable to the width of the outflow plume underneath the ice shelf (Figure 3).
Due to the confinement of high velocities to a narrow (
$\sim 2 $ km) outflow plume, the area of high melt is small compared to the non-rotating 3D simulations. While the majority of the ice shelf exhibits melt rates below 4 m yr−1, the melt rate in the outflow region reaches maximum values of around 20 m yr−1 for the 10 and 5 km wide simulations. Similar to the non-rotating 3D simulations, the melt rate maxima decrease with decreasing fjord width (Figure 4). Overall the averaged and shelf-integrated melt from 3D simulations is lower than from 2D circulation, calculated assuming the same width (Table 1 and Figure 5). The average melt rate decreases with increasing fjord width (Figure 5a) initially but the decrease saturates upon increasing the width from 5 to 10 km. While the area of very high melt does not increase significantly when comparing these two simulations (Figure 3e,f), the area of low to intermediate melt rate (2–6 m yr−1) increases, which is enough to even increase the shelf averaged melt rate from 3.1 m yr−1 (5 km width) to 3.3 m yr−1 (10 km width, Table 1).
Shelf averaged (a) and shelf integrated (b) melt rate from 2D and 3D simulation of varying fjord widths. In panel (a) the dashed line shows the predicted melt rate expected from the presence of a constant area of frictional wall layers (with a total width W 0), with lower melt rates. This relationship for the mean basal melt is given by 
$M= M_{0}(1-W_{0}/W)$; where 
$M_{0} = 15.2$ m yr-1 and 
$W_{0}=0.35$ km are fitted to match the 5 and 10 km results. The fitted W 0 corresponds approximately to the width of four grid cells. In panel (b), the dashed lines give order-of-magnitude estimates of rotating integrated melt rates based on non-rotating ones adjusted by the factor 
$R1/W$, where 
$R1 \sim 2$ km is the estimated Rossby radius in the ice cavity; see the text for details.
Figure 5b shows that the integrated basal melt increases approximately linearly with the fjord width for both the non-rotating and the rotating simulations, and that the integrated melt rates are roughly comparable for the smallest fjord widths. However, the increase in integrated melt rates with fjord width is significantly lower for the rotating simulations. This indicates that rotational effects act to decrease the basal melt, either by restricting the fjord circulation delivering the heat flux towards the glacier, and/or by affecting the local processes determining basal melt rates. The fact that the highest basal melt rates in the rotating simulations are confined to a limited area on one side of the ice shelf allow us to crudely estimate the integrated melt rates in the rotating simulations from the ones in the non-rotating ones as follows: assuming that the integrated melt in the rotating simulations comes from a fixed area (controlled by the ice cavity Rossby radius 
$R_1\sim 2$ km), then the rotating integrated melt rate will be a factor of 
$R_1/W$ smaller that the non-rotating ones (where melt rates are assumed to be constant over ice shelf). As shown in Figure 5b, this procedure brings rotating and non-rotating melt rates reasonably close, and can serve as a leading-order correction of rotational effects on the integrated basal melt rates. Additional difference between the shelf-integrated melt from rotating 3D simulations and from non-rotating 3D simulations could be due to the Ekman boundary layer developing below the gently sloping ice shelf in rotating 3D simulations (Jenkins, Reference Jenkins2016, Reference Jenkins2021). However, a detailed investigation of how these are realized and how they impact the melt rates as parameterized in the MITgcm configuration is left to future work.
4. Discussion and conclusions
Using the MITgcm we evaluate to what extent submarine glacier melt rate estimates from 2D simulations (effectively non-rotating) represent melt rates in rotating 3D simulations in long and narrow Greenlandic silled fjords with floating ice shelves. It is commonly assumed that for fjords narrower than the baroclinic Rossby radius of deformation (
$R_1 \approx 10$ km), rotational effects are small and the fjord circulation and the ice–ocean interactions can be represented by two-dimensional ocean-circulation model simulations (Xu and others, Reference Xu, Rignot, Menemenlis and Koppes2012; Sciascia and others, Reference Sciascia, Straneo, Cenedese and Heimbach2013; Cai and others, Reference Cai, Rignot, Menemenlis and Nakayama2017; Reinert and others, Reference Reinert, Lorenz, Klingbeil, Büchmann and Burchard2023; Wiskandt and others, Reference Wiskandt, Koszalka and Nilsson2023). Our simulations show that even for fjords with widths less than 10 km, rotational effects lead to strong across-fjord variability in circulation and melt rates below the ice shelf, creating a fjord-circulation pattern similar to that found in previous 3D model simulations in narrow fjords with ice–ocean interactions at a quasi-vertical ice cliff rather than at an ice shelf (Zhao and others, Reference Zhao, Stewart and McWilliams2022). We note that Millgate and others (Reference Millgate, Holland, Jenkins and Johnson2013), in ocean-circulation model simulations of Petermann Glacier with 20 km wide and 70 km long idealised ice shelf, obtained flow and melt patterns that are qualitatively similar to the one obtained here for the 10 km wide rotating Ryder simulation. This suggests that the rotational suppression in the area-averaged submarine melt seen in the present study (Figure 4a) will also be similar in fjords with widths on the order of 20 km.
We found that a contributing factor to the strong impact of rotation on the melt rates is due to the small Rossby radius inside the ice cavity (
$R_1 \approx2$ km), where the upper ocean halocline is absent and the stratification is weaker and the water column depth is smaller. We are able to link this circulation pattern to the distribution of submarine melt and the resulting shelf averaged melt rate drops fivefold compared to shelf averaged melt rates from 2D simulations. Despite investigating a fjord with a marine terminating ice shelf in this study, we also expect rotational effects on marine melt on quasi-vertical ice cliffs lacking floating ice shelves (Zhao and others, Reference Zhao, Stewart and McWilliams2021, Reference Zhao, Stewart and McWilliams2022).
The present idealized model setup is based on Ryder Glacier in Sherard Osborn Fjord in North Greenland, using a rather crude approximation of the fjord–glacier geometry (ice shelf base and sea bed). Marine melt rate estimates from Ryder Glacier based on satellite data and radar measurements suggest a total melt of 60 m3 s−1 (Wilson and others, Reference Wilson, Straneo and Heimbach2017). This value lies between the estimates of the integrated melt rates from the present simulation of a 10 km wide fjord in 2D (122.5 m3 s−1) and 3D (25.3 m3 s−1). For 79∘ North Glacier, a comparison of 2D and 3D simulations also shows a decrease due to rotational effects, albeit less drastic (Kanzow and others, Reference Kanzow2024).
In this study, we focus on the effect of a 3D circulation on the submarine melt rate distribution underneath an ice shelf, but there are several outstanding questions regarding the fjord circulation left to future work, for example, 3D hydraulics (Whitehead, Reference Whitehead1998; Nilsson and others, Reference Nilsson, Van Dongen, Jakobsson, O’Regan and Stranne2023) and the response of the 3D circulation to different sill depths and geometry. Note that we used a quadratic drag with a drag coefficient value as in Holland and Jenkins (Reference Holland and Jenkins1999) and Wiskandt and others (Reference Wiskandt, Koszalka and Nilsson2023). The dependence of submarine melt rates on drag coefficient has been done before (e.g. Dansereau and others, Reference Dansereau, Heimbach and Losch2013). However, the effects of friction will depend on model resolution. A study investigating friction and resolution is left to the future, as here we focus on the effects of rotation. A realistic, more complex bathymetry will further affect the circulation strength and pattern, similar to the variability of the geometry of the ice base and presence of channels that can generate cross-fjord variability in melt rate (Millgate and others, Reference Millgate, Holland, Jenkins and Johnson2013). Seasonal and interannual variability and external drivers like tides, local and non-local winds or effects from icebergs and sea ice are also not considered here. Another very important factor is the SGD flux during the summer (Slater and Straneo, Reference Slater and Straneo2022). In this work we aim at isolating the effects of rotation, so we do not apply SGD. Moreover, we do not have any observational estimates of SGD for Ryder Glacier to make our simulations relevant to summer conditions. In Wiskandt and others (Reference Wiskandt, Koszalka and Nilsson2023), we investigated the impact of SGD by a suite of 2D model experiments setting the SGD flux as a fraction of the average submarine melt flux for winter simulations. We found that the melt rates increase in a fractional manner with the SGD, consistent with previous studies (Slater and coauthors, Reference Slater, Goldberg, Nienow and Cowton2016), and the maxima in plume velocity and melt rates move closer to the grounding line. Based on these results, we expect that adding the SGD flux will move the velocity and melt rate maxima in Figure 3 closer to the grounding line and increase average melt rates. Our expectations are consistent with 3D model simulations of an Antarctic ice shelf of a fixed width and varying SGD flux by Vankova and others (Reference Vaňková2024). But they also found that the effect of SGD can change substantially depending on how SGD flux is distributed (several SGD channels instead of an uniform distribution, Slater and coauthors, Reference Slater, Nienow, Cowton, Goldberg and Sole2015).
We expect that the inclusion of all aforementioned effects will make the three-dimensional aspect of the circulation even more pronounced and complex, and will require extensive modelling studies in the future. We note however that in order to verify the future modelling studies investigating the 3D fjord circulation and its effects on the submarine melt, and to make them more realistic, it is imperative to collect more observations of the circulation in ice cavities of the Greenland glaciers and their submarine melt rate patterns.
To conclude, the presented results provide an evidence that adding a third dimension and rotation into glacial fjord simulations leads to considerable differences in circulation and melt rates compared to 2D simulations. These differences should be considered when interpreting results from 2D simulations and when developing submarine melt parameterizations for climate ice sheet modelling.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.10073.
Acknowledgements
JW would like to acknowledge the funding from the Faculty of Science, Stockholm University (grant no. SUFV-1.2.1-0124-17). The computations and data analysis were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at NSC partially funded by the Swedish Research Council through grant agreement no. 2022-06725.
Open access
