Methanesulfonic acid (MSA) migration in polar ice: data synthesis and theory

Abstract

Methanesulfonic acid (MSA; CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$H) in polar ice is a unique proxy of marine primary productivity, synoptic atmospheric transport, and regional sea-ice behavior. However, MSA can be mobile within the firn and ice matrix, a post-depositional process that is well known but poorly understood and documented, leading to uncertainties in the integrity of the MSA paleoclimatic signal. Here, we use a compilation of 22 ice core MSA records from Greenland and Antarctica and a model of soluble impurity transport in order to comprehensively investigate the vertical migration of MSA from summer layers, where MSA is originally deposited, to adjacent winter layers in polar ice.

We find that the shallowest depth of MSA migration in our compilation varies over a wide range ($\sim$ 2 to 400 m) and is positively correlated with snow accumulation rate and negatively correlated with ice concentration of Na${}^{+}$ (typically the most abundant marine cation). Although the considered soluble impurity transport model provides a useful mechanistic framework for studying MSA migration, it remains limited by inadequate constraints on key physico-chemical parameters – most notably, the diffusion coefficient of MSA in cold ice ($D_{\mathrm{MS}}$). We derive a simplified version of the model, which includes $D_{\mathrm{MS}}$ as the sole parameter, in order to illuminate aspects of the migration process. Using this model, we show that the progressive phase alignment of MSA and Na${}^{+}$ concentration peaks observed along a high-resolution West Antarctic core is most consistent with 10${}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ < $D_{\mathrm{MS}}$ < 10${}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$, which is 1 order of magnitude greater than the $D_{\mathrm{MS}}$ values previously estimated from laboratory studies. More generally, our data synthesis and model results suggest that (i) MSA migration may be fairly ubiquitous, particularly at coastal and (or) high-accumulation regions across Greenland and Antarctica; and (ii) can significantly change annual and multiyear MSA concentration averages. Thus, in most cases, caution should be exercised when interpreting polar ice core MSA records, although records that have undergone severe migration could still be useful for inferring decadal and lower-frequency climate variability.

Introduction

Measurements of soluble impurity species in polar ice cores provide important high-resolution proxies of past climatic phenomena, including past changes in sea-ice extent, marine and terrestrial productivity, volcanism, biomass burning, atmospheric cycling, and anthropogenic pollution (e.g., Legrand and Mayewski, 1997). A foundational premise, however, is that these species undergo negligible post-depositional redistribution in the ice column, an assumption unsupported by numerous ice core records from Greenland and Antarctica. Processes acting within the upper firn layer – including wind pumping, diffusion, photolysis, volatility, sublimation, and melt (Wolff, 1996) – can affect the stability of chemical species soon after deposition and undermine the climatic interpretation of these records down-core (Wagnon et al., 1999; Weller et al., 2004). Deeper in the ice column, both observation (e.g., Barnes et al., 2003b) and theory (Nye, 1991; Rempel et al., 2001, 2002) indicate the potential for solid and liquid-state chemical migration, impacting the stability of chemical species much later after deposition. In this study, we focus on one species particularly susceptible to vertical migration in polar ice, methanesulfonic acid, or MSA (CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$H; Pasteur and Mulvaney, 2000).

The processes leading to the production, transport, and deposition of MSA onto an ice sheet are complex (e.g., Abram et al., 2013). The progenitor compound of MSA, dimethylsulfoniopropionate, or DMSP ((CH${}_{\mathrm{3}}{)}_{\mathrm{2}}$S${}^{+}$CH${}_{\mathrm{2}}$CH${}_{\mathrm{2}}$COO${}^{-})$, is produced by certain marine algae as an osmotic regulator (Dickson and Kirst, 1986). Planktonic life cycle processes ultimately release DMSP to the water column, whereupon the ensuing bacterial-mediated cleavage of the compound promotes the formation of dimethylsulfide, or DMS ((CH${}_{\mathrm{3}}{)}_{\mathrm{2}}$S), a highly insoluble gaseous compound (Yoch, 2002). Once freed to the atmosphere, DMS is rapidly photo-oxidized (Saltzman et al., 1983), branching to form either non-sea-salt sulfate (nssSO${}_{\mathrm{4}}^{-\mathrm{2}})$ or, to a lesser extent, MSA. Unlike nssSO${}_{\mathrm{4}}^{-\mathrm{2}}$, however, DMS production appears to be the exclusive source of MSA (Abram et al., 2013).

Early studies on DMS and DMSP production in the Southern Ocean (Curran and Jones, 2000) and in Arctic waters (Leck and Persson, 1996) reported the greatest fluxes of DMS near the marginal sea-ice zone at the onset of spring–summer decay (Turner et al., 1995; Curran et al., 1998). Since concentrations of atmospheric MSA rapidly decrease with increasing altitude and distance from their marine source (the mean atmospheric lifetime of MSA is estimated to about 7 days; Hezel et al., 2011), deposition of MSA at coastal ice sheet localities near-ubiquitously exhibits well-defined annual peaks during the late-spring to summer months. Conversely, the relatively low wintertime MSA deposition may be jointly attributed to (1) limited marine productivity during polar darkness; (2) increased wintertime sea-ice extent and, accordingly, atmospheric transport distances; and (or) (3) diminished atmospheric OH radical concentrations, the primary oxidant of airborne DMS (Jourdain and Legrand, 2001). The strong seasonality and unique marine source of MSA in ice cores have led to its predominant use as a high-resolution proxy for past sea-ice cover (e.g., Welch et al., 1993; Abram et al., 2013).

Mulvaney et al. (1992) appear to be the first authors to report migration of MSA in polar ice, using data from Dolleman Island, Antarctic Peninsula. In the shallow portions of the Dolleman Island core, concentrations of MSA exhibited well-defined summer maxima, as expected. With increasing depth in the ice column, however, a distinctive shift to predominantly winter [MS${}^{-}$] maxima was found (hereafter, MSA is denoted as MS${}^{-}$ when referring to its anionic form, CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}^{-}$, as measured by ion chromatography). Concluding that a change in the seasonality of peak MSA production and (or) deposition was unlikely, the authors postulated this shift to result from post-depositional vertical migration. Since then, numerous ice core studies from both Antarctica and Greenland have reported the migration phenomenon over a wide range of depths, temperatures, and ionic concentrations in the ice column (Table 1 and references therein).

Ice core sites where MSA migration (or lack thereof) has been reported (“$\sim$” denotes unavailable data and “n/a” denotes a non-applicable field).

Ice core name	Lat/long	Year	Depth reached/	Elevation,	Distance to	SAT,	$\dot{b},$	MSA	$z_{\mathrm{fo}}$, m	${\mathit{\rho }}_{\mathrm{fo}}$,	[$\stackrel{\mathrm{‾}}{{\mathrm{MS}}^{-}}$],	[$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$],	Reference
		collected	reported,	m	coast, km	${}^{\circ }$C	m w.e. yr${}^{-\mathrm{1}}$	migration	(year)${}^{\mathrm{c}}$	kg m${}^{-\mathrm{3}}$	$\mathrm{\mu }$g L${}^{-\mathrm{1}}$	$\mathrm{\mu }$g L${}^{-\mathrm{1}}$
			m (years)					reported?${}^{\mathrm{a}}$
West Antarctic MSA records
Dyer Plateau	70.65${}^{\circ }$ S, 65.02${}^{\circ }$ W	1988–1989	56 m w.e. (103)	1943	190	$-$21.7	0.48	U	n/a	$\sim$	$\sim$	$\sim$	Mulvaney et al. (1992),
													Pasteur and Mulvaney (2000)
Dolleman Island (1)	70.58${}^{\circ }$ S, 60.93${}^{\circ }$ W	1985–1986	96 m w.e. (297)	398	20	$-$16.8	0.34	Y	10 (10)	$\sim$	16.6	382	Mulvaney et al. (1992)
Dolleman Island (2)	70.58${}^{\circ }$ S, 60.93${}^{\circ }$ W	1992–1993	18.5 m w.e. (45)	398	20	$-$16.8	0.34	Y	8.5 m w.e.	$\sim$	$\sim$	$\sim$	Pasteur and Mulvaney (2000)
									(19)
Gomez Nunatak	74.02${}^{\circ }$ S, 70.63${}^{\circ }$ W	1980–1981	37 m w.e. (42)	1130	135	$-$17	0.88	N	n/a	n/a	$\sim$	$\sim$	Mulvaney and Peel (1988),
													Pasteur and Mulvaney (2000)
Beethoven Peninsula	71.88${}^{\circ }$ S, 74.57${}^{\circ }$ W	1992–1993	30 m w.e. (28)	580	16	$-$12.5	1.20	N	n/a	n/a	$\sim$	$\sim$	Pasteur and Mulvaney (2000)
Berkner Island north (1)	78.30${}^{\circ }$ S, 46.28${}^{\circ }$ W	1989–1990	11 (21)	730	50	$-$22.5	0.22	Y	9 (16)	560	19.0	583${}^{\mathrm{d}}$	Wagenbach et al. (1994)
Berkner Island north (2)	78.30${}^{\circ }$ S, 46.28${}^{\circ }$ W	1994–1995	39 m w.e. (174)	730	50	$-$22.5	0.20	Y	$\sim$	$\sim$	$\sim$	$\sim$	Pasteur and Mulvaney (2000)
Berkner Island south	79.60${}^{\circ }$ S, 45.62${}^{\circ }$ W	1989–1990	11 (28)	940	150	$-$24.5	0.17	Y	6 (16)	520	16.0	317${}^{\mathrm{d}}$	Wagenbach et al. (1994)
Siple 94-1	81.65${}^{\circ }$ S, 148.8${}^{\circ }$ W	1994–1995	148 (1150)	621	400	$-$25	0.13	Y	2 (6)	$\sim$	21.7	147	Kreutz et al. (1998)
DIV2010	76.80${}^{\circ }$ S, 101.7${}^{\circ }$ W	2010–2011	112 (216)	1329	180	$-$24	0.41	Y	26 (38)	640	7.7	37.6	Criscitiello et al. (2014),
													Criscitiello (2014)
THW2010	77.0${}^{\circ }$ S, 121.2${}^{\circ }$ W	2010–2011	62 (145)	2020	340	$-$28	0.28	Y	17 (32)	610	9.0	27.5	Criscitiello et al. (2014),
													Criscitiello (2014)
Bruce Plateau	66.03${}^{\circ }$ S, 64.07${}^{\circ }$ W	2009–2010	448 ($\sim$)	1976	30	$-$14.8	1.98${}^{\mathrm{b}}$	Y	395 (560)	$\sim$	$\sim$	$\sim$	Goodwin (2013),
													Porter et al. (2016)
Byrd NBY-2	80.02${}^{\circ }$ S, 119.5${}^{\circ }$ W	1989–1990	164 (629)	1530	650	$-$28	0.11	Y	2.6 (13)	$\sim$	6.7	37	Langway et al. (1994)
Filchner-Ronne D235	77${}^{\circ }$ S, 64${}^{\circ }$ W	1987–1988	4.3 m w.e. (20)	$\sim$	125	$\sim$	0.18	Y	1.8 m w.e.	$\sim$	14.4	361${}^{\mathrm{c}}$	Minikin et al. (1994)
									(9)
Ferrigno	74.57${}^{\circ }$ S, 86.90${}^{\circ }$ W	2010–2011	136 (309)	1354	475	$-$24.7	0.35	Y	25 (30)	620	5.9	88	Thomas and Abram (2016)
East Antarctic MSA records
Law Dome: W20k	66.77${}^{\circ }$ S, 112.35${}^{\circ }$ E	1997	45 (191)	1370	110	$-$22	0.15	Y	$\sim$	$\sim$	$\sim$	$\sim$	Curran et al. (2002)
Law Dome: DSS	66.77${}^{\circ }$ S, 112.42${}^{\circ }$ E	1997; 2000	124 (156)	1370	120	$-$21.8	0.64	U	n/a	n/a	7.0	85.5	Curran et al. (2002, 2003)
Law Dome: DE08	66.72${}^{\circ }$ S, 113.18${}^{\circ }$ E	1986	196 (145)	$\sim$	100	$-$19	1.27	N	n/a	n/a	$\sim$	$\sim$	Curran et al. (2002)
WHG – Victoria Land	72.90${}^{\circ }$ S, 169.08${}^{\circ }$ E	2006	105 (130)	400	12	$-$15	0.61	N	n/a	n/a	22	1901	Sinclair et al. (2012, 2014)
Greenland MSA records
Summit2010	72.33${}^{\circ }$ N, 38.28${}^{\circ }$ W	2010	87 (268)	3213	360	$-$29.5	0.22	Y	43.7–46.9	710–720	3.4	5.1${}^{\mathrm{e}}$	Maselli et al. (2017)
									(120–130)
D4	71.4${}^{\circ }$ N, 44.0${}^{\circ }$ W	2004	145 (270)	2710	300	$\sim$	0.42	Y	54.3–59.4	760–770	2.5	5.3${}^{\mathrm{d}}$	MSA data unpublished
									(90–100)
2Barrell	76.94${}^{\circ }$ N, 63.15${}^{\circ }$ W	2011	21.3 (21)	1685	100	$\sim$	0.51	N	n/a	n/a	$\sim$	$\sim$	Osterberg et al. (2015)

${}^{\mathrm{a}}$ Y $=$ yes; N $=$ no; U $=$ unclear.${}^{\mathrm{b}}$ Bruce Plateau annual mean accumulation rate appears to be highly variable over the core-depth; listed value represents the AD 1750–2010 estimate (Goodwin, 2013). ${}^{\mathrm{c}}$ See Supplement Sect. S1 and Table S1 for more information on how $z_{\mathrm{fo}}$ is defined at each site. ${}^{\mathrm{d}}$ [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$] estimated using the reported site value for [$\stackrel{\mathrm{‾}}{{\mathrm{Cl}}^{-}}$] and converted assuming a site-ratio [$\stackrel{\mathrm{‾}}{{\mathrm{Cl}}^{-}}]$ $/$ $\left[\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}\right]=\mathrm{1.8}$, the mean sea-water ratio (e.g., Chesselet et. al, 1972).${}^{\mathrm{e}}$ [$\stackrel{\mathrm{‾}}{\mathrm{Na}}$] measurements made using continuous-flow inductively coupled plasma mass spectrometry (ICP-MS), as opposed to [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$] measurements made using ion chromatography (IC). As ICP-MS measures both the soluble andinsoluble mass content, [$\stackrel{\mathrm{‾}}{\mathrm{Na}}$] values are likely slightly higher than IC-based estimates of [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$].

