Subglacial‐Discharge Plumes Drive Widespread Subsurface Warming in Northwest Greenland's Fjords

Greenland's glacial fjords modulate the exchange between the ice sheet and ocean. Subglacial‐discharge‐driven plumes adjacent to glaciers may exert an important influence on fjord water properties, submarine glacier melting and the export of glacially‐modified waters to the shelf. Here we use a numerical plume model in conjunction with observations from proximal to 14 glaciers in northwest Greenland to assess the impact of these plumes on near‐glacier water properties. We find that in late summer, waters emanating from glacial plumes often make up >50% of the fjord water composition at intermediate depths. These plume waters are comprised largely of upwelled Atlantic Water, warming the near‐glacier water profile and likely increasing submarine melting. Our findings demonstrate the key role played by plumes in driving water modification in Greenland's fjords, and the potential for simple models to capture these impacts across a range of settings.

Observational studies have frequently identified temperature and salinity anomalies in Greenland's fjords, attributed to the effects of glacial modification (e.g., Inall et al., 2014;Mortensen et al., 2020;Straneo et al., 2011). It is however difficult to constrain the specific impacts of subglacial-discharge plumes and submarine ice melting based on hydrographic observations alone, reflecting challenges including scarce measurements of plume properties and the multiple source waters present in fjords (Beaird et al., 2015(Beaird et al., , 2018. Combining observations with numerical modeling of plume properties has proven advantageous in this respect, providing valuable insight into plume mechanics (Bendtsen et al., 2015;Jackson et al., 2017) and aiding with the interpretation of plume modification in individual fjords (Muilwijk et al., 2022;Stevens et al., 2016).
Nevertheless, if the effects of glacial modification on ice sheet-ocean interaction are to be more fully accounted for, an assessment of these impacts across a diverse range of fjords is required. Here, we combine numerical plume modeling with hydrographic observations from NASA's Ocean's Melting Greenland (OMG) mission to characterize and quantify the impacts of subglacial-discharge plumes on near-glacier hydrography at 14 tidewater glaciers in northwest Greenland over a period of five years. This provides a novel and wide-ranging assessment of the impact of glacial modification across one of Greenland's most heavily-glaciated coastlines.

Methods
To evaluate the form and causes of near-glacier water modification, we consider three sets of water properties: (a) from close to the terminus of each study glacier; (b) from the continental shelf adjacent to each study glacier; and (c) of the plume modified water (PMW) generated by subglacial-discharge plumes at each study glacier.

