Rapid ice loss from the Greenland ice sheet since 1992 is due in equal parts to increased surface melting and accelerated ice flow. The latter is conventionally attributed to ocean warming, which has enhanced submarine melting of the fronts of Greenland's marine-terminating glaciers. Yet, through the release of ice sheet surface meltwater into the ocean, which excites near-glacier ocean circulation and in turn the transfer of heat from ocean to ice, a warming atmosphere can increase submarine melting even in the absence of ocean warming. The relative importance of atmospheric and oceanic warming in driving increased submarine melting has, however, not been quantified. Here, we reconstruct the rate of submarine melting at Greenland's marine-terminating glaciers from 1979 to 2018 and estimate the resulting dynamic mass loss. We show that in south Greenland, variability in submarine melting was indeed governed by the ocean, but, in contrast, the atmosphere dominated in the northwest. At the ice sheet scale, the atmosphere plays a first-order role in controlling submarine melting and the subsequent dynamic mass loss. Our results challenge the attribution of dynamic mass loss to ocean warming alone and show that a warming atmosphere has amplified the impact of the ocean on the Greenland ice sheet.

Ocean forcing is understood to have driven half of the mass loss from the Greenland ice sheet during 1992-2018 [1][2][3] , leading to the ice sheet's contribution to global sea level rise increasing from 5% in 1993 to 25% in 2014 4 . Ocean forcing manifests as melting of the submerged fronts of Greenland's marine-terminating glaciers (submarine melting) and can result in glacier retreat and acceleration either directly or by increasing the rate of iceberg production [5][6][7] .

The rate of submarine melting is controlled both by the heat content of the ocean near the glacier and by the turbulent processes that flux heat from the ocean into the ice 8,9 . These turbulent processes are strongly influenced by vigorous upwelling of ocean waters in buoyant plumes adjacent to the ice, driven by the release of fresh subglacial discharge from beneath the glacier (Fig. 1) 10,11 . In turn, subglacial discharge is sourced largely from atmospheric-driven melting of the ice sheet surface 12 . Thus, submarine melting experiences both oceanic and atmospheric influences (Fig. 1), with oceanic variability controlling the heat available for melting and atmospheric variability determining to a large extent how much of the heat is actually used.

The relative importance of atmospheric and oceanic sources of variability to submarine melting and the subsequent ice sheet dynamic mass loss has, however, not been quantified. The prevailing view is that submarine melt rate variability and the consequent impact on marine-terminating glacier dynamics result predominantly from changes in the ocean 1,3,5,13 , but previous work has either been limited to individual glaciers, or has not accounted for or not isolated the atmospheric influence on submarine melting. This limits our understanding of the fundamental sensitivity of the ice sheet to climate warming and restricts our ability to project future change.

Here, we use observations and reanalysis of the atmosphere and ocean to reconstruct submarine melt rates at Greenland's marine-terminating glaciers over a 40-year period from 1979 to 2018. We isolate variability driven by the atmosphere and ocean respectively subglacial discharge (Extended Data Fig. 1) and ocean thermal forcing (Extended Data Fig. 2), we estimate monthly submarine melt rate using monthly subglacial discharge and monthly ocean thermal forcing before averaging over each year from 1979 to 2018. When splitting out variability in submarine melting driven by the atmosphere or ocean we allow one factor to vary while holding the other constant at its mean 1979-1993 value (Supplementary Information).

The 40-year mean annual submarine melt estimate shows that at the ice sheet scale, spatial variability in submarine melting is controlled to first order by the distribution of ocean heat governed by the large-scale ocean circulation (Fig. 2a). The highest melt rates are found in south Greenland where glaciers and fjords have ready access to warm Atlantic waters from the Irminger Sea. There is a marked latitudinal gradient in submarine melting along west Greenland north of the Davis Strait, with melt rates declining with increasing distance from the warm source waters in the south. Even more pronounced is the gradient across Denmark Strait, meaning that glaciers in northeast Greenland experience much colder ocean waters and relatively little submarine melting. Glaciers in north Greenland, farthest from the source of warm water, experience the coldest ocean and least submarine melting of all (Fig. 2a).

Grouping of the glaciers into five regions shows that submarine melt rate varies within a region due to order-of-magnitude differences in subglacial discharge and smaller differences in ocean thermal forcing (Fig. 2b). Variability in subglacial discharge is linked to hydrological catchment size, and there exist statistically significant relationships between ocean thermal forcing and grounding line depth (Extended Data Fig. 5), and submarine melt rate and grounding line depth, in each region (Fig. 2c). As a rule, annual submarine melt rate increases by 0.2 to 0.4 m d -1 for every 100 m of grounding line depth. Therefore, although grounding line depth does not appear explicitly in the submarine melt rate parameterization, it plays an important role in determining glacier-to-glacier variability in submarine melt rate within a region.