Despite the importance of [MS${}^{-}$] records in paleoclimatic studies, there have been few systematic evaluations of the environmental and (or) chemical conditions promoting MSA migration from summer to winter layers in polar ice. Furthermore, aspects of this migration process exhibit behavior distinct from other post-depositional processes and are important to understand mechanistically. First, the [MS${}^{-}$] maxima that are formed in winter are converse to what would be expected from typical (Fickian) diffusion, which would instead weaken the amplitude of summertime [MS${}^{-}$] peaks. Second, MSA movement has been reported to occur in the up-core direction (Curran et al., 2002), ruling out gravitational forcing as the sole mechanism for migration. These observations, corroborated by evidence for highly concentrated regions of sulfuric acid at the interface of individual ice crystals (Mulvaney et al., 1988), led to speculation that liquid migration of soluble impurities could occur along the boundaries of individual ice crystals, likely along concentration gradients (Mulvaney et al., 1992). However, critical questions remain. For example, why should MSA in particular exhibit migration while associated soluble impurities and acids (e.g., nssSO${}_{\mathrm{4}}^{\mathrm{2}-})$ do not (Pasteur and Mulvaney, 1999)? How could diffusion result in clearly defined concentration maxima in winter layers? Could the “trapping” of migrating MSA in the adjacent winter layer stem any subsequent “spillover” of MSA into an adjacent annual layer, as assumed in prior studies (e.g., Kreutz et al., 1998; Pasteur and Mulvaney, 2000; Becagli et al., 2009; Thomas and Abram, 2016)? Clearly, investigating the mechanisms(s) responsible for MSA migration in polar ice is necessary if we are to answer these questions and build confidence in the use of [MS${}^{-}$] in polar ice as a paleoclimatic indicator.

The overarching goal of this study is to develop a better understanding of the environmental and physico-chemical processes that are conducive to MSA migration in cold, polar ice, in order to improve the interpretation of ice-core [MS${}^{-}$] records. This paper is organized as follows. In Sect. 2, we consider published [MS${}^{-}$] and ancillary measurements from a variety of Greenlandic and Antarctic ice cores in an effort to identify site-specific factors that influence MSA migration. In Sect. 3, we present observations from the high-resolution, precisely dated DIV2010 ice core (West Antarctica) as a case study for MSA migration. In Sect. 4, we summarize our current understanding of the physico-chemical processes leading to MSA migration, utilizing an existing model describing soluble impurity transport along ice grain boundaries (Rempel et al., 2002). We derive a simplified version of this model to illuminate these processes and to test the ability of different values of MSA diffusivity to reproduce the down-core change in phase relationship between the concentrations of MSA (deposited primarily in summer) and Na${}^{+}$ (deposited primarily in winter) observed along the DIV2010 core. In Sect. 5, we assess the integrity of the DIV2010 [MS${}^{-}$] record and discuss the broader implications of our results for the interpretation of [MS${}^{-}$] records across a wide range of polar conditions. Section 6 concludes with an overview of our results and suggestions for future research.

Compilation of data for Antarctic (circles) and Greenland (diamonds) ice cores for which [MS${}^{-}$] records meet the criteria of this study. (a) Annual mean accumulation rate versus annual mean surface air temperature. (b) Core-averaged [MS${}^{-}$] versus distance to the coast. (c) Core-averaged [Na${}^{+}$] versus distance to the coast. The colors indicate whether MSA migration has been reported, deemed as unclear, or not reported in the original publications. The numbers indicate ice core sites with $\dot{b}$ > 0.45 m w.e. yr${}^{-\mathrm{1}}$: 1. Bruce Plateau (length of record: 448 m), 2. DS08 – Law Dome (196 m, corresponding to a time span of 145 years), 3. Beethoven Peninsula (47 m, 28 years), 4. Gomez Nunatak (56 m, 42 years), 5. WHG – Victoria Land (106 m, 130 years), 6. DSS – Law Dome (124 m, 156 years), and 7. Dyer Plateau (80 m, 103 years).

[Figure omitted. See PDF]

MSA migration in ice cores

In this section, we use observations to evaluate the relative importance of site-specific variables for MSA migration in polar ice. We compiled 22 ice core [MS${}^{-}$] records originating from 20 sites in Greenland and Antarctica (see Supplement), taken from both the literature and unpublished datasets (Table 1). The following criteria are adopted for selection of the records: (1) high temporal resolution and dating accuracy (annually resolved and dating uncertainty < 3 years); (2) multi-decadal record length (20 years minimum); and (3) documented changes in seasonality in the [MS${}^{-}$] record and (or) an explicit mention of MSA migration from summer to winter layers. Post-depositional surficial losses of MSA may occur via gaseous diffusion in the top 1–2 m of the firn at low-accumulation sites (Wagnon et al., 1999; Delmas et al., 2003; Weller et al., 2004). As a result, we exclude records from sites where annual mean accumulation rate is less than 100 kg m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$ and assume that vertical redistribution of MSA via gas-phase diffusion (Kuo et al., 2011) is negligible at all considered sites and depths. Also excluded are records from sites subject to moderate to severe melt, which, through percolation, can also rapidly redistribute MS${}^{-}$ along the firn (Moore et al., 2005).

Of the 20 sites, 13 indicate MSA migration, 5 do not indicate migration, and at 2 sites migration is deemed unclear at the greatest depths sampled (Fig. 1). The sites and records cover a wide range of climatologic and glaciological conditions, as represented by annual mean surface air temperature (SAT), annual mean accumulation rate, distance inland, and impurity concentrations in the ice (Table 1 and Fig. 1). Figure 1a illustrates the range of annual mean accumulation rates ($\dot{b})$ and annual mean surface air temperatures of all sites. A strong, nonlinear relationship is apparent between both variables, as expected (Cuffey and Patterson, 2010). Similarly, for the 14 sites where glaciochemical data are available (Table 1), Fig. 1b and c show, respectively, the depth-averaged concentrations of Na${}^{+}$ ([$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$]; see Table 1 for proper symbolic representation of Na${}^{+}$ concentration), an ionic species whose relevance we expound upon in Sect. 2.2, and MS${}^{-}$ ([$\stackrel{\mathrm{‾}}{{\mathrm{MS}}^{-}}$]). Both [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$] and [$\stackrel{\mathrm{‾}}{{\mathrm{MS}}^{-}}$] are taken here as rough measures of winter and summer deposition, respectively, across sites. While impurity deposition largely covaries with distance inland (Fig. 1b and c), with higher [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$] and [$\stackrel{\mathrm{‾}}{{\mathrm{MS}}^{-}}$] at sites nearer to the coast, this pattern does not always hold true, likely due to effects associated with seasonal atmospheric transport and (or) local differences in the production of the progenitor aerosols of both ions (Iizuka et al., 2016). There is also a notable geographic difference: Greenland values of [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$] and [$\stackrel{\mathrm{‾}}{{\mathrm{MS}}^{-}}$] (observed at the high-elevation, inland Summit2010 and D4 sites) are lower than those reported for all Antarctic ice cores. In the following subsections, we explore in detail the effects of snow accumulation, ice impurity concentration, surface air temperature, and ice density on MSA migration.

Snow accumulation

Previous work has suggested that high snow accumulation acts as the primary deterrent for MSA migration (Pasteur and Mulvaney, 2000), which is a result of (i) the longer distance required for an ion to travel from the summer to winter layer, (ii) the corresponding suppression of summer–winter layer ionic concentration gradients at higher accumulation sites, or (iii) both (Curran et al., 2002). General support for this suggestion is illustrated in Fig. 1a, showing that sites with $\dot{b}>$ 0.45 m w.e. yr${}^{-\mathrm{1}}$ (meter of water equivalent per year) appear less likely to exhibit clear signs of MSA migration than sites with lower accumulation rates. This consideration, however, does not take into account the depth required for migration to occur. As such, we consider below whether MSA migration invariably occurs across the full range of accumulation rates, given sufficient depth within the firn column.

At low-to-moderate accumulation rate sites ($\dot{b}=$ 0.1–0.45 m w.e. yr${}^{-\mathrm{1}})$, MSA migration seems to universally occur, with the shallowest reported depth of migration showing a positive relationship with accumulation rate (Table 1 and Fig. 2). Data from sites characterized by moderately high accumulation rate ($\dot{b}=$ 0.45–0.65 m w.e. yr${}^{-\mathrm{1}})$ are not as straightforward to interpret. For example, in the Dome Summit South (DSS; Law Dome, Antarctica) and Dyer Plateau ice cores, signs of MSA migration were noted in the original studies. At DSS ($\dot{b}=$ 0.64 m w.e. yr${}^{-\mathrm{1}})$, moderate indications of [MS${}^{-}$] annual maxima extending into adjacent autumn layers were reported at depths greater than about 40–50 m (Curran et al., 2002). Similarly, at Dyer Plateau ($\dot{b}=$ 0.48 m w.e. yr${}^{-\mathrm{1}})$ [MS${}^{-}$] showed suppressed, localized maxima in both winter and summer layers throughout a core section from $\sim$ 51–54 m depth, indicating, perhaps, initial or transitory stages of migration (Pasteur and Mulvaney, 2002). On the other hand, two [MS${}^{-}$] records at sites characterized by annual mean accumulation rates similar to DSS and Dyer Plateau, the 2Barrell ($\dot{b}=$ 0.51 m w.e. yr${}^{-\mathrm{1}})$ and WHG ice cores ($\dot{b}=$ 0.61 m w.e. yr${}^{-\mathrm{1}})$, respectively, showed no signs of migration. However, neither WHG nor the 2Barrell records extend as deep as their counterparts, suggesting WHG and 2Barrell may be of insufficient length and thus may also show migration at similar depths. Finally, of the four sites characterized by high accumulation rates (> 0.65 m w.e. yr${}^{-\mathrm{1}})$, three [MS${}^{-}$] records – the Gomez Nunatak ($\dot{b}=$ 0.88 m w.e. yr${}^{-\mathrm{1}})$, Beethoven Peninsula ($\dot{b}=$ 1.2 m w.e. yr${}^{-\mathrm{1}})$, and DE08 ($\dot{b}=$ 1.27 m w.e. yr${}^{-\mathrm{1}})$ ice cores – did not present signs of migration. The exception is the highest accumulation site, Bruce Plateau ($\dot{b}=$ 1.98 m w.e. yr${}^{-\mathrm{1}})$, where clear evidence of summer-to-winter migration was reported at $\sim$ 395 m depth (Porter et al., 2016). Notably, Bruce Plateau is also the deepest drilled record among the high accumulation rate sites of our compilation.

(a) Shallowest depth of MSA migration, $z_{\mathrm{fo}}$, versus annual mean accumulation rate ($\dot{b}$) for Antarctic (circles) and Greenland sites (diamonds). The dashed line is the least-squares fit $z_{\mathrm{fo}}$ $=$ 125 $\cdot$ ${\dot{b}}^{\mathrm{1.77}}$ through the Antarctic data with $\dot{b}$ < 0.45 m. w.e. yr${}^{-\mathrm{1}}$ ($n$ $=$ 10; $r^{\mathrm{2}}$ $=$ 0.82) and used to predict $z_{\mathrm{fo}}$ at sites where MSA migration was not reported (grey and blue). At two of these sites (WHG and DE08; grey), maximum sampling depth exceeds the predicted $z_{\mathrm{fo}}$. The solid black line is the fit when including the Bruce Plateau data ($\dot{b}$ $=$ 1.98 m. w.e. yr${}^{-\mathrm{1}}$). (b) Expanded view of panel (a) (rectangle bounded with dashed lines). Error estimates for $z_{\mathrm{fo}}$ are generally crude (Supplement Sect. S1). Some are based on unit conversion from m w.e. to m using a firn densification model (Harron and Langway, 1980) constrained by site diagnostic observations (Table 1) and an assumed surface snow density range of 300–400 kg m${}^{-\mathrm{3}}$. The dashed line is the same as in panel (a).

[Figure omitted. See PDF]

We use these observations to determine the relationship between $\dot{b}$ and the depth of first occurrence of MSA migration, $z_{\mathrm{fo}}$, defined as the shallowest reported depth where [MS${}^{-}$] consistently shows its annual maximum in winter (see Supplement for details about the estimation of $z_{\mathrm{fo}}$ for each site). We first consider just Antarctic sites reporting MSA migration (Fig. 2a). A least-squares fit of a power law of $z_{\mathrm{fo}}$ against $\dot{b}$, with the intercept fixed to the origin, yields $r^{\mathrm{2}}=\mathrm{0.99}$ (Fig. 2a). To ensure that the fit is not dominated by the data pair with the highest $\dot{b}$ value (Bruce Plateau), $z_{\mathrm{fo}}$ is regressed against $\dot{b}$ using only those records with accumulation rates < 0.45 m w.e. yr${}^{-\mathrm{1}}$, which yields $r^{\mathrm{2}}=\mathrm{0.82}$ (Fig. 2b). Remarkably, extension of the latter fit to the Bruce Plateau data point produces a value of $z_{\mathrm{fo}}$ that is within $\sim$ 15 m of the reported value (Fig. 2a). The ability of the power law to describe the data in both scenarios suggests that, at least at Antarctic sites, the annual mean rate of snow accumulation has a strong (nonlinear) influence on the shallowest depth of MSA migration in the ice column.

Applying this power law to the five sites where MSA migration was not observed or reported (Fig. 1), we find that three of these sites did not reach the predicted $z_{\mathrm{fo}}$. The exceptions are (i) the WHG core, where little information on MSA migration in the deepest portions of the core is provided in the original study (Sinclair et al., 2014); and (ii) the Law Dome DE08 core, where the largest depth reported by Curran et al. (2002) remains within $\sim$ 15 m of the predicted $z_{\mathrm{fo}}$. These results suggest that MSA migration may be more general than commonly thought and that records for which MSA migration was not reported either may have not reached the requisite depth or they did not have the necessary resolution required to observe the phenomenon. The two high-elevation records from inland Greenland, which have low overall impurity concentration relative to the coastal Antarctic records, appear anomalous (Fig. 2a). This anomaly may result from the effect of ice impurities – in particular cationic sea salts – on MSA migration, a point we investigate in more detail below.

