Decadal Evolution of Ice‐Ocean Interactions at a Large East Greenland Glacier Resolved at Fjord Scale With Downscaled Ocean Models and Observations

In recent decades, the Greenland ice sheet has been losing mass through glacier retreat and ice flow acceleration. This mass loss is linked with variations in submarine melt, yet existing ocean models are either coarse global simulations focused on decadal‐scale variability or fine‐scale simulations for process‐based investigations. Here, we unite these scales with a framework to downscale from a global state estimate (15 km) into a regional model (3 km) that resolves circulation on the continental shelf. We further downscale into a fjord‐scale model (500 m) that resolves circulation inside fjords and quantifies melt. We demonstrate this approach in Scoresby Sund, East Greenland, and find that interannual variations in submarine melt at Daugaard‐Jensen glacier induced by ocean temperature variability are consistent with the decadal changes in glacier ice dynamics. This study provides a framework by which coarse‐resolution models can be refined to quantify glacier submarine melt for future ice sheet projections.


Introduction
In recent decades, the Greenland ice sheet has been losing mass and is presently contributing approximately 0.5 mm/yr to global sea-level rise (Mankoff et al., 2021).This mass loss is approximately evenly split between enhanced melt on the surface of the ice sheet that runs off to the ocean and increased discharge of ice from Greenland's marine-terminating glaciers.The latter process, termed dynamic mass loss, is instigated by submarine glacier melt which causes glaciers to retreat (Wood et al., 2021) and subsequently accelerate (King et al., 2020), increasing the flux of ice into the ocean (Mankoff et al., 2021).The magnitude of melting depends on the volume and distribution of subglacial discharge at glacier termini, and on the thermal forcing of ocean water (Slater et al., 2015;Xu et al., 2013).Thus, the decadal evolution of submarine melt rates, driven by a combination of variations in regional air and ocean temperatures (Slater & Straneo, 2022) is directly linked to sea-level rise contributions from the ice sheet.
Despite recent advances in our knowledge of the factors contributing to Greenland's ice loss, ice models predicting sea-level rise contributions from the ice sheet systematically underestimate historical ice loss rates, casting doubt on the future projections (Aschwanden et al., 2021;Goelzer et al., 2020).These model-observation differences are largely attributed to uncertainties in the magnitude and variability of submarine melt rates and the subsequent processes that induce dynamic mass loss (Fox-Kemper et al., 2021)-in effect, the oceanic boundary conditions that drive dynamic change in ice sheet models.One of the main sources of this uncertainty arises from a lack of observations and model simulations of ocean properties within Greenland's fjords.From an observational perspective, there are few in situ measurements near Greenland's glaciers, and the majority of those which exist are from recent years after dynamic mass loss had initiated (Mouginot et al., 2019).From an ocean modeling standpoint, most ocean models that resolve fjord circulation and melt are constructed for idealized experiments, are typically forced at their boundary by a single observation, and are run for a period of a few years or less (e.g., Slater et al., 2015;Xu et al., 2013;Zhao et al., 2021).To our knowledge, there are no models, idealized or otherwise, which have been constructed to investigate temperature variability within the fjord network on a decadal scale.As a result, the glaciological modeling community has been left to estimate submarine melt on Greenland's ice from coarse resolution models and sparse ocean observations outside of the fjords.This limits our ability to understand how variations in submarine melt have induced glacier change over the past several decades and it hinders our ability to estimate how glaciers will continue to lose ice and contribute to global sea-level rise in the coming decades.
To address this gap, we present a novel ocean model framework to resolve circulation, temperature variability, and ice-ocean interactions at the fjord scale around Greenland.The numerical ocean model framework links a global-scale ocean-sea ice state estimate with two downscaled, regional-scale ocean model configurations that increase in resolution but decrease in spatial extent to make efficient use of computational time and resources.Using this approach, we leverage decades of remotely-sensed and in situ ocean observations via the global state estimate while also resolving the fjord-scale circulation that modulates submarine melt rates.The framework is released as open-source code and designed to be adaptable to other parts of the polar oceans to investigate iceocean interactions and other local-scale processes which are not resolved by global-or regional-scale dataassimilative models.
We illustrate our new framework on the Scoresby Sund fjord network-the longest fjord system in Greenland.Several large glaciers terminate in this region including Daugaard-Jensen Glacier (DJG, Figure 1), which accounts for approximately 9 Gt per year of ice discharge from the ice sheet (Mouginot et al., 2019).For more than two decades prior to 2010, DJG was dynamically stable: its ice velocity remained within 5% of its annual mean flow speed of 2,500 m/yr, and the ice front location varied seasonally within ±1 km of its annual mean position (Bevan et al., 2012).In 2010, the glacier underwent an abrupt change, starting a multi-kilometer retreat of the ice front and accelerating by more than 10% (Mouginot et al., 2019;Wood et al., 2021).The acceleration and retreat continued until 2018 when it began to re-advance and slow down.The onset of these dynamic changes on DJG have been hypothesized to be linked with an increase in subsurface ocean melt, which caused the glacier to retreat (Wood et al., 2021), yielding subsequent acceleration and thinning (King et al., 2020).We test this hypothesis using our numerical framework to simulate the variability of submarine melt on DJG during the progression of these observed dynamic glacier changes.