Near-Glacier and Shelf Water Properties
We select 14 glaciers in northwest Greenland, based on the availability of contemporaneous OMG Conductivity-Temperature-Depth (CTD) data from near to the glacier and on the adjacent shelf ( Figure 1a) (Fenty et al., 2016;OMG, 2020). Data is available from five glaciers in 2016, one in 2017, five in 2018, eleven in 2019 and nine in 2020, and is available from August to October only ( Figure 1b).
To obtain near-glacier water properties (water type 1, above), we select the closest CTD to each of the named glaciers ( Figure 1, Table S1 in Supporting Information S1). We do not consider across-fjord variability, which we assume to be small compared to fjord-to-shelf differences (e.g., Muilwijk et al., 2022). The distance from the named glacier to the near-glacier CTD ranges from <1 to 16 km (median 5 km) ( Table S1 in Supporting Information S1). To obtain the corresponding shelf water properties (water type 2, above), we select a second CTD cast for each glacier, lying close to the fjord mouth or on the cross-shelf trough ( Figure 1a). In most cases the shelf and the near-glacier cast were recorded within 1-2 days (maximum difference of 6 days; Table S1 in Supporting Information S1), and we attribute differences between these to be due to spatial rather than temporal variability in water properties. We convert all values to Absolute Salinity (S, g kg − 1 ) and Conservative Temperature (Θ, °C), henceforth referred to as salinity and temperature, using the TEOS10 Gibbs-SeaWater Oceanographic Toolbox (McDougall & Barker, 2011).
For context, we additionally extract monthly temperature data for this region from Hadley Centre EN.4.2.2 temperature objective analyses (including .g10 bias corrections) (Good et al., 2013(Good et al., , 2023Gouretski & Cheng, 2020;Gouretski & Reseghetti, 2010) (Figure 1b and Figure S1 in Supporting Information S1). We exclude months where the stated observation weight is <0.6 at the sample location (where ∼1 shows observations have strongly influenced the analysis, and 0 indicates no observations were available).

Plume Modified Water
As direct observations of PMW properties are scarce, we model these using a buoyant plume model , forced with subglacial discharge, grounding line depth and ambient water properties. Similar models have been used extensively for tidewater glaciers (e.g., Carroll et al., 2016;Jenkins, 2011), and perform well when tested against observations (Jackson et al., 2017;Mankoff et al., 2016).
We use data from BedMachine v5 (Morlighem et al., 2017(Morlighem et al., , 2023 to approximate the depths of glacier grounding lines and the shallowest sills between each of the paired shelf and near-glacier CTD casts, as well as to delimit glacier catchments based on the hydraulic potential gradient from ice surface and bed topography (Shreve, 1972). We estimate timeseries of subglacial discharge ( Figure S2 in Supporting Information S1) based on daily ice sheet surface runoff outputs from the regional climate model RACMO2.3, statistically downscaled to 1 km resolution (Noël et al., 2016), summed over the catchment of each glacier and assumed to enter the ocean with minimal delay.
Following Jackson et al. (2017), we use a "truncated line plume" geometry of 200 m width, though a point source of buoyancy  gives similar results (not shown). We assume that all runoff enters fjords subglacially at the maximum grounding line depth, and allow the plume to overshoot the point of neutral buoyancy, terminating when the momentum decreases to zero or the plume reaches the fjord surface. We then calculate the depth of neutral buoyancy based on the plume's final density relative to the equivalent near-glacier profile. We run the plume model only for the named glaciers, though we recognize that near-glacier CTD profiles may also be influenced by other glaciers terminating nearby.
To estimate PMW properties, we run the plume model for each glacier based on the subglacial discharge at daily intervals through each melt season, using both the near-glacier and shelf temperature and salinity profiles as ambient conditions. Although the former is closest to the location of the modeled plumes, the latter allows us generate PMW properties independent of PMW already present in the near-glacier observations. In practice, these approaches produce similar results and we focus on those generated using shelf profiles, with equivalent results based on near-glacier profiles shown in Supporting Information S1. To account for obstruction of deep water ingress into fjords, we replace below-sill shelf temperature and salinity profiles with sill-level values (e.g., Morlighem et al., 2019), and extrapolate these down to the grounding line depth when necessary.  (Bjørk et al., 2015) of study glaciers (green) and shelf CTDs (red) ( Table S1 in Supporting Information S1). Near-glacier CTD locations overlap with glacier locations and are not shown. (b) EN.4.2.2 shelf subsurface water temperature (Good et al., 2013) (1 × 1° cell marked by dashed rectangle in (a)). Vertical bars in (b) show times during which CTD observations were obtained.
For a given subglacial discharge, the plume model allows us to derive the depth at which the plume is predicted to reach neutral buoyancy, and the volume flux, temperature and salinity of the PMW input to the fjord from the plume at this depth. As subglacial discharge varies over the melt season, this traces out a vertical profile of PMW input volume and properties, which can be compared directly to the equivalent profiles of shelf and near-glacier water properties. Because we are restricted to a single survey per year, we are unable to incorporate the evolution of ambient conditions on sub-annual timescales. The accuracy of modeled PMW properties, and the impact of this PMW input on the recorded near-glacier water properties, thus both likely increase toward the point of acquisition of observations in late summer/fall.

Mixing Model
As subglacial-discharge plumes find neutral buoyancy and spread out into the fjord, PMW will be intruded into and mixed with ambient waters of similar density to itself. Where these plumes are the dominant mechanism of water modification, near-glacier waters in a given density class should therefore primarily comprise waters of shelf origin mixed with intruded PMW of similar density. We can therefore use a simple linear mixing model, based on temperature and salinity as conservative tracers, to derive the ratio of shelf waters and PMW required to generate the observed near-glacier water properties.
To do this we solve (in a least-squares sense) for the fraction F of PMW (F p (ρ)) and shelf water (F s (ρ)) giving the best fit to Ѳ g (ρ) and S g (ρ), where subscripts g, s, and p denote near-glacier, shelf, and PMW respectively, and (ρ) denotes water masses of equal density: We normalize Ѳ and S values with respect to the range encountered across the three water types in each density class (Supporting Information S1), and discard results where the sum of the squared residuals is >10 −2 (50% of cases), indicating that the near-glacier water properties cannot be satisfactorily reproduced within the assumptions of the mixing model. The majority of the retained results are located between ∼100 and 300 m depth ( Figure  S3 in Supporting Information S1).

Results
Shelf water properties align with the main shelf water masses typically observed in this region (e.g., Muilwijk et al., 2022)  Near-glacier water properties are distinct from those on the shelf (Figure 2). The upper layer is typically cooler and fresher in the near-glacier profiles relative to those from the shelf, with mean differences (near-glacier minus shelf) of ΔѲ = −2.84°C and ΔS = −0.38 g kg −1 across all years (Figures 2c and 2d). The intermediate layer is conversely typically warmer (ΔѲ = 0.34°C) and saltier (ΔS = 0.03 g kg −1 ) in the near-glacier profiles during years 2016-2019, switching to slightly cooler (ΔѲ = −0.14°C) and fresher (ΔS = −0.01 g kg −1 ) in 2020. In the deep layer, near-glacier waters are slightly cooler (ΔѲ = −0.18°C) and fresher (ΔS = −0.11 g kg −1 ), particularly where sills block the flow of warm, salty water into fjords at depth.
There is a close agreement between PMW profiles modeled using the shelf and near-glacier profiles as ambient conditions in the plume model (Figure 3), due to the similarities between shelf and near-glacier properties across the deep layer in which the plume is predominantly entraining (Figures 2 and 3). Most PMW is predicted by the plume model to find neutral buoyancy in the intermediate and upper layers (Figure 4), and within these layers the similarity between the modeled PMW temperature profiles and the near-glacier profiles is in many cases striking (Figure 3 and Figure S9-18 in Supporting Information S1).
PMW fractions, calculated for all near-glacier profiles using the mixing model, are shown in Figure 4 and Figure S4 and S5 in Supporting Information S1. Values of PMW >0.5 (where 1 is pure PMW) are obtained near to many glaciers, with maximum concentrations in the intermediate layer ( Figure 4). The fraction of PMW decreases to close to zero below ∼500 m and at the surface, depths at which the modeled input of PMW is also limited (Figure 4). Sensitivity experiments show that uncertainty in grounding line depth, plume width, the plume entrainment coefficient and subglacial discharge predominantly affect the vertical distribution of PMW, whilst uncertainty in water properties predominantly affects the calculated PMW fractions (shading in Figure 4; Supporting Information S1).

Interpretation
PMW properties represent an integration of water properties between the grounding line and terminal depth of the plume. Plumes reaching the upper layer contain a mix of PW and AW, and can partly explain the cooling of this layer relative to the seasonally-warmed surface waters on the shelf (e.g., Figures 3b and 4l). The greatest impact is however observed in the intermediate layer, where plumes contain a high concentration of AW and are thus warm relative to the PW that otherwise dominates this layer (Figures 2 and 3). This plume-driven warming is seen in all years except 2020 and is most marked for glaciers at the southern end of the transect, where near-glacier waters are as much as ∼2°C warmer than shelf waters at equivalent depth ( Figure 3). In 2020, near-glacier temperatures are contrastingly cooler than those on the shelf through much of the depth profile (Figure 2c and Figure S17 in Supporting Information S1). This may be explained by a sharp warming of shelf waters shown in the EN.4.2.2 data in fall 2020 (Figure 1b and Figure S1 in Supporting Information S1), which may not have been transmitted into fjords by the time of the survey. This results in comparatively little PMW being identified by the mixing model in 2020 ( Figure S5 in Supporting Information S1), highlighting one of the challenges of relying on annual snapshots of shelf and near-glacier conditions.
Integrated over the summer, the greatest input of PMW is predicted by the plume model to occur at a shallower depth (∼50 m) than the peak concentrations of PMW identified by the mixing model (∼175 m) (Figure 4o). This may be in part due to the timing of observations in late summer/fall, when discharge decreases and plumes reach neutral buoyancy at greater depths ( Figure S8 in Supporting Information S1). It is thus possible that the impact of PMW on the intermediate layer may be reduced during peak summer, when PMW input is typically concentrated around the upper-intermediate layer interface (Figure S8 in Supporting Information S1), though we note that observations from fjords adjacent to glaciers J-K (Figure 1) in July-August 2014 also show a marked warming of the intermediate layer relative to the shelf (Bartholomaus et al., 2016, their Figure 2). While PMW input may typically be greatest around the upper-intermediate layer interface, it is likely that this also drives a relatively rapid fjord-shelf exchange at this depth (e.g., Cowton et al., 2016), which may prevent high concentrations of PMW being maintained. Compared to within the intermediate layer, the temperature of PMW is also relatively similar to that of shelf waters around the upper-intermediate layer interface (Figure 3 and Figures S9-18 in Supporting Information S1), reducing its impact on fjord water properties at this depth. Additionally, the presence of PMW may be masked by processes (such as ice melt, terrestrial runoff input and atmospheric exchanges) occurring predominantly in the upper layer. This is reflected in the spike in the number of results failing the mixing model residual threshold above ∼100 m depth ( Figure S3 in Supporting Information S1), indicating that the two component mixing model does not adequately reflect the processes contributing to water modification at these depths.
The concentration and distribution of PMW also varies between locations (Figure 4). To examine controls on this, we tested for Spearman's rank correlation between PMW depth (calculated as the mean depth of the PMW fraction distribution) and abundance (calculated as the depth-integrated PMW fraction) and key variables in the glacier/fjord system (Table S3 in Supporting Information S1). We find a positive correlation between PMW depth and grounding line depth (r = 0.58, p < 0.05), with shallower grounding lines helping plumes to reach shallower depths (e.g., Carroll et al., 2016). We do not however find a correlation between PMW abundance and distance of the observations from the glacier terminus, total annual subglacial discharge, or measures of fjord length, width or sill depth. While these factors may be expected to influence the input and dispersal of PMW (e.g., Carroll et al., 2017), it may be that this effect is masked by the relatively small sample population, uncertainties inherent in the models, and challenges in accurately quantifying PMW presence and complex fjord geometry as simple metrics .

Discussion
We find that intermediate layer water properties can frequently be explained as an approximately equal mixture of PMW and shelf PW, whilst alternative processes (such as submarine ice melt) become increasingly important toward the fjord surface. This agrees with studies from Upernavik Fjord (H in Figure 1; Muilwijk et al., 2022) and Sermilik Fjord, east Greenland (Beaird et al., 2018), both of which found high concentrations (>50%) of upwelled AW and low concentrations (<0.3%) of submarine melt water at intermediate depths. Our findings demonstrate that this high concentration of relatively warm PMW at intermediate depths is characteristic of a wide range of locations along Greenland's northwest coast, highlighting a key driver of summer water modification proximal to Greenland's tidewater glaciers.
Of particular interest is the extent to which water temperatures adjacent to tidewater glaciers, key to determining submarine melt rates, differ from those on the shelf (Catania et al., 2020;Slater et al., 2020). We find that in most cases the temperature of the deep layer is similar between shelf and near-glacial profiles, except where sills serve to block in the ingress of water (Figures 2 and 3 and Figures S9-S18 in Supporting Information S1). Contrastingly, we find that at the time of the surveys in late summer/fall, plume-driven upwelling of AW has substantially warmed the intermediate layer adjacent to glaciers (at times by as much as 2°C relative to the shelf), whilst the surface layer is cooled relative to the shelf (Figures 2 and 3). These modifications will affect both the total submarine melting and the vertical distribution of melt (which may in turn affect calving rates (O'Leary & Christoffersen, 2013;Slater et al., 2021)) compared to melt rates calculated based on shelf temperatures alone, and should be considered when modeling the oceanic forcing of tidewater glaciers (e.g., Slater et al., 2020).
Meltwater from the Greenland Ice Sheet may have a profound impact on ocean circulation if it serves to freshen surface waters in deep water formation sites, weakening the AMOC (Golledge et al., 2019;Swingedouw et al., 2022). Our results illustrate however that subglacial discharge from major tidewater glaciers usually finds neutral buoyancy below the fjord surface (with PMW input rates peaking around 50 m depth), with the addition of freshwater counteracted by the upwelling of and mixing with deep, salty AW. As such, while the efficacy of plume entrainment amplifies the volume of PMW by two orders of magnitude relative to the subglacial discharge , it limits the freshening effect of this input. While there is a substantial freshening of the upper layer of the near-glacier profiles relative to the shelf (Figure 2d), this is not explained by the major subglacial discharge inputs modeled in this study, which produce sub-surface terminating plumes of similar salinity to ambient waters (Figures S9-S18 in Supporting Information S1). Instead it is likely to be driven by a combination of submarine melting of glacier termini and icebergs (e.g., Beaird et al., 2018), shallower surface-terminating subglacial-discharge plumes (at smaller glaciers or the periphery of larger glaciers) and surface runoff inputs.

Summary
We have analyzed hydrographic observations from proximal to 14 tidewater glaciers in northwest Greenland over the period 2016-2020. By comparing paired profiles from the shelf and near-glacier environment with modeled subglacial-discharge plume properties, we have characterized and quantified the impact of plume upwelling on near-glacier water properties. Whilst the greatest input of PMW occurs around 50 m depth, we find the strongest influence of plumes on fjord water properties within the intermediate layer (∼50-250 m), where PW is mixed with and displaced by PMW. Within this layer, PMW can form the largest component of near-glacier waters. Warm, salty AW entrained in the PMW warms near-glacier waters at intermediate depths, but acts counter to the initial input of freshwater. As such, subglacial-discharge plumes will increase submarine melting and the export of glacially modified waters to the shelf, but result in comparatively little change to the salinity stratification. The implications of this export for wider ocean circulation are thus likely very different to that of surface freshwater inputs, and should be assessed as a matter of priority.