Annual submarine melt rate increased in all regions from the 1990s to peak in the 2000s at 13-50% above 1979-1993 mean values (Fig. 3a-d, black). Since the 2000s, melt rates have decreased but have not returned to pre-1990s values. Note that here we highlight only the three regions that dominate Greenland's dynamic sea level contribution (south, central-west and northwest Greenland); northeast and north Greenland are described in the Extended Data Figs. 6 and 7.

The dominant driver of change in submarine melt rate varies by region and is revealed by allowing only the atmosphere or ocean to vary while holding the other constant. In south and central-west Greenland, ocean variability alone explains most of the change in submarine melt rate (Fig. 3a,b). The rapid increase in submarine melting from the mid-1990s in these regions was driven almost exclusively by increased ocean thermal forcing, and even high subglacial discharge during the 2010s has played only a minor role (Fig. 3e,f). By contrast, in northwest Greenland atmospheric and oceanic variability seem to be approximately equally important to the evolution of submarine melt rate (Fig. 3c). The importance of the atmosphere in this region results from the dramatic increase in subglacial discharge beginning in 2000 relative to the much more muted increase in ocean thermal forcing (Fig. 3g).

For Greenland as a whole, the ocean seems to be more important in controlling submarine melting, especially during the 1990s, but the atmosphere has played a substantial role, especially since 2000 (Fig. 3d,h). For example, the ten-year smoothed submarine melt rate in 2010 was 22% above the 1979-1993 value, but only three-fifths of this elevated value can be accounted for by the ocean alone. and force a two-stage glacier model with the resulting time series to estimate the impact on ice sheet dynamic mass loss. Our findings suggest that atmospheric variability is nearly as important as oceanic variability to submarine melting at the ice sheet scale and is dominant in some regions, challenging the paradigm that equates Greenland's dynamic mass loss with ocean variability alone.

We quantify plume-driven submarine melt rate, ṁ, for 123 marine-terminating glaciers around the ice sheet as ṁ = 0.142 Q 0.31 TF 1.19 , in which Q is the subglacial discharge arising from both surface melting and frictional melting at the bed and TF is the ocean thermal forcing, defined as the difference between the in situ ocean temperature and freezing point (Methods and Extended Data Fig. 3). While there exists substantial uncertainty regarding the absolute values predicted by the parameterization 9 , the form of the parameterization, which dictates the relative change in melt rates used in this study, has extensive support from theory 14,15 , observations 10,[16][17][18] and modelling 19,20 . Subglacial discharge arising from surface melting is estimated by summing the surface meltwater volume, as simulated by the regional climate model RACMO 21 , over the subglacial hydrological catchment of the glacier 22 (Methods). The assumption that the majority of surface meltwater enters the fjord subglacially is supported by a high moulin and crevasse density at the margin of the ice sheet 23 and by studies that find agreement between RACMO-estimated subglacial discharge and plume freshwater content 10,17 . Discharge from melting at the ice sheet bed is assumed constant and is estimated following ref.

. Ocean thermal forcing is estimated using a multi-model mean of the ocean reanalyses ASTE_R1, ORAS5 and CHORE_RL [25][26][27] and the observational dataset EN.4.2.1 28 , with properties extrapolated into fjords from the nearest available values on the continental shelf, taking into account fjord bathymetry (Methods). Given the seasonality of both Freshwater emerging from beneath marine-terminating glaciers (subglacial discharge) drives vigorous upwelling plumes. Subglacial discharge can be sourced either from surface melting or frictional melting at the bed. Plumes entrain warm fjord waters, which are replenished by the ocean and fjord circulation processes. The combination of fast-moving and warm water in the plume drives melting of submerged ice at the glacier terminus (submarine melting). By the described chain of processes, submarine melting is influenced by both atmospheric and oceanic variability.

We assess the impact of submarine melting on Greenland's marine-terminating glaciers using a simple two-stage model that has been shown to reliably capture a glacier's dynamic response to external forcing when compared to more complex models 29,30 . The model accounts for terminus ice loss due to submarine melting and the subsequent potential amplification of calving. For each of the 123 glaciers, we run the model forwards from 1979 holding surface mass balance constant and assuming a linear prograde bed slope, but imposing the three submarine melt rate time series (atmosphere only, ocean only, both vary). We generate a large ensemble of simulations by varying all remaining free parameters (Methods), which control the timescale and magnitude of the glacier's dynamic response. Lastly, we group the simulations by ice sheet region. We strongly emphasize that due to their simplicity, particularly regarding geometry and their neglect of forcings beyond submarine melting, these simulations are not intended to accurately simulate any individual glacier. Instead, these simulations quantify the dynamic response to submarine melt rate of an idealized glacier located in a given region of the ice sheet, essentially providing a mapping from submarine melt rate (Fig. 3a-d) to retreat and dynamic sea level contribution.

