Greenland Subglacial Discharge as a Driver of Hotspots of Increasing Coastal Chlorophyll Since the Early 2000s

Subglacial discharge emerging from the base of Greenland's marine‐terminating glaciers drives upwelling of nutrient‐rich bottom waters to the euphotic zone, which can fuel nitrate‐limited phytoplankton growth. Here, we use buoyant plume theory to quantify this subglacial discharge‐driven nutrient supply on a pan‐Greenland scale. The modeled nitrate fluxes were concentrated in a few critical systems, with half of the total modeled nitrate flux anomaly occurring at just 14% of marine‐terminating glaciers. Increasing subglacial discharge fluxes results in elevated nitrate fluxes, with the largest flux occurring at Jakobshavn Isbræ in Disko Bay, where subglacial discharge is largest. Subglacial discharge and nitrate flux anomaly also account for significant temporal variability in summer satellite chlorophyll a (Chl) within 50 km of Greenland's coast, particularly in some regions in central west and northwest Greenland.

• We use buoyant plume theory to model subglacial discharge-driven nitrate fluxes across Greenland's largest marine-terminating glaciers • The largest positive nitrate flux anomalies are concentrated in a few major systems with the largest subglacial discharge fluxes • Runoff and modeled nitrate upwelling can explain temporal variability in surface cholorophyll in some coastal areas in west Greenland

Supporting Information:
Supporting Information may be found in the online version of this article.

Citation:
Oliver, H., Slater, D., Carroll, D., Wood, M., Morlighem, M., & Hopwood, M. J. (2023). Greenland subglacial discharge as a driver of hotspots of increasing coastal chlorophyll since the early 2000s. Geophysical Research Letters, 50, e2022GL102689. https://doi. org/10.1029/2022GL102689 et al., 2017. The estimated upwelling flux of nitrate, a nutrient that limits primary production in most Arctic systems (Codispoti et al., 2013), though not southeast and south of Greenland (Garcia et al., 2014), theoretically increases nonlinearly with a deeper grounding line and increasing subglacial discharge (Hopwood et al., 2018). This upwelling mechanism has been linked to higher productivity and Halibut catches within fjords containing marine-terminating glaciers along west Greenland (Meire et al., 2017). Additionally, some proportion of these upwelled nutrients could be exported to the shelf (Oliver et al., 2020), potentially leading to broader-scale fertilization effects.
Though the physical oceanography of only around a dozen glacier fjord systems has been characterized around Greenland, hundreds of marine-terminating glaciers surround Greenland's coastline (Morlighem et al., 2017) and therefore significant buoyant upwelling is likely to be occurring on a pan-Greenland scale. The spatial extent of the downstream biogeochemical signature of subglacial discharge in the ocean is challenging to assess on a pan-Greenlandic scale, however, due to the limited number of case studies and the broad variability in glacier and fjord properties (Morlighem et al., 2017;Straneo & Cenedese, 2015). Further difficulties arise when distinguishing local processes driven by runoff from broader-scale changes, such as sea-ice loss, which also drive shifts in phytoplankton growth across the high-latitude oceans (Arrigo & van Dijken, 2015;. Satellite work in the Labrador Sea and Baffin Bay has shown increasing primary production in September across most of the region, a shift towards a later spring bloom, and correlation between bloom timing and modelled runoff arrival in the southeast Labrador Sea (Arrigo et al., 2017;York et al., 2020). Concerning the potential impact of increasing Greenland subglacial discharge on marine primary production, several critical questions remain: (a) how do subglacial discharge-driven nutrient fluxes vary around Greenland; (b) to what extent does increasing subglacial discharge impact upwelled nutrient fluxes; and (c) on what scale does subglacial discharge have a detectable influence on productivity?
To address these questions, we use buoyant plume theory to model nitrate upwelling from marine-terminating glaciers on a pan-Greenland scale. We compile a suite of nitrate profiles from the Greenland shelf and pair them with newly available data products, including high-resolution Greenland bathymetry, modeled runoff estimates, and fjord properties from the NASA Oceans Melting Greenland project (OMG; Fenty et al., 2016). We then compare model results with satellite-based estimates of Chl to assess the degree to which subglacial discharge corresponds to surface-ocean productivity across sectors of Greenland's shelf.