Ice impurities

Earlier studies hypothesized that the presence of well-defined [MS${}^{-}$] peaks in winter layers reflects interaction with sea salts, which are preferentially deposited during winter months at most coastal locations (Legrand and Mayewski, 1997). The attributed mechanism, first hypothesized by Mulvaney et al. (1992) and echoed by subsequent studies (Wolff, 1996; Kreutz et al., 1998; Pasteur and Mulvaney, 2000; Curran et al., 2002), posits that an MSA molecule rejected to and dissolved in the undercooled liquid veins present at the interface of ice grains (Mulvaney et al., 1988) would be transported by diffusion along a concentration gradient. This mechanism could transport MS${}^{-}$ in either the down or up gradient direction and would not necessarily conform to the concentration gradient measured on bulk ice samples (Sect. 4). Mulvaney et al. (1992) further proposed that, upon contact with a cation, the precipitation of a stable cation salt through a metathesis reaction could effectively remove MS${}^{-}$ from solution, thus sequestering the resultant precipitate in the winter layer.

Here we evaluate whether observations support this mechanism, by comparing $z_{\mathrm{fo}}$ to the bulk concentrations of winter-deposited cations at sites showing evidence of MSA migration. Our evaluation is premised by the notion that higher ionic concentrations have a greater ability to regulate the chemical composition of liquid veins, when in situ temperatures exceed the eutectic temperatures of the dominant progenitor salt species of those ions (Nye, 1991; Sect. 4). We constrain our analysis to the relationship between MSA and Na${}^{+}$ concentrations for the following reasons. First, at coastal localities, the most abundant winter-maximum cation is typically Na${}^{+}$, though the relative abundance of Mg${}^{\mathrm{2}+}$, Ca${}^{\mathrm{2}+}$, and K${}^{+}$ may vary from site to site (Legrand and Mayewski, 1997). Second, past studies have noted that MSA concentration peaks tend to coincide with Na${}^{+}$ concentration peaks down-core (Pasteur and Mulvaney, 2000; Kreutz et al., 1998). Finally, Na${}^{+}$ is considered to be relatively nonreactive within the ice column (Barnes et al., 2003a, b; Legrand and Mayewski, 1997), while less abundant cations, such as Mg${}^{\mathrm{2}+}$, appear to be more susceptible to post-depositional effects (Kreutz et al., 1998; Wolff, 1996).

In Fig. 3, $z_{\mathrm{fo}}$ is plotted against the core-averaged concentration of Na${}^{+}$ ([$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$]) for the 11 sites where [Na${}^{+}$] data are available. These sites have comparable annual mean accumulation rates, within 0.13–0.42 m w.e. yr${}^{-\mathrm{1}}$, so that the effect of [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$] on MSA migration could be isolated with some confidence. We find that as [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$] decreases, MSA migration tends to be observed at greater depths in the firn or ice column (Fig. 3). The depth $z_{\mathrm{fo}}$ appears to be particularly sensitive to [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$] for low concentrations of $\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$, as can be seen most prominently for the two Greenland cores showing migration, D4 and Summit2010. This result conforms to the suggestion that, as the concentration of Na${}^{+}$ is reduced, the concentration gradient of MS${}^{-}$ existing in the liquid vein network between the winter and summer layers is also reduced, thereby decreasing the rate at which MSA migration occurs. Another possibility is that the reduction of [Na${}^{+}$] closes off the liquid vein network within winter layers, inhibiting the vertical transport of ions along grain boundaries, although data are lacking to explore this further.

Temperature

Temperature may influence MSA migration in at least two ways. First, the rate of diffusion of MS${}^{-}$ along grain boundaries may depend on temperature, with higher rates occurring at warmer sites, as suggested by Pasteur and Mulvaney (1999). Such dependence would imply that sites characterized by lower in situ temperatures would tend to exhibit larger values of $z_{\mathrm{fo}}$. Second, the precipitation of an MS${}^{-}$ sea salt from solution may be required to sequester MS${}^{-}$ in the winter layer, as proposed by Mulvaney et al. (1992). If effective, the formation of winter-layer [MS${}^{-}$] annual maxima would be inhibited at sites where in situ temperature exceeds the eutectic temperature(s) of the dominant precipitated MS${}^{-}$ sea salt(s). We can test both of the foregoing thermal influences with the available observations. However, due to a lack of firn temperature profiles, we substitute the more widely available annual mean surface air temperature as a surrogate for in situ temperature. This approach appears justified, as vertical thermal gradients in polar ice sheets are typically small below 10 m and their magnitude is dictated primarily by local annual mean SAT (Cuffey and Patterson, 2010). We note that, deeper in the ice, thermal gradients could be important (Rempel et al., 2001), particularly near the bedrock due to geothermal fluxes and (or) frictional heat dissipation, but we do not consider those cases here.

Shallowest depth of MSA migration versus core-averaged [Na${}^{+}$] for different core sites in Antarctica (circles) and Greenland (diamonds). The dashed line is the least-squares fit $z_{\mathrm{fo}}$ $=$ 142 $\cdot$ [Na${}^{+}{]}^{-\mathrm{0.58}}$ through all the data ($n$ $=$ 11, $r^{\mathrm{2}}$ $=$ 0.85).

[Figure omitted. See PDF]

As a rough test of a thermal influence on MSA migration, we regress $z_{\mathrm{fo}}$ against SAT (not shown). Unlike the relationships between $z_{\mathrm{fo}}$ and $\dot{b}$ or [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$], no significant relationship between $z_{\mathrm{fo}}$ and SAT is observed. Likewise, no significant relationship is found between the $z_{\mathrm{fo}}$ residuals of the relationship between $z_{\mathrm{fo}}$ and $\dot{b}$ (arguably, the best predictor for $z_{\mathrm{fo}})$ and SAT (also not shown). These results suggest that temperature has a negligible influence on MSA migration, at least for the sites considered (Table 1).

Similarly, available observations do not indicate that sites reporting MSA migration can be discriminated on the basis of SAT. Although sites presenting no evidence for MSA migration tend to experience high SAT ($-$19 to $-$12.5 ${}^{\circ }$C, Fig. 1a), two sites where MSA migration has been reported are also characterized by equal or higher SAT. We interpret this observation as an artifact of the strong relationship between SAT and accumulation rate in polar regions (Fig. 1), whereupon $\dot{b}$ is deemed to be the most relevant variable in driving $z_{\mathrm{fo}}$ (Sect. 2.1). This conclusion, in turn, indicates that either precipitating MS${}^{-}$ salts may not drive MSA migration (Mulvaney et al., 1992) or the eutectic temperature(s) of the dominant MS${}^{-}$ sea salt(s) is (are) less than $-$29.5 ${}^{\circ }$C, the lowest SAT reported in Table 1. We return to this discussion in Sect. 5.1, after incorporating recent estimates of the eutectic temperatures of the dominant MS${}^{-}$ sea salts.

Firn or ice density

The mechanism of MSA migration described above, involving diffusion along liquid grain boundaries, relies on the assumption that the ice grains are sufficiently well compacted to form interconnected premelted veins between the summer and winter layers. This mechanism seems to imply there should be a relationship between the onset of MSA migration and ice density ($\mathit{\rho })$, yet to our knowledge such possible relationship has not been investigated. Here, we consider the firn or ice densities (${\mathit{\rho }}_{\mathrm{fo}})$ observed at the shallowest depth of MSA migration again using our compilation of Greenland and Antarctic data.

In Fig. 4, $z_{\mathrm{fo}}$ is plotted against ${\mathit{\rho }}_{\mathrm{fo}}$ for the seven sites where available data permit (Table 1). A weighted least-squares fit to the data, with the weighting provided by the errors in the $z_{\mathrm{fo}}$ estimates, shows a positive relationship between $z_{\mathrm{fo}}$ and ${\mathit{\rho }}_{\mathrm{fo}}$ ($p$ < 0.01), as expected. Data for low-accumulation Berkner Island sites, where annual mean $\dot{b}=\mathrm{0.18}$–0.22 m w.e. yr${}^{-\mathrm{1}}$, indicate the occurrence of MSA migration at densities exceeding about 515–560 kg m${}^{-\mathrm{3}}$ (Wagenbach et al., 1994), while higher accumulation sites, including THW2010, Ferrigno, and DIV2010 ($\dot{b}=\mathrm{0.28}$–0.41 m w.e. yr${}^{-\mathrm{1}})$, show ${\mathit{\rho }}_{\mathrm{fo}}$ values in the range of 610–650 kg m${}^{-\mathrm{3}}$. The highest ${\mathit{\rho }}_{\mathrm{fo}}$ values, in the range 710–760 kg m${}^{-\mathrm{3}}$, are found for the Greenland cores, D4 and Summit2010.

Interestingly, the firn or ice densities at the shallowest depths of MSA migration ($\sim$ 500 to $\sim$ 800 kg m${}^{-\mathrm{3}}$, Fig. 4), are all comparable to or higher than the critical value of 550 kg m${}^{-\mathrm{3}}$ that corresponds to the theoretical closest random packing of spherical ice grains (Benson, 1962). Densities greater than 550 kg m${}^{-\mathrm{3}}$ can be achieved only by bond formation, or sintering, at the contact of individual grains (Cuffey and Patterson, 2010). Persistent vertical transport of undercooled liquid at the grain boundaries (Mulvaney et al., 1992) requires reaching or exceeding this effective critical density. The compiled data also indicate that the critical density for MSA migration may vary between sites (Fig. 4), which might be due to variability in grain shape, grain size, and (or) impurity content. For example, firn samples with densities of 350–400 kg m${}^{-\mathrm{3}}$ can typically be cut into blocks (Cuffey and Patterson, 2010), which suggests that bond formation in the upper firn pack can begin at densities below 550 kg m${}^{-\mathrm{3}}$. Furthermore, at the lowest accumulation sites – Byrd Station, Siple Dome, and the Filchner-Ronne Ice Shelf – MSA migration was reported in the shallowest $\sim$ 5 m of the firn ($z_{\mathrm{fo}}=\mathrm{2.6}$, 2, and 3–4 m, respectively). Although density data are not available for these specific records, this observation suggests that MSA migration may begin at bulk densities substantially lower than 550 kg m${}^{-\mathrm{3}}$, perhaps through high-density microlayers, such as wind-blown crusts, which are features commonly found at lower accumulation sites (Cuffey and Patterson, 2010). While more data are needed, available evidence presented here supports the notion that a density of at least 500 kg m${}^{-\mathrm{3}}$, comparable to the theoretical value of 550 kg m${}^{-\mathrm{3}}$, is needed for MSA migration.

Shallowest depth of MSA migration versus firn or ice density for different core sites in Antarctica (circles) and Greenland (diamonds). The vertical bars are error estimates (Supplement Sect. S1); the dashed line is the weighted least-squares fit $z_{\mathrm{fo}}=-\mathrm{956}+\mathrm{0192}\mathit{\rho }$ through the data ($n$ $=$ 7; $r^{\mathrm{2}}$ $=$ 0.89), with the weighting provided by the $z_{\mathrm{fo}}$ error estimates; and the shaded region is the region of 95 % confidence.

[Figure omitted. See PDF]

Synthesis

Our data analysis suggests that the annual mean rate of snow accumulation appears to have a strong influence on the shallowest depth at which MSA migration is observed. The concentration of Na${}^{+}$ in the ice or firn appears to also be an important factor, especially at sites characterized by low [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$]. Annual mean SAT appears to play a less important role in determining $z_{\mathrm{fo}}$, at least in the SAT range where data are available ($-$29.5 to $-$12.5 ${}^{\circ }$C). Lastly, the onset of migration appears to be associated with a critical density near 450–550 kg m${}^{-\mathrm{3}}$, though this may not hold true for all low-accumulation sites. Overall, our analysis suggests that MSA migration may be more common than usually thought and that existing [MS${}^{-}$] records not exhibiting evidence for migration may not be deep enough, or have the necessary (sub-annual) resolution, to observe the phenomenon.

The shallowest depth of MSA migration ($z_{\mathrm{fo}})$ appears most readily predictable from $\dot{b}$ in coastal Antarctica. This result appears in part attributable to the comparably higher concentrations of Na${}^{+}$ in the coastal Antarctic cores than in the two Greenland cores showing evidence of migration (Table 1). The joint effect of annual mean accumulation rate and depth-averaged Na${}^{+}$ concentration on the shallowest depth of MSA migration is illustrated in Fig. 5, which shows MSA migration tends to occur deeper in the firn or ice column for larger values of $\dot{b}$ and lower values of [$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$]. Notably, Antarctic sites with comparably low concentrations of Na${}^{+}$ tend to occur further inland and at higher elevations, and thus also tend to have lower accumulation rates than considered in this study (i.e., $\dot{b}$ < 0.1 m w.e. yr${}^{-\mathrm{1}})$. Greenland, in contrast, experiences comparably large accumulation rates even at its highest-altitude, inland locations (e.g., Summit2010). For a recent overview of Greenland-wide [Na${}^{+}$] deposition, see Rhodes et al. (2017).

A case study: the DIV2010 MSA record

In this section, we use a well-dated multi-century ice core record from DIV2010, an intermediate accumulation rate site in coastal West Antarctica (Table 1), to document in detail the phenomenon of MSA migration. Prior studies have investigated aspects of the DIV2010 core, including variability in accumulation (Medley et al., 2013, 2014) and chemical composition (Criscitiello et al., 2013, 2014; Pasteris et al., 2014). As this is the first report of the DIV2010 [MS${}^{-}$] record below the zone of MSA migration, however, we describe the record in more detail here. Inorganic salt ions (e.g., Na${}^{+})$ and MS${}^{-}$ concentrations were measured on discrete samples at a constant sampling interval $\mathrm{\Delta }z=\mathrm{5}$ cm using standard suppressed ion chromatography methods (Curran and Palmer, 2001). The age-depth relationship of the core was established independently by identifying summer maxima in three parameters – [nssS], [H${}_{\mathrm{2}}$O${}_{\mathrm{2}}$], and ${\mathit{\delta }}^{\mathrm{18}}$O – measured at $\sim$ 2 cm resolution (Pasteris et al., 2014). Here, we examine data from the top 60.4 m of DIV2010 (AD 1905–2010) covering the zone of progressive migration of MSA. Dating uncertainty over this depth interval is estimated to be less than 1 year based on tie points to well-known volcanic events (Pasteris et al., 2014), and the sampling frequency ($\sim$ 10–15 samples per year) allows for seasonally resolved records throughout this depth interval (Criscitiello et al., 2013; Criscitiello, 2014). At a depth of 60.4 m in the core, annual layer thinning estimated using a thinning model assuming frozen basal conditions and a linear vertical strain rate (Nye, 1963) is considerably less than the distance required for a chemical species to migrate from the summer to the winter layer ($\sim$ 30 cm). As a result, no correction for thinning is applied.