All simulations show retreat and a positive dynamic sea level contribution in response to increased submarine melting (Fig. 4). In south and central-west Greenland, simulations in which submarine melt rate by the black ticks. Other labels are the Irminger Sea (Irm), Davis Strait (Dav)

and Denmark Strait (Den). Bathymetry is from ref. 42 and ref. 43 . b, Subglacial discharge (x axis, note logarithmic scale) and ocean thermal forcing (y axis) for each marine-terminating glacier. The background shading shows the resulting submarine melt rate. Glaciers are coloured by their regional grouping. The larger squares and error bars show the median and interquartile range for each region, respectively. c, Submarine melt rate versus grounding line depth by region with fitted linear trends (all significant at the 5% level) as dashed lines. All the results shown in these plots are annually averaged over 1979-2018.

https://doi.org/10.1038/s41561-022-01035-9

fluctuates only due to oceanic variability give ~75% of the retreat and dynamic sea level contribution compared with simulations in which both atmospheric and oceanic variability are considered (Fig. 4a,b,e,f), consistent with the limited importance of atmospheric variability to the submarine melt rate (Fig. 3a,b). Conversely, in northwest Greenland, atmospheric-sourced variability in submarine melting alone explains ~60% of retreat and dynamic sea level contribution, with oceanic variability playing a lesser role (Fig. 4c,g), once more consistent with the greater importance of the atmosphere to submarine melting in northwest Greenland (Fig. 3c). There is only little sensitivity of these percentages to the choice of a particular ocean reanalysis product, except in northwest Greenland where use of only the ORAS5 dataset would suggest an even greater role for the atmosphere (Extended Data Fig. 10). At the ice sheet scale, oceanic-sourced variability in submarine melting alone accounts for ~70% of retreat and dynamic sea level contribution (Fig. 4d,h), such that atmospheric-sourced variability in submarine melting (which accounts for the remainder) has played a substantial amplifying role in glacier retreat and ice sheet dynamic mass loss.

The relative importance of atmospheric and oceanic sources of variability in submarine melting to ice sheet dynamic mass loss (Fig. 1) is determined firstly by how these factors combine in the submarine melt rate parameterization, secondly by their observed relative variability (Figs. 2 and 3) and thirdly by the sensitivity of glaciers to submarine melting (Fig. 4). Owing to plume dynamics and the physics of heat transfer from ocean to ice, the melt rate parameterization contains the subglacial discharge and ocean thermal forcing raised to the power 0.31 and 1.18, respectively. Thus, submarine melt rate is approximately four times more sensitive to equivalent variability in ocean thermal forcing than subglacial discharge. Yet, in all regions, the relative increase in subglacial discharge from the 1990s into the 2000s is larger than for the thermal forcing, dramatically so in northwest Greenland (Fig. 3e-h). Once the glacier model is forced by these submarine melt rate time series, seasonality also modulates the relative importance of atmosphere and ocean. The increased subglacial discharge from a warming atmosphere manifests as rapid submarine melting during summer but makes little difference during winter, while the impact of a warming ocean is felt year-round. In our simple simulations the glaciers seem to be particularly sensitive to rapid summer melting, lending greater importance to atmospheric-sourced variability than might be expected from examining submarine melt rates alone (Figs. 3 and 4).

In sum, our results show that the retreat and speed-up of Greenland's marine-terminating glaciers over past decades was a result of a combination of atmospheric and oceanic variability whose relative role varied by sector: in south and central-west Greenland it was largely a result of ocean warming, while in northwest Greenland the atmosphere played the greater role (Fig. 4). We suggest that these regional differences arise from the same large-scale atmospheric and oceanic circulations that explain spatial variability in submarine melting (Fig. 2). South Greenland is directly exposed to North Atlantic Ocean variability 1 , while northwest Greenland sees only a muted signal because waters lose much of their heat as they travel up the west coast of Greenland (Fig. 2) 3,13,31,32 . Similarly, owing to Arctic amplification and a differing cloud response to large-scale atmospheric changes, northern Greenland has seen greater relative increases in air temperature and subglacial discharge than southern Greenland, although changes in the south remain larger in an absolute sense 33 . Our conclusions thermal forcing (red) relative to the 1979-1993 mean for SO (e), CW (f), NW (g)

