Marine melt in three dimensional Greenlandic sill fjord simulations

Abstract

Submarine glacier melt rates of the Greenland Ice Sheet remain a major uncertainty in climate model projections of future sea level rise. Development of submarine melt parameterizations has to a high degree relied on ocean circulation modelling of glacial fjords, designed to quantify effects such as ocean thermal forcing and fjord–glacier geometry. Greenlandic fjords are relatively narrow, and it is frequently assumed that across-fjord flow variations are small enough to allow marine melt to be quantified with two-dimensional ocean-circulation models. Here, we present three-dimensional model simulations showing that the interplay between fjord–glacier geometry, side wall friction, and Earth’s rotation makes the circulation in ice-shelf cavities three-dimensional even in narrow fjords. Remarkably, we find that Earth’s rotation changes the flow pattern in the cavity below the ice shelf, leading to a decrease in the marine melt on a 10 km wide ice shelf by a factor of five compared to a non-rotating simulation. Our study prompts using three-dimensional model configurations of Greenlandic fjords.

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, 1989; Jenkins, 1991, 2016, 2021, Losch, 2008; Lazeroms and others, 2018, 2019). 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, 2021). The development of ocean ice-melt parameterizations for climate ice sheet models is thus a priority (Nowicki and Seroussi, 2018). 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, 2013; Straneo and Cenedese, 2015). Many of the calving glaciers typically develop a floating shelf-like ‘tongue’ in the winter (e.g. Moyer and others, 2019).

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, 2020; Schaffer and others, 2020; Nilsson and others, 2023). 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, 2015). 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, 2015; Wiskandt and others, 2023 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, 2022; Wiskandt and others, 2023).

Figure 1. Schematic illustration of the Ryder Glacier and fjord, based on Figure 6 in Jakobsson and others (2020), 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, 2014; Cai and others, 2017; Reinert and others, 2023; Wiskandt and others, 2023; Bao and Moffat, 2024).

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 (2013) and Prakash and others (2022) have investigated Petermann Glacier and Wekerle and others (2024) and Kanzow and others (2024) studied the 79^∘N Glacier and Slater and coauthors (2018) 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, 2016; Fraser and others, 2018; Zhao and others, 2021, 2022). 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, 2018).

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, 2018). Notably, Sherard Osborne Fjord has a prominent sill that partly blocks the inflow of warm Atlantic Water to Ryder Glacier (Jakobsson and others, 2020). 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 (2023) 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 (2020).

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, 1997; Adcroft and others, 2004). MITgcm 2D configurations have been used to study ice-sheet ocean interactions in Greenlandic fjords (Xu and others, 2012; Sciascia and others, 2013; Carroll and others, 2015; Cai and others, 2017; Jordan and others, 2018). 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\) m² s⁻¹ and \(\nu_v=10^{-3}\) m² s⁻¹ 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 (1999), 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 ρ: (1)\[ \rho = \rho_0\bigl[1-\alpha(T - T_0) + \beta(S - S_0)\bigr], \]

where the values of coefficients are given in SI.

The melt parametrization in MITgcm (SHELFICE package, Losch (2008)) solves the three equation melt parametrization (Holland and Jenkins, 1999; Jenkins and others, 2001), as in Hellmer and Olbers (1989), Holland and Jenkins (1999) and Wiskandt and others (2023): (2) \begin{align} T_b &= \lambda_1 S_b + \lambda_2 + \lambda_3 P_b \end{align} (3) \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} (4) \begin{align} \rho_i \gamma_S (S_w - S_b) &= -S_b q \end{align}

The interface boundary layer temperature (T_b) is the in-situ freezing point temperature obtained from the boundary layer pressure and salinity (P_b and S_b, 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 S_b and q, where q is the upward freshwater flux (negative melt rate, in units of freshwater mass per time) and L_i is the latent heat of fusion. Upward heat flux implies submarine melting (a downward freshwater flux), hence the minus sign (Losch, 2008). 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 (1999)), 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, 2024). 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 (2014): (5)\[ R_1 = \frac{1}{\pi f} \int_{-H}^0 \ N \ dz, \]

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⁻¹ and the rotating 3D simulations suggest 3.1–4.8 m yr⁻¹ depending on the fjord width.

Table 1.

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⁻¹)

	Average melt [m yr−1]
	1 km	2 km	5 km	10 km
2D	15.8	15.8	15.8	15.8
3D f=0	8.8	12.2	14.1	14.6
3D f≠0	4.8	4.2	3.1	3.3
	Integrated melt [m3 s−1]
	1 km	2 km	5 km	10 km
2D	12.5	25.0	62.5	122.5
3D f=0	6.9	19.2	55.9	113.3
3D f≠0	3.8	6.7	12.4	25.3

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 82^o 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, 2015). 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, 2023). 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.