Core-averaged [Na${}^{+}$] versus annual mean accumulation rate for different core sites in Antarctica and Greenland. The different colors correspond to different values of the shallowest depth of MSA migration.

[Figure omitted. See PDF]

The migration of MSA from summer to winter layers in the DIV2010 ice core appears to be progressive. To document this progression and contrast the behavior of Na${}^{+}$ and MS${}^{-}$, the month of the annual maximum of [MS${}^{-}$] ($m_{{\mathrm{MS}}^{-}})$ and the month of the annual maximum of [Na${}^{+}$] ($m_{{\mathrm{Na}}^{+}})$ are each plotted versus age (Fig. 6). The month of the annual maximum of [MS${}^{-}$] tends to change down-core from summer to predominantly winter, as revealed by the significant linear trend of $m_{{\mathrm{MS}}^{-}}$ between AD 1905 and 1999 ($r=\mathrm{0.75}$; $p$ < 0.001). In contrast, [Na${}^{+}$] consistently shows maxima in the winter layers; i.e., $m_{{\mathrm{Na}}^{+}}$ portrays no significant trend over this period ($r=\mathrm{0.09}$, $p$ $=$ 0.41). Similarly, [nssSO${}_{\mathrm{4}}^{\mathrm{2}-}$] displays its annual maximum consistently in the summer layers during this period (not shown). Note that no significant decrease in annual mean [MS${}^{-}$] is found deeper than 10 m, suggesting post-depositional losses are negligible at the DIV2010 core site.

Timing of annual maximum of [MS${}^{-}$] (red) and annual maximum of [Na${}^{+}$] (blue) versus calendar age (year AD) for the DIV2010 core. No obvious sign of migration in the up or down-core direction is observed in DIV2010. As such, the timing of annual maximum [MS${}^{-}$] and annual maximum [Na${}^{+}$] is defined in terms of the number of degrees out of phase from 1 January (0${}^{\circ }$), where 180${}^{\circ }$ indicates an annual maximum on 1 July. The two dashed lines are the least-squares fits for [MS${}^{-}$] and [Na${}^{+}$] over the period AD 1999–1905, corresponding to depths below which $\mathit{\rho }$ $=$ 550 kg m${}^{-\mathrm{3}}$. The fit is highly significant for [MS${}^{-}$] ($n$ $=$ 95, $r$ $=$ 0.75, $p$ < 0.001) and not significant for [Na${}^{+}$] ($n$ $=$ 95, $r$ $=$ 0.09, $p$ $=$ 0.41).

[Figure omitted. See PDF]

Records of [MS${}^{-}$] (red) and [Na${}^{+}$] (blue) from the DIV2010 ice core. (a) The entire record considered in this study, with a depth scale at the top and a timescale (year AD) at the bottom (raw data). (b–d) The 11-year long portions of the records within the shallow zone (b), the transition zone (c), and the deep zone (d). Three-point running averages are displayed for each zone. Dashed vertical lines denote 1 January of each year.

[Figure omitted. See PDF]

Monthly mean concentration of Na${}^{+}$ and MS${}^{-}$ in the shallow (a), transition (b), and deep (c) zones of DIV2010. The shaded regions indicate $\pm$1 standard error of the mean.

[Figure omitted. See PDF]

(a–c) Correlograms of [MS${}^{-}$] and [Na${}^{+}$] in the three zones of the DIV2010 ice core. The horizontal dashed lines show the 95 % confidence interval for a time series with no autocorrelation (white sequence). (d–f) Cross-correlation between [MS${}^{-}$] and [Na${}^{+}$] in the three zones of the DIV2010 ice core. The horizontal dashed lines show the 95 % confidence interval for two uncorrelated time series.

[Figure omitted. See PDF]

Identification of DIV2010 migration zones

Notwithstanding the progressive nature of MSA migration in the DIV2010 core, we identify three distinct zones in the [MS${}^{-}$] record for this core (Fig. 7). First, a shallow zone is defined as the depth interval where density is less than 550 kg m${}^{-\mathrm{3}}$; at DIV2010, this depth occurs at approximately 9.1 m (Medley et al., 2014). The upper 9.1 m of the core appears to contain the original (unaltered by migration) [MS${}^{-}$] variations, to the extent that annual maxima of [MS${}^{-}$] are found in the summer layers and are out of phase with winter [Na${}^{+}$] maxima. Second, a transition zone is defined further down-core, where the [MS${}^{-}$] record exhibits no consistent seasonality. Finally, a deep zone is defined as the deepest portion of the record considered, where [MS${}^{-}$] and [Na${}^{+}$] annual maxima appear to be broadly in phase (Fig. 7). In order to facilitate our analysis, we linearly interpolated the [MS${}^{-}$] record onto a monthly scale, with the three zones defined to entail an equal number of data points ($n=\mathrm{132}$ months, or 11 years).

To better reveal the variable phase relationship between [MS${}^{-}$] and [Na${}^{+}$] in the shallow, transition, and deep zones, monthly mean values of [MS${}^{-}$] and [Na${}^{+}$] are calculated and the standard error of the monthly means are computed for each zone (Fig. 8). In all three zones, Na${}^{+}$ shows concentration maxima in winter, albeit with broad variation between April and September. In contrast, MS${}^{-}$ exhibits concentration maxima during summer (January) in the shallow zone but during winter (broadly, May to August) in the deep zone. In the transition zone, the monthly mean [MS${}^{-}$] values are not significantly different between summer (DJF) and winter (JJA), though local maxima in both spring (MAM) and fall (SON) are apparent.

Cross-correlation analysis

Cross-correlation coefficients are calculated to quantify the degree of linear relationship between the monthly mean [MS${}^{-}$] and [Na${}^{+}$] at different lags in the three zones. Due to the finite length of the records, only the coefficients from lag 1 to lag 24 (2 years) are calculated (Fig. 9). In the shallow zone, [MS${}^{-}$] and [Na${}^{+}$] show a negative correlation at lag 0 and a positive correlation at lag 6 (i.e., 6 months); both correlations appear significant at the 5 % level compared to two uncorrelated series (Chatfield, 1996). These results reflect MSA and Na${}^{+}$ being deposited primarily during summer and winter, respectively. In contrast, in the transition zone, and more prominently in the deep zone, the [MS${}^{-}$] and [Na${}^{+}$] records show a positive correlation at lag 0 and a negative correlation at lag 6; in the deep zone the correlations at lags 0 and 6 are both significant at the 5 % level. The positive correlation at zero lag in the transition and deep zones indicates that positive deviations in [MS${}^{-}$] tend to coincide with positive deviations in [Na${}^{+}$]. The distinct phase relationships between the two species in the shallow and deep zones are consistent with MSA migration.

Towards a mechanistic understanding of MSA migration

In this section, we discuss the physico-chemical processes that may be responsible for MSA migration using an existing model of soluble impurity transport. We then derive a linearized version of the model in order to further illuminate the processes that lead to the movement of MSA from summer to winter layers.

The impurity transport model of Rempel et al. (2002)

The physical mechanisms responsible for the presence of a liquid phase in ice cores, even well below the freezing point of pure water, are well documented (e.g., Nye and Frank, 1973; Nye, 1991) and attributed to two distinct physical processes. The first relates to the fact that the atomic radii of most impurity species (not including requisitely small ionic species, e.g., F${}^{-}$, Cl${}^{-}$, and NH${}_{\mathrm{4}}^{+}$; Wolff, 1996) possess a misfit strain energy that inhibits their incorporation into the tightly packed, crystalline ice lattice. During densification, these impurities are thus preferentially expelled to the boundaries of individual grains of ice. This process has been observed using optical measurements (Mulvaney et al., 1988; Bartels-Rauch et al., 2014). As the concentration of impurities at the grain boundary is increased, the local equilibrium temperature is decreased, depressing the freezing point of the water-impurity mixture. At temperatures greater than a system's eutectic point, premelted aqueous solutions are assumed to exist at equilibrium as interconnected, submicron veins at the grain boundaries. The second process responsible for the presence of a liquid phase in polar ice pertains to the interstitial curvature occurring at the interface of three (i.e., triple junctures) or four (i.e., nodes) ice grains. Known as the Gibbs–Thompson effect, this thermodynamic phenomenon is related to the deviation in chemical potential of a vapor surrounding a curved surface from that of the same vapor at equilibrium with a flat liquid surface. In effect, it allows smaller (more curved) ice grains of a given composition to melt at lower temperatures than larger (less curved) ice grains of the same composition (Wettlaufer and Worster, 2006).

Importantly, both of the foregoing processes should respond collaterally to a temperature change, although they may not have the same importance in the maintenance of a liquid phase: scaling arguments suggest that the curvature effect required to reach a given volume of premelted liquid is negligible in ice sheets in comparison to the effect of impurity-driven undercooling (Rempel et al., 2001). Focusing on the latter process, Rempel et al. (2002) developed an elegant model describing the movement of soluble impurities along crystal grain boundaries. They considered impurity migration due to (1) temperature gradients (Rempel et al., 2001), which vary gradually down-core (typically < 1 ${}^{\circ }$C per 100 m in ice sheet interiors; Cuffey and Patterson, 2010); and (2) impurity concentration gradients, which are characterized by length scales typically on the order of centimeters (Rempel et al., 2002). They derived the following impurity migration equation, referred to below as the RWW model: $${\displaystyle \frac{\partial }{\partial t}}{c}_{\mathrm{B},k}=-\mathrm{\nabla }\cdot (v+v_k)c_{\mathrm{B},k},$$ where, $$v_k=D_k{\displaystyle \frac{\mathrm{\nabla }T}{\left(T_{\mathrm{m}},-,T\right)}}+D_k{\displaystyle \frac{\sum {\mathrm{\Gamma }}_i\mathrm{\nabla }{c}_{\mathrm{B},i}-\frac{\mathrm{\nabla }{c}_{\mathrm{B},k}}{c_{\mathrm{B},k}}\sum {\mathrm{\Gamma }}_i{c}_{\mathrm{B},i}}{\sum {\mathrm{\Gamma }}_i{c}_{\mathrm{B},i}}}.$$ Here $c_{\mathrm{B},k}$ is the bulk concentration of the $k$th impurity species, i.e., the mass of the $k$th impurity species per unit ice volume as measured during standard chemical analyses; $t$ is time; and $v_k$ is an effective velocity of the $k$th impurity species relative to the surrounding ice. The ice velocity, $v$, would arise, for example, from a vertical strain rate and appears to have a small influence on MSA transport on timescales of years to decades (Nye, 1963). As a result, it is systematically neglected in this paper (see also, Rempel et al., 2002).

As indicated by Eq. (2), the relative velocity $v_k$ has two distinct contributions. The first is due to the motion of molecules in the liquid along a temperature gradient and leads to solute transport, even in the absence of concentration gradients. It is proportional to the diffusivity of the $k$th impurity in the liquid, $D_k$, and inversely proportional to the difference $T_{\mathrm{m}}-T$ between the melting point of pure ice ($T_{\mathrm{m}})$ and the temperature of the liquid ($T)$. The second contribution to $v_k$ arises from the concentration gradients of all solutes present in the liquid veins, including the $k$th impurity. As with the first contribution, it is proportional to the diffusivity of the $k$th impurity in the liquid, but it would occur in the absence of a temperature gradient. It also depends on the slope of the liquidus curve, ${\mathrm{\Gamma }}_i$, of the various solutes that are present in the liquid veins. Scaling arguments suggest that $v_k$ is often well approximated by the second contribution in Eq. (2), at least in portions of the ice column where temperature gradients are small (Rempel et al., 2002). The small effect of the first contribution on MSA migration is further supported by observations made in Sect. 2.3.