and all Greenland (h). In all plots the lines show the median value over all marineterminating glaciers in the region, while the shading shows the interquartile range. The fact that this shading is relatively narrow shows that all glaciers in a region experience similar variability in forcing. All time series have been smoothed using a ten-year window. Equivalent plots for the two ice sheet regions not highlighted here are given in Extended Data Fig. 6.

https://doi.org/10.1038/s41561-022-01035-9

for northwest Greenland are comparable to ref. 34 , which found that marine-terminating glaciers on the other side of Baffin Bay in Arctic Canada have primarily responded to atmospheric temperature. The primary implication of our results is that a full understanding of submarine melting and its impacts requires careful accounting for both atmospheric and oceanic change, particularly as large-scale circulation ensures that these changes may not be coupled either in space or in time. For example, global climate models project nuanced and regionally variable ocean warming over the coming century [35][36][37][38] . Yet the same models, once downscaled using regional climate models, also project vast, ice-sheet-wide increases in subglacial discharge [38][39][40] that may amplify or overwhelm ocean-sourced variability in submarine melting, leading to the possibility that submarine melt rates could continue to increase and drive glacier retreat even if for some regions or time periods the ocean has not warmed. The sea level implication is clearly most relevant to regions of the ice sheet that are exposed to the ocean through many tidewater glaciers, such as southeast and northwest Greenland 41 .

The dynamics of plumes at Greenland's tidewater glaciers ensure that atmospheric variability is a first-order control on submarine melting and the associated glacier retreat and dynamic mass loss. A simple glacier model suggests that, without concurrent atmospheric warming, the retreat and dynamic sea level contribution of Greenland's marine-terminating glaciers would have been reduced by around a third (and by nearly a half in northwest Greenland). If we are to attribute the recent dramatic retreat and acceleration of Greenland's marine-terminating glaciers to increased submarine melting, then this is due to a warming of both the atmosphere and ocean. As such, the atmosphere has substantially amplified the response of the ice sheet to ocean warming and will likely continue to do so over the coming century.

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-022-01035-9.

Only a brief description of the methods is given here. For full details, see the Supplementary Information.

To define the set of marine-terminating glaciers considered, we start from a list of Greenland's 243 fastest-flowing glaciers provided in Table S1of ref. 43 . To focus on marine-terminating glaciers, we remove 69 glaciers that are either land-terminating or effectively land-terminating (maximum grounding line depth <50 m). We also remove 41 glaciers for which the deepest path to the ocean crosses a point shallower than 25 m. These glaciers are either relatively isolated from the ocean, or if the bathymetry is inaccurate then we cannot reliably infer ocean properties for these glaciers. Lastly, we remove 10 glaciers with a mean annual subglacial discharge <2.5 m 3 s -1 . The final dataset of 123 marine-terminating glaciers (Fig. 2) contains 76 of Greenland's fastest 100 glaciers 43 .

We estimate monthly subglacial discharge for each marine-terminating glacier (Extended Data Fig. 1) by delineating its hydrological catchment and estimating surface runoff and basal melt rates following the method of ref. 22 . The hydrological catchment is defined by performing flow routing 44 on the subglacial hydropotential defined in terms of the bed topography and ice thickness according to BedMachine v3 43 , and assuming a subglacial water pressure equal to the ice overburden pressure. Monthly surface runoff is obtained from the regional climate model RACMO2.3p2 21 , while a constant basal melt rate is assumed following ref. 24 . Subglacial discharge is then the sum of the surface and basal meltwater runoff over the hydrological catchment of the glacier. We assume that the catchments are stable throughout the time period and that little supraglacial redirection into neighbouring catchments occurs.

In the absence of continuous, long-term oceanographic observations from a substantial number of glacial fjords around Greenland, we have used the ocean reanalysis and observational products ASTE_R1, ORAS5, CHORE_RL and EN.4.2.1 [25][26][27][28] to quantify properties on Greenland's continental shelf (for example, Fig. 2). For each glacier, we apply an algorithm similar to ref. 45 and ref. 38 to identify the deepest, closest water on the continental shelf that is able to access the calving front without being impeded by bathymetry (Extended Data Fig. 2). Note that we assign only a single grounding line temperature and salinity to the glacier rather than a full profile-we assume that owing to upwelling in the plume, this grounding line value is the most relevant for submarine melting of the glacier. Grounding line water properties are then converted to thermal forcing, defined as the difference between the in situ temperature and in situ freezing point, using the linearized freezing point equation of ref.