Numerical Ocean Modeling
We use the MIT General Circulation Model (MITgcm, Marshall et al., 1997) to simulate circulation around and through the Scoresby Sund fjord network via three models of increasing resolution and decreasing spatial extent (L0, L1, and L2).The domains of the models are shown in Figure 1.Details of the downscaled models are described here and summarized in Table S1 in Supporting Information S1.The methodology to generate downscaled model fields for initial, boundary, and external forcing conditions is described in Text S1 in Supporting Information S1.All model code required to rerun and reproduce the model results is provided in an opensource repository.
For our L0 model, we begin with an existing global coupled ocean-sea ice model ("LLC270" a.k.a."ECCOv5 Alpha") with a nominal 1/3 degree resolution constructed as part of the ECCO consortium (Zhang et al., 2018b).For convenience, we refer to this model as "L0" in this manuscript.This model was constructed using the adjoint of the MITgcm to minimize the misfit between the model solution and observations, resulting in a physically-consistent simulation which assimilates available observations.The model is constructed on the "Latitude-Longitude-Cap" (LLC) grid with 270 grid cells on each of 13 model faces that cover the globe, resulting in a horizontal resolution of approximately 15 km near Scoresby Sund.The model time domain over which the model fit to observations was achieved spans from 1992 through 2017.We extend the time duration of this model to the end of 2021 using external forcing conditions from recent atmospheric reanalysis, but it should be noted that the years 2018-2021 are simulated with a forward extension which does not assimilate observations.In L0, the runoff fields are climatological (Fekete et al., 2002).All other external forcing conditions are derived from the adjusted 6-hourly ERA-interim reanalysis product (Dee et al., 2011).Note that a spin-up is not required for the L0 model because it employs the MITgcm adjoint to assimilate observations, which operates backwards in time.
The first downscaled model (L1) of our study is constructed on a subset of the "LLC1080" grid with a four-fold higher resolution (Table S1 in Supporting Information S1).This configuration has a nominal resolution of ∼3 km near Scoresby Sund which is the approximate length of the mode 1 and mode 2 Rossby radii of deformation in this region of the sub-Arctic ocean (Nurser & Bacon, 2014); further off the shelf, the grid resolution increases to ∼4 km.The model is similar to the L0 model except for two key differences.First, since the resolution of the L1 model is eddy-permitting, temperature, salinity, and other tracers are advected and diffused in the model by eddies that form during the simulation.This contrasts with the L0 model which is too coarse to resolve eddies; instead, the advection and diffusion of tracers by eddies in L0 is parameterized using the "GMRedi" scheme (Gent et al., 1995;Redi, 1982).Second, as a regional model, L1 has boundary conditions provided by output from the L0 model.Boundary conditions are prescribed hourly for potential temperature, salinity, and velocity, as well as sea ice parameters for area, thickness, snow thickness, and ice velocity.In addition to the prescribed values, a "sponge" layer is added that relaxes values within 10 grid cells of the boundary toward the prescribed values to compensate for inconsistencies with the model solution inside the domain.External forcing conditions are identical to those in L0 and updated every 6 hr.The L1 model is initialized 1 Jan 1992 using the interpolated instantaneous ocean state from the L0 model.
The second downscaled model (L2) has a resolution of 500 m and is similar to the L1 model except for three key differences.First, additional vertical layers are added to enhance the vertical resolution-the surface grid cell is 1 m thick, telescoping to thicker cells at depth.Boundary conditions are balanced a priori so that no ocean volume is lost or accumulated in the domain as a result of the change in the model grid.Second, subglacial discharge, which is emitted at depth from underneath glacier termini, along with the resulting submarine melt along the glacier face, is added using the MITgcm iceplume package (Cowton et al., 2015).Subglacial discharge rates and locations are generated via hydrologic routing of runoff from the MAR regional climate model (Mankoff et al., 2020) and melt rates are estimated with recent updates in the melt parameterization for vertical fronts (Schulz et al., 2022).The melt plume is modeled as a line plume with a width of 500 m.To compensate for a computationally-necessary simplification in the melt rate parameterization, that is, not calculating the temperature and salinity properties of the small plumes formed by ambient melting, we modify the temperature applied in the thermal balance at the ice-ocean interface as a linear function of fjord temperature.Details on this modification can be found in Text S2 in Supporting Information S1.Third, the runoff fields in L2 are taken from Mankoff et al. (2020) and are distributed into the closest wet model grid cell to the discharge location.Mass which enters the model domain through subglacial discharge or the surface is removed via an adjustment calculated online during the model run and applied to all wet boundary grid cells.