Materials and Methods
We used a buoyant plume model combined with runoff estimates from Mankoff et al. (2020), temperature-salinity profiles from OMG, bathymetry from BedMachinev4 (Morlighem, 2021), and shelf nitrate profiles to estimate the vertical flux of nitrate adjacent to marine-terminating glaciers around Greenland. The physics of the plume model and the physical datasets used are the same as in Slater et al. (2022), which enables us to consider 136 of Greenland's largest marine-terminating glaciers.
The relevance of the nitrate flux to changing nitrate distribution within the water column depends on the depth of the nitracline and the distribution of nitrate in fjord source waters (Cape et al., 2019). We compiled a suite of nitrate observations into seven mean regional spring-summer shelf profiles (April-September) ( Figure 1): north (NO), northwest (NW), northeast (NE), central west (CW); central east (CE), southwest (SW), and southeast (SE) (Text S1 in Supporting Information S1). There is a notable large spread in NE and CE profiles which may reflect a combination of factors, including the broad width of these shelf regions and variability in the influence of sea-ice cover which influences near-surface nutrient drawdown (Michel et al., 2015).
The plume model solves for plume volume, momentum, heat and salt fluxes as the plume rises due to its buoyancy while being diluted by entrainment of fjord waters(for a full description of the model see Slater et al., 2016). We add an evolution equation for nitrate concentration (supporting information), enabling us in this study to focus additionally on nitrate fluxes. Two critical inputs to the plume model are the subglacial discharge flux and the plume width. We input mean, minimum, and maximum monthly subglacial discharge estimates averaged over the melting season (year days 146-255), and over years 2010-2019; our results should be interpreted as broadly representative of recent summers. Plume volume flux is sensitive to the assumed plume width for a line plume-idealized theory suggests flux scales as plume width raised to the power 2/3 (e.g., Jenkins, 2011)-yet, few observations are available to constrain plume width. Here we assume a line plume width of 200 m following the observational constraints of Jackson et al., 2017, but also probe the sensitivity of our results by considering plume widths of 100 and 300 m. Note that for simplicity and given the scale of the study, we treat any glaciers with floating ice shelves as approximately vertical front tidewater glaciers.
As a plume rises from the grounding line of a glacier, its density increases due to entrainment of ambient fjord waters. At some depth the density of the plume will be equal to the density of the fjord; this is the neutral buoyancy depth. We assume the plume intrudes into the fjord at this depth, transporting a flux equal to the upwelling flux at this depth. In a homogenous water column, the vertical nitrate flux would not impact nitrate distribution; to account for this, we use the nitrate flux anomaly (NFA). For each glacier, we calculate the NFA (in mol s −1 ) from the plume model at each glacier as: where the nitrate anomaly is computed as the difference of NB (modeled nitrate concentration at the neutral buoyancy depth) and SP (nitrate concentration at the neutral buoyancy depth in the regional shelf profile), and NB is the volume flux of the plume at neutral buoyancy.
We also assess whether glacier-driven nutrient fluxes may be sufficient to drive temporal variability in local phytoplankton biomass by examining how well runoff (which includes both surface runoff and subglacial discharge from Mankoff et al. (2020)), subglacial discharge, and modeled NFAs correspond to local surface Chl.
As surface runoff can also affect phytoplankton growth potential (e.g., Oliver et al., 2018), we include it in our analysis to explore whether there is a Chl response in proximity to land-terminating glaciers without turbulent buoyant plumes. Greenland marine-terminating glaciers where upwelled nitrate flux is modeled, colored by their nitrate region. Xs show locations of shelf (depth <600 m) nitrate profiles used to construct regional smoothed means. (b) Regional nitrate profiles. Colored lines show smoothed regional mean values, with the concentrations used to create the regional means shown by gray circles.
For the Chl analysis, we used published estimates of Arctic surface-ocean Chl concentrations from 2003-2018 . While satellite-derived estimates of coastal Chl concentrations can be subject to large uncertainties (Zheng & DiGiacomo, 2017), high turbidity has not been found to impact Greenland coastal fluorescence line height over the Greenland shelf offshore of fjord waters, so estimates of Chl a using the MODIS algorithm are likely accurate (Arrigo et al., 2017). Chl concentrations were calculated over two shelf masks: a near-coastal mask (defined as 50 km from the nearest land pixel) and over the entire shelf (defined at the 600-m isobath).
We calculated the early summer (June-July), late summer (August-September), and seasonal (June-September) mean Chl within 1-degree latitude bands for both E (from 60° to 80°N) and W (from 60° to 76°N) Greenland (36 masks total). We then paired each resulting Chl estimate with mean runoff, subglacial discharge, and NFA from all glaciers within its corresponding latitude band. We use mixed layer depths from ECCO LLC270 for comparison with Chl (Text S3 in the Supporting Information).

