RESEARCH LETTER
10.1029/2019GL086913
Key Points:• Internal variability plays a dominantrole in modulating summertime
Greenland ice sheet melt conditions
at high frequencies (<~20 years)• The North Atlantic Oscillation is thedominant contributor to thismodulation and may act as a weaknegative feedback on future ice melt
• Climate change will affect internalvariability impacts on Greenland by
decreasing its relative importancebut increasing its area of impact
Supporting Information:• Supporting Information S1

Correspondence to:
M. McElroy,
mbm@seas.harvard.edu

Citation:
Sherman, P., Tziperman, E., Deser, C.,& McElroy, M. (2020). Historical andfuture roles of internal atmosphericvariability in modulating summertime
Greenland Ice Sheet melt. GeophysicalResearch Letters, 47, e2019GL086913.
https://doi.org/10.1029/2019GL086913
Received 2 JAN 2020
Accepted 8 MAR 2020
Accepted article online 12 MAR 2020

Historical and Future Roles of Internal Atmospheric
Variability in Modulating Summertime
Greenland Ice Sheet MeltPeter Sherman1, Eli Tziperman1,2, Clara Deser3, and Michael McElroy1,21Department of Earth and Planetary Sciences, Harvard University, Cambridge, MA, USA, 2School of Engineering andApplied Sciences, Harvard University, Cambridge, MA, USA, 3Climate and Global Dynamics Division, National Center
for Atmospheric Research, Boulder, CO, USA

Abstract Understanding how internal atmospheric variability affects Greenland ice sheet (GrIS)
summertime melting would improve understanding of future sea level rise. We analyze the Community
Earth System Model Large Ensemble (CESM‐LE) over 1951–2000 and 2051–2100. We ﬁnd that internalvariability dominates the forced response on short timescales (~20 years) and that the area impacted byinternal variability grows in the future, connecting internal variability and climate change. Unlike priorstudies, we do not assume speciﬁc patterns of internal variability to affect GrIS melting but derive them frommaximum covariance analysis. We ﬁnd that the North Atlantic Oscillation (NAO) is the major source ofinternal atmospheric variability associated with GrIS melt conditions in CESM‐LE and reanalysis, with thepositive phase (NAO+) linked to widespread cooling over the ice sheet. CESM‐LE and CMIP5 project anincrease in the frequency of NAO+ events, suggesting a negative feedback to the GrIS under futureclimate change.
Plain Language Summary Melting of the Greenland ice sheet is an important anticipatedconsequence of human‐induced climate change. The human‐induced signal on the Greenland ice sheet
can be modulated by internal variability (inherent atmospheric “noise” that arises due to nonlinear
processes associated with atmosphere, ocean, and atmosphere‐ocean dynamics). Here, we study historicaland future climate simulations to assess the role of summertime internal variability on changes in melting
conditions in Greenland. We ﬁnd that the impact from internal variability changes in the future, whichsuggests a connection between this atmospheric noise and anthropogenic climate change. We then identify
using a statistical technique the patterns of summertime internal variability that are most connected toGreenland ice sheet melting. We ﬁnd that the North Atlantic Oscillation is the dominant mode of internal
atmospheric circulation variability that affects Greenland ice melt conditions, both in the models andobservations. Changes in Earth's radiative balance attributable to anthropogenic climate change are found toincrease the frequency of positive summertime NAO events in climate models, with signiﬁcant implications
for Greenland's future mass balance. Understanding how these modes of internal variability inﬂuence
Greenland conditions is an important and necessary step in projecting future changes in the Greenland
ice sheet.
1. Introduction©2020. The Authors.
This is an open access article under theterms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

SHERMAN ET AL.

Land‐ice melt is an important expected consequence of climate change and has serious societal implications.
Increased ice melt can weaken the Atlantic Ocean overturning circulation, which has the potential to affecttemperatures around the North Atlantic (Vellinga & Wood, 2002). This has bearing on future changes inboth global temperature and sea level. The growth in atmospheric greenhouse gases has coincided largely
with changes in Greenland's surface temperatures, with an increase of roughly 2 °C over the period 1990
to 2014 (Reeves Eyre & Zeng, 2017). As a consequence of the rising temperatures, land‐ice melt from theGreenland ice sheet (GrIS) has been an important contributor to global sea level rise over the past three decades and is responsible for an increase of approximately 0.33 mm/year of the global 3.2 mm/year from 1993
to 2010 (Vaughan et al., 2013). Despite the recent net melting of the GrIS over the past few decades, previous
literature has noted a pause in mass loss from 2013 to 2014 (Bevis et al., 2019). Understanding what drives
such pauses is an important and necessary step to project future changes in the GrIS.

1 of 10

Geophysical Research Letters10.1029/2019GL086913