In Situ Ocean Observations
We compare our model results with a collection of temperature and salinity profiles in the Scoresby Sund fjord network and in the regional ocean on and off of the shelf.The majority of measurements we use for comparison were accessed via the EN4.2.2 database (Good et al., 2013).This data set includes individual measurements from the global Argo array, PI-led cruises, instrumented pinnipeds, and other instruments.These measurements are interpolated in space and time to generate a 1°gridded objective analysis product, which we utilize off the shelf to assess regional averages.
Within the fjord network, we use additional conductivity-temperature-depth (CTD) measurements from a few PIled efforts that have not yet been included in the EN4.2.2 data set.In 1990, a cruise on the RV Polarstern took CTD measurements at several locations within the fjord including in close proximity to DJG (Dowdeswell et al., 1992).To our knowledge, these measurements were not repeated until 2016 when NASA's Oceans Melting Greenland (OMG, Fenty et al., 2016) mission surveyed the fjord each year 2016-2021.

Remotely-Sensed Observations of Ice Dynamics
To reconstruct the dynamic evolution of DJG, we quantify changes in glacier ice velocity, elevation, and front position.For ice velocity, we rely on annual data from NASA MEaSUREs products including the ITS_LIVE mosaics (Gardner et al., 2019) and mosaic compilations derived from Interferometric Synthetic Aperature Radar data (Joughin, 2020;Joughin et al., 2010).We denote these products as MEaSUREs-ITS_LIVE and MEaSURES-InSAR, respectively.We sample velocity by averaging available data within a 500 m radius around 71.875°N, 28.682°W, located approximately 5 km from the mean terminus position of DJG prior to the onset of retreat.For ice elevation, we use a variety of data sets from airborne and satellite platforms.For historical data prior to 2010, we rely on reconstructions of ice elevation from ortho-rectified historical photographs taken in Scoresby Sund in 1985 (Korsgaard et al., 2016) and the GIMP DEM (Howat et al., 2014) for 2005.From 2010 to 2018, we use digital elevation models from the ArcticDEM dataset constructed from orthorectified pairs of WorldView images (Porter, 2018) and for 2016-2019, we use ice elevation measurements from the GLISTIN-A radar collected during NASA's OMG campaign (OMG, 2020).Finally, for years 2018-2021, we use laser altimetry observations from the ICESat-2 satellite (Smith et al., 2021).We sample the elevation data in the same manner as velocity.For ice front positions, we use a compilation of manually-digitized front positions provided in the TermPicks data set (Goliber et al., 2022).The digitized ice fronts in this data set have a median error of approximately 100 m.