Comparison between the RWW model and the linearized model. (a) Profiles of [MS${}^{-}$] and [Na${}^{+}$] simulated by the RWW model with $D_{\mathrm{MS}}=D_{\mathrm{Na}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ and ${\mathrm{\Gamma }}_{\mathrm{MS}}={\mathrm{\Gamma }}_{\mathrm{MS}}=\mathrm{6.5}$ K M${}^{-\mathrm{1}}$. (b) Profiles of [MS${}^{-}$] simulated by the linearized model with an effective velocity based on the [Na${}^{+}$] profile shown in the panel and with $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ (see text and Supplement for details).

[Figure omitted. See PDF]

A binary mixture provides the simplest context to discuss the mechanisms of MSA migration in the RWW model. In a mixture comprising MS${}^{-}$ and Na${}^{+}$, Eqs. (1)–(2) reduce to (for negligible $v)$ $$\begin{array}{}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}{c}_{\mathrm{MSA}}=-\mathrm{\nabla }\cdot \left(v_{\mathrm{MS}},c_{\mathrm{MS}}\right),& {\displaystyle }& {\displaystyle \frac{\partial }{\partial t}}{c}_{\mathrm{Na}}=-\mathrm{\nabla }\cdot \left(v_{\mathrm{Na}}{c}_{\mathrm{Na}}\right),\end{array}$$ where, $$\begin{array}{}& {\displaystyle }& {\displaystyle }{v}_{\mathrm{MSA}}={\displaystyle \frac{D_{\mathrm{MS}}{\mathrm{\Gamma }}_{\mathrm{Na}}}{{\mathrm{\Gamma }}_{\mathrm{MS}}{c}_{\mathrm{MS}}+{\mathrm{\Gamma }}_{\mathrm{Na}}{c}_{\mathrm{Na}}}}\left(\mathrm{\nabla },c_{\mathrm{Na}},-,{\displaystyle \frac{c_{\mathrm{Na}}}{c_{\mathrm{MS}}}},\mathrm{\nabla },c_{\mathrm{MS}}\right),& {\displaystyle }& {\displaystyle }{v}_{\mathrm{Na}}={\displaystyle \frac{D_{\mathrm{Na}}{\mathrm{\Gamma }}_{\mathrm{MS}}}{{\mathrm{\Gamma }}_{\mathrm{MS}}{c}_{\mathrm{MS}}+{\mathrm{\Gamma }}_{\mathrm{Na}}{c}_{\mathrm{Na}}}}\left(\mathrm{\nabla },c_{\mathrm{MS}},-,{\displaystyle \frac{c_{\mathrm{MS}}}{c_{\mathrm{Na}}}},\mathrm{\nabla },c_{\mathrm{Na}}\right).\end{array}$$ In this case, the model implicitly assumes that an appreciable amount of the MSA and Na${}^{+}$-containing impurities are rejected from the crystalline lattice of individual ice grains during firnification. The postulated expulsion would concentrate impurities at the grain boundaries, depressing the freezing point of the intergranular medium to form submicron, undercooled liquid veins. In general, higher impurity concentrations would lead to higher abundance of premelt liquid. Thus, more premelt liquid is predicted to occur in winter layers, where [Na${}^{+}$] is typically maximum, than in summer layers, where the comparatively low (bulk) concentration of MSA shows a maximum. The ensuing network of liquid veins would allow the ionic impurities to diffuse under their own concentration gradients, such that a large proportion of the MS${}^{-}$ from the MS${}^{-}$-rich summer layer migrates to the MS${}^{-}$-poor winter layer. Conversely, a comparatively small proportion of Na${}^{+}$ migrates to the summer layer, because the Na${}^{+}$ concentration difference between the summer and winter layers is reduced by the larger amount of premelt liquid in the winter layer than in the summer layer. The net result of the different transport rates of MS${}^{-}$ and Na${}^{+}$ is that variations in [MS${}^{-}$] ultimately become in phase with variations of [Na${}^{+}$] (Fig. 10a).

Physico-chemical parameters of MSA migration

The RWW model as applied to the system containing MS${}^{-}$ and Na${}^{+}$ (Eqs. 3–4) includes four parameters: the slopes of the liquidus curve for relevant MS${}^{-}$- and Na${}^{+}$-containing soluble impurity species (${\mathrm{\Gamma }}_{\mathrm{MS}}$ and ${\mathrm{\Gamma }}_{\mathrm{Na}})$, and the grain-boundary diffusivities of MS${}^{-}$ and Na${}^{+}$ ($D_{\mathrm{MS}}$ and $D_{\mathrm{Na}})$. Below, we review the existing literature on each of these quantities.

Grain boundary diffusivity of MS${}^{-}$ and Na${}^{+}$

We first consider the grain boundary diffusion coefficient of MS${}^{-}$, $D_{\mathrm{MS}}$. Lacking empirical constraints, Rempel et al. (2002) approximated the diffusivity for ionic constituents in Eq. (2) as one-third the molecular diffusivity of a bulk liquid (i.e., $D_{\mathrm{MS}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{10}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}})$, scaled so as to account for the random orientation of premelted liquid veins in the ice (Lemlich, 1978). Smith et al. (2004) reported a value of $D_{\mathrm{MS}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{13}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ for solid ice, estimated by measuring variations in [MS${}^{-}$] across horizontal sections of an ice core from Law Dome, Antarctica, following nearly 15 years of freezer storage at $-$20 ${}^{\circ }$C. Using a similar ice substrate and experimental setup, Roberts et al. (2009) revised this estimate to ($\mathrm{4.1}\times {\mathrm{10}}^{-\mathrm{13}}\pm \mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{14}})$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ at $-$20 ${}^{\circ }$C, which was interpreted by the authors to represent diffusive losses of volatile MSA occurring during extended periods of freezer storage. Notably, this estimate is 1–3 orders of magnitude larger than that reported for solid-state diffusion of HCl (Thibert and Dominé, 1997), HNO${}_{\mathrm{3}}$ (Thibert and Dominé, 1998), HCHO (Barret et al., 2011), and deuteriorated water (Lu et al., 2009) determined in single ice crystals, despite the molecular radius of MSA greatly exceeding that of each of these species (Roberts et al., 2009). In fact, subsequent studies have contended that the $D_{\mathrm{MS}}$ estimate of Roberts et al. (2009) is unlikely to represent pure solid-state diffusion of MSA in firn or ice and suggested that at least some of the storage-related losses of MSA occurred via liquid transport along grain boundaries (McNeil et al., 2012; Bartels-Rausch et al., 2014). We thus consider the two values of $D_{\mathrm{MS}}$ as suggested by Rempel et al. (2002) and Roberts et al. (2009) as potential extrema.

Although the grain boundary diffusion coefficient for Na${}^{+}$, $D_{\mathrm{Na}}$, is also under-constrained, empirical evidence supports a relative immobility of Na${}^{+}$ in polar ice. For example, Barnes et al. (2003a) noted no detectable changes in the amplitude of [Na${}^{+}$] peaks over the past $\sim$ 11 000 years (top 350 m) in the low-accumulation Dome C ice core record (East Antarctica), while the amplitudes of both [Cl${}^{-}$] and [SO${}_{\mathrm{4}}^{\mathrm{2}-}$] peaks were found to change over the same period. Furthermore, optical measurements at Dome C suggest a predisposition for Na${}^{+}$ to be situated at grain boundaries and for Cl${}^{-}$ to be located preferentially within the crystalline structure (Barnes et al., 2003b). The importance of these findings is twofold: (i) Na${}^{+}$ appears to be situated in the requisite location to favor the presence of premelt liquid at the grain boundaries, enabling MS${}^{-}$ migration as envisioned in the RWW model to occur; and (ii) Na${}^{+}$ shows greatly reduced mobility relative to MS${}^{-}$ (or similar sulfur-based acidic species).

Liquidus relationships for relevant sea-salt species

The slopes of the liquidus curves, ${\mathrm{\Gamma }}_{\mathrm{MS}}$ and ${\mathrm{\Gamma }}_{\mathrm{Na}}$, represent linear approximations of the undercooling as a function of impurity concentration in the liquid phase present near the grain boundaries, $c$, i.e., $T_{\mathrm{m}}-T=\mathrm{\Gamma }c$, where $T$ is the in situ temperature and $T_{\mathrm{m}}$ is the melting point for pure ice (see Supplement, Fig. S2). Knowledge of $\mathrm{\Gamma }$ requires knowledge of the dominant precursor (bonded) molecular state(s) of the MS${}^{-}$ and Na${}^{+}$ ions present in the ice (thus, ${\mathrm{\Gamma }}_{\mathrm{MS}}$ and ${\mathrm{\Gamma }}_{\mathrm{Na}}$ should be viewed as shorthand notations for ${\mathrm{\Gamma }}_{\mathrm{MS}\ast }$ and ${\mathrm{\Gamma }}_{\mathrm{Na}\ast }$, where $\ast$ represents some unknown cationic–anionic pair). Unfortunately, such data remain sparse beyond those reported in a few notable studies (e.g., Barnes et al., 2003b; Sakurai et al., 2010; Iizuka et al., 2016).

It is generally assumed that all measured MS${}^{-}$ present in polar firn or ice samples derives solely from MSA (Sakurai et al., 2010). In the binary system MSA–H${}_{\mathrm{2}}$O, MSA reaches its eutectic temperature at $-$75 ${}^{\circ }$C (Stephen and Stephen, 1963). Thus, any MSA molecules expelled to and concentrated at grain boundaries are expected to exist in liquid solution with H${}_{\mathrm{2}}$O. By contrast, Na${}^{+}$ in polar ice may have a number of precursors. For coastal ice cores, however, it seems reasonable to expect that the majority of Na${}^{+}$ is deposited either as NaCl derived primarily from sea spray during storm activity (Legrand and Mayewski, 1997) or as sodium-sulfate salts such as mirabilite, Na${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ $\cdot$ 10H${}_{\mathrm{2}}$O, derived from brine rejection in sea ice or from atmospheric sea-salt sulfatization (Rankin et al., 2002; Iizuka et al., 2016). While the NaCl–H${}_{\mathrm{2}}$O system reaches its eutectic at $-$21.3 ${}^{\circ }$C (Stephen and Stephen, 1963), the eutectic of the Na${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$–H${}_{\mathrm{2}}$O system is $-$1.6 ${}^{\circ }$C (Hougen et al., 1954), suggesting that Na${}^{+}$ deposited as Na${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ should be relatively immobile at most polar ice core sites. Consequentially, the majority of Na${}^{+}$ relevant to grain boundary migration is likely derived from NaCl. At DIV2010, for example, the molar ratio Cl:Na in the top 60.4 m of the core averages 1.806, similar to the mean molar ratio Cl : Na $=$ 1.8 of seawater (e.g., Chesselet et al., 1972). This indicates a primary marine source of NaCl aerosols at DIV2010, as both brine rejection and sea salt sulfatization in the atmosphere would tend to produce an offset in the amount of Cl${}^{-}$ deposited (Iizuka et al., 2016).

Slopes of the liquidus curves ($\mathrm{\Gamma }$) estimated for various binary mixtures composed of impurity species of likely relevance and water. All values listed in columns prior to the $\mathrm{\Gamma }$ column are required for calculation of $\mathrm{\Gamma }$ (see Supplement Sect. S2).

Species	Mol. mass	Eutectic	Eutectic	Density at	$\mathrm{\Gamma }$	Reference
	(g mol${}^{-\mathrm{1}}$)	temperature (${}^{\circ }$C;	composition (wt %;	eutectic (g mL${}^{-\mathrm{1}}$;	(K M${}^{-\mathrm{1}})$
		binary with H${}_{\mathrm{2}}$O)	binary with H${}_{\mathrm{2}}$O)	${}^{\ast }$ $=$ approximated)
CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$H	96.11	$-$75.0	51.1	1.20	11.7	Stephen and Stephen (1963)
Na(CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}})$	119.11	$-$29.3	47	1.15${}^{\ast }$	6.5	Sakurai et al. (2010)
Mg(CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}{)}_{\mathrm{2}}$	214.50	$-$5.0	14.2	1.15${}^{\ast }$	6.6	Sakurai et al. (2010)
Ca(CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}{)}_{\mathrm{2}}$	230.27	$-$32.6	47	1.20${}^{\ast }$	12.9	Sakurai et al. (2010)
NaCl	58.44	$-$21.3	23.3	1.16	4.6	Stephen and Stephen (1963)
MgCl${}_{\mathrm{2}}$	95.21	$-$33.0	21.6	1.13	12.9	Stephen and Stephen (1963)
CaCl${}_{\mathrm{2}}$	110.98	$-$51.0	30	1.19	15.9	Stephen and Stephen (1963)
H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$	98.08	$-$62.0	35.6	1.19	14.3	Hornung et al. (1956)
Na${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$	142.04	$-$1.6	4.0	1.12	5.5	Hougen et al. (1954)
MgSO${}_{\mathrm{4}}$	120.37	$-$3.6	17.3	1.22	2.3	Marion and Farren (1999)
CaSO${}_{\mathrm{4}}$	136.14	$-$0.7	18.0	1.23	0.4	Rolnick (1954)

Although the annual mean surface air temperature at most sites showing MSA migration is observed to be well below $-$21.3 ${}^{\circ }$C (Sect. 2.3), the upper firn (< 10–20 m) undergoes seasonal temperature fluctuations that may exceed this value during summer months (Cuffey and Paterson, 2010), thereby providing a potential mechanism to temporarily free Na${}^{+}$ from its bonded state with Cl${}^{-}$. When premelted liquid solutions containing Na${}^{+}$ and Cl${}^{-}$ refreeze, Cl${}^{-}$ may thus be preferentially allocated within the ice structure (Tokumasu et al., 2016; Barnes et al., 2003b).

With Na${}^{+}$ and MS${}^{-}$ both situated at the grain boundaries, the resulting binary system is the sodium-salt of MS${}^{-}$, CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na $\cdot$ $n$H${}_{\mathrm{2}}$O, and water (Mulvaney et al., 1992). Recent experimental data indicate that the eutectic temperature for the CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O system occurs at approximately $-$29.3 ${}^{\circ }$C (Sakurai et al., 2010). For comparison, the eutectics for the binary systems Ca(CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}{)}_{\mathrm{2}}$ $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O and Mg(CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}{)}_{\mathrm{2}}$ $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O amount to $-$32.6 and $-$5.0 ${}^{\circ }$C, respectively (Sakurai et al., 2010). Table 2 lists the slopes of the liquidus curves for these MS${}^{-}$ salts in addition to those for alternative relevant sea salts containing Ca${}^{\mathrm{2}+}$, Mg${}^{\mathrm{2}+}$, and (or) SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (see also Supplement Sect. S2). For the system CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O the slope amounts to 6.5 K M${}^{-\mathrm{1}}$.

A simplified model of MSA migration

Although the RWW model provides significant insight into the mechanisms of MSA migration, the system of nonlinear partial different Eqs. (3)–(4) does not permit a straightforward analysis (e.g., no closed form solution of these equations with general initial and boundary conditions is available to our knowledge). In this section, we develop a linearized version of the model for the binary system comprising CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na and H${}_{\mathrm{2}}$O (Eqs. 3–4) in order to further our understanding of MSA migration. Of course, the insight to be gained is only as reliable as the assumptions upon which the linearized model relies.