Regional climate internal variability arises from processes associated with the nonlinear atmosphere, ocean,
and atmosphere‐ocean dynamics. The North Atlantic Oscillation (NAO), for example, is the dominant modeof large‐scale interannual variability in sea level pressure (SLP) over the North Atlantic (Hurrell, 1995). The
NAO is associated with changes in the strength and position of the jet stream (Woollings & Blackburn, 2012),
affecting temperatures over Greenland and north and central Europe (Chylek et al., 2004). These NAO
effects are often discussed in the context of winter climate; however, the NAO has also been shown to havean important role also in summer, when conditions are conducive to Greenland ice melt (Bevis et al., 2019;
Fettweis et al., 2013; Folland et al., 2009; van Angelen et al., 2014). This connection is reﬂected in both observations (Bevis et al., 2019; Box et al., 2012) and models (Ding et al., 2014; Hahn et al., 2018). Notably, therecent 2013–2014 pause in the west coast GrIS melting coincided with a strong and positive NAO (Bevis
et al., 2019). More generally, Box et al. (2009) attributed the NAO‐GrIS connection over the period 1840–
2007 to an anomalous anticyclone that forms during the negative phase of the NAO, which promotes advection of warm air into the region, with an increased downward ﬂux of radiation and reductions in snowfall
rates. Two important questions are addressed in this study: (1) Which mode of internal variability is mostdirectly implicated in GrIS changes? and (2) Will climate change exacerbate the NAO's impact on GrIS melt?
Here, we analyze summertime (June–July–August) atmospheric conditions from the 35‐member
Community Earth System Model “Large Ensemble” (CESM‐LE; Kay et al., 2015) over the period
1951–2100 to assess how internal variability has contributed to melting of the GrIS in the simulations. We
focus speciﬁcally on summer months because that is when melting conditions are persistent overGreenland. While the NAO may play a more dominant role on climate variability in other seasons, surface
air temperatures seldom reach above freezing in boreal spring and fall even in the future under a high emissions scenario (see http://webext.cgd.ucar.edu/Multi‐Case/CVDP_repository/cesm1.lens_2018‐2100/tas.
mean.mam.png). Rather than assuming that a certain atmospheric variability mode (e.g., the NAO) is
responsible for the response of GrIS melting and studying their temperature composites (Hahn et al., 2018),
we use maximum covariance analysis (MCA) to identify the spatial structure of SLP variability which moststrongly correlates with GrIS melting patterns. We ﬁrst compare the magnitudes of the ensemble mean (i.e.,
forced) response and atmospheric internal variability on the GrIS melting rate and spatial patterns characterized by positive degree days on a monthly basis (PDDs; the cumulative number of degree days above freezing over a month using daily surface air temperatures) as a surrogate for ice melt. How these contributions
could change in the future under a high emissions scenario (RCP8.5) is also assessed. We then determine thedominant variability modes of atmospheric circulation connected to this internal variability and propose
physical mechanisms for how these patterns are manifest in PDD over the GrIS. We address also how theNAO and its impacts on the GrIS are expected to change in the future. Section 2 describes the methods
and data used. Section 3 presents the results. Section 4 summarizes the study and provides a
general discussion.

2. Method
2.1. Data Overview—CESM‐LE, 20CR, and CMIP5
CESM‐LE is an ensemble of fully coupled CESM1 simulations at an approximate spatial resolution of 1° × 1°,
covering the period 1920–2100. Simulations are initialized from the same model state on 1 January 1920,
except for minor perturbations in initial air temperatures (on the order of 10−14 K). These slight differences
between simulations impact the global climate system by altering the sequence of unpredictable internalvariability (Deser et al., 2012). While the sequence of internal variability differs from simulation to simulation, external forcing is consistent for each member of the ensemble: CMIP5 historical radiative forcingsfrom 1920 to 2005 and RCP8.5 from 2006 to 2100. This fact can be exploited to assess the strength of the externally forced signal relative to the uncertainty associated with internal variability within CESM‐LE (i.e.,
Deser et al., 2012; Vega‐Westhoff & Sriver, 2017). Further elaboration of the ensemble is discussed in Kay
et al. (2015). Positive degree days (PDDs) represent a conventional metric for evaluating atmospheric conditions conducive to melting (Reeh, 1989). To study the monthly variability in GrIS changes, we aggregate
PDD to a monthly resolution. Monthly PDD, SLP, temperature at 500 hPa (T500), and winds at 500 hPa(u500 and v500) are considered here. Data at 500 hPa are considered because they represent the atmospheric
dynamics better than surface variables, being less inﬂuenced by the orographic discontinuity between the
GrIS and surrounding waters.
SHERMAN ET AL.

2 of 10

Geophysical Research Letters10.1029/2019GL086913