Figure 2. Across fjord averages of temperature (background colour) and streamfunction (in m² s⁻¹, 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 m² s⁻¹ in the 10 km wide simulation (Figure 2b) to a maximum of 1.5 m² s⁻¹ for the 1 km wide simulation (Figure 2e). The exchange flow across the sill decreases from 2.0 m² s⁻¹ at 10 km width to 1.0 m² s⁻¹ 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⁻¹ for 10 km wide simulation to around 20 m yr⁻¹ 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).

Figure 3. Melt rate (background colour) and depth integrated streamfunction (in m³ s⁻¹, 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).

Figure 4. Across fjord averages of temperature (background colour) and streamfunction of the across fjord averaged flow per unit width (in m² s⁻¹, 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, 2023; Reinert and others, 2023), 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⁻¹, the melt rate in the outflow region reaches maximum values of around 20 m yr⁻¹ 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⁻¹) increases, which is enough to even increase the shelf averaged melt rate from 3.1 m yr⁻¹ (5 km width) to 3.3 m yr⁻¹ (10 km width, Table 1).

Figure 5. 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₀), 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₀ 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, 2016, 2021). 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, 2012; Sciascia and others, 2013; Cai and others, 2017; Reinert and others, 2023; Wiskandt and others, 2023). 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, 2022). We note that Millgate and others (2013), 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, 2021, 2022).

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 m³ s⁻¹ (Wilson and others, 2017). 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 m³ s⁻¹) and 3D (25.3 m³ s⁻¹). 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, 2024).

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, 1998; Nilsson and others, 2023) 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 (1999) and Wiskandt and others (2023). The dependence of submarine melt rates on drag coefficient has been done before (e.g. Dansereau and others, 2013). 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, 2013). 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, 2022). 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 (2023), 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, 2016), 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 (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, 2015).

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.

Author Affiliation

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.

References

Adcroft AJ, Hill C, Campin JM, Marshall J and Heimbach P (2004) Overview of the formulation and numerics of the MIT GCM. In Proceedings of the ECMWF seminar series on Numerical Methods, Recent developments in numerical methods for atmosphere and ocean modelling, 139–149, ECMWF https://www.ecmwf.int/en/elibrary/7642-overview-formulation-and-numerics-mit-gcm

Bao W and Moffat C (2024) Impact of shallow sills on circulation regimes and submarine melting in glacial fjords. The Cryosphere 18(1), 187–203. doi: 10.5194/tc-18-187-2024

Cai C, Rignot E, Menemenlis D and Nakayama Y (2017) Observations and modeling of ocean-induced melt beneath Petermann Glacier Ice Shelf in northwestern Greenland. Geophysical Research Letters 44(16), 8396–8403. doi: 10.1002/2017GL073711

Carroll D, Sutherland DA, Shroyer EL, Nash JD, Catania GA and Stearns LA (2015) Modeling turbulent subglacial meltwater plumes: Implications for fjord-scale buoyancy-driven circulation. Journal of Physical Oceanography 45(8), 2169–2185. doi: 10.1175/JPO-D-15-0033.1

Cowton T, Sole A, Nienow P, Slater D, Wilton D and Hanna E (2016) Controls on the transport of oceanic heat to Kangerdlugssuaq Glacier, East Greenland. Journal of Glaciology 62(236), 1167–1180. doi: 10.1017/jog.2016.117

Dansereau V, Heimbach P and Losch M (2013) Journal of Geophysical Research: Oceans. Journal of Geophysical Research: Oceans 2001(3805), 1–26. doi: 10.1002/2013JC008846.Received

Fox-Kemper B and 17 others (2021) Ocean, Cryosphere and Sea Level Change. In V Masson-Delmotte and 18 others (eds.), International Panel for Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, 1211–1362. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA. doi: 10.1017/9781009157896.011

Fraser NJ, Inall ME, Magaldi MG, Haine TWN and Jones SC (2018) Wintertime Fjord–Shelf Interaction and Ice Sheet Melting in Southeast Greenland. Journal of Geophysical Research: Oceans 123(12), 9156–9177. doi: 10.1029/2018JC014435

Hellmer H and Olbers D (1989) A two-dimensional model for the thermohaline circulation under an ice shelf. Antarctic Science 1(4), 325–336. doi: 10.1017/S0954102089000490

Hill EA, Rachel Carr J, Stokes CR and Hilmar Gudmundsson G (2018) Dynamic changes in outlet glaciers in northern Greenland from 1948 to 2015. Cryosphere 12(10), 3243–3263. doi: 10.5194/tc-12-3243-2018

