The recent acceleration of outflowing marine-terminating glaciers at the margins of the Greenland Ice Sheet and Antarctic Ice Sheet has received widespread attention (van den Broeke et al. 2016).
In the Greenland Ice Sheet, the accelerated melting is postulated to result from warming of deep ocean currents that come into contact with the termini of tidewater glaciers (Holland et al. 2008;Straneo and Heimbach 2013;Wood et al. 2018;Cowton et al. 2018) as well as a growing surface melt contribution (e.g., Hofer et al. 2020). This submarine melt at the sides of marine-terminating glaciers drives glacial retreat and also amplifies iceberg calving, depending on the properties of the glacier and fjord (Chauché et al. 2014;Rignot et al. 2015;Wagner et al. 2016;Morlighem et al. 2016;Fried et al. 2018;Slater et al. 2021;Wood et al. 2021). The submarine melt rate consists of ambient face-wide melt and discharge plume-driven melt (Straneo and Cenedese 2015;Jackson et al. 2019). Although subglacial discharge plumes have the potential to drive a melt rate of more than a meter per day in the glacial area near the plume, they only occupy a small fraction of the glacial face (Cowton et al. 2015;Slater et al. 2018). By comparison, face-wide melting can occur along the entire glacial face as a result of either convection (Magorrian and Wells 2016) or fjord circulation (Bartholomaus et al. 2013).
The focus of previous 2D and 3D simulations of the shelf-to-fjord system has been to understand the sensitivity of glacial melt and the overturning circulation/fjord renewal to various fjord characteristics and atmospheric/oceanic drivers (e.g., Gladish et al. 2015b, Sciascia et al. 2013, Xu et al. 2012, and Jackson et al. 2018). Along with the relative scarcity of ocean observations near marine-terminating glaciers, only recently has the horizontal recirculation within fjords and their sensitivity to fjord and forcing parameters received attention in models (Zhao et al. 2019(Zhao et al. , 2021)), which has been suggested to have an influence on the face-wide melt rates (Carroll et al. 2017;Slater et al. 2018;Jackson et al. 2019;Zhao et al. 2021). Existing melt parameterizations either do not take into account horizontal near-glacier velocities (e.g., Xu et al. 2012, Sciascia et al. 2013) or do not resolve the horizontal flows necessary for accurate melt rate predictions (e.g., Cowton et al. 2015, Carroll et al. 2017).
To better understand these processes in the context of realistic fjord geometries, we conduct highresolution fjord simulations of three major Greenland deep-water fjords and compare the emergent dynamics to results from previous studies. We will use simple dynamical theories of overturning circulation, horizontal recirculation in the fjord interior, and glacial melt rate from previous studies (Zhao et al. 2021;Zhao 2021) to understand the model behavior. Using these results, we address a gap in understanding of how 3D fjord circulation drives melt in realistic fjord geometries, which has important implications for glacial retreat at the oceanic margins of ice sheets.
In Sect. 2, we present our model setup, configuration, and design philosophy. We also present an overview of the phenomenology and dynamics of the three representative regional fjord simulations used in our study: Ilulissat, Sermilik, and Kangerdlugssuaq. In Sect. 3, we present the phenomenlogy of the long-lived standing eddies for each fjord simulated. As a representative example with important melt implications, we discuss the properties of the standing eddy near the Ilulissat glacial face, its spinup, transport, mergers with other eddies, and its influence on Lagrangian and Eulerian residence times within the fjord. We also briefly discuss the properties of the other eddies in our three regional simulations and their preferred locations relative to fjord bathymetry. In Sect. 4, we discuss the vorticity balance within Ilulissat fjord, which demonstrates the importance of the standing eddies to the integrated vorticity budget and use this to develop a scaling theory for the vertical profile of the near glacial horizontal velocity. We use this theory to predict the vertical profile of melt rate and discuss how standing eddies influence and in some cases, potentially amplify glacial melt. In Sect. 5, we summarize and discuss the major implications and caveats of our findings and suggest future avenues of research.
The design of our model setup is primarily motivated by the need for an improved understanding of the 3D circulation within warm, deep-water fjords with realistic geometries. Various characteristics of the 3D circulation were previously shown using fjord idealized geometries to be an important factor in determining glacial melt rates (Zhao 2021).
Fig. 1 shows the bathymetry around Greenland using the Bedmachine V3 dataset (Morlighem et al. 2017), which is a compilation of sonar measurements, depth soundings, and gravity inversions. Along the perimeter, glacial fjords connect marine-terminating glaciers to the ocean on the continental shelf. The zoomed-in panels show the three Greenlandic fjord-shelf domains (Ilulissat, Kangerdlugssuaq, and Sermilik) selected for this study, which are some of the widest and longest of Greenland's warm, deep-water fjords. We selected these fjords because they generally have a larger flux of solid ice and freshwater into the ocean and can be more easily resolved. For each of these regional simulations, our primary aim is to capture the key drivers of the steady-state summertime fjord circulation and to ignore time-varying and secondary effects or those that we cannot currently adequately represent. We anticipate that these models will greatly benefit from the inclusion of parameterizations of unresolved ice-ocean processes and have significant room for improvement in the future. The model configuration specifics are described in the following two subsections.
The model used in the study is the Massachusetts Institute of Technology General Circulation Model (MITgcm) (Marshall et al. 1997; see data availability statement). Using this model, we solve the hydrostatic, Boussinesq primitive equations with a nonlinear equation of state based on Jackett and McDougall (1995) in three high-resolution configurations of Ilulissat, Sermilik, and Kangerdlugssuaq fjords and a small area of the adjacent continental shelf for each simulation.
The model bathymetry for each of the three regions uses the Bedmachine V3 data (Morlighem et al. 2017), which has a 150 m horizontal resolution (see Fig. 1). Any gridpoints with a bathymetry shallower than 20 m were modified to dry cells and the entire model grid bathymetry field was then slightly smoothed with a 5-gridpoint Gaussian filter to reduce spurious sources of vorticity.
The model domain dimensions 𝐿 × 𝑊 × 𝐻 vary by region and are presented for each region separately in Sect. 2c. The model horizontal resolution is 150 m (the same as bathymetry data) and the vertical resolution varies slightly between the different regions. We use a Smagorinsky biharmonic horizontal viscosity and the K-Profile Parameterization (KPP) of the vertical viscosity and diffusivity (Smagorinsky, 1963;Large et al., 1994), in addition to a background vertical diffusivity of 10 -6 m 2 s -1 . There is quadratic bottom drag with a coefficient of 2 × 10 -3 . We use an 𝑓 -plane approximation with a representative Coriolis parameter of 𝑓 = 1.31 × 10 -4 s -1 , which approximately corresponds to the latitude of the fjords in this study. The model experiments are run for 1 year because the fjord recirculation adjusts slowly and requires multiple months of spinup.
Our boundary conditions allow us to achieve a fully spun-up steady-state 3D circulation within each fjord in as simplified a way as possible. The buoyancy drivers in our simulations are supplied by the open ocean and (vertical) glacial face boundary conditions, which makes their representation critical in our simplified forcing. In this subsection, we discuss the following: first, the open-ocean forcing, second, the glacial face forcing, and finally, the exclusion of surface forcing.
Each of the domains has three open-ocean boundaries on the shelf region with one inflow boundary on the shelf. These inflow boundaries are forced by an inflow of a time-invariant vertical profiles of temperature, salinity, and velocity. The inflow boundaries are the southern, eastern, and northern boundaries of Ilulissat, Sermilik, and Kangerdlugssuaq, respectively (see Fig. 1 for the inflow boundaries in each domain). All other boundaries used the Orlanski radiation boundary conditions (Orlanski 1976). The inflow temperature and salinity uses the Oceans Melting Greenland (OMG) AXCTD (Airborne eXpendable Conductivity Temperature Depth) data (Fenty et al. 2016; see data availability statement for individual AXCTDs used). For each fjord, the AXCTDs used correspond to those that were positioned closest to the shelf region within each domain during the 2020 summer season (July to September).
The boundary-normal velocity is imposed as a constant over the full cross-sectional area at the inflowing open-ocean boundary (see subsection c for the imposed boundary condition total barotropic transport for each region) and is derived from approximations of the along-coast integrated transport. This coastal transport has broadband temporal variability, but is assumed to constant during the summer months in our idealized setup. The integrated transport is motivated by a combination of Estimating the Climate and Circulation of the Ocean (ECCO) Version 5, Release alpha (≈ 11 km horizontal resolution, Zhang et al. 2018) coastal transport, existing data (Gladish et al. 2015b, Sutherland et al. 2014, Straneo and Cenedese 2015, and references therein), and near-coast (≈ 10 km offshore in these domains) sea surface height gradients from the Making Earth System Data Records for Use in Research Environments (MEaSUREs) dataset (Zlotnicki et al. 2019). None of these methods provide high accuracy of the near-coast integrated transport near these fjords, which are specified for each domain in the following subsection. We did not find strong sensitivity of the fjord's circulation to changes in the inflow velocity. However, future iterations would benefit from improved realistic coastal variability, which is likely to lead to fjord flushing events and is not explored in this study (e.g., Gladish et al. 2015a).
Subglacial discharge exits at the base of the glacier and is typically 200 to 400 m 3 /s in the summer for these fjords in our simulations (but may vary from 0 to 2000 m 3 /s throughout the melt season) and nearly zero in winter (Straneo and Cenedese 2015;Chu 2014). Areas of elevated mixing within the fjord-to-shelf region are primarily forced by subglacial and ambient melt plumes as they are a dominant mechanism of mixing for the majority of Greenland's fjords (Carroll et al. 2017;Gladish et al. 2015b;Magorrian and Wells 2016). We therefore use plume parameterizations for both the discharge plume and the melt plume across all ocean-glacial boundaries. The glacial geometry used Bedmachine V3 data with a submerged glacier-ocean interface that is assumed to be vertical. This is because the slope of the ice-ocean boundaries at the submerged glacier wall is not well measured and is likely changing substantially on subseasonal or longer timescales (see e.g., Sutherland et al. 2019).
The plume parametrizations used in our model setup are based on buoyant plume theory for a point plume, which solves 1-dimensional equations (vertical profiles) for mass and momentum conservation within the plume. The adjacent temperature/salinity profiles evolve in response to advection, entrainment of ambient waters into the plume and outflow from the plume, and the turbulent transfer of heat and salt between the plume and the ice (Hellmer and Olbers 1989).
The plume is coupled to the circulation and stratification in the MITgcm model configuration and is a slightly modified version of that proposed by Cowton et al. (2015), optimized to work efficiently in high resolution simulations (see data availability statement). This is identical to the parameterization package detailed in Cowton et al. (2015) for point plumes, except that we redistribute the buoyancy anomalies from the solutions to the discharge plume equations over a 5-gridpoint-radius semi-circle in the horizontal and apply a 3-gridpoint smoothing in the vertical while conserving the overall buoyancy anomaly and entrainment (similarly to Zhao et al. 2022, which uses a similar smoothing process). This prevents prohibitive restrictions set by the Courant-Friedrichs-Lewy (CFL) condition on the model timestep in our high resolution simulations as well as spurious mixing caused by sharp gradients in the forcing at the gridscale. Tests conducted with a 10-gridpoint-radius semi-circle in the horizontal did not significantly alter the vorticity balance presented in Section 4. The subglacial discharge uses a steady 2-month time average of summer discharge averaged over the years 2017 to 2019 based on the outflow locations and discharge magnitudes from the Mankoff et al. (2020) dataset. Lastly, we discuss briefly the exclusion of surface boundary forcing. We note that these simplified fjord-shelf regional configurations are not intended to fully represent the dynamics of Greenland's fjords, but rather to capture a few salient features that include more realism not present in previous studies (Gladish et al. 2015b;Carroll et al. 2017;Zhao et al. 2021;Zhao 2021). However, we did test the sensitivity of our regional simulations to steady winds (stresses of up to 0.15 N/m 2 ), a thin layer of static sea ice throughout the domain (using the sea ice model from (Losch et al. 2010)), and mean summer atmosphere forcing (temperature and freshwater fluxes), which did not have a noticeable influence on the fjord circulation below 100 m depth. Intermittent, strong katabatic winds are likely important for fjord dynamics as they may lead to flushing events (e.g., Spall et al. 2017), but the temporal variability of fjord dynamics is not investigated in the present study. One reason for this is that the effect of coastal Greenland air-sea interaction (atmospheric temperature, air-sea freshwater fluxes, winds, floating ice) in the abutting shelf seas likely influence the interior stratification while the air-sea interaction within fjords have only been observed to impact the near-surface fjord waters. Thus, we do not anticipate these factors to qualitatively change our findings.
In this subsection, a phenomenological description of the hydrography, circulation, and melt is presented for three major Greenlandic fjords (Ilulissat, Sermilik, and Kangerdlugssuaq; see locations in Fig. 1). In these regional case studies, we quantify the fjord overturning circulation via the overturning streamfunction, which is calculated as
Here, 𝑢 is the time-averaged velocity in the 𝑥-direction (and defined to be 0 below bathymetry),
𝑊 is the width of the domain in the 𝑦-direction, and 𝑧 𝐵 (𝑥, 𝑦) is the bathymetric elevation. To quantify the horizontal recirculation, we use the depth-integrated (for
Here, the horizontal quasi-streamfunction is approximately a streamfunction for the nearly nondivergent horizontal flow.
Ilulissat fjord in central-west Greenland has been discussed in many previous studies (e.g., Gladish et al. 2015b;Beaird et al. 2017;Khazendar et al. 2019). In Fig. 2, we present diagnostic Fig. 2h,i show the melt rate at the glacial face and the across fjord-averaged melt rate (comparing the discharge plume only and total melt), respectively. Although this fjord has one of the largest discharge rates in Greenland, over half of the total melt occurs outside of the discharge plume (and the other fjord locations, discussed below).
Due to the shallow sill at the fjord mouth, the sill overflow is hydraulically controlled i.e., the Froude number of the 𝜎 = 28.5 kg/m 3 density layer is approximately critical (not shown). As a result, the melt and subglacial plumes drive an overturning (primarily the deeper melt-plume driven overturning) that is limited to density classes up to this threshold. This lower warm water availability within the fjord due to the hydraulically-controlled sill overflow results in lower glacial melt rates.
However, a smaller range of density variation leads to a stronger overturning and recirculation strength for a fixed buoyancy flux forcing (acting on a weaker overall stratification), which leads to higher melt rates overall due to the stronger, primarily horizontal velocities at the ice face. See Zhao et al. (2019) for additional discussion on hydraulically-controlled fjord overturning and Pratt and Whitehead (2007) for background on hydraulically-controlled flows. Note in this simulation, the hydraulically-controlled sill overflow is not well represented using a hydrostatic model and future work using either a nonhydrostatic model or parameterized sill overflow turbulence (e.g., Jackson et al. 2008) will likely improve the representation of the sill-induced mixing and overturning reflux (Hager et al. 2022) at the fjord mouth.
Sermilik fjord in southeast Greenland has been discussed in multiple previous studies (e.g., Straneo et al. 2011;Straneo and Cenedese 2015). In Fig. 3, we present a series of diagnostic fields for Sermilik fjord and Helheim glacier melt rates similarly to Fig. 2. Table 1 lists the fjord parameters and bulk diagnostics from this simulation. This fjord domain has been rotated 81 degrees clockwise in this figure for ease in visualizing the overturning.
Finally, we discuss Kangerdlugssuaq fjord in central-east Greenland. In Fig. 4, we present a series of diagnostic fields for Kangerdlugssuaq fjord and the corresponding glacier melt rates similarly to Fig. 2. Table 1 lists the fjord parameters and bulk diagnostics from this simulation. The model domain has been rotated 45 degrees counterclockwise in this figure for ease in visualizing the overturning. This weakens the access of warm waters, which then weakens the overturning and recirculation below 𝑧 = -500. This is partly the reason this fjord has comparatively weaker melt near the grounding line (see Fig. 4g,h), with melt rates that peak at mid-depth.
surface in Fig. 5 (labeled Eddies A, B, and C)). These three eddies have different sizes and vertical vorticity profiles, but all emerge and begin spinning up around day 50 into the simulation around mid-depth (𝑧 = -200 to -400 m). We find that the existence of these eddies requires realistic fjord geometries (with multiple bowl-like depressions that can confine these eddies), an adequate horizontal resolution (a simulation at half this resolution was not adequate to produce these eddies)
and a favorable choice of horizontal viscosity parameterization.
The emergence of these eddies occurs as an apparent detachment of boundary currents that hug the sloping sidewalls of the fjord. This detachment first leads to a series of along-isobath gyre-like recirculation cells with alternating vorticity sign along the fjord (not shown but similar to the surface recirculation cells in Figs. 2a). This boundary current reorganization is primarily due to the nonuniformity of fjord width along the fjord. After detachment, the recirculation cells evolve to become more axi-symmetric over time and emerge as standing eddies. In Ilulissat fjord, the axisymmetrization or evolution towards more radially-symmetric circulation patterns lead to three mid-depth long-lived cyclonic eddies. Near the surface, the recirculation is much more variable in time and leads to a mix of cyclonic and anticycylonic eddies that advect out of the fjord. Following this early axisymmetrization, the eddies grow in their vertical extent, primarily downwards to the seafloor over a period of 150 days -in Eddy C, this coincides with the spinup of the along-glacier face velocity discussed below.
Fig. 5f illustrates the 3D cyclonic vorticity surface 𝜁/ 𝑓 = 0.75, which shows that in addition to these eddies, there is a positive vorticity source from the sill overflow region due to vorticity generation from water mass transformation, which is connected to and advects vorticity into Eddy A, and a positive vorticity source at the glacier face near Eddy C, which intermittently interacts with Eddy C, but primarily flows out of the fjord in the surface 125 m. Note that these eddies do not extend to the surface region where the vorticity field is dominated by the strong outflow interacting with bathymetry along the fjord sides (e.g., in Fig. 5a). However, there are occasional instances of vertical alignment of the outflowing near-surface eddies in the surface 150 m with the deep-water eddies that are the focus of this study (not shown). and circulation spin up, the pool of warm water in the deep fjord is consumed via mixing and export, which over time contributes to a lower melt rate. However, the circulation (overturning and recirculation) strengthens during this time, which compensates the cooler waters such that the overall melt rate does not change significantly -however, the total melt does increase slightly (by ≈ 10%) over this time. In Ilulissat fjord, the spinup process of the standing eddies and fjord recirculation takes approx. 3-6 months based on the near-glacier circulation and melt rate. The spinup time depends on the renewal rate of deep fjord waters, which depend on the fjord volume below a given depth and the overturning circulation (driven primarily by the buoyancy flux from the glacier at depth). However, there is also a spinup time associated with the tendency of the vorticity balance, which is discussed further in the next section.
Note that the melt and near-glacier velocity increases at depth over this time period, which is important for accurate predictions of rates of glacial undercutting (see Zhao 2021 for further discussion on the melt-circulation feedback in fjords). Due to this long duration, the temporal response/adjustment of the fjord circulation forcing variability on seasonal timescales (such as 1-2 month peak in summer subglacial freshwater discharge and seasonal winds) are important to take into account in an observational context because circulation may take months to spin up or spin down following a warm water renewal or summer subglacial discharge. the glacier face. The total integrated melt rates when the eddy is closest to the ice face (defined here as periods where eddy core location > 83 km) is approximately 10% higher than the remainder of the time series. This effect on the melt rate is primarily due to the eddy-induced velocity field rather than the eddy-influenced temperature field since the horizontal velocities also increase by approximately 10% during these time periods.
In Fig. 7, we show the azimuthally-averaged properties of the near-glacier eddy (Eddy C) at day 200. This eddy is located near the Ilulissat (formerly Jakobshavn) glacial face and is of particular importance due to its influence on the near-glacier velocity field and melt rate. Fig. 7a shows a snapshot of depth-averaged vorticity in the near-glacier region, which shows a radially-symmetric cyclonic eddy core centered at approx. 𝑥 = 82 km in addition to a positive vorticity region near the glacier face. In the panels of Fig. 7b-e, we observe that the azimuthally-averaged vorticity, azimuthal velocity, and isopycnal structure are consistent with an cyclonic submesoscale coherent vortex, as discussed in previous literature (e.g., McWilliams 1990). In particular, the stronger stratification at the core of the eddy drives a geostrophic circulation and a cold, fresh anomaly at the top of the eddy, which is consistent with the downward vertical velocity, and a warm, salty anomaly at the base of the eddy (in panels f through k). Both of these anomalies extend to but are much weaker at the glacier face in this snapshot.
Although eddy mergers do not have a significant effect on fjord overturning or heat transport (eddy momentum and heat transport terms are weak compared to the mean transport terms, which are not shown), they do have a significant effect on the maintenance of these standing eddies via vorticity advection, which is further discussed in Section 4.
In Fig. 8a-c, we show a Hovmöller diagram of the across fjord-averaged vorticity, which highlights the preferred locations of the standing eddies in Ilulissat fjord and merger activity at various depths.
Compared to Eddies A and C, Eddy B has a larger range and participates in more mergers (with C).
Eddies A and C have strong deep vorticity signatures while Eddies A and B also have near-surface vorticity signatures and exhibit more clear instances of vertical alignment between the surface and deeper eddies. For the deeper eddy dynamics (below 125 m, panels b and c), Eddy B periodically propagates eastward and merges with Eddy C, and we can see that the signal of Eddy B is weak in time periods following these propagation events. In addition, there is an eastward propagation of vorticity from the region near Eddy A to Eddy B, although this signature is partially due to the vorticity advection from the near sill region (see Fig. 5f). The eastward propagation velocities are approximately 1 km/day and are consistent with the across fjord-and depth-averaged horizontal mean flow.
For the surface eddies (above 125 m, panel a), Eddy B periodically propagates westward with a propagation velocity of 1.25 km/day and merges with Eddy A, while Eddy A occasionally exits westward (out of the domain) where the vorticity is advected/diffused by the strong exiting flow.
The periodical eddy mergers occur approximately every 30 days for the deeper eddies and it occurs approximately every 70 days for the surface eddies.
In the remainder of this study, we ignore the surface eddies, which are not standing eddies and do not strongly influence the fjord circulation, vorticity balance, or melt rates as much as the deeper standing eddies. In addition, these surface eddies are likely to be strongly influenced by the surface forcing (e.g., drag against sea ice/mélange), which is not represented in these model configurations.
In contrast, the deeper eddies are by comparison are much less likely to be influenced by the surface forcing.
Fig. 8d-f shows the 3D vorticity surface (𝜁/ 𝑓 = 0.75) at days 108, 112, and 116, highlighting an eddy merger event between Eddy B (cyan) and C (pink) and their merger (lime green). The maximum vorticity over the course of the merger is shown in Fig. 8g that results in a peak vorticity for Eddy C that is approx. 60% greater than the pre-merger peak vorticity. Note the existence of an outflowing surface eddy in these panels (centered at 𝑥 = 70 km, 𝑧 = -100 m).
Lagrangian tracer experiments have been informative in fjord and estuaries in both observations and numerical simulations (Pawlowicz et al. 2019). However, near glacier fjords, these simulations have only tested the fjord outflow on shelves instead of within the fjords themselves, e.g., in the West Antarctic Peninsula (Pinones et al. 2011) and Kangerdlugssuaq regions (Gelderloos et al. 2017). To better understand the influence of these eddies on tracers and residence times within these fjords, we use a tracer release experiment within the Ilulissat fjord simulation.
We deploy 100,000 tracer particles over a period of 20 days (5000 per day) starting at day 100 at even spacings in the x and z directions at the inflow boundary condition (50 tracers in the vertical direction and 100 tracer in the horizontal direction). The particle trajectories were not sensitive to the deployment rate because those that stay within the fjord spend a much longer period of time trapped within the fjord than the deployment duration. Based on our tracer release experiment, approximately 7.5% of shelf tracers enter the fjord (consistent with the Eulerian streamfunction). After 100 days, half of the tracers still remain in the fjord, where a majority of the tracers are near or within one of the three eddies. A key takeaway from this experiment is that the along-fjord and vertical deviations from the mean fjord residence time suggest that melt rates inferred from tracer concentration measurements within fjords such as noble gas and oxygen isotope measurements might be biased higher than actual melt rates. This is especially the case if they are within the near-glacier eddies where they may be trapped for a substantially longer time than the average Eulerian residence time.
In addition to the near-glacier Ilulissat Eddy C (discussed in subsection a), we also investigate and track the position of each eddy in all three fjords over time. To do so, we use the Angular Momentum for Eddy Detection and tracking Algorithm (AMEDA) (Vu et al. 2018). This detection and tracking algorithm has been effectively used in both numerical simulations and satellite data (Morvan et al. 2020). The algorithm uses gridded velocity, deformation radius, and a few userdefined tuning parameters as inputs and tracks individual eddy locations and radii over time. The algorithm uses this to find eddy centers that correspond to an extremum of the local normalized angular momentum that is contained within a closed streamline. The streamlines surrounding this center are then computed and the eddy mean radius is defined as the equivalent radius of a disc with the same area as one delimited by the closed streamline with the maximum area.
In our implementation of this algorithm, we use the depth-averaged horizontal velocity field below the discharge plume neutral buoyancy depth (approximately 𝑧 = -150 m for all three fjords), a deformation radius of 4 km that is approximately representative of all three fjord interiors, and the default tuning parameters from the algorithm (our results were not sensitive to these parameters).
In addition to the cyclonic eddies, the algorithm also detected smaller anticyclonic eddies, but these were short-lived eddies with weaker cores and are therefore not included in the discussion in the remaining sections. Fig. 10 shows the time-mean eddy radii (over days 200 to 300) and box plots of eddy locations compared to the width of the three fjords at selected depths. For all three fjords, there is a suggestive visual correlation of eddy locations with the widest parts of the fjords. Most of the eddy radii are contained within the 𝑧 = -500 m isobath and the eddies are fairly evenly spaced with no two preferred eddy locations within 10 km (a few eddy radii) of each other. A few of the eddy positions have a large horizontal extent (e.g., Ilulissat Eddy B), which seem to occur when the fjord width does not change much in the along-fjord direction. These same eddies also undergo eddy mergers more frequently (as discussed in Sect. 3b). Importantly, there are eddies (with varying properties) positioned near the glacier face in all three of the fjords tested, which has implications for glacial melt rate. However, the near-glacier eddies in the Sermilik and Kangerdlugssuaq fjords are much weaker than the one in Ilulissat because these fjords have a weaker subglacial discharge.
Recent modeling results show that the horizontal recirculation plays an important and potentially dominant role in glacial melt rates in deep-water fjords (Zhao 2021). Specifically, the near-glacier horizontal velocity, which owes its magnitude to the horizontal recirculation within the fjord, plays an important role in driving ambient front-wide glacial melt and may be comparable to the subglacial discharge-driven melt (Slater et al. 2018, Jackson et al. 2019). While the horizontal velocity can be locally complicated (e.g., within the standing eddies and near topographic features),
we take an overall view that focuses on the fjord-scale recirculation, for which a dynamical analysis of the circulation (i.e., the horizontally-integrated vorticity balance) is the appropriate diagnostic.
In order to develop scaling theories for the recirculation and near-glacier horizontal velocities, we present a vorticity balance analysis using Ilulissat fjord as an example. We then use this to develop a simple theory to predict the near-glacier horizontal velocities and glacial melt-rate, which builds on the theory from Zhao (2021).
To provide a theoretical scaling prediction for the near-glacier horizontal velocity, we first diagnose the fjord vorticity balance. We start with the horizontal inviscid momentum equations using the Boussinesq approximation on an 𝑓 -plane,
for bottom friction parameterized as F fric = 𝜕 𝑧 τ for a stress
where 𝐶 𝑑 = 2 × 10 -3 is the bottom drag coefficient. The stress in the interior includes viscous terms parameterized by the K-Profile Parameterization (Large et al. 1994), but this is negligible in our simulations. Note the difference in notation used here for the 2D velocity vector and Laplacian operator u ℎ = (𝑢, 𝑣), ∇ ℎ = (𝜕 𝑥 , 𝜕 𝑦 ), and the 3D velocity vector and Laplacian operator u = (𝑢, 𝑣, 𝑤), ∇ = (𝜕 𝑥 , 𝜕 𝑦 , 𝜕 𝑧 ).
Taking the horizontal curl of Eq. (3), we obtain the equation for the vertical component of vorticity
Fig. 11b-e shows the time-averaged across fjord-integrated along fjord-cumulative integrals of each term in Eq. ( 5) vertically-integrated over four distinct depth bands. We use a along fjordcumulative integral, which starts at zero at the glacial face and is integrated westward, to visually reduce noisy vorticity sources generated by bathymetry. The terms in this balance shown in Fig. 11 include the vorticity generation, bottom stress curl, and total vorticity advection (including horizontal and vertical vorticity advection, vorticity stretching, and tilting). The vortex tilting terms are comparatively much smaller than the other terms and grouped with this total vorticity advection term. The four depth bands are chosen based on the sign of 𝜕 𝑧 𝑤 (which is consistent with the two overturning cells in Fig. 2b) and correspond to the inflow/outflow of the melt plume-driven overturning (𝑧 < -380 m, -380 < 𝑧 < -235 m) and the inflow/outflow of the discharge plumedriven overturning (-235 < 𝑧 < -125 m, -125 < 𝑧 m), which is more apparent in the near-glacier area integral discussed below.
In the bottom layer, we start at the glacier face and integrate the terms in Eq. ( 5) westward cumulatively. In the near-glacier region, vorticity generation is primarily balanced by the total vorticity advection term, which is dominated by boundary currents. The along-fjord changes in the total advection are primarily due to horizontal vorticity advection (as a result of vorticity generation near rough bathymetry) as well as vortex stretching. Throughout the fjord interior, vorticity advection balances vorticity generation, but for the cumulative integral over the entire fjord (𝑥 > 35 km), the vorticity generation is primarily balanced by the integrated bottom stress curl. We use the approximate balance between the vorticity generation and bottom stress curl in the bottom layer to develop our scaling prediction for the melt rate in the next subsection.
In the other three layers, the vorticity generation is predominantly balanced by total vorticity advection over most of the fjord, which is primarily due to horizontal vorticity advection within the boundary current. See Zhao et al. (2022) for a similar fjord vorticity balance with a simpler bathtub bathymetry.
We now calculate vertical profiles of the terms in Eq. ( 5) over a defined region of the fjord by taking an area-integral and applying the Stokes' and divergence theorems
where n is the unit vector normal to and t is the unit vector tangential to the boundary of area 𝐴 and the vortex tilting refers to the last component of the total vorticity advection term in Eq. ( 5).
Fig. 11e shows the vertical profiles of each of the terms in Eq. ( 6) integrated over two regions:
the near-glacier recirculation area (bounded by 𝑥 = 75 km and the glacial face) in the solid lines and the entire fjord region excluding the sill (bounded by 𝑥 = 35 km and the glacial face) in the dotted lines. Note that over the near-glacier region, the vorticity generation in the bottom layer is balanced by total vorticity advection while bottom stress curl is weak. However, the bottom stress curl is the dominant term in the bottom 200 m when integrated over the full fjord domain.
In order to understand the sensitivity of glacial melt rates to fjord circulation, we extend previous theories (Zhao et al. 2021(Zhao et al. , 2022) ) to relate fjord vorticity balance to glacial melt. These previous theories showed that the vertically-integrated melt rate could be accurately predicted in a similar model setup with a simplified bathtub fjord geometry. These results were validated over a number of key geometric and forcing parameter dependencies. The discussion here is therefore framed around how these previous theories can explain the vertical structure of circulation and melt rate when realistic bathymetry is used instead.
Based on the vorticity balance in the bottom layer in Eq. ( 5) and Fig. 11, we calculate a prediction for the along-isobath velocity by balancing the area-integrated bottom stress curl and vorticity generation terms from Eq. ( 6) following on Zhao et al. (2022). An underlying assumption is that this along-isobath velocity is assumed to be approximately constant around the isobath, which includes the near-glacier region as part of its circuit. This assumption allows us to predict vertical profiles of horizontal velocity at the glacier face. The vorticity generation term is defined as 𝑓 𝑄(𝑧), where 𝑄(𝑧) ≡ ∬ 𝑤 𝑧 d𝐴 .
The bottom stress curl term
can be approximated by using a mean along-isobath velocity, 𝑣 bdy , which approximates the path integral with a boundary perimeter length 𝒞 circumscribing region A using an along-path averaged velocity scale 𝑣 bdy = 𝒞 -1 ∮ 𝜕 𝐴 u ℎ • t d𝑠 within a bottom boundary layer scale height, 𝐻 eff . In the simplified boundary layer parameterization using a bulk drag coefficient (in MITgcm with an unresolved bottom boundary layer), 𝐶 𝜏 ≡ 𝐶 𝑑 (Δ𝑧) -1 and 𝑣 bdy is evaluated at the deepest wet grid cell just above bathymetry. This is used as the prediction of the mean along-perimeter horizontal velocity at each depth in the discretized vertical grid, which we denote as 𝑣 theory .
Setting the terms from Eqs. ( 7) and ( 8) equal, we have a prediction for the near-glacier velocity
We compare the prediction of the vertical profile of horizontal velocity to the simulated results of Ilulissat fjord. Fig. 12a,b show the near glacier circulation, which is consistent with the near-glacier eddy discussed in previous sections. However, the along-face horizontal velocity in Fig. 12b,c exhibit complex across-fjord and vertical structure. Fig. 12d shows the depth-averaged velocity based on Eq. ( 9), compared to the across fjord-averaged along-glacier velocities diagnosed from simulations. The differences between the theoretical depth-averaged along-glacier velocity and simulated velocity profiles are substantial at most depths. Specifically, in the bottom layer, the small recirculation region in the southeast corner leads to a southward along-glacier velocity, which is not captured in our theory. This negative (clockwise) recirculation region accounts for the difference between our theory over -620 < 𝑧 < -400 m. Below these depths, the specific pathways of currents guided by bottom bathymetry dominates the near-glacier velocity. However, in the bottom layer the theoretical prediction of the depth-averaged along-glacier velocity is at most 10% larger than the simulation-diagnosed value. This prediction is not expected to be accurate in the other depth bands because the vorticity generation is balanced by the total vorticity advection term instead of bottom stress curl.
Using this prediction of the near-glacier velocity magnitude, we develop predictions for the glacial melt rate. Assuming that the melt is primarily driven by horizontal velocities external to the discharge plume and vertical velocities within the discharge plume, we use the 3-equation thermodynamics (using e.g., Hellmer andOlbers 1989, Holland andJenkins 1999) and assume ice temperatures that are approximately freezing. This allows us to simplify this relationship to a linear melt rate 𝑀 (in m/s) that is approximately proportional to 𝑣 theory for the melt rate external to the plume as
where 𝐿 𝑖 = 3.35 × 10 5 J kg -1 is the latent heat of fusion of ice, 𝑐 𝑤 = 3.974 × 10 3 J kg -1 K -1 is the specific heat capacity of water, Γ 𝑇 = 2.2 × 10 -2 is the thermal transfer constant, and 𝑇 𝑏 = 𝑇 𝑓 and 𝑇 𝑎 are the boundary layer (assumed to be at freezing temperature) and ambient temperature, respectively. The vertical ambient temperature profile is diagnosed from the model over near-glacier region, 𝑥 > 85 km.
The total melt rate is the sum of the region outside of the plume (from Eq. ( 10)) and the region within the plume
where the vertical velocities are predicted used plume dynamics (Morton et al. 1956), which is used in the parameterization of plume entrainment in our model (Cowton et al. 2015) and can be diagnosed directly (or explicitly included in the theory).
The melt rate predictions, 𝑀 v,theory and 𝑀 total,theory , are shown in Fig. 12e compared to the simulation-diagnosed across fjord-averaged melt rate. Within the bottom layer, the melt rate is well predicted by the theory in Eq. ( 11), with the vertical and horizontal velocity components each driving approximately half of the melt rate in both the theory and simulations. This improved accuracy compared to the near-glacier velocity is due to the fact that the temperature profile is diagnosed from the model results. Note that this theory only accounts for thermal contribution to the elevated melt rates near the grounding line (leading to glacial undercutting) and not the increased along-glacier velocity at the grounding line since it uses a layer-averaged velocity. However, this only demonstrates a partial success of the theory as shown in Fig. 12d,e. For the profiles shapes in Fig. 12d, the theory for the third layer provides an accurate prediction of near-glacier velocity, and in Fig. 12e the second layer is an inaccurate prediction for melt rate due to advection playing a significant role in the vorticity balance in Fig. 11. On the other hand, the general magnitudes of velocity and melt rate are still reasonable.
In summary, we find that vertical profiles of velocity and melt at the glacial face are dependent on complex bathymetric features and restricts/guides access of dense warm water near the grounding line. The melt rate is strongly influenced by a complicated flow that is not easily captured in our simple theories. However, a simple vorticity balance and melt rate prediction can aid with the interpretation of realistic simulations to better understand how circulation drives glacial melt.
In this study, we use a high-resolution numerical model with realistic geometry to simulate the circulation within three major Greenlandic fjords (Section 2, Figs. 234). These simulation results reveal multiple standing eddies in each fjord (see Fig. 5). We discuss the properties of these eddies and their role in fjord circulation and tracer advection, and the role of bathymetry in determining their preferred locations (Section 3,Figs. 7,9,10). To understand the influence of eddies and the resulting circulation within realistic fjord geometries, we analyze the fjord vorticity balance, which allows us to extend previous theories for the glacial melt rate (Section 4).
We find that eddies within glacial fjords are generated by a combination of two vorticity sources (see Fig. 11): (1) vorticity generation from the subglacial discharge and meltwater plumes at the glacial face and (2) vorticity advection into the fjord through the fjord mouth at mid-depths.
These eddies take months to spin up (Fig. 6) and eventually reach a steady state with bottom stress curl balancing/dissipating the vorticity input at depth. The eddies are large perturbations on the horizontal streamfunction within the fjord, undergo mergers with other eddies (Fig. 8), and significantly increase the Lagrangian and Eulerian residence times within the fjord (see Fig. 9).
These eddies prefer deep and wide regions within the fjords (see Fig. 10).
Most important, when these eddies are close to the ice face, they serve to amplify glacial melt rates (as shown in Fig. 12). To develop a scaling theory for the vertical melt rate, we used the vorticity balance in Section 4a. In this balance, the eddies play an important role in the vorticity budget at deeper depths, where the near-eddy regions dissipate the majority of the vorticity through bottom stress curl. Although the near-eddy circulation is not entirely set by the eddy dynamics, the aggregation of most eddies near the deepest and widest portions of the fjord occurs where the along-isobath circulation dissipates most of the vorticity beneath these eddies. Our results show the existence of long-lived eddies for all three fjords studied. However, standing eddies are also likely to exist in many of the other deep-water fjords around Greenland.
Observations have likely missed these dynamically-significant features due to their small scale and temperature/salinity anomalies being less apparent (but still observable) compared to anomalies in the less-sampled velocity and vorticity fields because the eddies appear to exist primarily in the well-mixed fjord interior. In addition, although these are standing eddies, they do move periodically over distances larger than their radii, making them difficult to observe.
In real fjord systems, these eddies may take months to spin up or spin down following a warm water renewal or summer subglacial discharge. It is possible then that these eddies do not fully spin up during a melt season. However, although the timescale associated with a full eddy spinup process is longer than the melt season (as Zhao 2021 suggests) the subglacial discharge-driven circulation acts to trigger the melt-circulation feedback. This feedback then dominates the deep circulation (below the neutral buoyancy depth of the subglacial discharge plume) and lasts for a much longer time period. In the absence of subglacial discharge, the melt-driven circulation would take longer to spin up the eddies, but they would likely still exist.
There are numerous caveats in this study due to the limitations of our model configuration.
These include the absence of atmospheric fluxes, simple vertical mixing representation, the lack of sea ice, mélange, and icebergs, which can supply substantial buoyancy input (Enderlin et al. 2016). Another caveat is the prescription of a time-invariant open-ocean boundary, which limits the shelf variability within our simulations; there can be a shelf current-induced increase/ decrease in the exchange between the fjord and shelf (Zhao et al. 2021). We also do not consider the effect of winds, which likely exhibits a larger effect on the shelf region via fjord overturning driven by coastal upwelling (not included in our domain), but may also directly drive fjord circulation/renewal for strong enough katabatic wind events (Zhao et al. 2021;Spall et al. 2017). Also, in our glacial boundary parameterization, the melt rates are calculated using the closest grid point of horizontal and vertical fjord velocities, which is an imperfect representation; in general, a better understanding and representation of the ice-ocean boundary layer needed to improve glacial melt rate estimates. Another caveat is that much of the analysis in this study is specific to Illulisat, which can be extended to more fjords in the future.
Following this study, there are a number of open questions that require further attention. Additional work is needed to investigate other fjords at higher resolution as well as conduct this analysis over a larger sample of Greenlandic fjords. In particular, analyzing the measures of 3D circulation, vorticity balance, and melt rate for more fjords may help us understand the range of circulation-melt interaction across the fjord population. Another future avenue is to investigate boundary layer parameterizations at the glacial face, which is currently not well-supported by observations (Jackson et al. 2019). Improved representations of the boundary layer may also influence the interaction of submesoscale-microscale dynamics. A final avenue is to investigate the interaction between multiple neighboring fjords, which is likely important when many fjords are closely packed along the shelf.