Results
Here, we present the results of our model simulations by comparing simulated and observed in situ ocean properties.We work inward toward Scoresby Sund, starting from the regional ocean in the off-shelf region of the L1 model (Section 3.1); then the continental shelf outside Scoresby Sund (Section 3.2); and finally within the fjord itself (Section 3.3).We then describe our modeled melt rates (Section 3.4) and the reconstruction of glacier dynamics at DJG (Section 3.5).

Temperature in the Regional Ocean
We assess our model results in the regional ocean by comparing the L0 and L1 model (Figure 2a).Since the aim of this study is on melt rates on glaciers in Scoresby Sund, we focus on depths between 200 m, the approximate upper-bound on the layer of warm, salty, Atlantic-sourced water, and 500 m, which is the approximate depth of the continental shelf near the entrance to the sound.We sample the models off the continental shelf in a regional polygon (Figure 1b) which was used for the creation of ocean boundary conditions for ice sheet models (Slater et al., 2020) in recent projections of sea-level rise (Goelzer et al., 2020).This polygon encompasses the East Greenland Current and the upstream source waters for Scoresby Sund.In both the L0 and L1 model simulations, we find regional ocean temperature remained between 0.5 and 1.3°C in the years 1992-2011 with warmer temperature in the period 1995-2000 and in 2005 (Figure 2a).After 2011, temperatures increased in both models by nearly 1°C through the end of the simulation in 2021.For the majority of the simulation, the L1 model is colder than the L0 by 0.2-0.4°Calthough both models start at the same temperature of 0.8°C in 1992 and end at approximately the same temperature around 1.8°C in 2021 for an overall warming of +1°C.Observations in this

Geophysical Research Letters
10.1029/2023GL107983 region of the global ocean, even off the shelf, are quite sparse so we rely on the Hadley EN4.2.2 Objective Analyses (OA) product.Relative to the OA estimate, the L0 and L1 models had overall biases of +0.6 and +0.5°C , respectively, although we note that there is significant noise in the OA product owing to the sparsity of observations in this region (Figure S1 in Supporting Information S1).

Temperature on the Continental Shelf
On the continental shelf where all three models overlap (Figure 2b), we find a similar progression of temperature as in the regional ocean with temperatures increasing between 2011 and 2021.However, we observe a marked difference between the coarse L0 model and the finer-resolution L1 and L2 models.First, the L0 model exhibits stronger seasonal variability compared with L1 and L2.Second, between 2005 and 2010 the higher resolution simulations feature an increase in temperature by 0.4°C while the L0 model declines.Third, while temperature increases in all three models between 2011 and 2020, the L0 model exhibits a larger increase compared to the higher resolution models, but we note that the temperatures are approximately equal by the end of the simulation in 2021.Since the L0 and L1 models are subject to the same external forcing from the adjusted ERA-interim reanalysis, we attribute these model differences to the mechanisms by which heat is advected onto the continental shelf.In L0, the heat advection onto the shelf is parameterized, leading to a larger magnitude of both seasonal and interannual change.In L1 and L2, heat is advected via eddies that are resolved by the simulation.These model comparisons highlight the role that eddies play in heat transport across the continental shelf.
To compare with observations, we calculate model errors relative to CTD measurements from the OMG campaign.We find that the L0, L1, and L2 models had an overall bias of 0.2°C, 0.42°C and 0.40°C relative to observations (Figure S1 in Supporting Information S1).For L1 in particular, this contrasts with the warm bias of L1 relative to the Hadley OA product, suggesting this product may have a regional bias relative to the OMG measurements.