Modeling Pan-Greenland Nitrate Upwelling Variability
The  (Figures 2c and 2d), and had a weaker relationship with sill depth (R 2 = 0.33, p < 0.001) and grounding line nitrate concentration (R 2 = 0.18, p < 0.001). Deep grounding line depths did not always correspond to large NFAs (R 2 = 0.16, p < 0.001), with several deeply grounded glaciers in SE Greenland having low NFAs (Figure 2e). Modeled NFAs were also less sensitive to plume width than subglacial discharge rate; increasing the plume width by 100 m resulted in a 4% increase in the totaled pan-Greenland NFA; decreasing the width by 100 m resulted in a 10% decrease ( Figure S1 in Supporting Information S1).
Systems with large nitrate upwelling were present across most coastal sectors, except SE Greenland, which had very low NFAs relative to all other sectors (Figure 2a). Low SE NFAs result from both comparatively low subglacial discharge fluxes and low nitrate anomalies at the neutral buoyancy depth, as near-surface nitrate in the SE profile is not as depleted compared to other regions (Figure 1). Model estimates of nitrate, nitrate anomaly, and nitrate fluxes are in general agreement with estimates calculated from data collected from well-sampled Greenland fjord systems (Text S4 and Table S1 in Supporting Information S1). Notably, while basal melt is significant at the 79°N Glacier, a floating ice shelf, our entrainment estimates are similar to those made with observations (Krisch et al., 2021).
The 19 largest upwelling systems (14% of the studied glaciers) accounted for over 50% of the total NFA summed over all 136 glaciers ( Figure S2 in Supporting Information S1), with the largest NFAs corresponding to glaciers with large subglacial discharge fluxes. Jakobshavn Isbrae, which flows into Ilulissat Icefjord and Disko Bay in CW Greenland, is the largest upwelling system and is responsible for 6%-8% of the total subglacial discharge and 6%-14% of the total modeled NFA, with the proportion decreasing with increasing subglacial discharge. At Jakobshavn Isbrae, increasing discharge by a factor of 2.2 increases the modeled NFA at the glacier by 50%, from 297 mol s −1 to 445 mol s −1 . Most of the pan-Greenland NFA is concentrated in W Greenland (SW + CW + NW, 69%-76%). In contrast, only 13%-17% of subglacial discharge is released from the 27 SE Greenland glaciers combined, which constitute 11%-12% of the total pan-Greenland NFA.
We also consider the modeled plume neutral buoyancy depth (Figure 2e), which affects (a) the bioaccessibility of upwelled nitrate, which depends on the plume position relative to the euphotic zone depth, and (b) the nitrate anomaly at neutral buoyancy, which depends on the plume position relative to the nitracline. With the mean subglacial discharge, only 26% (most in W Greenland) of plumes reach a modeled neutral buoyancy depth shallower than 40 m, which is roughly the maximum euphotic depth estimated for Greenland fjord systems (Holding et al., 2019;Holinde & Zielinski, 2016;Mascarenhas & Zielinski, 2019;Murray et al., 2015). This result does not discount the possibility that short periods of high subglacial discharge fluxes could still drive plumes into the photic zone, however: with the maximum monthly subglacial discharge rates, 59 plumes (43%) achieved a modeled neutral buoyancy depth shallower than 40 m.