Holland DM and Jenkins A (1999) Modeling thermodynamic ice-ocean interactions at the base of an ice shelf. Journal of Physical Oceanography, 29(8 PART 1), 1787–1800. doi: 10.1175/1520-0485(1999)029<1787:mtioia>2.0.co;2

Jackson RH, Lentz SJ and Straneo F (2018) Dynamics of Shelf Forcing in Greenlandic Fjords. Journal of Physical Oceanography 48, 2799–2827. doi: 10.1175/JPO-D-18-0057.1

Jakobsson M and 30 others (2020) Ryder Glacier in northwest Greenland is shielded from warm Atlantic water by a bathymetric sill. Communications Earth &amp; Environment 1(1), 1–10. doi: 10.1038/s43247-020-00043-0

Jenkins A (1991) A one-dimensional model of ice shelf-ocean interaction. Journal of Geophysical Research 96(C11), 20671. doi: 10.1029/91JC01842

Jenkins A (2016) A Simple Model of the Ice Shelf-Ocean Boundary Layer and Current. Journal of Physical Oceanography 46(6), 1785–1803. doi: 10.1175/JPO-D-15-0194.1

Jenkins A (2021) Shear, Stability and Mixing within the Ice-Shelf-Ocean Boundary Current. Journal of Physical Oceanography 51, 2129–2148. doi: 10.1175/JPO-D-20-0096.1

Jenkins A, Hellmer HH and Holland DM (2001) The Role of Meltwater Advection in the Formulation of Conservative Boundary Conditions at an Ice-Ocean Interface. Journal of Physical Oceanography 31(1), 285–296. doi: 10.1175/1520-0485(2001)031<0285:TROMAI>2.0.CO;2

Jordan JR and 7 others (2018) Ocean-Forced Ice-Shelf Thinning in a Synchronously Coupled Ice-Ocean Model. Journal of Geophysical Research: Oceans 123(2), 864–882. doi: 10.1002/2017JC013251

Kanzow T and 28 others (2024) The atmosphere-land/ice-ocean system in the region near the 79N Glacier in Northeast Greenland: Synthesis and key findings from GROCE. doi: 10.5194/egusphere-2024-757

Lazeroms WM, Jenkins A, Hilmar Gudmundsson G and Van De Wal RS (2018) Modelling present-day basal melt rates for Antarctic ice shelves using a parametrization of buoyant meltwater plumes. Cryosphere 12(1), 49–70. doi: 10.5194/tc-12-49-2018

Lazeroms WM, Jenkins A, Rienstra SW and Van De Wal RS (2019) An analytical derivation of ice-shelf basal melt based on the dynamics of meltwater plumes. Journal of Physical Oceanography 49(4), 917–939. doi: 10.1175/JPO-D-18-0131.1

Losch M (2008) Modeling ice shelf cavities in a z coordinate ocean general circulation model. Journal of Geophysical Research: Oceans 113(8), 1–15. doi: 10.1029/2007JC004368

Marshall J, Hill C, Perelman L and Adcroft A (1997) Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. Journal of Geophysical Research C: Oceans 102(C3), 5733–5752. doi: 10.1029/96JC02776

Millgate T, Holland PR, Jenkins A and Johnson HL (2013) The effect of basal channels on oceanic ice-shelf melting. Journal of Geophysical Research: Oceans 118(12), doi: 10.1002/2013JC009402

Moyer AN, Nienow PW, Gourmelen N, Sole AJ, Slater DA, Truffer M and Fahnestock M (2019) Spatio-temporal variations in seasonal ice shelf submarine melt rate at a tidewater glacier in southwest Greenland. Journal of Glaciology 65(252), 523–530. doi: 10.1017/jog.2019.27

Nilsson J, Van Dongen E, Jakobsson M, O’Regan M and Stranne C (2023) Hydraulic suppression of basal glacier melt in sill fjords. The Cryosphere 17(6), 2455–2476. doi: 10.5194/tc-17-2455-2023

Nowicki S and Seroussi H (2018) Projections of future sea level contributions from the Greenland and Antarctic ice sheets: Challenges beyond dynamical ice sheet modeling. Oceanography 31(2), 109–117. doi: 10.5670/oceanog.2018.216

Nurser AJG and Bacon S (2014) The Rossby radius in the Arctic Ocean. Ocean Science 10(6), 967–975. doi: 10.5194/os-10-967-2014

Prakash A, Zhou Q, Hattermann T, Bao W, Graversen R and Kirchner N (2022) A nested high-resolution unstructured grid 3-D ocean-sea ice-ice shelf setup for numerical investigations of the Petermann ice shelf and fjord. MethodsX 9+, 101668. doi: 10.1016/j.mex.2022.101668