Temperature Inside Scoresby Sund
The Scoresby Sund fjord network features a large embayment approximately 40 km in width and 130-180 km in length from the mouth of the sound.From the end of the embayment, there is a network of three long fjords where several glaciers terminate.The northwestern fjord-Nordvestfjord, the longest fjord in Greenland-is approximately 150 km in length and terminates at the face of DJG.Since the focus of this manuscript is on ice-ocean interactions at DJG, we assess temperature at a location mid-sound near the mouth of Nordvestfjord (24.888°W , 71.094°N, Figure 2c) and a location near the front of DJG (28.350°W, 71.938°N, Figure 2d).At the location mid-sound, our simulated temperature shows a similar progression to that on the shelf: an increase in temperature of approximately 0.3°C, a decline through 2012, and then a subsequent increase from 2014 through the end of the model run.We note that the magnitude of variations and warming trend in the L1 model at the mid-sound point are smaller compared with the L2 model, owing to its resolution which is 6 times coarser in Scoresby Sund.The midsound progression of temperature is also observed at the near-glacier point nearly 150 km upstream, although the magnitudes are subdued: there is a 0.2°C increase between 2005 and 2010 and then a further increase' 2 years later in 2017.In comparison with available observations the L1 and L2 models have overall (cold) biases of 0.76°C and 0.62°C, respectively, relative to OMG CTDs at the mid-sound location; at the near-glacier site, which is only resolved by the L2 model, there is an overall cold bias of 0.69°C (Figure S1 in Supporting Information S1).
To investigate variations of heat transport through Scoresby Sund and Nordvestfjord, we examine the relative depths of the ρ PA = 1,028.5kg/m 3 isopycnal (the approximate level delineating Polar water vs. Atlantic Water (AW) on the East Greenland Shelf) and the relative temperature at ρ AW = 1,029 kg/m 3 , which is within the AW layer (Figure S2 in Supporting Information S1).In the periods leading up to the 2005, 2010, and 2015 warm temperature anomalies near DJG, two conditions in the fjord contributed to enhanced heat transport to the glacier: (a) ρ PA was anomalously shallow by approximately 20 m and (b) the temperature at ρ AW was anomalously high by approximately +0.2°C.In effect, a larger proportion of dense AW was advected into the sound and the temperature of the dense water is anomalously warm.The Hovmöller diagram (Figure S2 in Supporting Information S1) for temperature anomalies through the sound and fjord reveals that the 6-7 km mouth of Nordvestfjord acts as a "filter" for heat transport, subduing temperature anomalies from the wider Scoresby Sund and leading to smaller variations in temperature near the glacier relative to those in the sound.

Model Melt Rates on Daugaard-Jensen Glacier
The modeled melt rate timeseries at the terminus of DJG has two notable components reflecting the atmospheric and oceanic processes that dictate the melt rate magnitude (Figure 3a).First, the timeseries has a distinct seasonal cycle, increasing 1.5-2 m/day above ambient background rates in mid-summer.The timing of this increase results from the increase in subglacial discharge from inland ice surface melt.The highest peaks are found in 2010, 2012, and 2019 when summertime ice sheet surface melt was anomalously high.The seasonal cycle is overlain on a background melt rate whose magnitude is determined by the near-ice fjord temperature in the model run (Figure 2d).The background melt rate reaches peaks of approximately 4.4 m/day in 2005, 4.6 m/day in 2010, and increases to 4.7 m/day by 2020 near the end of the model simulation.The elevated subglacial discharge rates in 2010, 2012, and 2019 occur simultaneously with elevated fjord temperatures so that the largest overall melt rates occurred in those years.

Glacier Dynamics Reconstruction
The reconstruction of glacier dynamics at DJG reveals a period of dynamic stability from 1985 until around 2005 at which point the ice started to accelerate (Figure 3c).Prior to 2005, the annual mean ice velocity was relatively stable around 2,500 m/yr.Subsequently, DJG accelerated by 200 m/yr 2 between 2005 and 2010, and then another 200 m/yr 2 between 2012 and 2020, peaking around 2,800 m/yr (∼12% increase in speed).For elevation, a combination of data sources show a thinning rate of more than 2.5 m/yr over a 10 years span during 2010-2020 (Figure 3d).Due to the sparsity of historical elevation data prior to 2010, it is difficult to determine the progression of thinning.However, Airborne Topographic Mapper data from NASA's Operation Ice Bridge measured a surface elevation in 2010 similar to that derived from airborne ortho-photographs in 1987 while the GIMP DEM, with a nominal date between 2002 and 2009, estimates a surface 10 m lower.For ice front retreat, the progression is similar-the ice front position varied seasonally around an approximately constant location until around 2011 when the mean front position began to retreat inland by approximately 2 km (Figure 3e).Overall, the years 2010-2011 marked a dynamic inflection point in which DJG began to accelerate, retreat, and thin, losing mass dynamically.