To verify the historical results derived from CESM‐LE, we use the 20th Century Reanalysis (20CR) version 2(Compo et al., 2006; Compo et al., 2011; Whitaker et al., 2004). 20CR is a global assimilated meteorologicalproduct covering the period 1871–2012 at a spatial resolution of 1° × 1°. The data are generated by assimilating only surface pressures (and using monthly SSTs and sea ice distributions as boundary conditions)through the use of an Ensemble Kalman Filter. Ensemble mean 20CR SLP and PDD data are retrieved ata monthly resolution for comparison with CESM‐LE. Plots of the mean and standard deviation of PDD
for both CESM‐LE and 20CR are shown in supporting information Figures S1 and S2: Both quantities show
general spatial consistency across the two data sets.
To explore model uncertainty, we also use historical (1951–2000) and future (2051–2100) scenario experiments from the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2012). We use theﬁrst ensemble member from the 31 different models that conducted historical RCP4.5 and RCP8.5 simulations. We make use of output from the Climate Variability Diagnostics Package (CVDP; Phillips et al., 2014),
which computes diagnostics of the leading modes of climate variability (data can be found at http://www.
cesm.ucar.edu/working_groups/CVC/cvdp/data‐repository.html). Speciﬁcally, we obtain the NAO index
for each of the 31 CMIP5 models and the CESM‐LE, calculated as the ﬁrst EOF of SLP variability over the
North Atlantic‐European domain (20–80°N, 90°W to 40°E) following Hurrell and Deser (2009), to study
how the frequency of positive NAO events (NAO+) is expected to change in the future under radiative
forcing scenarios.
2.2. Maximum Covariance Analysis (MCA) Implementation
MCA is a linear technique, which has been employed in numerous climate studies (e.g., Bretherton
et al., 1992; Wallace et al., 1992) to assess how two spatio‐temporal variables covary in time. Here, we use
MCA to analyze the covariance between the SLP and PDD to study how internal atmospheric variability
is connected with changes in melting patterns over the ice sheet in CESM‐LE. The model's forced SLP
and PDD components, estimated by averaging over the 35 ensemble members, are removed from the twodata sets, respectively. SLP values are considered only in the North Atlantic‐European domain (90°W to
40°E and 20–80 N), and PDD values are considered only over the GrIS, deﬁned here as the land mass boundary. While temperatures in ice‐covered grid cells will be constrained close to the freezing point and ice‐freecells will not, there should not be much sensitivity in our results to this deﬁnition given that the GrIS does
not differ substantially from the land margin (Lipscomb & Sacks, 2010). The data are area‐weighted to avoid
overestimating the role of densely gridded high latitudes. For each ensemble member, SLP and PDD data
matrices are constructed such that each column corresponds to all grid points of a given ﬁeld, and different
columns correspond to different months. All SLP ensemble members are then concatenated together along
the temporal dimension, forming a matrix with n = 50 years × 3 months × 35 ensemble members = 5,250
columns. The same process is implemented for the PDD matrices. These updated SLP and PDD matrices
are the X and Y matrices used for the MCA. Monthly summer SLP and PDD values are considered fortwo periods: 1951–2000 and 2051–2100.
We also perform MCA on SLP and PDD in 20CR over the period 1951–2000. Given that we do not have multiple realizations of the reanalysis product, we remove a linear temporal trend from the SLP and PDD,
assuming that it approximates the forced response. These data are considered over the same domains as speciﬁed for CESM‐LE. As with the CESM‐LE MCA calculation, the data are area‐weighted. Because of the lackof ensemble members, there are n = 50 years × 3 months = 150 columns for the 20CR analysis. The difference in spatial resolutions between CESM‐LE and 20CR affects the magnitude of the MCA patterns. In orderto account for the higher resolution in CESM‐LE, we normalize the MCA vectors in both CESM‐LE and20CR by dividing by the standard deviation of the vector.
To assess the statistical signiﬁcance of the MCA patterns in CESM‐LE, a bootstrap reshufﬂing technique isused following Wang et al. (2014). With this technique, the X and Y matrices used in the MCA are reconstructed by resampling of the original SLP and PDD matrices at every grid point, that is, by randomly shufﬂing the time steps independently at each grid point (such that the shufﬂing differs from point to point), butusing the same time‐reshufﬂing for both data sets. A “resampled” covariance matrix is then calculated using
the shufﬂed X and Y matrices. The singular values are recalculated with the new covariance matrix. This
process is repeated 1,000 times, and the modiﬁed singular values are compared with the true values to evaluate statistical signiﬁcance. If the magnitude of the true singular value at a given grid cell is larger than the
SHERMAN ET AL.

3 of 10

Geophysical Research Letters10.1029/2019GL086913

Figure 1. Spatial distribution of the time of emergence [years] of PDD for (a) 1951–2000 and (b) 2051–2100, and (c) theirdifference based on CESM‐LE. The time of emergence is deﬁned as the time when the forced PDD signal ﬁrst exceeds twostandard deviations of the internal variability of PDD. Dark purple values in panel (a) indicate regions where the PDD
signal has not emerged in the historical period. The gray shading in panel (c) indicates where this signal becomes emergent in the future period.