Reinert M, Lorenz M, Klingbeil K, Büchmann B and Burchard H (2023) High–Resolution Simulations of the Plume Dynamics in an Idealized \(79^\circ\mathrm{N}\) Glacier Cavity Using Adaptive Vertical Coordinates. Journal of Advances in Modeling Earth Systems 15(10), e2023MS003721. doi: 10.1029/2023MS003721

Schaffer J, Kanzow T, von Appen WJ, von Albedyll L, Arndt JE and Roberts DH (2020) Bathymetry constrains ocean heat supply to Greenland’s largest glacier shelf. Nature Geoscience 13(3), 227–231. doi: 10.1038/s41561-019-0529-x

Sciascia R, Cenedese C, Nicolì D, Heimbach P and Straneo F (2014) Impact of periodic intermediary flows on submarine melting of a Greenland glacier. Journal of Geophysical Research: Oceans 119(10), 7078–7098. doi: 10.1002/2014JC009953

Sciascia R, Straneo F, Cenedese C and Heimbach P (2013) Seasonal variability of submarine melt rate and circulation in an East Greenland fjord. Journal of Geophysical Research: Oceans 118(5), 2492–2506. doi: 10.1002/jgrc.20142

Slater DA, Goldberg DN, Nienow PW and Cowton TR (2016) Scalings for Submarine Melting at Tidewater Glaciers from Buoyant Plume. Theory. J. Phys. Oceanogr. 46, 1839–1855. doi: 10.1175/JPO-D-15-0132.1

Slater DA, Nienow PW, Cowton TR, Goldberg DN, Sole AJ (2015) ffect of near-terminus subglacial hydrology on tidewater glacier submarine melt rates. Geophys. Res. Lett. 42, 2861–2868. doi: 10.1002/2014GL062494

Slater DA and Straneo F (2022) Submarine melting of glaciers in Greenland amplified by atmospheric warming. Nature Geoscience 15(10), 794–799. doi: 10.1038/s41561-022-01035-9

Slater DA, Straneo F, Das SB, Richards CG, Wagner TJ, Nienow PW (2018) Localized plumes drive front-wide ocean melting of a Greenlandic tidewater glacier. Geophys. Res. Lett. 45, 12350–12358. doi: 10.1029/2018GL080763

Straneo F and 14 others (2013) Challenges to understanding the dynamic response of Greenland’s marine terminating glaciers to oceanic and atmospheric forcing. BAMS 94, 1131–1144. doi: 10.1175/BAMS-D-12-00100

Straneo F and Cenedese C (2015) The Dynamics of Greenland’s Glacial Fjords and Their Role in Climate. Annual Review of Marine Science 7(1), 89–112. doi: 10.1146/annurev-marine-010213-135133

Vaňková I and 7 others (2024) Subglacial discharge effects on basal melting of a rotating, idealized ice shelf. EGUsphere [preprint]. doi: 10.5194/egusphere-2024-2297

Wekerle C and 8 others (2024) Atlantic Water warming increases melt below Northeast Greenland’s last floating ice tongue. Nature Communications 15(1), 1336. doi: 10.1038/s41467-024-45650-z

Whitehead JA (1998) Topographic control of oceanic flows in deep passages and straits. Reviews of Geophysics 36(3), 423–440. doi: 10.1029/98RG01014

Wilson N, Straneo F and Heimbach P (2017) Satellite-derived submarine melt rates and mass balance (2011–2015) for Greenland’s largest remaining ice tongues. The Cryosphere 11(6), 2773–2782. doi: 10.5194/tc-11-2773-2017

Wiskandt J and N Jourdain (2024) Brief Communication: Representation of heat conduction into the ice in marine ice shelf melt modeling. EGUsphere, Copernicus GmbH, publisher, pp.1–10. doi: 10.5194/egusphere-2024-2239

Wiskandt J, Koszalka IM and Nilsson J (2023) Basal melt rates and ocean circulation under the Ryder Glacier ice tongue and their response to climate warming: a high-resolution modelling study. The Cryosphere 17(7), 2755–2777. doi: 10.5194/tc-17-2755-2023

Xu Y, Rignot E, Menemenlis D and Koppes M (2012) Numerical experiments on subaqueous melting of Greenland tidewater glaciers in response to ocean warming and enhanced subglacial discharge. Annals of Glaciology 53(60), 229–234. doi: 10.3189/2012AoG60A139

Zhao KX, Stewart AL and McWilliams JC (2021) Geometric Constraints on Glacial Fjord–Shelf Exchange. Journal of Physical Oceanography 51(4), 1223–1246. doi: 10.1175/JPO-D-20-0091.1

Zhao KX, Stewart AL and McWilliams JC (2022) Linking Overturning, Recirculation, and Melt in Glacial Fjords. Geophysical Research Letters 49(15), e2021GL095706. doi: 10.1029/2021GL095706