Discussion
In this study, we leveraged one global and two downscaled models to simulate circulation in the regional ocean and within the Scoresby Sund fjord network.This approach allows for the resolution of circulation within the fjord network while also leveraging the data-assimilative capabilities of the coarse L0 model-a model which would be too computationally expensive to run at the timescale and resolution required to resolve fjord-scale circulation.The resolution of the downscaled models reflects the scales for heat transport from the regional ocean to the sound and the near-glacier environment.The L1 model, with a resolution of approximately 3 km at the shelf break, permits the resolution of small eddies and associated heat transport onto the shelf, which is not well-resolved by the L0 model.The L2 model, with a resolution of 500 m, permits the resolution of some submesoscale dynamical features as well as circulation inside the sound, into the narrow mouth of Nordvestfjord, and near the glacier where melt occurs, which is not well resolved in the L1 model.At the mid-sound location, the temperature variations in the L1 model are damped relative to those in the higher-resolution L2 model; near the glacier, the L1 model does not resolve any temporal variability at all.One of the most encouraging results from our downscaled model simulations is the comparison of the timing and relative magnitude of modeled melt rates with the progression of glacier dynamics.The preliminary increase in simulated melt rate around 2005 is approximately coincident with an initial speed-up of the ice front.These results are consistent with a range of previous studies that demonstrate that variations in ice velocity are triggered by changes that occur at the glacier terminus (e.g., King et al., 2020;Nick et al., 2009).The second increase in melt rate around 2010 is coincident with ice front retreat and thinning, and subsequently acceleration in 2013, which are the three main signals of dynamic glacier mass loss.We note that these interannual variations in melt are not observed in two existing estimates of ice front melt from Wood et al. (2021) and Slater and Straneo (2022) which both sample ocean properties outside of the fjords from coarse-resolution models or observational products (Figure S3 in Supporting Information S1).
To estimate the effect of our modeled melt variations on glacier retreat, we calculate the integrated anomaly in ice front melt relative to the first 5 years of the L2 simulation (2000)(2001)(2002)(2003)(2004)-a metric of how much extra ice was melted as a result of increased ocean melt.We find that the ablation anomaly totals more than 1 km by the end of the simulation in 2021, which is of similar magnitude as the 2 km of observed retreat.While the approximate timing and magnitude of the integrated anomaly is expected to be similar to retreat, the two values are not expected to match perfectly.On one hand, destabilization via enhanced ocean melt induces acceleration and increased calving rates through the melt-induced calving multiplier effect (Slater et al., 2021), resulting in retreat that exceeds total melt by the ocean.Second, there is uncertainty inherent in the absolute magnitude of melt because the melt rates are parameterized in the model run.We expect further studies using our melt rate estimates as a forcing in dynamical ice models will elucidate the complex links between ocean forcing and the dynamical ice response.
As indicated above, the downscaled models L1 and L2 are systematically colder relative to observations taken at the mouth of Scoresby Sund and inside the fjords.Since the melt rate magnitude is a near-linear function of temperature (Rignot et al., 2016;Slater et al., 2020;Xu et al., 2013), we suspect our estimated melt rate may also have a consistent and systematic low bias.To examine potential biases in our assessment of melt rate variations, we examine modeled differences in temperature and salinity profiles on the continental shelf and compare with similar differences from CTD observations (Figure 4).Despite the systematic low temperature bias in the models, remarkably the modeled fjord-shelf temperature differences match quite well with observations (Figures 4c and  4f).The effect is similar for salinity at depths lower than approximately 200 m where subglacial melt occurs.These results have two implications.First, for this study, the close comparison indicates that while the absolute magnitude of melt rates may be too low, the relative magnitude and timing of the melt rate variations are robust, supporting the interpretation of melt rate variations as drivers of dynamic glacier change.Second, these results indicate that the high-resolution modeling framework described here is a viable method to quantify shelf-toglacier variations in ocean properties induced by the fjord geometry and circulation, even in the longest and most complex fjord systems.In effect, the high-resolution model defines a fjord-specific transfer function that can be used to reliably estimate near-glacier properties given ocean information near the continental shelf.We propose that similarly-derived fjord transfer functions for all of Greenland's fjords, in combination with coarse resolution ocean projections from Atmosphere-Ocean General Circulation Models, can be used for projections of ocean temperature within the fjords.These projections could be used as a basis for generating ocean conditions on the boundaries of ice sheet models used for future ice loss projections, such as those contributing to IPCC Comparison between model output and in situ observations for temperature (top panels) and salinity (bottom panels).The left column corresponds to profiles on the continental shelf, the middle column is for profiles next to Daugaard-Jensen glacier, and the right column shows the difference of the profiles in the two locations.For comparison below 400 m, the depth of the shelf at the sample location, the shelf profiles are extrapolated downward to 1,000 m from their deepest point.sea-level rise assessments.For recent ice models, ocean conditions from coarse-resolution ocean projections are applied within fjords unless shallow bathymetry obtrudes the profile below a particular depth (Slater et al., 2020), while the observations and model results from this study demonstrate that significant variations occur between the shelf and fjord environment.
We point out three caveats of our downscale model framework as recommendations for improvements that may bring model simulations closer in line with existing observations if the ocean modeling community were to take up similar efforts around Greenland, Antarctica, or elsewhere in the polar oceans.First, our downscaled L1 and L2 models use atmospheric external forcing fields interpolated from those used in the nominal 1/3°resolution L0 global model.Since fjord circulation is known to depend on local katabatic winds, modeled circulation may be improved by higher-resolution, region-specific wind fields (Spall et al., 2017).Second, our models lack icebergs within the fjords.Recent work to simulate iceberg drift and melt within a fjord in MITgcm has revealed that icebergs have the potential to promote up-fjord heat fluxes toward glaciers, which would enhance variations in ice melt (Davison et al., 2020;Hager et al., 2023).Third, our model lacks the influence of barotropic tides on fjord circulation.In other fjord networks, such as Godthåbsfjord in West Greenland, tides have been shown to alter circulation regimes and provide additional heat toward glacier termini, which would increase melt rates (Mortensen et al., 2011).The inclusion of each of the above processes would likely bring model simulations in closer agreement with observations since they encapsulate processes known to affect fjord-scale circulation.