90th percentile at the corresponding grid location in the reshufﬂed data sets, the pattern is deemed to bestatistically signiﬁcant at the 90th percentile. Statistical signiﬁcance of the MCA patterns is not evaluatedfor 20CR, given the sparsity of data relative to CESM‐LE.
After performing MCA between SLP and PDD, winds and temperature patterns at 500 hPa associated withthe ﬁrst and second SLP dominant modes of covariability are used to understand the physical mechanisms
governing variability in ice melt conditions. To ﬁnd the wind and temperature anomalies that are connectedto the dominant covariability SVD modes, we compute a “similarity score” by calculating the area‐weighted
Euclidean distance between the vector representing SLP spatial structure at every time step for each ensemble member, and each of the ﬁrst two SLP singular vectors being analyzed. We take the 500 months of the
CESM‐LE out of the total 5,250 months with the smallest Euclidean distances to the SLP pattern as representative of that mode of covariability and use the mean wind and temperature patterns over these 500 months
as representative of the winds and temperatures characteristic of the singular vectors. Wind and temperature
anomalies are found then by subtracting out the mean state. Given the sparsity of data for 20CR relative to
CESM‐LE, we take the 20 time steps with the smallest Euclidean distance out of the total 150 months to thedominant SLP singular vector patterns for 20CR.

3. Results
3.1. Comparing the GrIS Impacts From Internal Variability and the Forced Response
The forced PDD component at a given location is deﬁned as the 10‐year running mean of PDD monthly
anomalies (K days) averaged over all ensemble members over the studied period, while the amplitude of
the internal variability is calculated as two standard deviations (K days) over the 35 monthly anomalies
across the ensemble, as a function of time. The time of emergence (ToE) is deﬁned here as the time in years
since the start of the period examined (1951 or 2051) for the forced response to exceed the amplitude of theinternal variability (Lehner et al., 2017). The ToEs for the periods 1951–2000 and 2051–2100 are shown in
Figure 1. In this historical period, the contribution from internal variability plays an important role in modulating ice melt conditions along the coastlines where the ToE is roughly 10–15 years in southwest
Greenland and up to 30 years in parts of northwest Greenland. Much of the high frequency variability in
PDD can be attributed therefore to unforced components of the climate system. The ToE decreases slightly
to around 5–10 years on Greenland's periphery in the future scenario. The relative importance of internalvariability along the coasts does not have a clear spatial change in the future under a high emissions scenario, but there is a general decline further inland. That being said, the area of impact for this variability
is much greater and extends deeper into the country, as indicated by the dark gray in Figure 1c. This suggests
a connection between the ice sheet internal variability and anthropogenic changes; that is, as the GrISSHERMAN ET AL.

4 of 10

Geophysical Research Letters10.1029/2019GL086913

Figure 2. MCA results for the covariability of summer (JJA) SLP and PDD. The ﬁrst singular vectors in the 20th Century Reanalysis (20CR) are shown for (a) SLP(hPa) and (d) PDD (K days) for the period 1951–2000. The ﬁrst singular vectors in the CESM‐LE are shown in (b) and (e) for 1951–2000, and in (c) and (f) for
2051–2100. The percent variance explained for each singular vector is displayed at the top of each panel. The percentage of the covariability explained between SLP
and PDD (“COVAR”), and the pattern correlation between the MCA mode and the North Atlantic Oscillation (NAO; Figure S3), “R (MCA1,NAO)” are given
between the rows. For example, for 20CR the covariance explained is 56% and the pattern correlation is 0.86. Stippled regions denote locations where the value issigniﬁcant at a 90% conﬁdence level according to a bootstrap reshufﬂing.

recedes, it becomes more susceptible to internal variability in its interior because the future period has days
where the surface temperatures are above freezing, while the historical period does not. The expanded roleof internal variability in the future period can be seen also in the difference of ToE between the two periods(Figure 1c) and in the standard deviation of the PDD (Figure S2). The relative importance of internalvariability in parts of the south Greenland over the high‐elevation ice sheet actually increases in the
future, likely because as temperatures rise closer to 0 °C, temperature ﬂuctuations attributable to internalvariability can determine whether the atmospheric conditions are conducive to melting or freezing. Thechanging role of internal variability on the GrIS in a warming climate is analogous to the situation foundfor Arctic sea ice, where models have shown that as sea ice thins, it becomes more susceptible toatmospheric variability (Kay et al., 2011). Such a result may be intuitive, but it emphasizes the importanceof understanding the connection between forced and unforced components to constrain future GrIS
changes.
3.2. Dominant Modes of Variability
To identify the spatial structure of SLP internal variability that leads to strong GrIS melting, we employ MCA
between SLP and PDD as described in section 2.2. The ﬁrst singular vectors for 1951–2000 are shown for
20CR (Figures 2a and 2d) and from 1951–2000 and 2051–2100 for CESM‐LE (Figures 2b, 2c, 2e, and 2f).
In the historical period, the ﬁrst singular vector captures 34% of the SLP variance in 20CR and 46% in
CESM‐LE; in the future scenario, it accounts for 44% of the SLP variance in CESM‐LE. This SLP patternis optimally correlated with a cooling (PDD) pattern along southwest Greenland in the historic record,
which extends along the entire coastline in the future scenario, as shown by the ﬁrst PDD singular vectors(Figures 2d–2f). In 20CR, this cooling pattern explains 70% of the PDD variability, and 63% and 54% in
CESM‐LE for the historical and future periods, respectively. The MCA‐derived melt patterns optimally
SHERMAN ET AL.