Consider the governing Eq. (3a–b), making it explicit that concentration gradients are strictly vertical, $$\begin{array}{}& {\displaystyle }& {\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial t}}=-{\displaystyle \frac{\partial }{\partial z}}\left(w_{\mathrm{MS}},c_{\mathrm{MS}}\right),& {\displaystyle }& {\displaystyle \frac{\partial c_{\mathrm{Na}}}{\partial t}}=-{\displaystyle \frac{\partial }{\partial z}}\left(w_{\mathrm{Na}},c_{\mathrm{Na}}\right),\end{array}$$ where, $$\begin{array}{}& {\displaystyle }& {\displaystyle }{w}_{\mathrm{MS}}=D_{\mathrm{MS}}{\displaystyle \frac{{\mathrm{\Gamma }}_{\mathrm{Na}}}{{\mathrm{\Gamma }}_{\mathrm{Na}}{c}_{\mathrm{Na}}+{\mathrm{\Gamma }}_{\mathrm{MS}}{c}_{\mathrm{MS}}}}\left({\displaystyle \frac{\partial c_{\mathrm{Na}}}{\partial z}},-,{\displaystyle \frac{c_{\mathrm{Na}}}{c_{\mathrm{MS}}}},{\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}}\right),& {\displaystyle }& {\displaystyle }{w}_{\mathrm{Na}}=D_{\mathrm{Na}}{\displaystyle \frac{{\mathrm{\Gamma }}_{\mathrm{MS}}}{{\mathrm{\Gamma }}_{\mathrm{Na}}{c}_{\mathrm{Na}}+{\mathrm{\Gamma }}_{\mathrm{MS}}{c}_{\mathrm{MS}}}}\left({\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}},-,{\displaystyle \frac{c_{\mathrm{MS}}}{c_{\mathrm{Na}}}},{\displaystyle \frac{\partial c_{\mathrm{Na}}}{\partial z}}\right).\end{array}$$ Here, $w_{\mathrm{MS}}$ and $w_{\mathrm{Na}}$ are the vertical components of MS${}^{-}$ and Na${}^{+}$ migration velocity, respectively, and $z$ is depth. Three assumptions are made (see also Rempel et al., 2002). First, the slope of the liquidus curve is taken to be the same for the two ionic species, i.e., ${\mathrm{\Gamma }}_{\mathrm{Na}}={\mathrm{\Gamma }}_{\mathrm{MS}}$ (Sect. 4.2.2). This assumption appears plausible if the MS${}^{-}$ salt species CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na $\cdot$ $n$H${}_{\mathrm{2}}$O dominates in the premelt liquid present near the grain boundaries. With this assumption, the slopes of the liquidus curves cancel out in the defining relationships for $w_{\mathrm{MS}}$ and $w_{\mathrm{Na}}$ (Eq. 6a–b). Second, the concentration of MS${}^{-}$ is taken to be much smaller than the concentration of Na${}^{+}$ in the liquid veins, i.e., $c_{\mathrm{MS}}\ll c_{\mathrm{Na}}$. This assumption is generally supported by [MS${}^{-}$] and [Na${}^{+}$] measurements on ice core samples originating from most coastal sites (Sect. 2.3; Table 1). Under the two assumptions above, relation (6a) becomes $$w_{\mathrm{MS}}=D_{\mathrm{MS}}\left(\mathrm{1},+,\mathrm{0},\left[{\displaystyle \frac{c_{\mathrm{MS}}}{c_{\mathrm{Na}}}}\right]\right)\left({\displaystyle \frac{\mathrm{1}}{c_{\mathrm{Na}}}},{\displaystyle \frac{\partial c_{\mathrm{Na}}}{\partial z}},-,{\displaystyle \frac{\mathrm{1}}{c_{\mathrm{MS}}}},{\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}}\right),$$ upon expansion of the denominator in a Taylor series. Thus, to the first order in $c_{\mathrm{MS}}/c_{\mathrm{Na}}$, the speed of MS${}^{-}$ migration can be approximated as $${\displaystyle }{w}_{\mathrm{MS}}=D_{\mathrm{MS}}\left({\displaystyle \frac{\mathrm{1}}{c_{\mathrm{Na}}}},{\displaystyle \frac{\partial c_{\mathrm{Na}}}{\partial z}},-,{\displaystyle \frac{\mathrm{1}}{c_{\mathrm{MS}}}},{\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}}\right).$$ A similar development for $w_{\mathrm{Na}}$ leads to $${\displaystyle }{w}_{\mathrm{Na}}=D_{\mathrm{Na}}\left({\displaystyle \frac{\mathrm{1}}{c_{\mathrm{MS}}}},{\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}},-,{\displaystyle \frac{\mathrm{1}}{c_{\mathrm{Na}}}},{\displaystyle \frac{\partial c_{\mathrm{Na}}}{\partial z}}\right){\displaystyle \frac{c_{\mathrm{MS}}}{c_{\mathrm{Na}}}},$$ which is also first order in $c_{\mathrm{MS}}/c_{\mathrm{Na}}$. The ratio of the migration speeds for the two ionic species is thus $${\displaystyle \frac{w_{\mathrm{Na}}}{w_{\mathrm{MS}}}}=\left({\displaystyle \frac{D_{\mathrm{Na}}}{D_{\mathrm{MS}}}}\right){\displaystyle \frac{c_{\mathrm{MS}}}{c_{\mathrm{Na}}}}.$$ If $D_{\mathrm{MS}}$ is comparable to or higher than $D_{\mathrm{Na}}$, then Na${}^{+}$ would migrate much more slowly than MS${}^{-}$. In this case, Na${}^{+}$ would be quasi-immobile and its concentration at a given depth would vary only slowly with time (compared to MS${}^{-})$. This consideration suggests the following third assumption. The concentration of Na${}^{+}$ at a given depth in the ice column and at a given time is decomposed into a mean value, ${\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}\left(z\right)$, and a fluctuation, $c_{\mathrm{Na}}^{\prime }(z,t)$: $$c_{\mathrm{Na}}\left(z,,,t\right)={\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}\left(z\right)+c_{\mathrm{Na}}^{\prime }\left(z,,,t\right).$$ Assuming that $c_{\mathrm{Na}}^{\prime }(z,t)\ll {\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}\left(z\right)$, as suggested by the relatively small mobility of Na${}^{+}$, the vertical speed of MS${}^{-}$ migration along the ice column can be further approximated as $$w_{\mathrm{MS}}=D_{\mathrm{MS}}\left({\displaystyle \frac{\mathrm{1}}{{\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}}},{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}}{\partial z}},-,{\displaystyle \frac{\mathrm{1}}{c_{\mathrm{MS}}}},{\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}}\right).$$ The insertion of Eq. (11) into Eq. (5a) yields $${\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial t}}+{\displaystyle \frac{\partial }{\partial z}}\left(w_{\ast },c_{\mathrm{MS}}\right)={\displaystyle \frac{\partial }{\partial z}}\left(D_{\mathrm{MS}},{\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}}\right),$$ where $w_{\ast }$ is an effective velocity of MS${}^{-}$ induced by vertical gradients in [Na${}^{+}$], $$w_{\ast }={\displaystyle \frac{D_{\mathrm{MS}}}{{\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}}}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}}{\partial z}}.$$ Thus, under the three stated assumptions, MSA migration can be described by a single, linear partial differential equation (Eq. 12) with the MS${}^{-}$ diffusivity as a single parameter. In this model (Eq. 12), MSA migration arises from two fundamental processes: (1) the convergence or divergence of MS${}^{-}$ driven by Na${}^{+}$ concentration gradients and (2) the diffusion of MS${}^{-}$ along its own concentration gradient. Albeit physically distinct, both processes depend on the diffusivity of MS${}^{-}$ in the intergranular liquid.

It is instructive to consider the character of the steady state distribution of [MS${}^{-}$] according to the linearized model. With the tendency term $\partial c_{\mathrm{MS}}/\partial t$ set to zero, Eq. (12) reduces to $${\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{c_{\mathrm{MS}}}{{\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}}},D_{\mathrm{MS}},{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}}{\partial z}}\right)={\displaystyle \frac{\partial }{\partial z}}\left(D_{\mathrm{MS}},{\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}}\right),$$ given the defining relation for $w_{\ast }$ (Eq. 13). If $D_{\mathrm{MS}}$ is uniform along the ice column ($\partial D_{\mathrm{MS}}/\partial z=\mathrm{0}$), and at depths where Na${}^{+}$ shows an extremum ($\partial {\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}/\partial z=\mathrm{0}$), Eq. (14) becomes $${\displaystyle \frac{c_{\mathrm{MS}}}{{\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}}{\partial z^{\mathrm{2}}}}={\displaystyle \frac{{\partial }^{\mathrm{2}}{c}_{\mathrm{MS}}}{\partial z^{\mathrm{2}}}}.$$ Since concentrations are positive quantities, the concentration ratio on the left-hand side of Eq. (15) is always positive, implying that the two second-order derivatives should always have the same sign. Thus, minima (maxima) of MS${}^{-}$ concentration will coincide with minima (maxima) of Na${}^{+}$ concentration. The [MS${}^{-}$] profile, regardless of its initial (i.e., unaltered) character, will evolve so as to become eventually in phase with the [Na${}^{+}$] profile. Figure 10 illustrates this evolution of the [MS${}^{-}$] profile to steady state as simulated by the linearized model and compares it with the evolution simulated with the RWW model (see Supplement Sect. S3 for details about the numerical solutions of these two models).

The following example illuminates the respective roles of the effective velocity, $w_{\ast }$, and of the diffusivity, $D_{\mathrm{MS}}$, in the MSA migration process. Consider a locally Gaussian profile of [Na${}^{+}$], $${\stackrel{\mathrm{‾}}{c}}_{\mathrm{Na}}\propto \exp \left(-,{\left[z,-,z_{\mathrm{0}}\right]}^{\mathrm{2}},/,\mathrm{2},{\mathit{\sigma }}^{\mathrm{2}}\right),$$ where $z_{\mathrm{0}}$ is the depth at which [Na${}^{+}$] is maximum and $\mathit{\sigma }$ describes the spread of Na${}^{+}$ on each side of the maximum (Fig. 11). In this case, the effective velocity $w_{\ast }=-D_{\mathrm{MS}-}(z-z_{\mathrm{0}})/{\mathit{\sigma }}^{\mathrm{2}}$ is positive above $z_{\mathrm{0}}$ and negative below $z_{\mathrm{0}}$, and the migration Eq. (12) becomes $${\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial t}}-{\displaystyle \frac{D_{\mathrm{MS}}}{{\mathit{\sigma }}^{\mathrm{2}}}}\left(z,-,z_{\mathrm{0}}\right){\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}}=D_{\mathrm{MS}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{c}_{\mathrm{MS}}}{\partial z^{\mathrm{2}}}}+{\displaystyle \frac{D_{\mathrm{MS}}}{{\mathit{\sigma }}^{\mathrm{2}}}}{c}_{\mathrm{MS}},$$ where it has been again assumed that $D_{\mathrm{MS}}$ is vertically uniform. Interestingly, the migration Eq. (17) has the familiar form of an advection–diffusion–reaction equation. The second term on the left-hand side corresponds to downward advection of MS${}^{-}$ above $z_{\mathrm{0}}$ and to upward advection of MS${}^{-}$ below $z_{\mathrm{0}}$, i.e., it tends to accumulate MS${}^{-}$ at $z_{\mathrm{0}}$. On the right-hand side, the first term represents Fickian diffusion, and the second term is a “reaction” term that stems from the vertical variation of $w_{\ast }$. The second term is always positive, effectively leading to MS${}^{-}$ production throughout the ice column at a rate proportional to the amount of MS${}^{-}$ initially present. As time progresses, all [MS${}^{-}$] maxima that may be present in the ice section where [Na${}^{+}$] is distributed according to Eq. (16) will be gradually shifted toward $z=z_{\mathrm{0}}$. At steady state, the [MS${}^{-}$] profile will be maintained by a balance between MS${}^{-}$ advection to the [MS${}^{-}$] maximum and effective production on the one hand, and the diffusion of MS${}^{-}$ away from the [MS${}^{-}$] maximum on the other hand (Fig. 11).

Evolution of [MS${}^{-}$] simulated by the linearized model assuming (i) an effective velocity based on a Gaussian profile of [Na${}^{+}$] (blue) and (ii) $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. The dashed, dotted-dashed, and solid red lines show, respectively, the initial, transient, and steady state profiles of [MS${}^{-}$]. The different symbols and arrows show the two transport processes affecting the steady state profile of [MS${}^{-}$]: MS${}^{-}$ convergence toward the [Na${}^{+}$] maximum and MS${}^{-}$ diffusion away from the [Na${}^{+}$] maximum (see text and Supplement for details).

[Figure omitted. See PDF]

Assessment of MS${}^{-}$ diffusivity

We next aim to constrain a range of values for MS${}^{-}$ diffusivity along ice grain boundaries consistent with observed MS${}^{-}$ concentrations in polar ice. To this end, the simplified model of MSA migration (Eq. 12), which includes $D_{\mathrm{MS}}$ as the sole parameter, is solved for different values of $D_{\mathrm{MS}}$ and model results are compared with data from the DIV2010 ice core.

The model is solved numerically using finite differences (Supplement Sect. S3; Figs. S5 and S6). The model domain, with an upper boundary situated at $z=\mathrm{9.1}$ m and a lower boundary at $z=\mathrm{60.4}$ m, is intended to represent the present-day shallow zone at the DIV2010 site. The model grid has a uniform spacing ($\mathrm{\Delta }z=\mathrm{0.05}$ m), with grid points coinciding with the sampling depths of the DIV2010 core. The grid cell interfaces coincide with the upper and lower boundaries of the domain. With this configuration of the grid, the boundary conditions of the model consist of a vanishing flux of MS${}^{-}$ prescribed at the upper and lower boundaries of the domain: $$w_{\ast }{c}_{\mathrm{MS}}-D_{\mathrm{MS}}{\displaystyle \frac{\partial c_{\mathrm{MS}}}{\partial z}}=\mathrm{0},\text{at}z=\mathrm{9.1}\text{ and }\mathrm{60.4}\mathrm{m}.$$ The initial conditions of the model consist of an idealized [MS${}^{-}$] profile obtained by linearly interpolating, at the model grid points, the (unaltered) monthly mean [MS${}^{-}$] values for the shallow zone (Fig. 5). The vertical profile of [Na${}^{+}$], which determines the effective velocity of MS${}^{-}$ ($w_{\ast })$, is directly derived from the measured profile of [Na${}^{+}$] in the shallow zone (since model grid points coincide with sampling depths, no interpolation is necessary).