Conclusion
In this study, we developed a numerical modeling framework to investigate heat transport through the Scoresby Sund fjord network and estimate melt rates at the terminus of DJG.We find that variations in modeled melt rates are consistent with dynamics changes on the glacier as observed from remote sensing platforms.The downscaled model framework provides a robust estimate of the fjord-to-shelf ocean environment, providing a basis for estimating ocean boundary conditions for ice model projections of ice loss and sea-level rise contributions from the Greenland ice sheet.We expect the framework outlined here to be applicable in other regions in Greenland as well as in Antarctica and other regions where glaciers interact with the ocean.

Figure 1 .
Figure 1.Ocean bathymetry in the domains of the downscaled configurations used to investigate circulation within Scoresby Sund and quantify melt rates at the terminus of Daugaard-Jensen Glacier.(a) Global domain (Level 0) with the extent of the Level 1 domain outlined in red.(b) Level 1 domain with the extent of the Level 2 domain outlined in orange and the regional sampling polygon outlined in the thin dashed white line.(c) Level 2 domain.Background imagery is from MODIS "Natural Color" imagery.

Figure 2 .
Figure 2. Comparison of temperature timeseries in and around Scoresby Sund for (a) regional comparison between the L0 and L1 models, (b) point-wise temperature comparison on the continental shelf between the L0, L1, and L2 models, (c) pointwise temperature comparison for L1 and L2 at a point mid-way through Scoresby Sund, and (d) L2 temperature near Daugaard-Jensen glacier.

Figure 3 .
Figure 3. (a) Modeled maximal melt rate at the melt plume location of Daugaard-Jensen Glacier (DJG) estimated from the L2 model.(b) Time-integrated modeled melt rate relative to the 2000-2004 mean melt rate.(c) Ice speed near the ice front of DJF.(d) Ice elevation near the ice front of DJG.(e) Ice front position relative to the mean position 1985-2010.Positive values indicate inland retreat.

Figure 4 .
Figure 4. Comparison between model output and in situ observations for temperature (top panels) and salinity (bottom panels).The left column corresponds to profiles on the continental shelf, the middle column is for profiles next to Daugaard-Jensen glacier, and the right column shows the difference of the profiles in the two locations.For comparison below 400 m, the depth of the shelf at the sample location, the shelf profiles are extrapolated downward to 1,000 m from their deepest point.