. Because each of the four ocean products considered have different resolutions and bathymetry, the process is repeated separately for each product, and the resulting four time series are combined into a single time series by averaging over the data products to create a multi-model mean (Extended Data Fig. 2).

Submarine melting induced by a plume can be parameterized in the form ṁ = κ Q α TF β , where κ, α and β are constants, Q is subglacial discharge and TF is ocean thermal forcing (for example, refs. 14,15,20 ). Note that within a discharge plume, submarine melt rate close to the grounding line is independent of terminus shape 46 and so provided we consider melt rate close to the grounding line (as in this study), we can equally apply this parameterization at vertical or undercut fronts and at glaciers having an ice tongue or shelf. We take ṁ to be the maximum submarine melt rate in the plume and fix the constants κ, α and β by comparison to a standard buoyant plume model 14,15 with a line plume geometry of width w (ref. 17 ), so that the subglacial discharge per unit width is then q = Q/w. We perform a parameter sweep over q, TF and grounding line depth and fit the parameterization to these simulations by least squares (Extended Data Fig. 3), obtaining ṁ = 0.732 q 0.31 TF 1.19 . Note that beyond its effect on thermal forcing, grounding line depth only weakly impacts maximum submarine melt rate and so does not appear explicitly in the parameterization. Observations support a line plume width w ≈ 200 m (ref. 17 ) so that the final parameterization may be written in terms of the total subglacial discharge Q and thermal forcing TF as ṁ = 0.142 Q 0.31 TF 1. 19 . When splitting out variability due to atmosphere and ocean, we apply the same parameterization but hold one of the factors Q or TF constant at its mean 1979-1993 value.

To estimate the impact of variability in submarine melt rate on glacier retreat, ice discharge and mass loss we use a simple two-stage glacier model adapted from ref. 29 , with the only difference to that study being the addition of submarine melt rate and calving amplification to the calving flux. The model consists of a glacier interior of length L, ice thickness H and surface mass balance P. The glacier interior delivers an ice flux Q to the terminus, which has ice thickness h g , while the grounding line/calving flux is Q g . Following ref. 29 , mass conservation for the interior and terminus regions gives

where t is time. The interior flux is estimated as

L n , in which n = 3 is the Glen's flow law exponent and ν is a constant related to the basal friction. The grounding line/calving flux is written Q g = Ωh γ g + (1 + θ) ṁh g , in which Ω is a constant related to the basal friction and ice softness, γ = 19/4 (ref. 47 ), ṁ is the submarine melt rate and θ is a calving multiplier. Finally, the terminus is assumed to be at flotation h g = -(ρ w /ρ i )b L , in which b L is the water depth at the terminus, ρ w = 1,028 kg m -3 is the density of ocean water and ρ i = 917 kg m -3 is the density of ice. The model applies equally to glaciers with an ice shelf, provided that the ice shelf provides little buttressing to flow, which is a simplifying assumption but affects very few of the glaciers in the dataset.

To initialize the glacier model we assume that during 1979-1993 the glacier is close to steady state (that is, the mean value of dH dt and dL dt during this period is 0). The model then has six remaining free parameters, which can be taken as initial ice velocity v 0 , initial calving front ice thickness h g0 , initial inland ice thickness H 0 , initial glacier length L 0 , bed slope b x and calving rate multiplier θ. All constants in the glacier model can be expressed in terms of these parameters. Rather than trying to model a specific glacier, we instead test two representative values for each of these six free parameters (Supplementary Table 2), giving a set of 2 6 = 64 simulations for each submarine melt rate time series. With this initialization, the model is run forwards from 1979 to 2018, forced by monthly submarine melt rate from each glacier in the dataset. Surface mass balance is held constant throughout. We consider three submarine melt rate time series per glacier: atmosphere varies, ocean varies and both vary. We then have 64 simulations for each of 123 glaciers and for each of 3 submarine melt scenarios (total 23,616 simulations).

The terminus position, L, is a direct output of the simulations while the cumulative sea level contribution per unit width of glacier in year t, SLR, is calculated as SLR = ∫ , where SLR b is the total 1979-2018 sea level contribution in the simulation with both atmospheric and oceanic variability in submarine melting (Extended Data Fig. 4). By this normalization, all simulations that account for both atmospheric and ocean-sourced variability in submarine melting have a normalized retreat of -1 and a sea level contribution of 1 in 2018, while simulations accounting for only atmospheric or oceanic variability alone have a smaller normalized retreat and sea level contribution (Fig. 4). Lastly, simulations are grouped by ice sheet region.