Cross-correlation between [Na${}^{+}$] and [MS${}^{-}$] at different vertical spacings (a) and different time lags (b). In each panel, the red line shows the cross-correlations calculated from the DIV2010 data below the shallow zone ($z$ > 9.1 m), the blue lines show the cross-correlations calculated from the linearized model for different values of $D_{\mathrm{MS}}$, the solid black line shows the cross-correlations corresponding to the initial conditions of the model, and the horizontal dashed lines show the 95 % confidence interval for two uncorrelated time series (see text and Supplement for details).

[Figure omitted. See PDF]

The model is integrated over a time interval that approximates the time it would take for the shallow zone to be buried by a layer of equal thickness through surface accumulation (see Supplement Sect. S4 for details). This final time, denoted as $t_{\mathrm{f}}$, is taken as 95 years, based on the difference between the ages of the samples at $z=\mathrm{9.1}$ and 60.4 m. At the end of the model integration ($t=t_{\mathrm{f}})$, the cross-correlation between the [MS${}^{-}$] and [Na${}^{+}$] profiles simulated by the model between $z=\mathrm{9.1}$ and 60.4 m is calculated and compared to the cross-correlation between the measured [Na${}^{+}$] and [MS${}^{-}$] profiles over the same depth interval. This procedure is repeated for four different values of $D_{\mathrm{MS}}$ – 10${}^{-\mathrm{10}}$, 10${}^{-\mathrm{11}}$, 10${}^{-\mathrm{12}}$, and 10${}^{-\mathrm{13}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ – encompassing the values assumed or suggested in prior studies (e.g., Rempel et al., 2002; Roberts et al., 2009). A “good” value of $D_{\mathrm{MS}}$ is expected to lead to a “good” agreement between the simulated and observed cross-correlations, at least at small lags.

Consider first the model solution with $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{13}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. At $t=t_{\mathrm{f}}$, the simulated [MS${}^{-}$] profile has not significantly deviated from its initial profile, which is negatively correlated at zero lag with the [Na${}^{+}$] profile (Fig. 12a). This result suggests that the value of $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{13}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ is too small to account for the down-core change in phase relationship between [Na${}^{+}$] and [MS${}^{-}$] observed at DIV2010. For 10${}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ < $D_{\mathrm{MS}}$ < 10${}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$, the simulated cross-correlation at zero lag between [Na${}^{+}$] and [MS${}^{-}$] switches from negative to positive. For values of $D_{\mathrm{MS}}\ge {\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$, it is positive but much stronger than observed, suggesting that these values may be too large. Thus, $D_{\mathrm{MS}}$ values that best explain the DIV2010 data would be in the range from 10${}^{-\mathrm{12}}$ to 10${}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$, i.e., greater than the value of ($\mathrm{4.1}\times {\mathrm{10}}^{-\mathrm{13}}\pm \mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{14}})$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ reported by Roberts et al. (2009) and lower than the value of $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{10}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ assumed by Rempel et al. (2002). We stress that this result is immune to potential dating errors in the sense that the cross-correlation coefficients are calculated for different vertical spacings along the core, not for different time lags; calculating cross-correlations at different time lags leads to a similar result (Fig. 12b; see also Supplement Figs. S5 and S6).

While the $D_{\mathrm{MS}}$ range estimated by a comparison to DIV2010 data is instructive, we note it is not necessarily universal, as diffusivities in polar ice are expected to vary in response to multiple glaciological factors. For example, the experimental results of Kim and Yethiraj (2008) show that the diffusion coefficients of ions in undercooled mixtures are a function of both ionic concentration and temperature. Additionally, physical properties of the firn and ice, including porosity, grain-boundary density, and crystal size, may affect the partitioning of chemical impurities between the liquid premelt and the ice lattice (Dominé et al., 2008; Spaulding et al., 2011), thereby affecting the amount of impurities subjected to anomalous diffusion as well as the interconnectivity of the liquid premelt or vein network. While the RWW model can account for this partitioning (Rempel et al., 2002), the proportions of total MS${}^{-}$ and Na${}^{+}$ that are present in liquid form remain poorly constrained (Sakurai et al., 2010). Even at a given site, seasonal and interannual variations in impurity concentrations may lead to down-core changes in $D_{\mathrm{MS}}$. Finally, $D_{\mathrm{MS}}$ does not take into account whether MS${}^{-}$ migration is dominated by diffusion at two-grain boundaries, or at triple junctures and node networks (Wettlaufer and Worster, 2006; Riche et al., 2012). As a result of all these complicating factors, $D_{\mathrm{MS}}$, as defined in the RWW model and constrained here, should probably be viewed as an effective diffusivity.

Paleoclimatic implications

Revisiting the effect of temperature on grain boundary migration

In Sect. 2.3, we tested the hypothesis that post-depositional formation of winter [MS${}^{-}$] maxima occurs solely as a result of the precipitation of MS${}^{-}$ salts from their grain boundary solutions in sea-salt-rich winter layers (Mulvaney et al., 1992; Wolff, 1996; Kreutz et al., 1998; Pasteur and Mulvaney, 2000; Curran et al., 2002). This hypothesis, denoted below as the “Mulvaney model”, suggests that MS${}^{-}$ in undercooled solutions should migrate along its concentration gradient via Fickian diffusion until reaching Na${}^{+}$-rich layers, where crystallization of CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na removes MS${}^{-}$ from the premelt solution, thereby perpetuating an [MS${}^{-}$] gradient between summer and winter layers in the residual premelt. Importantly, it suggests MSA migration would be inhibited at sites where in situ temperatures are greater than the eutectic temperature of the binary system CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O ($-$29.3 ${}^{\circ }$C), since CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na in this case would not be precipitated from the premelt liquid. However, such an inhibition is not apparent in our data compilation (Sect. 2.3).

The RWW model (Sect. 4) is fundamentally different than the Mulvaney model. The RWW model does not represent crystallization and metathetic removal of constituents from the liquid phase (Sect. 2). Rather, in the RWW model, MS${}^{-}$ is implicitly assumed to remain dissolved in the premelt liquid following migration from the summer to winter layers, provided in situ temperatures exceed the eutectic temperature of the binary system CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O.

The foregoing assumption can be evaluated using our data compilation. At 19 of the 20 sites considered (Table 1), the ice temperature as estimated from the annual mean SAT exceeds the eutectic temperature of the binary system CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O, suggesting MSA migration could occur according to the RWW model at these sites. The only site where ice temperature is estimated to be less than $-$29.3 ${}^{\circ }$C is the Summit2010 record from Greenland (Maselli et al., 2017), where MSA migration is also observed but not predicted to occur based on this model. We offer a few explanations for this sole discrepancy. First, ice temperatures may depart significantly from annual mean SAT. The annual mean SAT at Summit, Greenland, is $-$29.5 ${}^{\circ }$C (Giese and Hawley, 2015), just below but very close to the eutectic temperature for the system CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O. It is therefore conceivable that ice temperatures may in fact slightly exceed the eutectic temperature of this system, at least along some portions of the ice column and (or) during some time intervals in the past. Another possibility is that the Ca${}^{\mathrm{2}+}$ salt of MS${}^{-}$, Ca(CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}{)}_{\mathrm{2}}$, which is a component of the system Ca(CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}{)}_{\mathrm{2}}$ $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O with a lower eutectic of $-$32.6 ${}^{\circ }$C, may be more efficient than the Na${}^{+}$ salt of MS${}^{-}$, CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na, in driving MSA migration at Summit2010, given the much higher abundances of Ca${}^{\mathrm{2}+}$ in inland Greenland ice compared to coastal Antarctica (Iizuka et al., 2008).