Comparing Plume Model to Satellite Chl
Variability in subglacial discharge and NFA can broadly account for interannual ( Figure 3) variability in near-coastal Chl. Over annual timescales, a statistically significant proportion of the variability in late summer coastal Chl along CW and NW Greenland can be accounted for by variability in late summer runoff ( Figure 3a) and modeled NFA (Figure 3b). Though runoff fluxes are largest in SW Greenland, where most runoff is from land-terminating glaciers, the relationships between total discharge and Chl in this region are weak. There were also few masks with statistically significant relationships between Chl and runoff/NFA along E Greenland, where total seasonal runoff is 23% lower, neutral buoyancy depths are 53% deeper (Slater et al., 2022), and nitrate anomalies at neutral buoyancy are 63% lower (Figure 2). While runoff and NFAs are greatest in early summer, stratification is also strongest during this period, and the early summer Chl-runoff and Chl-NFA relationships were statistically significant in only a few masks. Over the entire summer, strong interannual relationships between Chl and runoff are mainly found in CW Greenland and are limited to the Disko Bay/Uummannaq Fjord region when considering NFA. These relationships lose statistical significance when expanding from 50-km to over the entire continental shelf ( Figure S3 in Supporting Information S1).
With the large NFA predicted at Jakobshavn Isbrae, the plume model suggests that Disko Bay would be the coastal region most susceptible to subglacial discharge-fueled increases in primary production, particularly the northern part of the bay given the inferred pathway for export of waters from Ilulissat Icefjord (Beaird et al., 2017). There is some evidence of this already occurring in satellite Chl data from 2003 to 2018, with seasonal Chl across northern Disko Bay positively correlated with subglacial discharge, nitrate anomaly, and NFA over this period (Figure 4). Chl similarly correlates to subglacial discharge in Uummannaq Fjord, just north of Disko Bay, but the slope of the relationship is much steeper here. While large NFAs are also predicted for Kangerlussuaq Fjord and Sermilik Fjord in E Greenland, coastal Chl does not have statistically significant correlations with subglacial discharge in these systems. There is no interannual relationship exists between mixed layer depth and Chl in any of these systems (Figure 4b).

