Greenland has experienced large changes since the last glacial with its summit warming by approximately 21°C, average accumulation rates tripling, and annual amplitudes of temperature and accumulation seemingly declining. The altered seasonal cycle of accumulation has been attributed to a combination of the large-scale dynamical response of the North Atlantic storm track to surface boundary conditions and the modulation of moisture availability due to changes in winter sea ice cover. Using atmospheric simulations of preindustrial and glacial climate, the contributions of these two mechanisms are evaluated. Estimates of moisture source footprints make it possible to distinguish between long-range transport related to the storm track and regional transport from the ocean surface near Greenland. It is found that the contribution of both mechanisms varies significantly with the background climate. With greater ice cover and the North Atlantic storm track locked to the topographically enhanced stationary wave during the glacial, seasonal migration of the sea ice edge becomes relatively important in controlling moisture availability. In contrast, the preindustrial simulation has relatively greater transient eddy activity and is less moisture limited by sea ice extent, so accumulation is more strongly related to synoptic variability in the North Atlantic. These results highlight how changes in atmospheric circulation and sea ice together explain the shifts in annual mean and seasonal moisture supply to Greenland. Also discussed are some implications of the inferred narrow source distribution of accumulation during the glacial for the interpretation of stable isotopes derived from the central Greenland ice cores.
Greenland’s temperature and accumulation was greatly altered during the Last Glacial Maximum (LGM; ~21 kyr ago). In the modern epoch, near-surface air temperatures in the vicinity of Summit Camp in Greenland average approximately −29°C, with a winter-to-summer amplitude of approximately 24°C (Steffen and Box 2001). Inverse models constrained by measurements of in situ borehole temperature from the Greenland Ice Core Project (GRIP; Dansgaard et al. 1993; Dahl-Jensen et al. 1998) and the Greenland Ice Sheet Project 2 (GISP2; Cuffey and Clow 1997) indicate that LGM temperatures were 21°–23°C colder than modern values. These estimates are similar to a combined δ18Oice and nitrogen gas isotope estimate from the more northerly North Greenland Ice Core Project (NGRIP; Andersen et al. 2006) ice core of approximately 19°C cooling (Kindler et al. 2013). It is difficult to empirically constrain the annual amplitude of LGM surface temperature from observations because of temporal smoothing and variable interdependence between the seasonality of temperature and accumulation (Sime et al. 2008), but modeling studies indicate that the winter-to-summer temperature difference was approximately 39°C (Li et al. 2010), giving an annual cycle of temperature approximately 60% larger in amplitude than under modern conditions.
Modern observations of accumulation near Summit Camp, Greenland indicate an annual average of approximately 26 cm yr−1 liquid water equivalent (van der Veen and Bolzan 1999) and a modest seasonal amplitude with approximately 3.0 cm yr−1 greater accumulation during summer than winter (Shuman et al. 1995). Annual layer thickness measurements from GISP2 indicate 75% less annual accumulation during the LGM (Cuffey and Clow 1997), once corrected for flow thinning, and estimates from NGRIP indicate 67% less annual accumulation (Andersen et al. 2006). Differences in these estimates likely reflect both uncertainties in the estimates and regional variations in accumulation. As with the annual cycle of temperature, past seasonal accumulation rates have not been directly constrained, but simulations suggest a 2.7–10-fold increase in the ratio of summer-to-winter accumulation [defined as June–August (JJA) for summer and December–February (DJF) for winter for consistency with previous studies] (Werner et al. 2000; Krinner and Werner 2003; Li et al. 2005; Langen and Vinther 2009). Note that the reduction of accumulation during the LGM is greater than would be expected from the thermodynamic response to temperature alone, suggesting that changes in dynamics are also important (Kapsner et al. 1995).
Snow line changes throughout the North Atlantic sector also suggest a greater seasonality during the last glacial. On the basis of lateral moraine elevations in Milne Land and isostatic uplift, Denton et al. (2005) estimated that the winter season was approximately 24°C cooler during the Younger Dryas (12.8 kyr ago), in order to be consistent with reconstructed glacier equilibrium line altitudes. Denton et al. (2005) also used similar methods to show that glacial winters in Europe were colder than those of the Holocene by a far greater margin than summers, indicating that more extreme seasonality was not limited to Greenland.
Sea ice cover is a well-recognized influence on the seasonal cycle of accumulation in Greenland (Steig et al. 1994; Krinner et al. 1997; Werner et al. 2000, 2001; Krinner and Werner 2003; Li et al. 2005; Jouzel et al. 2007). Some earlier studies found diffuse moisture sources that were not directly controlled by surface conditions in the North Atlantic (Charles et al. 1994; Armengaud et al. 1998; Werner et al. 2001; Langen and Vinther 2009), but higher-resolution simulations indicate the importance of localized moisture sources (Li et al. 2005; Sodemann et al. 2008b,a; Lewis et al. 2013), especially from the Nordic seas. Consistent with the Greenland and Nordic seas having had greater seasonal sea ice cover during glacial periods, changes in sea ice cover have been proposed as a mechanism for modulating accumulation in central Greenland (Li et al. 2010). It has not been clear, however, whether anomalies in ice extent influence temperature and accumulation in Greenland directly via surface moisture and heat fluxes or indirectly through an interaction with the large-scale circulation.
A second mechanism for influencing seasonality is suggested by analyses of storm track activity in simulations of modern and LGM climate. Observed modern accumulation in central Greenland is relatively constant through the year, with the primary seasonal signal stemming from greater snowfall in late summer to early fall during the intensification and migration of the North Atlantic storm track. The transient eddy activity associated with baroclinic development is diminished in LGM simulations despite there being a larger meridional temperature gradient throughout the North Atlantic sector (Li and Battisti 2008; Donohoe and Battisti 2009; Li et al. 2010). An examination of eddy life cycles conditional on initial properties in the storm track entrance region shows that individual eddies are not behaving contrary to basic theory; there are simply fewer of them, their distribution is skewed toward smaller initial sizes when compared with the modern climatology, and they are advected more rapidly through the baroclinic zone (Donohoe and Battisti 2009). That more recent studies identify transient eddy activity as important, unlike earlier studies (e.g., Kageyama et al. 1999), may relate to the use of higher-resolution models or different surface boundary conditions. The more recent studies (Li and Battisti 2008; Donohoe and Battisti 2009; Li et al. 2010) use the ICE-5G reconstruction specified in the Paleoclimate Modelling Intercomparison Project phase 2 (PMIP2) 21k experiment, which features a rougher Laurentide Ice Sheet than the ICE-4G reconstruction used in the PMIP1 simulations (Peltier 2004; Justino et al. 2006).
The two mechanisms—regional sea ice cover and storm track structure—are not independent of one another. Open water in the Nordic seas introduces a strong zonal asymmetry in surface boundary conditions that may induce a nonlocal dynamic response, in addition to a localized surface flux enhancement. On the other hand, synoptic variability within the storm track may not directly feed Greenland with moisture, perhaps only creating the necessary conditions for localized fluxes to occur. Regardless of what provides the moisture, the two mechanisms are intrinsically linked by the relationship between surface gradients and the vertical structure of the midlatitude westerlies. Here we explore Greenland accumulation using a Lagrangian moisture analysis of modern and LGM simulations in order to characterize the extent to which the controls on accumulation events are rooted in large-scale dynamics or regional moisture availability.
2. Dynamical model simulations
We simulate modern and glacial climate conditions using the National Center for Atmospheric Research (NCAR) Community Atmosphere Model (CAM), version 3 (CAM3; Collins and Rasch 2004). The surface topography, land use, sea ice cover, and sea surface temperatures (SSTs) are all specified, though similar results for annual averages and the seasonal cycle would be expected from coupled simulations because the boundary conditions are derived from fully coupled simulations using the Community Climate System Model, version 3 (CCSM3), of which CAM3 is the atmospheric component. The model is integrated using the Eulerian spectral dynamical core at T42 resolution, while the surface boundary conditions are specified at 64 × 128 gridpoint resolution (approximately 2.8°). A spinup is performed for 3 years, and the model is then integrated for another 20 years, with 3-h output saved for winds, specific humidity, temperature, precipitation, and vertical fluxes of heat and momentum. The prescribed boundary conditions make a run of this length sufficient for obtaining stable results. Additional perturbation experiments and the use of nonclimatological SSTs and sea ice cover could be useful extensions in future work.
Simulations and boundary conditions
For the preindustrial (PI) simulation, the standard model configuration is used with the default set of modern climatological SSTs and sea ice cover. The climatology was derived from the Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) optimal interpolation dataset from 1950 to 1981 and the Reynolds EOF dataset between 1981 and 2001. Other parameters were set to match those of the PMIP2 control experiment, with insolation forcing for the year 1950 and concentrations of CO2 at 280 ppm, CH4 at 760 ppb, and N2O at 270 ppb (Braconnot et al. 2007a). For the LGM, PMIP2 21k boundary conditions were used: topography and ice sheet extent are from ICE-5G (Peltier 2004) for 21 kyr ago. Land surface properties, SSTs, and sea ice cover were derived from the fully coupled PMIP2 21k experiment integration of CCSM3 (Otto-Bliesner et al. 2006), which is a 300-yr integration to quasi equilibrium, with the last 30 yr of the output used to generate the climatological boundary conditions. In that simulation, changes in sea ice stem from a combination of oceanic variability, freshwater forcing, atmospheric variability, radiative forcing, and internal dynamics of the ice itself; however, the resulting sea ice climatology is held fixed here. Concentrations for CO2 are set to 185 ppm, CH4 to 350 ppb, and N2O to 200 ppb.
The climate simulated by CCSM for both the PI and LGM scenarios has previously been compared with other simulations from the PMIP2 ensemble and with proxy data (Braconnot et al. 2007a,b; Otto-Bliesner et al. 2009; Pausata et al. 2009). In most metrics, CCSM3 is near the model median, though it suffers from a well-known wet bias wherein Arctic precipitation is overestimated relative to other models and observations (Pan et al. 2010). The wet bias has been traced to an interaction between cooling via cloud radiative forcing and excess diabatic heating. Studies of synoptic activity nevertheless show that the frequency of a given synoptic pattern is well simulated (Schuenemann and Cassano 2009), implying that each weather system type systematically produces too much precipitation and suggesting that the model is suitable for our purposes of exploring relative changes in precipitation. The sensitivity of the bias to resolution has been shown to be negligible, further pointing to an issue that is isolated to the microphysics or convective parameterizations rather than the dynamical core (Hack et al. 2006).
The PI simulation produces magnitudes and patterns of accumulation across Greenland that are similar to observations. It should be noted, however, that comparisons between the PI simulation and modern observations are only approximate on account of recent changes in climate that have generally led to increased accumulation in recent decades (Burgess et al. 2010), the bias toward higher accumulation in CCSM3, and the need to compare gridbox values against observations that are essentially at a point. The grid cell that we focus on (centered at 73.9°N, 39.4°W) contains the location of Summit Camp (72.6°N, 38.5°W; elevation 3216 m) and has a simulated annual mean accumulation of 27.9 cm yr−1 (Figs. 1a,b), whereas observations of accumulation at Summit Camp between 1963 and 1987 yield annual averages of approximately 26 cm yr−1 (van der Veen and Bolzan 1999) (all accumulation amounts are reported in units of liquid water equivalent). The simulated seasonal cycle of accumulation has an amplitude of 6 cm yr−1 (Fig. 1), giving a seasonally averaged summer-to-winter ratio of accumulation of 1.3 (Fig. 2b). Instrumental data able to resolve the seasonal cycle of accumulation are not as widely available as for temperature because of difficulties in obtaining continuous accumulation measurements at automated weather stations, but monthly samples from snow pits taken between 1986 and 1991 (Shuman et al. 1995) indicate an amplitude of the annual cycle of 3.0 cm yr−1, giving a summer-to-winter accumulation ratio of approximately 1.2.
Monthly variability in the accumulation for the PI case is slightly lower than in observations, with modeled standard deviations of 0.73 cm month−1 and observed values of 0.83 cm month−1 (Shuman et al. 1995), although this is expected because the model is specified to have climatological SSTs and sea ice that do not have interannual variability. Though the annual cycle dominates the power spectrum of most dynamical variables, interannual variability associated with the North Atlantic Oscillation (NAO) also has a discernible signature in the modern accumulation and temperature signal in central Greenland (Sodemann et al. 2008b). A similar NAO-like mode is also present in the PMIP2 LGM simulations, although it generally explains less sea level pressure variance than the modern NAO (Pausata et al. 2009). Previous studies also suggest that these modes of variability may influence the isotopic composition of accumulation through precipitation intermittency (Casado et al. 2013) and weather regime variability (Ortega et al. 2014).
It is also possible to make a spatial comparison for which we rely upon the estimates for sites given in Steen-Larsen et al. (2011). Similar results would be obtained using the mapped observational estimates of Burgess et al. (2010). The model grid box containing the GRIP site (72.6°N, 37.6°W; elevation 3230 m) neighbors that containing Summit Camp and has an accumulation of 26.6 cm yr−1, compared with observations of 23 cm yr−1 (Fig. 1). The more distant and interior NGRIP site (75.1°N, 42.3°W; elevation 2919 m) has a simulated accumulation of 24.5 cm yr−1 and observed values of 19 cm yr−1. The North Greenland Eemian (NEEM) camp (77.5°N, 51.1°W; elevation 2484 m) and Camp Century (77.2°N, 61.2°W; elevation 1890 m) in the northwest have 34.6 and 40.3 cm yr−1 simulated accumulations versus observed values of 20 and 35 cm yr−1, respectively. Finally, the Dye-3 core site near the southern coast (65.2°N, 43.8°W; elevation 2490 m) has an observed accumulation of 50 cm yr−1 that is incommensurate with the simulated value of 87.7 cm yr−1, perhaps because the steep topography near that grid cell is poorly resolved in the model, which instead represents precipitation nearer to the coast. The effects of sublimation have not been incorporated into the above comparisons, but Box and Steffen (2001) estimated these to be approximately 5 cm yr−1 for NEEM camp, Camp Century, and the Dye-3 core site and approximately 2 cm yr−1 for Summit Camp, GRIP, and NGRIP, each of which would tend to reduce the model–data discrepancy. Loss from surface melt and runoff would also decrease model–data discrepancies.
As expected, the PI and LGM simulations feature dramatically different patterns of accumulation, with glacial boundary conditions inducing a change in the annual mean as well as the seasonal cycle (Figs. 1, 2). The grid cell containing Summit Camp has an annual mean accumulation that decreases from 27.9 cm in the PI simulation to 7.7 cm in the LGM simulation, generally consistent with the observed reduction from 24 cm yr−1 (modern) to 5.5–7 cm yr−1 (LGM) estimated by Cuffey and Clow (1997) for a range of margin retreat scenarios, after accounting for the wet bias in the model. Andersen et al. (2006) used annual layer thicknesses in the NGRIP ice core to estimate a slightly smaller decrease in accumulation from 19 cm yr−1 in the modern era to 6.2 cm yr−1 during the LGM. The simulated summer-to-winter ratio of accumulation at Summit Camp increases from 1.3 in the PI simulation to 4.0 in the LGM simulation (Fig. 1), an increase that is generally consistent with other model simulations (Braconnot et al. 2007a), though those models run using the more extreme changes in sea ice indicated by the CLIMAP (Climate: Long-Range Investigation, Mapping, and Prediction) boundary conditions show even larger ratios in the range of 11.9–18.3 (Werner et al. 2001; Li et al. 2005; Langen and Vinther 2009). The correspondence of our simulation with observations and other simulations, along with discrepancies that are of an anticipated sign, encourages our exploring accumulation source properties in more detail.
3. Source estimates
Explicit modeling of sources is computationally expensive, typically involving either adjoint methods (Rao 2007) or a set of numerical Eulerian tracers to estimate the footprint function S(x, y, t). An approximate estimate of moisture source footprints can be obtained using Lagrangian trajectories, which have been used extensively for air quality studies and in atmospheric chemistry applications (Arya 1998). A much larger time step can be used for propagating Lagrangian particles because they do not feed back on the dynamics (Draxler and Hess 1997), and we use 3-h time steps. Additionally, particle motions can be integrated in particular spatial regions or for temporal subsets of interest without reintegrating the dynamics.
Lagrangian trajectories have previously been used for a range of purposes, including the identification of dynamical regimes (e.g., Kahl et al. 1997), studying isotopic distillation using Rayleigh models (Ersek et al. 2010), examining the isotopic signatures of different air masses and their contributions to the composition of rainfall, accumulated snow, and near-surface water vapor (Bonne et al. 2013; Steen-Larsen et al. 2014), and estimating smooth moisture source footprints and their variability using large ensembles of trajectories (Sodemann et al. 2008b; Bonne et al. 2013). Here we use the Lagrangian technique to evaluate the relationship between moisture sources and the seasonal cycle. This approach is distinct from previous studies in that it uses large ensembles to estimate source footprints in paleoclimate simulations, as opposed to meteorological reanalyses.
a. Lagrangian footprints
A modified version of the National Oceanic and Atmospheric Administration (NOAA) Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model is used to estimate the spatiotemporal distribution of water vapor sources contributing to precipitation in Greenland (Draxler and Hess 1997). At each time step, 10 000 particles are released in the atmospheric column above the Summit Camp grid cell, distributed vertically in proportion to the moisture tendency due to cloud sedimentation so that the resulting footprint captures the source of accumulated snow. Full details of the modifications are provided in the appendix.
The Lagrangian source footprint estimates for both simulations have maxima centered on the Irminger Sea (the basin southwest of Iceland and east of the southern tip of Greenland) for all seasons (Fig. 3). In the PI simulation, the footprint covers broad regions of the Greenland and Nordic seas, North Atlantic, and parts of the northeastern Pacific. In contrast, the LGM simulation footprint is spatially limited by sea ice cover in the Atlantic sector and by the upstream blocking of flow over the Laurentide in the Pacific sector. The seasonality of the LGM footprint is greater than that of the PI simulation resulting from greater seasonal changes in the extent of LGM sea ice cover and a consequently greater modulation of surface fluxes in the North Atlantic and Nordic seas. Although the simulations lack realistic interannual variability as a result of the use of climatological SSTs and sea ice cover, the PI source footprint is broadly consistent with that derived by Sodemann et al. (2008a) for central Greenland under modern conditions using meteorological reanalysis.
b. Consistency with Eulerian tracers
As a consistency check, water vapor sources are explicitly represented at low resolution using a so-called “painted water” Eulerian tracer scheme (Werner et al. 2001) that was previously adapted for CAM3 (Noone and Sturm 2010). Eulerian tracers do not influence the dynamics of the simulation, but each must be advected and their tendencies computed individually by the moist physics parameterizations at every three-dimensional grid point for each time step. Although computationally intensive, Eulerian tracer schemes have the advantage of being approximately numerically conserved so that a closed moisture source budget can be constructed.
We defined 11 tracer source regions: 8 zones in the North Atlantic, Greenland, the Arctic Ocean, and the Nordic seas; 1 for the remaining Northern Hemisphere extratropical ocean surface (primarily the Pacific); 1 for the Mediterranean, Black, and Caspian Seas; and 1 for the global tropical oceans (15°S–15°N). The regions are automatically adjusted between the PI and LGM simulations to match the changing land–ocean mask. A threshold is applied such that only grid boxes with 75% ocean fraction or greater are included. No such distinction is made for sea ice cover, so the Eulerian surface fluxes capture zones of partial ice cover.
We define precipitation sources as follows. The spatial footprint of surface sources of precipitation occurring at time t in central Greenland, S(x, y, t), is given in units of m−2 such that the integral of S over the entire Earth’s surface equals one. In principle S(x, y, t) contains all source information, but it is typically unknown and difficult to estimate. The fraction, rj(t), of precipitation sourced from region j of area Aj is
where Pj(t) and P(t) are the source-region-specific precipitation and total precipitation in central Greenland, respectively. For the Eulerian estimates S(x, y, t) is estimated for each tracer region by simply dividing rj(t)—directly output from the model—by Aj.
The Eulerian tracers, which conserve the moisture budget partitioned between the various surface source regions, show that only 3.8% of central Greenland’s accumulation is sourced from the tropics in the PI case, compared with 3.3% in the LGM case. This is consistent with recent studies (Sodemann et al. 2008b; Langen and Vinther 2009), but in contrast to previous lower-resolution simulations that resolved the effects of topography and mixing less completely and predicted a 2–3 times greater fraction of tropical sources (Charles et al. 1994; Werner et al. 2001).
The Eulerian and Lagrangian footprint estimates show broad agreement in the annual mean for both the PI (Fig. 4) and LGM (Fig. 5) simulations. Both methods indicate that the primary source for annual mean accumulation for central Greenland is within the North Atlantic sector, with the highest concentrations occurring southeast of Greenland in the Irminger Sea. In the PI simulation, both methods capture a footprint that is less localized to the Greenland and Nordic seas (Fig. 4) than in the LGM simulation (Fig. 5). There are also some systematic differences. In both simulations, the Eulerian footprint estimates are more diffuse than the Lagrangian ones, owing to the presence of numerical diffusion in the model, which is used to represent subgrid-scale turbulent mixing from on the order of 200 km to the inertial subrange. This is expected, since the Lagrangian method does not attempt to resolve mixing that is approximated by numerical diffusion in the Eulerian model. The discrepancy is more pronounced in the PI simulation and is especially evident during the summer months (Figs. 4a,d and 5a,d). The presence of the intersimulation and interseasonal differences in footprint dispersion may both be related to the ratio of transient to stationary eddy fluxes. This ratio is greater during summer than winter in both simulations, and is greater for the PI simulation than the LGM simulation in all seasons. The Atlantic storm track is largely locked to a topographically generated stationary wave during modern winter and year-round in the LGM, presenting a barrier to cross-jet mixing. Under these conditions a greater fraction of the meridional fluxes are due to stationary waves, with the zonal mean equaling 2.3 PW at 53°N in the LGM simulation versus 1.1 PW in the PI simulation. Because baroclinic activity is suppressed by reduced seeding (Donohoe and Battisti 2009), meridional fluxes are dominated by the stationary wave. In contrast, the PI summer has meridional fluxes dominated by transient eddies that achieve proportionally more transport and mixing of tracers through numerical diffusion.
The inference that PI and summer discrepancies result from a greater fraction of the kinetic energy spectrum not being accounted for by the Lagrangian model is also consistent with previous Eulerian studies that used lower atmospheric resolution and, therefore, generally entailed more diffuse transport and obtained source footprints more heavily weighted at greater distances (Charles et al. 1994; Werner et al. 2001). Neither model is strictly correct: the Eulerian model uses a diffusive approximation in the vertical that depends on the vertical shear and the Richardson number and is influenced by the horizontal diffusion operator applied to enforce stability of the momentum equations, producing overly homogenized mixing, while the Lagrangian model ignores unresolved mixing outside of the boundary layer. Convergence of the Lagrangian and Eulerian footprint estimates would require accounting for all mixing processes.
4. Dynamics and seasonality of accumulation events
Because of stationary wave forcing by thermal and orographic variations in the Northern Hemisphere, the basic state of the LGM includes anomalously cold and dry conditions in the Atlantic sector, relative to other locations at the same latitude. The high topography of the Laurentide Ice Sheet induces a low-level anticyclone that increases poleward moisture and heat fluxes in the Pacific sector (Cook and Held 1988). In the lee of the Laurentide, the thermodynamic balance of the anticyclone is maintained primarily by equatorward advection of cold, dry air, effectively cutting off Greenland from long-range westerly or southerly sources of moisture under glacial conditions. The anticyclone is strongest in winter and maintains a similar low-level flow—albeit with reduced wind speeds—during the summer.
Thermal and mechanical forcing tends to lock the location of the jet stream to the stationary wave, displacing it to the south, and making it less variable on seasonal and synoptic time scales. Its southward displacement also means that the North Atlantic can maintain more significant sea ice cover. The vast majority of the poleward moisture flux in the Pacific is lost to precipitation prior to the flow entering the Arctic basin (Fig. 3). Though there is a source contribution of moisture from the Pacific (a total of 14% contribution over the whole basin in the PI case), it is even weaker (12% total contribution) in the LGM case, despite the relatively greater poleward moisture transport west of the Laurentide. Much of the moisture is lost orographically in the Cordilleran portion of the Laurentide in western North America, where it maintains the ice sheet’s surface mass balance against higher-than-zonal average temperatures (Roe and Lindzen 2001), and the remainder is further reduced in balance with radiative cooling over the Arctic. It has been proposed that some properties of the Greenland ice core records may be explained by the presence of a “split jet” that switches Greenland’s moisture source between the Pacific and Atlantic (Wunsch 2006; Langen and Vinther 2009). The lack of a significant moisture flux through the Arctic under glacial conditions suggests that switches to a Pacific source are unlikely, at least in this model, even if dry air flows northward around the Laurentide. This result differs from one of the simulations of Langen and Vinther (2009), which they note likely stems from that simulation’s use of CLIMAP surface boundary conditions featuring much more extensive perennial ice cover in the North Atlantic. The predominance of northerly cold air advection in the lee of the Laurentide raises the question of how substantial accumulation can occur in Greenland under glacial conditions at any time of the year. Is the strong seasonal cycle of accumulation that allows most of the accumulation there to occur only during summer a consequence of a weaker summer stationary wave, or are other mechanisms at work?
An examination of the synoptic conditions associated with accumulation events clarifies the relative importance of topography, sea ice cover, and seasonal dynamics of the jet stream in setting the structure of the source footprints. While eddy activity is large throughout the Atlantic storm track in the preindustrial winter, it is limited, except in the diffluent jet exit in the east Atlantic during the glacial winter (Fig. 6). Thus, despite the strong baroclinicity and meridional temperature gradient, the meridional moisture flux is almost completely eliminated over the ice-covered region of the northwestern Atlantic in the LGM simulation.
Lagged one-point regression maps of the 500-hPa geopotential height associated with Greenland accumulation events (Fig. 7) reveal the differing roles of the PI and LGM storm tracks in precipitation falling over central Greenland. Upper-level geopotential responses are not necessarily directly related to local surface conditions, reflecting instead a complex set of local and nonlocal thermodynamic and dynamic feedbacks (Cayan 1992; Kushnir et al. 2002). Nevertheless, the one-point correlations show that typical conditions leading to accumulation are not the same in the two simulations. In the PI simulation, central Greenland can rely upon regular synoptic events that draw moisture from throughout the region during the entire year. One-point maps with a broad range of lags support this interpretation, showing migrating perturbations with significant correlation in the 1–5-day range for the PI simulation (Figs. 7a,b for 1-day lag; longer lags not shown). The presence of localized sources of surface latent heat flux both sustains the synoptic events as they migrate poleward and provides them with additional moisture. This is consistent with observational studies of modern accumulation in Greenland (Sodemann et al. 2008b) that link the large interannual variability there to the North Atlantic Oscillation (i.e., displacement and modulation of the storm tracks).
One-point maps from the LGM simulation reveal a dependence on unusual southeasterly fluxes over the Irminger Sea (Figs. 7c,d), and one-point maps at longer lead times show rapid growth and decay of this structure. The short decorrelation time during summer—roughly 2 days before large-scale coherent patterns break down for the LGM simulation versus 5 days for the PI simulation—shows that the development does not follow the synoptic-scale norm of baroclinic life cycles. Again, this is consistent with previous studies that show that synoptic-scale transient eddy fluxes are greatly reduced in the Atlantic sector under glacial conditions (Donohoe and Battisti 2009). The peculiar structure of these events, coupled with the low accumulation rate in the LGM winter (7.7 versus 27.9 cm yr−1 in the PI), suggests that they are rare occurrences. Indeed, the distribution of daily accumulation rate shifts from a skewed distribution peaked at 0.41 mm day−1 (PI) to one peaked at 0.0 mm day−1 (LGM). The median preindustrial accumulation rate (0.56 mm day−1) is equal to that of the 90th percentile for the LGM.
5. Some implications for the interpretation of source properties
Changes in the dynamics of accumulation events between the PI and LGM simulations are indicative of a diminished role of the storm track under glacial conditions. Additionally, footprints suggest that source regions vary on interseasonal and glacial–interglacial time scales. Source properties are important for determining the isotopic composition of subsequent precipitation, and hence for the interpretation of stable isotope proxies of past temperature derived from the Greenland ice cores (Masson-Delmotte et al. 2005). Measurements of δ18O and δD in ice are used to infer past temperature changes, primarily recording the temperature difference between moisture source and sink (Jouzel et al. 2005). Consequently, the balance between moisture sourced regionally versus moisture sourced from greater distances by meridional fluxes in the storm track may influence the isotopic signal, warranting a more detailed examination of source properties.
The distribution of the distance between moisture sources and the sea ice edge can be thought of as a measure of the sensitivity of the source region geometry to sea ice cover. For each source point in the monthly average Lagrangian footprints, we find the distance to the nearest sea ice edge (defined by 70% area fraction) and then compute a weighted histogram of these distances, where the weights are equal to the product of the gridbox size and the source footprint at each location. Histograms are computed for monthly climatologies in order to illustrate the seasonal cycle (Fig. 8).
The difference between the two simulations is clear: the LGM distribution peaks near the sea ice edge, while the PI distribution shows that moisture is drawn from a broader region that is largely unconstrained by sea ice, especially during summer. PI conditions leave open many pathways for moisture to reach Greenland from the North Atlantic and Nordic seas via both storm track variability and localized fluxes; the different basic state during the LGM leads to Greenland being cut off from distant sources of moisture and more dependent on the seasonal retreat of sea ice in the Greenland and Nordic seas. This explanation is also consistent with both the spatial pattern of the summer-to-winter accumulation ratio (Fig. 2) and the one-point maps (Fig. 7), which together show that the seasonal growth and retreat of sea ice cover proximate to Greenland induces large accumulation seasonality in the LGM (Fig. 2d), but also for northwestern Greenland in the PI (Fig. 2b). This gives some indication that the spatial gradient of the seasonality of accumulation is controlled by the difference between summer and winter sea ice cover over the nearby ocean surface. The relatively sharp spatial dependence of this signal supports the notion that sources are localized and may help to explain proxy evidence of large differences in accumulation between different Greenland ice core sites during Dansgaard–Oeschger events (Guillevic et al. 2013).
The distance of the source region from central Greenland may also be important, in that it implies varying pathlengths under different climate conditions. In an ideal Rayleigh model of isotopic distillation, the distance traveled does not influence the process, as it is parameterized as a function of the remaining water vapor that is, in turn, driven by changes in pressure, temperature, and saturation specific humidity. However, in the more realistic context of an air parcel mixing turbulently and via diffusion with environmental air during advection, the length of the parcel’s trajectory becomes more important for setting the isotopic composition of its water vapor.
The migrating spatial footprint also has implications for the variability of source temperature. The PI distribution of source temperatures is in agreement with Sodemann et al. (2008a), who used high-resolution reanalysis data with a similar Lagrangian model. But more so than for the PI conditions, the simulated LGM conditions produce accumulation whose source is largely localized to the sea ice edge in space and near freezing in temperature (Fig. 8). This feature of LGM accumulation suggests that an approximation of a constant source temperature would be more appropriate under LGM conditions. Source relative humidity influences second-order isotopic parameters, such as deuterium excess and 17O excess (Jouzel et al. 2013; Pfahl and Sodemann 2014), and it would be useful in future work to compare simulated source relative humidity against these measures of excess isotopic variation. The lack of diffusivity in the Lagrangian approach may mean that the limited spatial extent of sources is exaggerated, though it is also the case that the ice edge rectifies surface fluxes so that the diffusivity effect may be less important than over open ocean.
Comparison of model results against deuterium excess and 17O excess (e.g., Uemura et al. 2008; Landais et al. 2012; Steen-Larsen et al. 2014) would provide further insights into the relationship between the localization of moisture source, close to the sea ice edge during the LGM, and isotopic fingerprints (Kurita 2011). Although beyond the scope of this study, it would be relevant to employ a model that explicitly simulates fractionation processes to better distinguish between the influences of source temperature, specific humidity, and wind speed on source properties, each of which is likely to be influenced by the proposed proximity of the sea ice edge (Lewis et al. 2013; Benetti et al. 2014). Also potentially important would be to account for in situ modification of snow and firn isotopic composition by equilibration with atmospheric water vapor (Steen-Larsen et al. 2014).
We apply a Lagrangian analysis for the source regions associated with accumulation in central Greenland to simulations of preindustrial (PI) and Last Glacial Maximum (LGM) climate. The Lagrangian estimate is consistent with a more standard Eulerian estimate up to differences attributable to a lower effective diffusivity in the Lagrangian approach. The Lagrangian approach is applied with a sufficiently large particle ensemble to provide for high-resolution estimates of source contributions. LGM accumulation in central Greenland is found to strongly depend on regional surface fluxes modulated by sea ice cover (Li et al. 2010) and low transient eddy activity in the North Atlantic (Donohoe and Battisti 2009). In contrast, general synoptic-scale variability is found to play a dominant role in providing accumulation across all seasons during the PI climate, with seasonal sea ice cover being the dominant factor only in northwestern Greenland.
A common assumption in ice core paleothermometry has been that accumulation is largely sourced from tropical or subtropical oceans (Dansgaard 1964; Johnsen et al. 1989), a conclusion not supported by the source footprint estimates presented here, nor by several previous studies (Werner et al. 2001; Sodemann et al. 2008a,b; Langen and Vinther 2009). Nor is a switch to a significant Pacific source during the LGM supported by this analysis. Instead, we find that moisture sources are primarily from the open North Atlantic during the PI climate and that they are confined to along the sea ice edge in the North Atlantic during the LGM.
This tracking of the source region along the sea ice edge leads to source temperatures being consistently near freezing. Such control by the sea ice edge suggests the possibility of a simplification whereby variations in source temperature may be assumed approximately constant when interpreting LGM isotopic composition in central Greenland ice cores, though variations in wind speed, water vapor saturation, and other environmental parameters could remain important variables. Further analysis comparing simulations of isotopic composition (e.g., Lewis et al. 2013) against observations, including triple oxygen isotope analysis (e.g., Uemura et al. 2008; Steen-Larsen et al. 2014), would be useful for better ascertaining the role of the sea ice edge in controlling the composition of Greenland’s accumulation.
The authors would like to acknowledge Camille Li, Ethan Butler, Valérie Masson-Delmotte, and two anonymous reviewers for helpful comments and discussions. This study was supported by NSF Grant 1304309.
Lagrangian Source Model
The modifications to the HYSPLIT model consist of a module to read netCDF-formatted meteorological output from the climate model and a back-deposition model that is designed to estimate moisture fluxes between the ocean and atmosphere. In the back-deposition model, virtual particles are released from all levels of the atmospheric column above the site of interest and are advected backward in time. Back trajectories are used partly for computational reasons, as the backward problem is more efficient than an explicit source–receptor forward integration in the Lagrangian framework. Each particle is assigned a dimensionless tracer quantity representing water vapor content with an initial value (i.e., at the time of precipitation) proportional to the instantaneous water vapor tendency due to cloud sedimentation at that point. The intrinsic tracer is then allowed to deposit over the surface of the ocean, with the rate of deposition controlled by relaxation in the boundary layer:
where m is the residual particle mass, hsl is the depth of the layer occupied by the particle (during deposition, the depth of the surface layer), and Δm is the loss of particle mass over a time step of length Δt. The so-called deposition velocity υd is parameterized using the resistance method:
where the Prandtl number Pr is 0.923, the von Kármán constant κ is 0.4, U* is the friction velocity, Z0 is the roughness length, ψh is the integrated stability function for heat and moisture, and Ts is the surface temperature (Draxler and Hess 1997). The effective Henry’s constant for water vapor Heff is 55/es, where es is the saturation vapor pressure. The dynamic boundary layer variables are estimated using surface momentum flux, sensible heat flux, surface roughness, and thermodynamic state variables that are output from CAM to compute profiles of vertical fluxes in the boundary layer using Monin–Obukhov similarity theory. In this way, vertical fluxes are tightly coupled to the CAM meteorology, despite the evaporation being estimated by a second layer of parameterization.
The resistance model does not explicitly account for moist convective processes, though the particle advection does include transport by bulk turbulent fluxes of constituent gases that result from the convective motions. The particles are advected in this way for 10 days so that virtually all of the particles cycle through to the surface (the residence time of atmospheric water vapor is approximately one week; Trenberth 1998), and we find consistent result for ranges of this period between 5 and 15 days. A “wall” boundary condition is imposed at 15°N to track particles that enter the tropics where vertical transport of moisture by deep convection cannot be adequately traced using large-scale velocity fields from CAM.