A currently poorly constrained situation arises for sites characterized by in situ temperatures less than $\sim$ $-$30 ${}^{\circ }$C and greater than $-$75 ${}^{\circ }$C (Table 2). In this temperature regime, MSA should remain in solution while Na${}^{+}$ is presumably immobile, either as solid-state NaCl, CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na, or Na${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ (Table 2). MSA migration as envisioned in the Mulvaney model, but not in the RWW model, may operate under such conditions. On the other hand, the Mulvaney model may not apply should summer concentrations of Na${}^{+}$ be high enough to sequester a large fraction of the [MS${}^{-}$] as CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na(s) in summer layers. This sequestration process appears supported by the lack of discernable MSA migration in the sub-annually resolved portion (i.e., down to $\sim$ 10.5 m) of the [MS${}^{-}$] record from the South Pole (SP-95), where annual mean SAT is $-$51 ${}^{\circ }$C (Meyerson et al., 2002). While SP-95 is not considered in our data compilation due to the site's low $\dot{b}$ (0.08 m w.e. yr${}^{-\mathrm{1}})$, the lack of clear MSA migration at SP-95 departs from the expected relationship found between $\dot{b}$ and $z_{\mathrm{fo}}$ in Antarctica (Sect. 2.1; Fig. 2). This observation leads us to speculate that MS${}^{-}$ at SP-95 may be immobilized in the summer layers through a metathesis reaction with Na${}^{+}$ allocated to the grain boundaries.

Extent of MSA migration in the DIV2010 ice core as estimated from the linearized model. (a) Top: DIV2010 [Na${}^{+}$] record (five-point smoothed). Middle: (i) annual means of [MS${}^{-}$] (black line) and standard deviations of [MS${}^{-}$], $\mathit{\sigma }$ (gray band, representing $\pm \mathrm{1}\mathit{\sigma }$ about the means) for the initial conditions (ICs) of the model experiments; and (ii) annual means of [MS-] (blue line) and standard deviations of [MS${}^{-}$] (blue band) for the terminal conditions (TCs, $\sim$ 95 years) of the model experiments for $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. Bottom: same as panel (a), middle (blue band), but for $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. (b) Probability density of the absolute difference between the annual mean [MS${}^{-}$] for the ICs and TCs, normalized to the annual mean [MS${}^{-}$] of the ICs, for all model experiments with $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ (blue) and ${\mathrm{10}}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ (red). (c) Absolute difference between the annual mean [MS${}^{-}$] for the ICs and TCs, normalized to the annual mean [MS${}^{-}$] of the ICs, for all model experiments with $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ (blue) and ${\mathrm{10}}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ (red), as a function of the data averaging interval. The shaded regions illustrate the dispersion ($\pm$1 standard deviation) of the normalized absolute differences among the model experiments.

[Figure omitted. See PDF]

In sum, we propose that MSA migration operates at most sites through anomalous diffusion as described in the RWW model. This migration process would occur if in situ temperatures exceed the relevant eutectic temperatures, so that MS${}^{-}$ can occur in solution in the presence of Na${}^{+}$ and Ca${}^{\mathrm{2}+}$ situated at the grain boundaries. The relevant eutectic temperatures are inferred here to range from $-$29.3 to $-$32.6 ${}^{\circ }$C, depending on the binary system considered (Sakurai et al., 2010). Importantly, most coastal ice core locations, despite their paleoclimatic significance, are susceptible to MSA migration given their relatively high SAT (Table 1), although high accumulation rates will mitigate this phenomenon to some degree (Sect. 2.1). One coastal region of exception may be northeastern Greenland, where relatively cold conditions (SAT of approximately $-$30 to 33 ${}^{\circ }$C; Weißbach et al., 2016) may keep the ice below the eutectic temperature of the primary MS${}^{-}$ salts.

Vertical extent of MSA migration

A reigning question in the use of [MS${}^{-}$] in polar ice as a paleoclimate proxy is the extent of MSA migration along the core. Past studies have circumvented this potential issue by assuming either that (1) MSA migration is confined within an annual layer (Kreutz et al., 1998; Curran et al., 2003; Thomas and Abram, 2016) or (2) multiyear averages of [MS${}^{-}$] are largely unaffected by migration (Wolff, 1996). While (2) appears a more conservative approach, the requisite averaging period, and thus the maximum resolution that can be achieved in a paleoclimatic reconstruction given MSA migration, remains unknown. In this section, we examine assumptions (1)–(2) using DIV2010 data and the linearized model of MSA migration.

In an effort to account for a range of initial (i.e., unperturbed) [MS${}^{-}$] profiles in the shallow zone of DIV2010 (Sect. 3.1), a large number (10 000) of numerical experiments of MSA migration are conducted. The initial [MS${}^{-}$] profile of a given experiment is obtained by adding, to the monthly mean [MS${}^{-}$] values observed in the shallow zone, a normal noise with a mean of zero and a variance equal to that of the mean monthly values in the shallow zone. If a negative concentration value arises in the initial profile, the procedure is repeated until all values in the profile are positive. Using this approach, interannual variability in the initial [MS${}^{-}$] profile is emulated, such that no 2 years should contain the same mean [MS${}^{-}$] in a given experiment, nor should a given year display the same mean [MS${}^{-}$] for different experiments. In contrast, all experiments rely on the same [Na${}^{+}$] profile measured at DIV2010. Given the uncertainties in MS${}^{-}$ diffusivity (Sect. 4.4), two sets of experiments are considered: a first set with $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ and a second with $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ (so that $\mathrm{2}\times \mathrm{10}\mathrm{000}=\mathrm{20}\mathrm{000}$ experiments are actually performed). For all experiments, the model is subjected to a condition of no MS${}^{-}$ flux at both the upper and lower boundaries of the domain and is integrated for $t_{\mathrm{f}}=\mathrm{95}$ years (Sect. 4.4).

Figure 13a shows (i) the [Na${}^{+}$] profile observed at DIV2010 and used to constrain the effective velocities $w_{\ast }$ in the model, and (ii) the annual mean [MS${}^{-}$] profiles simulated by the model at $t=t_{\mathrm{f}}$. For $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$, the changes in annual mean [MS${}^{-}$] relative to the initial [MS${}^{-}$] profile are much larger than for $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ (Fig. 13b), particularly prior to AD 1975 ($t=\mathrm{25}$ years). In some sections of the simulated profiles, dramatic positive or negative changes in annual mean [MS${}^{-}$] occur, depending on the magnitude of the local [Na${}^{+}$] gradients. In some individual years (e.g., AD 1954), relative changes in the annual mean [MS${}^{-}$] reach 60 to 100 %, clearly negating the assumption (1) above – that MSA migration is confined within an annual layer.

Given this finding, we next explore assumption (2) that multiyear averages of [MS${}^{-}$] data could be used to accurately reflect the original (i.e., unaltered) multiyear mean [MS${}^{-}$] signal. To this end, we average the simulated [MS${}^{-}$] profiles at $t=t_{\mathrm{f}}$ in different time intervals ranging from 3 to 15 years and compare the final (altered) averages to the initial (unaltered) averages (Fig. 13c). As expected, the difference between the final and initial averages of [MS${}^{-}$] decreases as the averaging period increases. Interestingly, the difference shows only modest reduction as the averaging period increases from 7 to 15 years.

In sum, while our results may pertain only to DIV2010 and rely on a series of modeling assumptions (Sect. 4.3), two points appear worthy of note. First, for [MS${}^{-}$] records showing evidence of MSA migration, the assumption that MSA remained confined within annual layers may not be generally valid, given, in particular, the high interannual variability in the concentrations of Na${}^{+}$ and other major impurities potentially conducive to MSA migration, which is typical of most ice cores originating from coastal sites (Legrand and Mayewski, 1997). High interannual variability in [Na${}^{+}$], for example, corresponds to large vertical [Na${}^{+}$] gradients along the core, which tend to enhance MSA migration. Second, at least for [MS${}^{-}$] records exhibiting severe MSA migration such as at DIV2010, averaging the data over a time period of approximately 10 years may constitute a reasonable compromise between accuracy and temporal resolution for paleoclimatic reconstruction.

Revisiting the combined influence of snow accumulation and [Na${}^{+}$] on MSA migration

In Sect. 2, we provided empirical evidence that two local factors appear to influence the shallowest depth of MSA migration in polar ice cores: annual mean accumulation rate ($\dot{b})$ and core-averaged Na${}^{+}$ concentration. Here, we assess whether the ability to predict $z_{\mathrm{fo}}$ from these two factors is also mechanistically grounded. Since the linearized model is not valid for small values of [Na${}^{+}$] $/$ [MS${}^{-}$], the original model of Rempel et al. (2002) is used in an effort to produce results of more general validity.

(a) Evolution of the phase difference between [MS${}^{-}$] and [Na${}^{+}$] maxima for (i) two different extents of the model domain, or “annual layer thickness” (0.1 and 1 m); and (ii) different layer averages of [Na${}^{+}$] (25–400 $\mathrm{\mu }$g L${}^{-\mathrm{1}}$), as calculated from the RWW model. Note the two different scales along the horizontal axis: the scale from 0.1 to 100 years applies to model results with $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$, and the scale from 1 to 1000 years applies to model results with $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. (b) Time required for approximate alignment of [MS${}^{-}$] and [Na${}^{+}$] maxima for different values of annual layer thickness and different layer averages of [Na${}^{+}$], as calculated from the RWW model. Conditions for the DIV2010 core site are indicated by the cross. (c) Shallowest depth of MSA migration for different values of annual mean accumulation rate and different core averages of [Na${}^{+}$] according to our data compilation, with DIV2010 denoted. Note the logarithmic scales in all panels.

[Figure omitted. See PDF]

The time required for approximate alignment of [MS${}^{-}$] and [Na${}^{+}$] maxima ($t_{\mathit{\phi }}$) as a function of annual layer thickness ($\mathit{\lambda }$) for different layer averages of [Na${}^{+}$] and $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$. The various curves are least-squares fits of the form $t_{\mathit{\phi }}=a{\mathit{\lambda }}^e$. The exponent $e$ is estimated to about 2 for all values of layer-averaged [Na${}^{+}$] (note that the fits for $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$, not shown, yield a value of $a$ that is a factor of 10 higher than for $D_{\mathrm{MS}}={\mathrm{10}}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$).

[Figure omitted. See PDF]

We first simulate, for a range of layer thicknesses and layer-averaged [Na${}^{+}$] values, the time it takes for an [MS${}^{-}$] maximum present in the annual layer, and initially out of phase with the [Na${}^{+}$] maximum in the layer ($\mathit{\phi }=\mathrm{180}$${}^{\circ }$), to align with the [Na${}^{+}$] maximum in the layer ($\mathit{\phi }=\mathrm{0}$${}^{\circ }$). Given the asymptotic nature of the concentration evolutions simulated by the model, we approximate this time as the time at which the phase difference between the [MS${}^{-}$] and [Na${}^{+}$] maxima drops to $\mathit{\phi }$ < 20${}^{\circ }$. Experiments are conducted for different extents of the model domain to represent different values of annual layer thicknesses ($\mathit{\lambda })$. For each experiment, the initial [MS${}^{-}$] and [Na${}^{+}$] profiles in the layer are sinusoidal functions of depth, with (i) a period set equal to the layer thickness, (ii) a [MS${}^{-}$] maximum present in the middle of the layer, and (iii) two [Na${}^{+}$] maxima present at the top and bottom of the layer (see Supplement Sect. S5 for details). The model is subjected to a condition of no flux both at the top and at the bottom of the layer. The model parameters are set to ${\mathrm{\Gamma }}_{\mathrm{MS}}={\mathrm{\Gamma }}_{\mathrm{Na}}=\mathrm{6.5}$ K M${}^{-\mathrm{1}}$ and $D_{\mathrm{MS}}=D_{\mathrm{Na}}={\mathrm{10}}^{-\mathrm{11}}$ or 10${}^{-\mathrm{12}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$.

The time required for approximate phase alignment ($\mathit{\phi }$ < 20${}^{\circ }$) of the simulated [MS${}^{-}$] and [Na${}^{+}$] profiles is shown for two different values of $\mathit{\lambda }$ and a range of layer-averaged [Na${}^{+}$] values (Fig. 14a); this time is referred to as $t_{\mathit{\phi }}$ below. It is seen that $t_{\mathit{\phi }}$ increases with $\mathit{\lambda }$ and decreases with layer-averaged [Na${}^{+}$]. Similar results are displayed in Fig. 14b in a form that is reminiscent of Figs. 5 (linear scales) and 14c (logarithmic scales), which both show the combined effect of annual mean accumulation and core-averaged [Na${}^{+}$] ([$\stackrel{\mathrm{‾}}{{\mathrm{Na}}^{+}}$]) on the shallowest depth of MSA migration ($z_{\mathrm{fo}})$ in our data compilation. To the extent that $\mathit{\lambda }$ increases with $\dot{b}$, the model results appear to be qualitatively consistent with the data, thereby providing a theoretical basis to the notion that $\dot{b}$ and [Na${}^{+}$] could be used to predict $z_{\mathrm{fo}}$.

The nonlinear relationship found between $z_{\mathrm{fo}}$ and $\dot{b}$ (Sect. 2) is also worth further exploration. Results from the RWW model can be well approximated by the power law $t_{\mathit{\phi }}\propto {\mathit{\lambda }}^{\mathrm{2}}$ (Fig. 15), which is reminiscent of the power law relationship $z_{\mathrm{fo}}\propto {\mathit{\lambda }}^{\mathrm{1.77}}$ derived from our data compilation (Fig. 2). The observed variation of the shallowest depth of MSA migration with accumulation rate would thus reflect the mere fact that, similar to Fickian diffusion with constant diffusivity, the timescale for anomalous diffusion ($t_{\mathit{\phi }})$ varies quadratically with the thickness over which the diffusion takes place ($\mathit{\lambda })$. Both current observations and the present set of experiments with the RWW model suggest that [MS${}^{-}$] records from ice cores characterized by high accumulation and low core-averaged [Na${}^{+}$] should undergo relatively small alteration by MSA migration.

Conclusions

Polar ice core records of methanesulfonic acid have been used to draw inferences about oceanic and atmospheric processes at polar latitudes on a range of timescales. However, both observation and theory suggest that MSA is mobile in the ice column, leading to uncertainties about its integrity as an indicator of past climatic conditions. Here, we synthesize existing data from a range of polar environments and consider an impurity transport model to study MSA migration in polar ice. Emphasis is placed on (i) the environmental conditions that favor MSA migration and (ii) the physico-chemical processes causing the movement of MSA in polar firn and ice.

Our analysis shows that the shallowest depth at which MSA migration occurs in coastal ice cores varies with annual mean accumulation rate. In Antarctica in particular, a power law characterizes this relationship accurately. It suggests that the absence of MSA migration observed in some ice cores from high-accumulation sites stems from the fact that chemical measurements for these cores have been conducted on samples that are not deep enough to have undergone migration. Thus, MSA migration in polar ice may be more general than commonly thought. Annual mean surface air temperature and the concentration of the dominant cation, Na${}^{+}$, appear to be less influential than accumulation rate under most circumstances, at least at most coastal Antarctic sites and in the temperature range from $-$29.5 to $-$12.5 ${}^{\circ }$C. A notable exception is for inland Greenland sites, where MSA migration tends to occur deeper in the core than would be predicted from surface accumulation alone – an offset hypothesized to stem from extremely low concentrations of marine-derived impurities relative to most coastal Antarctic sites. Our analysis further suggests that MSA migration generally takes place once firn or ice density reaches a critical value near 550 kg m${}^{-\mathrm{3}}$, which corresponds to the tightest packing of spherical ice grains in the firn and enabling the formation of premelted liquid veins at grain boundaries. However, at some low-accumulation sites ($\dot{b}=\mathrm{0.1}$–0.2 m w.e. yr${}^{-\mathrm{1}})$, MSA migration is observed at depths where bulk density is likely to be less than 550 kg m${}^{-\mathrm{3}}$. This result suggests that small-scale variability in ice density is important and (or) that other factors may also determine the onset of MSA migration along the firn or ice column.

New high-resolution data from the West Antarctic DIV2010 ice core show annual [MS${}^{-}$] maxima gradually shifting down-core from austral summer, when MSA deposition is high, to austral winter, when MSA deposition is low and Na${}^{+}$ deposition is high. As a result, a down-core change in the phase relationship between [MS${}^{-}$] and [Na${}^{+}$] is observed, whereby [MS${}^{-}$] and [Na${}^{+}$] are negatively correlated at zero lag in the upper part of the core and positively correlated at zero lag in the lower part of the core, providing evidence of the progressive nature of MSA migration.

A linearized version of the impurity transport model of Rempel et al. (2002) is derived for the binary system CH${}_{\mathrm{3}}$SO${}_{\mathrm{3}}$Na $\cdot$ $n$H${}_{\mathrm{2}}$O–H${}_{\mathrm{2}}$O in order to further understanding of the MSA migration phenomenon in polar ice. In this linearized model, MS${}^{-}$ transport is governed by a single linear partial differential equation with MS${}^{-}$ diffusivity in the undercooled liquid ($D_{\mathrm{MS}})$ as the sole parameter. In this model, MSA migration arises from two transport processes: (1) the convergence or divergence of MS${}^{-}$ driven by [Na${}^{+}$] gradients and (2) the diffusion of MS${}^{-}$ along its own concentration gradient. Analysis of this model shows that [MS${}^{-}$] maxima (minima) are bound to coincide with [Na${}^{+}$] maxima (minima) along the ice column, regardless of the timing of MSA deposition maxima. The model, therefore, provides a mechanistic explanation for the tendency for MSA, deposited mainly during summer, to present concentration peaks in winter layers in the deepest part of polar ice cores.

Finally, we use the linearized MSA migration model and the DIV2010 data to gain insight into two poorly constrained yet critically important aspects of MSA migration. First, we evaluate different values of MS${}^{-}$ diffusivities in polar ice. We find that $D_{\mathrm{MS}}$ values in the range from 10${}^{-\mathrm{12}}$ to 10${}^{-\mathrm{11}}$ m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$ lead to the most accurate simulations of the down-core change in the phase relationship between [MS${}^{-}$] and [Na${}^{+}$] observed at DIV2010. Second, using this range of values, we apply the model to determine the extent to which MSA migration has altered the original [MS${}^{-}$] record for DIV2010. We estimate the errors incurred by averaging [MS${}^{-}$] data over annual (and multiyear) intervals, an approach often adopted to reduce the effect of migration on the interpretation of [MS${}^{-}$] records. We find that MSA migration may have lead to significant changes in the annual and multiyear [MS${}^{-}$] averages at DIV2010. This result suggests that [MS${}^{-}$] records severely perturbed by MSA migration may best be used to infer decadal or lower-frequency climate variability, though a range of [MS${}^{-}$] records and a better constrained model are needed to investigate this further.

The migration of MSA in cold, polar ice is a fascinating but challenging phenomenon. This paper covers many, but not all, of its observational and theoretical aspects. For example, contentions of MSA migrating away from regions of high acidity in the core, as caused by the deposition of compounds of volcanic origin (Curran et al., 2002; Delmas et al., 2003), are not explored here. While the model of Rempel et al. (2002) provides an important mechanistic framework for understanding MSA migration in polar ice and perhaps for ultimately correcting its effects for paleoclimatic reconstruction, its usefulness remains limited by uncertainties about key physico-chemical parameters. These include most notably the diffusivities of the relevant migrating species in undercooled liquid, the slope of the liquidus curves for relevant, interacting species, and the partitioning of impurities between the ice lattice and the surface of the ice grains. Laboratory studies under a range of controlled conditions would help constrain these parameters, improve our understanding of MSA migration in polar ice, and make full use of the paleoclimatic potential of this compound.

The DIV2010 chemistry data used in the generation of this paper are available on request to M. Osman.

The Supplement related to this article is available online at https://doi.org/10.5194/tc-11-2439-2017-supplement.

MO, SBD, and OM designed the study. MO carried out the data synthesis. OM developed the linearized model and MO coded the RWW model and the linearized model. SBD collected the DIV2010 ice core. MJE measured the DIV2010 ice core ion data. MO, SBD, and OM interpreted the results. MO wrote the paper. All co-authors commented on the paper.

The authors declare that they have no conflict of interest.

Acknowledgements

We thank Elizabeth Thomas for contributing the Ferrigno data, Joseph McConnell for contributing the Summit2010 and D4 data, and Alison Criscitiello and Weifu Guo for early discussions about MSA migration in the DIV2010 record. Matthew Osman acknowledges government support awarded by DoD, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a. We thank Alan Rempel and one anonymous reviewer, whose comments and suggestions greatly improved the content of this paper. This work was supported by the US NSF (ANT-0632031 and PLR-1205196 to Sarah B. Das, and NSF-MRI-1126217 to Matthew J. Evans) and a Woods Hole Oceanographic Institution Interdisciplinary Research award to Sarah B. Das and Olivier Marchal. Edited by: Becky Alexander Reviewed by: Alan Rempel and one anonymous referee