Discussion
We demonstrate that (a) large NFAs estimated from buoyant plume theory are limited to a few major upwelling systems around Greenland and (b) that these subglacial discharge-driven anomalies explains can significant interannual variability in coastal surface-ocean Chl, in certain regions downstream of W Greenland glaciers. While the plume model is highly idealized and does not account for complex fjord hydrodynamics or higher-frequency variability in computed subglacial discharge fluxes, it generally agrees with estimates of plume nitrate concentration and nitrate fluxes made near well-sampled marine-terminating glacial fjord systems. Furthermore, a strong correlation between NFAs and Chl in specific cases increases confidence in the underlying plume theory because these two variables are determined independently. Given the significant relationships between Chl and NFA, the proposed buoyant plume mechanism appears to drive a significant proportion of interannual variability in late summer CW and NW Chl on a scale of ∼50 km from the coast and regionally within Disko Bay and Uummannaq Fjord. These regions in CW Greenland both host large inshore fisheries of Greenland halibut and northern shrimp (Pandalus borealis), Greenland's two most exported catches by volume and value (Berthelsen, 2014;Jørgensen & Arboe, 2013).
It is important to note that increased runoff can also drive increased stratification, shallower mixed layer depths, and increased turbidity . After particles sink, this could increase light availability for phytoplankton in late summer (Oliver et al., 2018), presenting another mechanism that could explain correlation between Chl and subglacial discharge on approximately the same spatial scale. However, modeled mixed layer depths do not explain variation in Chl near high NFA systems in W Greenland (Figure 4b). Such a stratifying effect would also be expected to occur broadly across most of the Greenland shelf, yet no evidence is found herein for a broad-scale Chl-runoff relationship. As the impact of runoff on shelf mixed layer depths occurs mainly after the peak meltwater season (Oliver et al., 2018), this stratifying mechanism appears to be at most a minor contributing factor to the observed interannual variability in seasonal Chl.
Furthermore, summertime phytoplankton growth is generally limited by nitrate availability rather than light, as evidenced by consistently low surface nitrate in all coastal sectors except the SE (Figure 1). Strong stratification in an inner-fjord environment has been demonstrated to lead to low summertime primary production via reduced nitrate supply (Holding et al., 2019). Increased stratification can suppress turbulent nitrate fluxes into the upper mixed layer (Randelhoff et al., 2020), which may explain why runoff and NFA explains little interannual variability in coastal Chl in early summer-the peak melting season. In our study, modeled NFAs correspond more strongly to nearby seasonal Chl along W Greenland than E Greenland. Differences in the response between W and E Greenland glaciers may be attributable to nitrate not being limiting in SE Greenland waters, and E Greenland glaciers generally having deeper modeled neutral buoyancy depths (Slater et al., 2022), which may have little immediate impact on surface-ocean Chl if the outflowing plume lies below the euphotic zone. E Greenland runoff may still have a downstream biological effect, however, because our method only assess the near-glacier effects on coastal Chl and does not account for the effects of advection.
Our choice of using a uniform coastal distance for Chl processing will also result in differing influences of boundary currents between systems and the coverage of fjord and coastal waters depending on local fjord geometry. For example, SE and CE Greenland glaciers are characterized by relatively long fjords, while Baffin Bay (NW) glaciers often terminate in shorter fjords (Slater et al., 2022). Subsurface subglacial plumes can also result in subsurface Chl maxima, as observed in Godthåbsfjord (Meire et al., 2017), and so subglacial discharge-driven subsurface maxima may not be detectable to the same extent on a pan-Greenland scale when using satellite-derived Chl.
If subglacial discharge continues to increase in coming decades, more fjords with deep glaciers may emerge from background variability as regional drivers of summertime phytoplankton blooms around Greenland. Totaled across all modeled glaciers, the modeled subglacial discharge-driven upwelling supplies 448 Gg additional N year −1 at the neutral buoyancy depth (using the mean 2010-2019 discharge), with the 10 largest systems out of the 136 considered herein accounting for 33% of this additional supply. 340 Gg of this additional yearly N supply is delivered shallower than 100 m. Using the maximum monthly subglacial discharge rates across all glaciers, the upwelled nitrate supply doubles (899 Gg additional N year −1 ).
It is possible to calculate the buoyant plume modeled nitrate flux at neutral buoyancy (R 2 = 0.99) using only the subglacial discharge rate, grounding line depth, and the nitrate concentration at the grounding line depth, with a power law function ( Figure S4 in Supporting Information S1). It could therefore be included in coarser-scale models (Carroll et al., 2020) without explicitly running an online plume model, to compute estimated changes in nutrient availability from ocean-glacier interactions, though these estimates would be sensitive to uncertainties in discharge rates and plume widths. This is also practical from a spatial perspective, as buoyant plumes alter nutrient distributions on sub-grid scales even in higher-resolution ocean simulations.

Conclusion
We used buoyant plume theory to model subglacial-discharge driven nutrient upwelling on a pan-Greenland scale and find that most of the total upwelled nitrate flux is concentrated in just a few key systems. Subglacial discharge and modeled NFA can account for significant temporal variability in seasonal satellite Chl near the studied glaciers in some areas, with the slopes of those relationships sensitive to the coastal sector and the grounding line depth. These results suggest that subglacial discharge may drive increased coastal primary production in certain high-runoff sectors of the Greenland Ice Sheet, including Disko Bay and Unmmannaq Fjord.