Historical and Future Roles of Internal Atmospheric Variability in Modulating Summertime Greenland Ice Sheet Melt

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 find that internal variability dominates the forced response on short timescales (~20 years) and that the area impacted by internal variability grows in the future, connecting internal variability and climate change. Unlike prior studies, we do not assume specific patterns of internal variability to affect GrIS melting but derive them from maximum covariance analysis. We find that the North Atlantic Oscillation (NAO) is the major source of internal atmospheric variability associated with GrIS melt conditions in CESM‐LE and reanalysis, with the positive phase (NAO+) linked to widespread cooling over the ice sheet. CESM‐LE and CMIP5 project an increase in the frequency of NAO+ events, suggesting a negative feedback to the GrIS under future climate change.


Introduction
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 affect temperatures around the North Atlantic (Vellinga & Wood, 2002). This has bearing on future changes in both 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 the Greenland 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.
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 mode of 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 have an 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 reflected in both observations (Bevis et al., 2019;Box et al., 2012) and models (Ding et al., 2014;Hahn et al., 2018). Notably, the recent 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  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 flux of radiation and reductions in snowfall rates. Two important questions are addressed in this study: (1) Which mode of internal variability is most directly implicated in GrIS changes? and (2) Will climate change exacerbate the NAO's impact on GrIS melt?

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 internal variability (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 forcings from 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 influenced by the orographic discontinuity between the GrIS and surrounding waters.
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 meteorological product 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 at a 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  and future (2051-2100) scenario experiments from the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2012). We use the first 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). Specifically, we obtain the NAO index for each of the 31 CMIP5 models and the CESM-LE, calculated as the first 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.

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 two data 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, defined here as the land mass boundary. While temperatures in ice-covered grid cells will be constrained close to the freezing point and ice-free cells will not, there should not be much sensitivity in our results to this definition 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 field, 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 for two 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 specified for CESM-LE. As with the CESM-LE MCA calculation, the data are area-weighted. Because of the lack of 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 order to account for the higher resolution in CESM-LE, we normalize the MCA vectors in both CESM-LE and 20CR by dividing by the standard deviation of the vector.
To assess the statistical significance of the MCA patterns in CESM-LE, a bootstrap reshuffling technique is used 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 shuffling the time steps independently at each grid point (such that the shuffling differs from point to point), but using the same time-reshuffling for both data sets. A "resampled" covariance matrix is then calculated using the shuffled X and Y matrices. The singular values are recalculated with the new covariance matrix. This process is repeated 1,000 times, and the modified singular values are compared with the true values to evaluate statistical significance. If the magnitude of the true singular value at a given grid cell is larger than the 10.1029/2019GL086913

Geophysical Research Letters
90th percentile at the corresponding grid location in the reshuffled data sets, the pattern is deemed to be statistically significant at the 90th percentile. Statistical significance of the MCA patterns is not evaluated for 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 with the first and second SLP dominant modes of covariability are used to understand the physical mechanisms governing variability in ice melt conditions. To find the wind and temperature anomalies that are connected to 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 first 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 the dominant SLP singular vector patterns for 20CR.

Comparing the GrIS Impacts From Internal Variability and the Forced Response
The forced PDD component at a given location is defined 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 defined 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 the internal 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 internal variability 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 GrIS 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 role of 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 internal variability 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 fluctuations attributable to internal variability can determine whether the atmospheric conditions are conducive to melting or freezing. The changing role of internal variability on the GrIS in a warming climate is analogous to the situation found for Arctic sea ice, where models have shown that as sea ice thins, it becomes more susceptible to atmospheric variability (Kay et al., 2011). Such a result may be intuitive, but it emphasizes the importance of understanding the connection between forced and unforced components to constrain future GrIS changes.

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 first 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 first 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 pattern is 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 first 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  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 is significant at a 90% confidence level according to a bootstrap reshuffling.

SHERMAN ET AL.
covarying with this SLP pattern, shown by the first PDD singular values, are consistent over both the historical 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 coefficients with the first EOF of SLP: 0.86 in 20CR, 0.85 in CESM-LE historical, and 0.83 in CESM-LE future. There are some discrepancies between the spatial patterns in the reanalysis and model products, specifically in the sign of the pressure anomalies over Africa (Figures 2a-2c). However, given the large distance from Greenland, it is unlikely that these pressure anomalies have a significant direct bearing on melt conditions over the GrIS. The second singular vector derived from the MCA explains roughly 15% of the SLP variability across the 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%).

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 the historical and future CMIP5 simulations from 31 models (Figure 4a)

Geophysical Research Letters
transitioning from the historical to future simulations, there is a noticeable shift in the PDF toward more frequent positive NAO events. Specifically, 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 for strong 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 positive summer 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

Geophysical Research Letters
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 frequency of 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 the summer NAO over the last century (Iles & Hegerl, 2017).

Discussion and Conclusions
The main purpose of this study was to assess how internal variability has influenced GrIS changes, and what modes of atmospheric variability matter most for GrIS melting patterns. We find 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 interior of the GrIS as temperatures increase in the future. These results suggest that internal variability plays a role in 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 covariance analysis, 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 first EOF of the SLP variability. Notably, we detect this pattern also using a simplified 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 mass loss from 2013 to 2014 coincides with a negative-to-positive shift of the summertime NAO ( Figure S5), somewhat anecdotally confirming our results.
Despite the interesting findings, they are limited to a single model Large Ensemble (the CESM-LE). As we have 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 have addressed this in part by studying the CMIP5 ensemble (Figure 4a). From this, we find that the frequency of 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.