5 of 10

Geophysical Research Letters10.1029/2019GL086913

covarying with this SLP pattern, shown by the ﬁrst PDD singular values, are consistent over both thehistorical periods and future scenario, but with a more widespread melt response in the future scenario(Figures 2d–2f). The SLP patterns indicate that an anomalous low‐pressure system centered over the
Irminger Sea coincides with an anomalous high‐pressure system to the south. This spatial pattern closely
resembles the NAO (Figure S3), and this is validated by the strong pattern correlation coefﬁcients withthe ﬁrst EOF of SLP: 0.86 in 20CR, 0.85 in CESM‐LE historical, and 0.83 in CESM‐LE future. There are somediscrepancies between the spatial patterns in the reanalysis and model products, speciﬁcally in the sign ofthe pressure anomalies over Africa (Figures 2a–2c). However, given the large distance from Greenland, it
is unlikely that these pressure anomalies have a signiﬁcant direct bearing on melt conditions over the
GrIS. The second singular vector derived from the MCA explains roughly 15% of the SLP variability acrossthe three data sets and resembles closely the East Atlantic (EA) pattern (Barnston & Livezey, 1987;
Figure S3). While the EA pattern is the second leading mode of SLP variability, it does not have a clear melting pattern associated with it across the reanalysis and model results and explains only a small portion of the
PDD (on the order of 5%).
3.3. Physical Mechanism
Now that a connection between GrIS melt (PDD) patterns and corresponding patterns of atmospheric variability (SLP) has been established, we will discuss the physical mechanisms connecting the speciﬁc SLP patterns represented by the ﬁrst MCA singular vectors to the corresponding PDD patterns. For this purpose, we
identify the months during which SLP is close to the ﬁrst MCA SLP singular vector, that is, using the similarity score explained in section 2.2. North Atlantic T500 anomalies are shown in Figures 3a–3c for the ﬁrst
MCA mode of variability. The T500 anomaly patterns in the 20CR and CESM‐LE historical periods are largely consistent with each other, featuring negative values over Greenland and the surrounding ocean. The
amplitudes, however, are larger in the model compared to reanalysis, consistent with the larger amplitude
SLP anomaly patterns (recall Figures 2a and 2b). Some of the discrepancy between CESM‐LE and 20CR
may be a consequence of the small sample size of NAO events in 20CR (roughly 30 months with an MCA
index >1), but there is also generally larger summertime air temperature variability over Greenland in
CESM‐LE relative to the observational record (see http://webext.cgd.ucar.edu/Multi‐Case/CVDP_repository/cesm1.lens_1920‐2018/tas.stddev.jja.png). The temperature pattern noted for the historical period is
also found in the CESM‐LE future scenario, and radiative balance provides a plausible mechanism for thestrong cooling response over Greenland in connection with the ﬁrst MCA mode in both the historical andfuture periods. The anomalous cyclone that forms over Greenland during the positive phase of this mode
encourages the ascension of air over the region, which in turn increases the cloud fraction(Figures 3d–3f). The increased cloudiness associated with this mode causes negative net radiation anomalies
at the surface (Figures 3g–3i). We suggest therefore that coupling between atmospheric circulation andradiative balance provides a mechanism that causes the anomalously low melt patterns associated withthe positive phase of the ﬁrst MCA mode found in both 20CR and CESM‐LE (Figure 2). However, other processes could also play important roles in local melt conditions. For example, anomalous southward advection of cold air over the Labrador Sea during the positive phase of this mode could explain theparticularly cool anomalies west of Greenland (Figures 3a–3c). As this mode has an SLP signal similar to
the NAO, this result agrees with previous literature, for example, Holland et al. (2008), who employed observational data to note a major change in subsurface ocean temperatures and wind patterns in winter thatcoincided with a shift in 1995/1996 from prolonged NAO+ conditions to prolonged NAO− conditions thereafter. Additionally, this result agrees with the low mass loss anomaly along Greenland's coastline in 2014
during a strong NAO+ (Bevis et al., 2019). Unlike these studies, we do not assume the NAO to be a player
in GrIS melting but derive it from MCA. The above results show that the NAO is the major northern hemisphere variability mode affecting net GrIS melting patterns.
3.4. Implications Under Future Climate Change
We next attempt to identify changes to the summer NAO under future climate change scenarios, which will,in turn, affect the GrIS melting patterns that are of interest in this work. For this purpose, we analyze thehistorical and future CMIP5 simulations from 31 models (Figure 4a) and CESM‐LE (Figure 4b) and construct probability density functions (PDFs) of the summer NAO index (in units of standard deviations) for
each model during 1951–2000 and 2051–2100 in both RCP4.5 and RCP8.5 scenarios (Figure 4a). In
SHERMAN ET AL.

6 of 10

Geophysical Research Letters10.1029/2019GL086913

Figure 3. North Atlantic JJA temperature anomalies (K) at 500 hPa (ΔT500) associated with the ﬁrst MCA mode of variability for (a) 1951–2000 in the 20th Century
Reanalysis (20CR) and (b) 1951–2000 and (c) 2051–2100 in the CESM‐LE; 500 hPa wind anomalies (m/s) associated with this mode are overlaid as vectors
2
in each panel (reference arrow of 5 m/s displayed at the top right). (d–f) Total cloud fraction anomalies (%) and (g–i) net surface radiation anomalies (W/m )associated with this mode.

transitioning from the historical to future simulations, there is a noticeable shift in the PDF toward more
frequent positive NAO events. Speciﬁcally, the frequency of a positive NAO event in the historical CMIP5
database is 51%: It is 55% in the future RCP4.5 scenario and 64% in RCP8.5. This trend persists even forstrong NAO events; the frequency of events with a summer NAO index greater than 1 increases from 17%
in the CMIP5 historical to 20% in RCP4.5 and 23% in RCP8.5. Increasing the frequency of positivesummer NAO index would have obvious implications for GrIS changes, as the positive phase coincides
with cooling conditions (section 3.2). This suggests that the NAO could act as a negative feedback to

SHERMAN ET AL.

7 of 10

Geophysical Research Letters10.1029/2019GL086913

Figure 4. (a) Probability density functions (PDF) of the summer NAO index for each member of the 31‐model CMIP5
ensemble during: (light green) 1951–2000, (light blue) 2051–2100 under RCP4.5, and (light red) 2051–2100 under
RCP8.5. (b) PDFs of the summer NAO index for each member of CESM‐LE during: (light green) 1951–2000 and (light red)
2051–2100 under RCP8.5. Thick dark curves indicate the ensemble means of the PDFs for each group.

climate change for the GrIS. There is notable inter‐model spread in these NAO changes; however, the NAO
trend dominates for high radiative forcing scenarios. Only one model has a mean negative summer NAO
index for RCP8.5. This can be seen within individual models as well; the mean summer NAO index is
positive for 33 of the 35 CESM‐LE ensemble members (Figure 4b). The possible increase in the frequencyof positive summer NAO events in the future could thus play an even greater role on future GrIS changes.
It is however noteworthy that the observational record does not show a discernible long‐term trend in thesummer NAO over the last century (Iles & Hegerl, 2017).

4. Discussion and Conclusions
The main purpose of this study was to assess how internal variability has inﬂuenced GrIS changes, and whatmodes of atmospheric variability matter most for GrIS melting patterns. We ﬁnd that internal variability
dominates the forced signal associated with GrIS melt on short timescales (roughly 15 years in southwest
Greenland and up to 30 years in the northwest) and that this internal variability expands into the interiorof the GrIS as temperatures increase in the future. These results suggest that internal variability plays a rolein ice sheet changes and that a greater area of the GrIS will be susceptible to internal variability as anthropogenic climate change induces further warming. Through the implementation of maximum covarianceanalysis, we determined the dominant modes of covariability between internal atmospheric circulation
variability and conditions conducive to ice melt over Greenland. From our analysis, we found that a large
portion of the melt variability (>54%) can be explained by the North Atlantic Oscillation, corresponding
to the ﬁrst EOF of the SLP variability. Notably, we detect this pattern also using a simpliﬁed analysis correlating PSL and PDD over Greenland (Figure S4). In its positive phase, this pattern is connected to cooler temperatures over most of Greenland during boreal summers, likely attributable to the coupling between
atmospheric circulation and radiation due to cloudiness changes—explaining the strong melt response
shown during this mode of variability. The increase in GrIS melting in 2011 and subsequent slowdown in

SHERMAN ET AL.

8 of 10

Geophysical Research Letters10.1029/2019GL086913

mass loss from 2013 to 2014 coincides with a negative‐to‐positive shift of the summertime NAO (Figure S5),
somewhat anecdotally conﬁrming our results.
Despite the interesting ﬁndings, they are limited to a single model Large Ensemble (the CESM‐LE). As wehave shown, this model simulates realistic patterns of internal atmospheric circulation variability and associated effects over the GrIS. However, future work should investigate the role of internal variability on the
GrIS in other models with large ensembles for intercomparison and enhanced understanding. We haveaddressed this in part by studying the CMIP5 ensemble (Figure 4a). From this, we ﬁnd that the frequencyof positive summer NAO events increases in most models in the future, suggesting that changes in the
NAO may partially oppose the human‐induced GrIS melting trend.
Acknowledgments
We thank the anonymous reviewers fortheir constructive comments andsuggestions. P. S. and M. M. aresupported by the Harvard Global
Institute. E. T. is supported by the NSF
Climate Dynamics program, joint
NSF/NERC Project Grant
AGS‐1924538, and by the Harvard
Global Institute, and the Weizmann
Institute is thanked for its hospitality
during parts of this work. The CESM
project is supported primarily by the
National Science Foundation (NSF).
This material is based on work
supported by the National Center for
Atmospheric Research (NCAR), which
is a major facility sponsored by the NSF
under Cooperative Agreement 1852977.
Support for the 20th Century
Reanalysis Project data set is provided
by the U.S. Department of Energy,
Ofﬁce of Science Innovative and Novel
Computational Impact on Theory andExperiment (DOE INCITE) program,
and Ofﬁce of Biological andEnvironmental Research (BER), and by
the National Oceanic and Atmospheric
Administration Climate Program
Ofﬁce. P. S. visited NCAR while
undertaking parts of this study. All data
used in this paper are available inpublic repositories. All meteorological
output from the CESM‐LE is publicly
available at the NCAR Climate Data
Gateway (https://www.
earthsystemgrid.org/dataset/ucar.cgd.
ccsm4.CESM_CAM5_BGC_LE.html).
20CR output is available from the
NOAA ESRL website (https://www.
esrl.noaa.gov/psd/data/20thC_Rean/).
The CMIP5 and CESM‐LE NAO indiceswere retrieved from the CVDP datarepository (http://www.cesm.ucar.edu/
working_groups/CVC/cvdp/
data‐repository.html).

SHERMAN ET AL.

References
Barnston, A., & Livezey, R. E. (1987). Classiﬁcation, seasonality, and persistence of low‐frequency circulation patterns. Monthly Weather
Review, 115, 1083–1126. https://doi.org/10.1175/1520‐0493(1987)115<1083:CSAPOL>2.0.CO;2
Bevis, M., Harig, C., Khan, S. A., Brown, A., Simons, F. J., Willis, M., et al. (2019). Accelerating changes in ice mass within Greenland, andthe ice sheet's sensitivity to atmospheric forcing. PNAS, 116(6), 1934–1939. https://doi.org/10.1073/pnas.1806562116
Box, J. E., Fettweis, X., Stroeve, J. C., Tedesco, M., Hall, D. K., & Steffen, K. (2012). Greenland ice sheet albedo feedback: Thermodynamics
and atmospheric drivers. The Cryosphere, 6, 821–839. https://doi.org/10.5194/tc‐6‐821‐2012
Box, J. E., Yang, L., Bromwich, D. H., & Bai, L.‐S. (2009). Greenland ice sheet surface air temperature variability: 1840–2007. Journal of
Climate, 22, 4029–4049.
Bretherton, C. S., Smith, C., & Wallace, J. M. (1992). An intercomparison of methods for ﬁnding coupled patterns in climate data. Journal of
Climate, 5, 541–560. https://doi.org/10.1175/1520‐0442(1992)005<0541:AIOMFF>2.0.CO;2
Chylek, P., Box, J. E., & Lesins, G. (2004). Global warming and the Greenland ice sheet. Climate Change, 63, 201–221. https://doi.org/
10.1023/B:CLIM.0000018509.74228.03
Compo, G. P., Whitaker, J. S., & Sardeshmukh, P. D. (2006). Feasibility of a 100 year reanalysis using only surface pressure data. Bulletin of
the American Meteorological Society, 87, 175–190.
Compo, G. P., Whitaker, J. S., Sardeshmukh, P. D., Matsui, N., Allan, R. J., Yin, X., et al. (2011). The twentieth century reanalysis project.
Quarterly Journal of the Royal Meteorological Society, 137(654), 1–28. https://doi.org/10.1002/qj.776
Deser, C., Phillips, A., Bourdette, V., & Teng, H. Y. (2012). Uncertainty in climate change projections: The role of internal variability.
Climate Dynamics, 38, 527–546. https://doi.org/10.1007/s00382‐010‐0977‐x
Ding, Q., Wallace, J. M., Battisti, D. S., Steig, E. J., Gallant, A. J., Kim, H. J., & Geng, L. (2014). Tropical forcing of the recent rapid Arctic
warming in northeastern Canada and Greenland. Nature, 509(7499), 209–212. https://doi.org/10.1038/nature13260
Fettweis, X., Hanna, E., Lang, C., Belleﬂamme, A., Erpicum, M., & Gallée, H. (2013). Brief communication “important role of the
mid‐tropospheric atmospheric circulation in the recent surface melt increase over the Greenland ice sheet”. The Cryosphere, 7, 241–248.
https://doi.org/10.5194/tc‐7‐241‐2013
Folland, C. K., Knight, J., Linderholm, H. W., Fereday, D., Ineson, S., & Hurrell, J. W. (2009). The summer North Atlantic Oscillation: Past,
present and future. Journal of Climate, 22, 1082–1103. https://doi.org/10.1175/2008JCLI2459.1
Hahn, L., Ummenhofer, C. C., & Kwon, Y.‐O. (2018). North Atlantic natural variability modulates emergence of widespread Greenland
melt in a warming climate. Geophysical Research Letters, 45, 9171–9178. https://doi.org/10.1029/2018GL079682
Holland, D. M., Thomas, R. H., De Young, B., Ribergaard, M. H., & Lyberth, B. (2008). Acceleration of Jakobshavn Isbrae triggered by warm
subsurface ocean waters. Nature Geoscience, 1, 659–664. https://doi.org/10.5670/oceanog.2016.98
Hurrell, J. W. (1995). Decadal trends in the North Atlantic Oscillation: Regional temperatures and precipitation. Science, 269(5224),
676–679. https://doi.org/10.1126/science.269.5224.676
Hurrell, J. W., & Deser, C. (2009). North Atlantic climate variability: The role of the North Atlantic Oscillation. Journal of Marine Systems,
78, 28–41. https://doi.org/10.1016/j.jmarsys.2008.11.026
Iles, C., & Hegerl, G. (2017). Role of the North Atlantic Oscillation in decadal temperature trends. Environmental Research Letters, 12,
114,010. https://doi.org/10.1088/1748‐9326/aa9152
Kay, J. E., Deser, C., Phillips, A., Mai, A., Hannay, C., Strand, G., et al. (2015). The Community Earth System Model (CESM) large ensemble
project: A community resource for studying climate change in the presence of internal climate variability. Bulletin of the American
Meteorological Society, 96(8), 1333–1349. https://doi.org/10.1175/BAMS‐D‐13‐00255.1
Kay, J. E., Holland, M. M., & Jahn, A. (2011). Inter‐annual to multi‐decadal Arctic Sea ice extent trends in a warming world. GeophysicalResearch Letters, 38, L15708. https://doi.org/10.1029/2011GL048008
Lehner, F., Deser, C., & Terray, L. (2017). Toward a new estimate of “time of emergence” of anthropogenic warming: Insights fromdynamical adjustment and a large initial‐condition model ensemble. Journal of Climate, 30, 7739–7756. https://doi.org/10.1175/
JCLI‐D‐16‐0792.1
Lipscomb, W., & Sacks, W. (2010). The CESM land ice model: Documentation and user's guide.
Phillips, A. S., Deser, C., & Fasullo, J. (2014). A new tool for evaluating modes of variability in climate models. Eos, 95, 453–455. https://doi.
org/10.1002/2014EO490002
Reeh, N. (1989). Parameterization of melt rate and surface temperature on the Greenland ice sheet. Polarforschung, 59, 113–128.
Reeves Eyre, J. E. J., & Zeng, X. (2017). Evaluation of Greenland near surface air temperature datasets. The Cryosphere, 11, 1591–1605.
https://doi.org/10.5194/tc‐11‐1591‐2017
Taylor, K. E., Stouffer, R. J., & Meehl, G. A. (2012). An overview of CMIP5 and the experiment design. Bulletin of the American
Meteorological Society, 93, 485–498. https://doi.org/10.1175/BAMS‐D‐11‐00094.1
van Angelen, J. H., Lenaerts, J. T. M., van den Broeke, M. R., & Wouters, B. (2014). Contemporary (1960–2012) evolution of theclimate and surface mass balance of the Greenland ice sheet. Surveys in Geophysics, 35, 1155–1174. https://doi.org/10.1007/
s10712‐013‐9261‐z
Vaughan, D. G., et al. (2013). Observations: Cryosphere. In T. F. Stocker, D. Qin, G.‐K. Plattner, M. Tignor, S. K. Allen, J. Boschung, A.
Nauels, Y. Xia, V. Bex, & P. M. Midgley (Eds.), Climate change 2013: The physical science basis. Contribution of working group I to the ﬁfth

9 of 10

Geophysical Research Letters10.1029/2019GL086913

assessment report of the intergovernmental panel on climate change (Chap. 4, pp. 317–382). Cambridge, UK and New York, NY:
Cambridge University Press.
Vega‐Westhoff, B., & Sriver, R. L. (2017). Analysis of ENSO's response to unforced variability and anthropogenic forcing using CESM.
Scientiﬁc Reports, 7(1), 1–10. https://doi.org/10.1038/s41598‐017‐18459‐8
Vellinga, M., & Wood, R. A. (2002). Global climatic impacts of a collapse of the Atlantic thermohaline circulation. Climate Change, 54,
251–267.
Wallace, J. M., Smith, C., & Bretherton, C. S. (1992). Singular value decomposition of wintertime sea surface temperature and 500‐mb
height anomalies. Journal of Climate, 5, 561–576. https://doi.org/10.1175/1520‐0442(1992)005<0561:SVDOWS>2.0.CO;2
Wang, Y.‐H., Magnusdottir, G., Stern, H., Tian, X., & Yu, Y. (2014). Uncertainty estimates of the EOF‐derived North Atlantic Oscillation.
Journal of Climate, 27, 1290–1301. https://doi.org/10.1175/JCLI‐D‐13‐00230.1
Whitaker, J. S., Compo, G. P., Wei, X., & Hamill, T. M. (2004). Reanalysis without radiosondes using ensemble data assimilation. Monthly
Weather Review, 132, 1190–1200.
Woollings, T., & Blackburn, M. (2012). The North Atlantic jet stream under climate change and its relation to the NAO and EA patterns.
Journal of Climate, 25, 886–902. https://doi.org/10.1175/JCLI‐D‐11‐00087.1

SHERMAN ET AL.

10 of 10