Pervasive diffusion of climate signals recorded in ice-vein ionic impurities

A theory of vein impurity transport conceived two decades ago predicts that signals in the bulk concentration of soluble ions in ice migrate under a temperature gradient. If valid, it would mean that some palaeoclimatic signals deep in ice cores (signals from vein impurities as opposed to matrix or grain-boundary impurities) suffer displacements that upset their dating and alignment with other proxies. We revisit the vein physical interactions to find that a strong diffusion acts on such signals. It arises because the Gibbs–Thomson effect, which the original theory neglected, perturbs the impurity concentration of the vein water wherever the bulk impurity concentration carries a signal. Thus, any migrating vein signals will not survive into deep ice where their displacement matters, and the palaeoclimatic concern posed by the original theory no longer stands. Simulations with signal peaks introduced in shallow ice at the GRIP and EPICA Dome C ice-core sites, ignoring spatial fluctuations of the ice grain size, confirm that rapid damping and broadening eradicates the peaks by two-thirds way down the ice column. Artificially reducing the solute diffusivity in water (to mimic partially connected veins) by 10³ times or more is necessary for signals to penetrate into the lowest several hundred metres with minimal amplitude loss. Simulations incorporating grain-size fluctuations on the decimetre scale show that these can cause the formation of new, non-migrating solute peaks. The deep solute peaks observed in ice cores can only be explained by widespread vein disconnection or a dominance of matrix or grain-boundary impurities at depth (including their recent transfer to veins) or signal formation induced by grain-size fluctuations; in all cases, the deep peaks would not have displaced far. Disentangling the different signal contributions – from veins, the ice matrix, grain boundaries, and grain-size fluctuations – will aid robust reconstruction from ion records.

1 Introduction

Chemical impurity concentrations in ice cores yield diverse palaeoclimatic information (e.g. Legrand and Mayewski, 1997; Wolff et al., 2006). As with other ice-core proxies, such as stable water isotopes, it is generally hoped that post-depositional modification of signals in such records, which may hamper their interpretation, is minimal. For signals of abrupt or discrete (e.g. volcanic) events, of interest is whether their position and shape record their timing and magnitude faithfully. If, for instance, signals diffuse in the ice, neighbouring peaks may merge as they descend towards the ice-sheet base. Conceivably, a range of physico-chemical processes may distort signals, limiting the resolution and accuracy of the retrievable information.

Figure 1 High-resolution records of dissolved ionic impurities in ice cores from Antarctica and Greenland, illustrating the occurrence and expressions of high peaks and signals at decimetre or shorter scales in deep ice. Vertical axes plot the total impurity concentration in bulk ice (this differs from our model variable c_B, which refers only to impurities in the water-vein network); 1 ng g⁻¹ is equivalent to 1 ppb or 1 µg L⁻¹. (a) The sulfate record from the EPICA Dome C core, Antarctica, measured by fast ion chromatography (Traversi et al., 2009); the arrow locates panel (b). (b) Zoomed portion of (a) showing two sulfate peaks. (c) A piece of the sodium record from the NEEM ice core, Greenland (Schüpbach et al., 2018), from continuous-flow analysis (CFA) measurements, showing large fluctuations on a stretch with high base level in 2187–2190 m as well as isolated peaks. (d) A piece of the sulfate record from the NGRIP ice core, Greenland (Svensson et al., 2013), measured by CFA, showing successive sub-decimetre-scale spikes. Gaps in (c) and (d) reflect missing data.

In a landmark paper, Rempel et al. (2001) proposed a theory to show that signals in the bulk concentration of dissolved ionic impurities – major ions such as SO${}_{\mathrm{4}}^{\mathrm{2}-}$, Cl⁻, and Na⁺ – in the water veins at polycrystalline grain junctions (Nye, 1989; Mader, 1992a, b) migrate relative to the ice when a temperature gradient is present. Driven by what they term “anomalous diffusion”, migration occurs in the direction of rising temperature and displaces signals with minimal distortion, so their apparent age deviates from their true age. Rempel et al. (2001) calculated cumulative signal displacements of ∼ 0.1–1 m in the lowest kilometre of the GRIP ice core. Their mechanism could decouple the ion records from other ice-core proxies and cause significant age errors in palaeoclimatic histories, especially in deep ice where temperature increases markedly towards the bed. This has prompted evaluation for signs of anomalous diffusion in some ice-core records (e.g. Tison et al., 2015) and simulation of the migration of specific species (e.g. methane sulfonic acid in firn; Osman et al., 2017). Note that it is non-trivial to infer the absence or amount of migration by comparing the records of signal peaks from different cores, due to uncertainty in depth–age scales, a lack of rigorous independent control of where individual peaks should lie, and spatial inhomogeneity in the atmospheric dispersal and deposition of impurities by environmental events. Some studies contend that signal migration may be limited – and thus the ion records dependable, because solute transport is hindered by disconnections in the vein network (Barnes et al., 2003; Barnes and Wolff, 2004) or because most ionic impurities are located at grain boundaries (Barnes and Wolff, 2004) or in salt micro-inclusions within ice crystals (Ohno et al., 2005), instead of in veins. How extensively ice-core chemical records have been altered by anomalous diffusion is unresolved, despite the relevance of this question to the synchronisation of ice-core age scales (e.g. Severi et al., 2007; Fujita et al., 2015) and dating of palaeoclimatic events.

Here we re-examine “Rempel's theory” of vein impurity transport (i.e. Rempel et al., 2001), discovering missing elements in it that change the predicted signal evolution. We find that the relevant processes cause diffusion of signals at rates that threaten their survival in ice cores, whether or not a temperature gradient drives their migration. Hence the theory's implications are radically revised. Our purpose is not to develop the theory of migrating signals to show how it can match observed ion records or be used in reconstructions, but to highlight the contrary: it struggles to explain the presence of distinct ionic peaks in deep ice. Figure 1 exemplifies such peaks in ice cores from Antarctica and Greenland. Much of our analytical and numerical work (Sects. 2 and 3) is spent on understanding the origin of the strong diffusion and showing its operation. We also correct the signal migration speed and show that grain-size fluctuations can cause new impurity signals to form. At the end (Sect. 4) we discuss what the revised theory means for the provenance of the deep ionic peaks in ice-core records.

Following Rempel's framework, we model processes below the firn–ice transition and one ion species only (SO${}_{\mathrm{4}}^{\mathrm{2}-}$ is used in our calculations). We let c_B denote its bulk concentration, using the unit mol L⁻¹ or M to refer to the amount of impurity in a unit volume of ice. Note that c_B accounts only for impurities dissolved in the water veins and excludes impurities at grain boundaries and in the ice matrix (inside grains), which are not modelled in this paper. Chemical alteration of signals via cation–anion associations (e.g. Iizuka et al., 2004; Traversi et al., 2009), reaction of vein impurities with dust (Barnes and Wolff, 2004), the segregation of impurities to locations outside the veins and thermodynamic coupling between multiple ion species (Rempel et al., 2002; Rempel and Wettlaufer, 2003) are ignored.

Figure 2 Interactions that cause the signal migration mechanism of Rempel et al. (2001) and our signal diffusion mechanism, and variables studied in this paper. (a) Conceptualisation of the vertical profiles of vein impurity concentration c and ice temperature T in an ice core; z is depth below the surface. The gradients shown here typically occur over the kilometre length scale and are most pronounced towards the base of an ice sheet; see Fig. 4 for real examples of T(z). Panels (b), (c), and (d) expand on the variations around a c_B signal, which occur at a much shorter scale of decimetres or centimetres. (b) Vein impurity concentration c. (c) Porosity ϕ. (d) Bulk impurity concentration c_B. (e) Vein cross-sectional geometry, showing the definition of the radius of curvature r_v. (f) Equilibrium vein conditions as the solution of model equations. An increase in c_B lifts the solution A to A', causing a perturbation on c (see panel b) that drives c_B diffusion. Neglect of the Gibbs–Thomson effect in the theory of Rempel et al. (2001) causes the bold curve in (f) to collapse onto the vertical dashed line; then no perturbation arises to diffuse the signal.

Before plunging into mathematics, we first outline our main finding; see Fig. 2a–d, where z denotes depth. In Rempel's theory, a centimetre- or decimetre-scale peak in the bulk solute concentration c_B, representing a climate signal, is mirrored by variations in the porosity ϕ (Fig. 2c, d), which represents the volume fraction of veins in the ice. With c_B encapsulating vein impurities, the relation c_B=cϕ holds¹, where c, the solute concentration of the vein water, is determined by temperature T through the liquidus. A background temperature gradient in the ice sets up a gentle gradient in c, driving downward solute diffusion through the vein network (Fig. 2a). Interestingly, the porosity modulates the solute diffusion flux so that on the trailing (upper) edge of the signal peak $\text{d}\mathit{\varphi }/\text{d}z>\mathrm{0}$ increases this flux with distance to draw down c_B, whereas on the leading (lower) edge of the peak, $\text{d}\mathit{\varphi }/\text{d}z<\mathrm{0}$ reduces the flux with distance to bump up c_B. Thus the peak signal in c_B translates (the same translation applies to trough signals). This mechanism is “anomalous” because solute diffusion causes the signal to move without signal diffusion, i.e. without changing its shape. In their calculations, Rempel et al. (2001) neglected the small effect of the vein-face curvature on the melting point (the Gibbs–Thomson effect), and a detailed justification of why this interfacial effect is small is given in their companion papers (Rempel et al., 2002; Rempel and Wettlaufer, 2003). Then c depends on T only, not on c_B. However, we find that when this seemingly reasonable approximation is not made, the c_B signal causes a perturbation on c(z) that drives non-negligible solute diffusion away from the peak (Fig. 2b), owing to the short length scale of the signal. Consequently, the c_B peak experiences pronounced broadening and amplitude reduction. In an ice core, these conspire with vertical compression to regulate the evolving peak shape.

Besides Rempel et al. (2001), Barnes et al. (2003) have modelled vein-mediated evolution of dissolved ion signals below the firn–ice transition. Their main interest was to explain the signal diffusion found in the top 350 m of the EPICA Dome C ice core (EPICA community members, 2004), which they inferred from observed trends of peak broadening and damping in the sulfate and chloride records with age. To explain the signal diffusivities estimated for these ions – respectively $\mathrm{4.7}\times {\mathrm{10}}^{-\mathrm{8}}$ and $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{7}}$ m² yr⁻¹, they conceived two models of vein solute transport driven by grain growth, motivated by the fact that the mean grain size in their stretch of the core increases with depth. One model invokes local gradients in c induced by porosity change during spatially non-uniform normal grain growth; this mechanism requires the presence of grain-size variations at the length scale of the c_B signals. The other, an analogue model, uses a cell-based simulation to demonstrate how the disconnection of veins by grain growth (modelled as random removal of cells and relocation of their impurity to neighbouring cells) causes diffusion even when the vein network is only partially connected. Both models predict no diffusion of c_B unless there is net grain growth.

We do not incorporate the Barnes et al. (2003) mechanisms into our model here, as the “Gibbs–Thomson diffusion” is much faster than their diffusion (by an order of magnitude at least), and adding the latter strengthens our conclusions. Our mechanism is independent of grain growth and occurs in the absence of grain-size variations. We do, however, account for the continual motion of veins during grain boundary migration, which causes a slow “residual diffusion” on c_B. This mechanism does not require a temperature gradient or grain-size variations to operate and is unaffected by vein disconnection.

Table 1 Constants used in this study.

^a One-third of the molecular diffusivity in water; see Rempel et al. (2001). ^b Value derived from Table 3.1 of Cuffey and Paterson (2010), after multiplying by $\mathrm{6}/\mathit{\pi }$ to correct for sectioning and stereological effects (Ng and Jacka, 2014).

2 The model

2.1 Key relationships

We treat polycrystalline ice in a continuum description with variables as functions of position x= (x, y, z) and time t. Ice with the mean grain size d_g (grain diameter) has the vein length density

$$\begin{array}{}\text{(1)}& l\approx {\displaystyle \frac{\mathrm{3}}{d_{\text{g}}^{\mathrm{2}}}};\end{array}$$

we adopt equality in this expression herein. If the vein faces have the radius of curvature r_v (Fig. 2e; Nye, 1989), then the ice porosity ϕ is

$$\begin{array}{}\text{(2)}& \mathit{\varphi }=\mathit{\alpha }{r}_{\text{v}}^{\mathrm{2}}l,\end{array}$$

where α=0.0725 is a geometrical factor (Nye, 1991; Mader, 1992a, b). Recalling the relation between the bulk solute concentration and the solute concentration of the vein water,

$$\begin{array}{}\text{(3)}& c_{\mathrm{B}}=\mathit{\varphi }c,\end{array}$$

and using Eqs. (1) and (2), leads to

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{d_{\text{g}}^{\mathrm{2}}{c}_{\mathrm{B}}}{\mathrm{3}\mathit{\alpha }}}=c{r}_{\text{v}}^{\mathrm{2}}.\end{array}$$

For ice at a temperature T below the reference temperature T₀, thermodynamic equilibrium between ice and vein water means that the liquidus relation is satisfied:

$$\begin{array}{}\text{(5)}& T_{\mathrm{0}}-T={\mathit{\theta }}_c+{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{{\mathit{\rho }}_{\text{i}}L{r}_{\text{v}}}}+{\mathit{\theta }}_p.\end{array}$$

This equation features in earlier studies of the subject (Nye, 1991; Mader, 1992b; Rempel et al., 2001, 2002; Rempel and Wettlaufer, 2003; Barnes et al., 2003; Dani et al., 2012). The three terms on its right-hand side describe the temperature depressions due to (i) solute, (ii) interfacial curvature (the Gibbs–Thomson effect), and (iii) pressure, respectively; ρ_i is ice density, L is latent heat, and γ is the interfacial energy. Table 1 lists the constants used in this paper, and Table 2 lists our model variables. For the first term in Eq. (5), a linear approximation θ_c=Γc is valid at temperatures not far below the melting point, e.g. Γ=4.53 K M⁻¹ for SO${}_{\mathrm{4}}^{\mathrm{2}-}$. The third term is typically small and may be absorbed into T₀ by accounting for glaciostatic overburden. Accordingly, the version of Eq. (5) that we use for the present analysis is

$$\begin{array}{}\text{(6)}& T_{\mathrm{0}}-T=\mathrm{\Gamma }c+{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{{\mathit{\rho }}_{\text{i}}L{r}_{\text{v}}}},\end{array}$$

but our simulations in Sect. 3 for specific ice-core sites will use Eq. (5), together with detailed nonlinear empirical formulas for θ_c and θ_p (Appendix A).

Table 2 Variables in our mathematical model and their units.

Given d_g, T, and c_B, Eqs. (4) and (6) can be solved for r_v and c (e.g. Barnes et al., 2003). Figure 2f illustrates the solution as the intersection point of two curves. In their theory, Rempel et al. (2001) employed the liquidus relation without including the Gibbs–Thomson term, and Rempel et al. (2002) and Rempel and Wettlaufer (2003) argued that under glaciological conditions, a large r_v makes the Gibbs–Thomson term negligible so that $\mathrm{\Gamma }c\approx T_{\mathrm{0}}-T$; then c(z) is dictated by the ice temperature profile (Fig. 2a). This approximation amounts to taking the intersection point in Fig. 2f to lie on the dashed line. However, as shown by Fig. 2f, the exact solution for c does depend on c_B – albeit weakly – for fixed T and d_g. This dependence lies at the heart of our “Gibbs–Thomson diffusion”. Specifically, when $\mathit{\gamma }{T}_{\mathrm{0}}/{\mathit{\rho }}_{\text{i}}L{r}_{\text{v}}\ll \mathrm{\Gamma }c$, a first-order approximate solution of Eqs. (4) and (6) is

$$\begin{array}{}\text{(7)}& \mathrm{\Gamma }c\approx (T_{\mathrm{0}}-T)-{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{{\mathit{\rho }}_{\text{i}}L{r}_{\text{v}}}},\end{array}$$

where $r_{\text{v}}\approx \sqrt{\frac{c_{\mathrm{B}}{d}_{\text{g}}^{\mathrm{2}}\mathrm{\Gamma }}{\mathrm{3}\mathit{\alpha }(T_{\mathrm{0}}-T)}}$, or

$$\begin{array}{}\text{(8)}& c\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{\Gamma }}}\left[(T_{\mathrm{0}}-T)-{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{{\mathit{\rho }}_{\text{i}}L}}\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\alpha }(T_{\mathrm{0}}-T)}{d_{\text{g}}^{\mathrm{2}}\mathrm{\Gamma }}}}{c}_{\mathrm{B}}^{-\mathrm{1}/\mathrm{2}}\right].\end{array}$$

This differs from Rempel's approximation by the γ (Gibbs–Thomson) term, which causes $\partial c/\partial c_{\text{B}}>\mathrm{0}$ and a perturbation to appear on c(z) where c_B(z) exhibits a signal (Fig. 2b, d). The rate of signal diffusion stemming from this minute perturbation will be quantified later. The result in Eq. (8) also shows that fluctuations in d_g and T could perturb c(z) to drive solute transport and influence c_B. We will explore the potential effect of grain-size fluctuations (Sect. 3.4), but not temperature fluctuations, as these decay quickly in ice.

If, instead of Eq. (6), the full liquidus relation Eq. (5) is analysed, with θ_c=f(c) being a mildly nonlinear function, then the above findings are qualitatively unchanged, and Eq. (8) would read $c\approx f^{-\mathrm{1}}$([ ]) with [ ] as given in Eq. (8) except that $f^{-\mathrm{1}}(T_{\mathrm{0}}-T)$ replaces ($T_{\mathrm{0}}-T)/\mathrm{\Gamma }$ in the square root.

2.2 Porosity, water, and solute conservation

The rest of the model is now formulated, for ice deforming with velocity u= (u, v, w). As the porosity is very small, $\mathit{\varphi }\sim {\mathrm{10}}^{-\mathrm{6}}\ll \mathrm{1}$, the incompressibility condition $\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0}$ holds. The total transport flux of porosity is uϕ−κ∇ϕ. The first term here describes advection by ice. The second (Fickian) term describes a net transport of porosity due to the random, unceasing vein motion that accompanies grain boundary migration. We detail its physical derivation in Appendix B. The diffusivity κ is given by

$$\begin{array}{}\text{(9)}& \mathit{\kappa }={\displaystyle \frac{K\left(T\right)}{\mathrm{3}{c}_{\mathrm{1}}}},\end{array}$$

where $K=K_{\mathrm{0}}\exp (-Q/RT)$ is the temperature-dependent grain growth rate, R is the gas constant, Q is activation energy, and c₁≈2 to 3 (we use 2.5 in our simulations).

Accordingly, porosity conservation is described by

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+\mathrm{\nabla }\cdot (\mathit{u}\mathit{\varphi }-\mathit{\kappa }\mathrm{\nabla }\mathit{\varphi })={\displaystyle \frac{m}{{\mathit{\rho }}_{\text{i}}}},\end{array}$$

where the melt rate m accounts for phase change occurring at the interfacial boundaries of the veins; m is the rate of mass melted per unit volume of ice. This result can also be derived from ice mass conservation.

If q is the water flux percolating through the vein network, then water mass conservation requires

$$\begin{array}{}\text{(11)}& {\mathit{\rho }}_{\text{w}}\left[{\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+\mathrm{\nabla }.(\mathit{q}+\mathit{u}\mathit{\varphi }-\mathit{\kappa }\mathrm{\nabla }\mathit{\varphi })\right]=m,\end{array}$$

and we deduce from this together with Eq. (10) that

$$\begin{array}{}\text{(12)}& \mathrm{\nabla }\cdot \mathit{q}=-(\mathrm{1}-r){\displaystyle \frac{m}{{\mathit{\rho }}_{\text{i}}}},\end{array}$$

where $r={\mathit{\rho }}_{\text{i}}/{\mathit{\rho }}_{\text{w}}\approx \mathrm{0.92}$ is the ratio of ice density to water density.

The local solute transport flux is $\mathit{u}{c}_{\text{B}}+\mathit{q}c-D\mathit{\varphi }\mathrm{\nabla }c-\mathit{\kappa }\mathrm{\nabla }{c}_{\text{B}}$. This describes the summed effects of advection by ice and water flow, molecular diffusion in vein water, and random vein motion (the same as that driving porosity diffusion; Appendix B); D is the diffusion coefficient for the ionic impurity in water, and the diffusivity κ was defined in Eq. (9). The corresponding solute conservation equation is

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+\mathrm{\nabla }\cdot (\mathit{u}{c}_{\mathrm{B}}+\mathit{q}c)=\mathrm{\nabla }\cdot (\mathit{\varphi }D\mathrm{\nabla }c+\mathit{\kappa }\mathrm{\nabla }{c}_{\mathrm{B}}).\end{array}$$

Following Rempel's theory, we define the “anomalous velocity”

$$\begin{array}{}\text{(14)}& {\mathit{u}}_c=-D{\displaystyle \frac{\mathrm{\nabla }c}{c}}\end{array}$$

and use c_B=ϕc to rewrite the D contribution in Eq. (13) to derive

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+\mathrm{\nabla }\cdot \left[\right(\mathit{u}+{\mathit{u}}_c)c_{\mathrm{B}}+\mathit{q}c]=\mathrm{\nabla }\cdot \left(\mathit{\kappa }\mathrm{\nabla }{c}_{\mathrm{B}}\right).\end{array}$$

If both water flux q and the κ-term (“residual diffusion” due to vein motion) are ignored, this simplifies, after using the incompressibility condition $\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0}$, to the Rempel et al. (2001) equation:

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+(\mathit{u}+{\mathit{u}}_c)\cdot \mathrm{\nabla }{c}_{\mathrm{B}}=-(\mathrm{\nabla }\cdot {\mathit{u}}_c)c_{\mathrm{B}}.\end{array}$$

As is well rehearsed in their theory, the advection term here predicts migration of chemical signals at the velocity u_c relative to the ice, with u_c controlled by the temperature profile via the liquidus and Eq. (14). Our analyses are to reveal departures from these predictions.

2.3 One-dimensional model

For an ice core beneath a flow divide or summit, where the ice motion is downward ($\mathit{u}=(\mathrm{0},\mathrm{0},w)$) and horizontal variations in ϕ and c_B are negligible, Eqs. (10) and (15) become (on using $\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0}$ again)

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+w{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}={\displaystyle \frac{m}{{\mathit{\rho }}_{\text{i}}}}+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\kappa }{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right),\end{array}$$

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+(w+w_c){\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(qc\right)={\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\kappa }{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z}}\right)-c_{\mathrm{B}}{\displaystyle \frac{\partial w_c}{\partial z}},\end{array}$$

where $w_c=-(D/c)\partial c/\partial z$ is the anomalous velocity downward. Equation (18) tracks the evolution of the bulk impurity profile c_B(z,t) and its signals, but the water flux q needs to be found via

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\partial q}{\partial z}}=-(\mathrm{1}-r){\displaystyle \frac{m}{{\mathit{\rho }}_{\text{i}}}},\end{array}$$

with m calculated from Eq. (17). This problem is supplemented by Eqs. (4) and (5) (or 6), which give the instantaneous distributions of c and ϕ from the depth profiles of c_B, T, and d_g. There are five equations (Eqs. 4, 5, 17–19) for the five unknowns c_B, c, ϕ, q, and m.

2.4 Correction to the signal migration velocity

Here we explain the first departure from Rempel's original predictions: the migration speed will be faster than w_c, because of water percolation induced by the migration. To see this, we need to put all advective parts of Eq. (18) (notably $\partial \left(qc\right)/\partial z$) in terms of c_B. Let c(z) be time-invariant as in Rempel's theory, and ignore the κ-term (vein motion), to write Eq. (18) as

$$c\left({\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+{\displaystyle \frac{\partial q}{\partial z}}\right)+(w+w_c){\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z}}=-c_{\mathrm{B}}{\displaystyle \frac{\partial w_c}{\partial z}}-q{\displaystyle \frac{\partial c}{\partial z}}.$$

Substituting for $\partial q/\partial z$ from Eq. (19), using Eq. (17) for m, converts this to

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& r\left({\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+w{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z}}\right)+w_c{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z}}=\\ & -c_{\mathrm{B}}{\displaystyle \frac{\partial w_c}{\partial z}}-q{\displaystyle \frac{\partial c}{\partial z}}-(\mathrm{1}-r)w\mathit{\varphi }{\displaystyle \frac{\partial c}{\partial z}}.\end{array}\end{array}$$

The z derivatives on the left-hand side describe signal advection by ice flow and anomalous diffusion, respectively. Following the scaling argument in the Supplement of Rempel et al. (2001), all terms on the right, based on the background gradient in c (Fig. 2a), are negligible on the length scale of climatic signals (≲10⁻¹ m). In Eq. (20), the prefactor r of the time derivative means that the migration speed is $w_c/r$, which exceeds w_c by ≈ 9 %.

What causes this correction? As introduced in Sect. 1, a migrating c_B peak is mirrored by a moving variation in porosity. This evolution implies freezing and melting on the peak's trailing edge and leading edge (Fig. 2c), causing water to be expelled from and absorbed into these regions, respectively, owing to the density change during phase change (see Eq. 19). Hence a water flow localised about the peak transports solute in the same direction to speed up its migration. Although modest in size, the correction applies at all depths and also in the three-dimensional model.

2.5 Diffusion of impurity signals

The second (more crucial) departure from Rempel's predictions is that signals in c_B will diffuse, as we emphasised at the outset. In Rempel's theory, c_B signals suffer no distortion as they migrate under anomalous diffusion except for a slight amplitude change due to the right-hand side of Eq. (16). But as anticipated in Sect. 2.1, they will perturb c(z) to drive solute diffusion (Fig. 2). The consequence can be studied via a stripped-down version of Eq. (15) (or 16) where ice flow, water flow, and vein motion (residual diffusion) are ignored:

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+\mathrm{\nabla }\cdot \left({\mathit{u}}_c{c}_{\mathrm{B}}\right)=\mathrm{0}.\end{array}$$

For the anomalous velocity, substituting for c from Eq. (8) into Eq. (14) yields, after some algebra,

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\mathit{u}}_c=& -D{\displaystyle \frac{\mathrm{\nabla }c}{c}}\\ \approx & {\displaystyle \frac{D}{(T_{\mathrm{0}}-T)}}\left[\mathrm{\nabla }T-{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{\mathrm{2}{\mathit{\rho }}_{\text{i}}L}}\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\alpha }(T_{\mathrm{0}}-T)}{d_{\text{g}}^{\mathrm{2}}\mathrm{\Gamma }}}}{c}_{\mathrm{B}}^{-\mathrm{3}/\mathrm{2}}\mathrm{\nabla }{c}_{\mathrm{B}}\right].\end{array}\end{array}$$

Note that this result assumes that d_g varies slowly. If grain size fluctuates strongly over short distances (e.g. decimetres), the evaluation of ∇c would introduce another term involving ∇d_g (Sect. 3.4). With the present result in Eq. (22), Eq. (21) becomes

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& {\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+\mathrm{\nabla }\cdot \left({\displaystyle \frac{D\mathrm{\nabla }T}{T_{\mathrm{0}}-T}}{c}_{\mathrm{B}}\right)=\\ & D{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{\mathrm{2}{\mathit{\rho }}_{\text{i}}L}}\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\alpha }}{d_{\text{g}}^{\mathrm{2}}\mathrm{\Gamma }}}}\mathrm{\nabla }\cdot \left({\displaystyle \frac{c_{\text{B}}^{-\mathrm{1}/\mathrm{2}}\mathrm{\nabla }{c}_{\text{B}}}{\sqrt{T_{\mathrm{0}}-T}}}\right).\end{array}\end{array}$$

This partial differential equation implies not only signal migration driven by ∇T, but also signal diffusion, due to the right-hand (Gibbs–Thomson) term, which encapsulates the c perturbation. The diffusion is nonlinear and independent of ∇T. Both the diffusion and advection terms in Eq. (23) are controlled by D because they originate from the molecular diffusion of solute in vein water. Setting γ→0 (which forces $\mathrm{\Gamma }c\equiv T_{\mathrm{0}}-T)$ recovers Rempel's theory.

How fast do signals diffuse? One way to gauge their rate of broadening and damping versus migration is by taking the magnitude ratio of the two terms in the square bracket in Eq. (22):

$$\begin{array}{}\text{(24)}& \mathit{\chi }={\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{\mathrm{2}{\mathit{\rho }}_{\text{i}}L}}\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\alpha }(T_{\mathrm{0}}-T)}{d_{\text{g}}^{\mathrm{2}}\mathrm{\Gamma }}}}{c}_{\mathrm{B}}^{-\mathrm{3}/\mathrm{2}}{\displaystyle \frac{\left|\mathrm{\nabla }{c}_{\mathrm{B}}\right|}{\left|\mathrm{\nabla }T\right|}}.\end{array}$$

The dimensionless number χ, which also sizes diffusion against advection in Eq. (23), resembles the inverse of the Péclet number in fluid mechanics. But the context here is unique as our model concerns the transport of bulk concentration signals, and the underlying physics involve the geometry of the vein network and two-phase (solute–ice) thermodynamic equilibrium.

As an example, taking $T_{\mathrm{0}}-T\sim \mathrm{25}$ K, |∇T|∼20 K km⁻¹ = 0.02 K m⁻¹, d_g∼5 mm, c_B∼1 µM, and |∇c_B|∼1 µM/10 cm = 10⁻⁵ M m⁻¹ (a decimetre-scale signal), which roughly approximate conditions in the lower third of the GRIP and EPICA ice cores, gives χ∼1.7. The effective diffusivity in the right-hand term of Eq. (23) is then $\mathrm{2.1}\times {\mathrm{10}}^{-\mathrm{6}}$ m² yr⁻¹, which is much higher than the diffusivities reported by Barnes et al. (2003) (Sect. 1). Although these ballpark estimates depend on the chosen values (including the signal magnitude) and will change if we consider elsewhere in the ice column, they show that diffusion can pervasively modify signals, whose short length scales play a key role in amplifying the perturbations on c to cause diffusion at a rate rivalling migration. We confirm this by numerical simulation next.

Figure 3 (a) Submergence velocity w(z) in an ice core and the corresponding (b, c) z–t and z^′–t plots of material trajectories. The bold line in (b) is the depth–age relation. In (c), z^′ is a depth coordinate that measures distance from ice located on the depth–age relation in (b). The z^′-reference frame thus moves with that ice as it descends towards the bed, and material trajectories that do not start at the surface at t=0 converge towards $z^{\prime }=\mathrm{0}$ due to vertical compression. An impurity signal centred initially at $z^{\prime }=\mathrm{0}$ can migrate by anomalous diffusion into $z^{\prime }>\mathrm{0}$ (dashed curve) and diffuse, while it experiences vertical compression. The z^′ coordinate system is used in the signal evolution experiments reported in Figs. 5–12.

3 Ice-core numerical experiments

Using the model above, which revises and extends Rempel's theory, we proceed to simulate the distortion of impurity signals in ice cores in both Greenlandic and Antarctic settings, by doping shallow ice with chemical peaks and seeing how they evolve. Most experiments explore predictions with a fully connected vein network, but we include some that simulate partial disconnections (“blocked veins”) by decreasing the molecular diffusivity D.

3.1 Material reference frame and set-up

A general simulation of Eqs. (4), (5), and (17)–(19) would couple them to time-varying velocity, temperature, and grain-size distributions in the ice column. But while the corresponding ice flow and thermal calculations are well established (e.g. Cuffey and Paterson, 2010), reliable grain-size modelling remains out of reach, especially for deep ice, where the way in which the mean grain size d_g is governed by strain-induced recrystallisation processes is poorly understood (Faria et al., 2014; Ng and Jacka, 2014). Consequently, we prescribe time-invariant “background” profiles of w(z), T(z), and d_g(z) in most of our experiments. In our final experiments in Sect. 3.4, we impose short-scale fluctuations on d_g to see what could happen to the impurity evolution, in a simple way without modelling recrystallisation processes.

We track signals in a reference frame moving with the ice, because their length scale (λ≲10⁻¹ m) is much shorter than the core (ice thickness H∼10³ m). The computational burden of modelling a short ice section rather than the whole core is also substantially less. To measure the separation distance of signals from ice material that was deposited on the surface at t=0 and now lies at depth at time t (so t represents age), we use the displacement variable $z^{\prime }=z-g\left(t\right)$, where g(t) is the core's depth–age scale, defined by

$$\begin{array}{}\text{(25)}& t=g^{-\mathrm{1}}\left(\mathit{\zeta }\right)=\underset{\mathrm{0}}{\overset{\mathit{\zeta }}{\int }}{\displaystyle \frac{\text{d}z}{w\left(z\right)}}.\end{array}$$

The change of variable from z to z^′ involves $\partial /\partial z\to \partial /\partial z^{\prime }$ and $\partial /\partial t\to \partial /\partial t-g^{\prime }\left(t\right)\partial /\partial z^{\prime }$, with $g^{\prime }\left(t\right)\equiv w\left(g\right(t\left)\right)$, so Eqs. (17) and (18) become

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+\stackrel{\mathrm{̃}}{w}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z^{\prime }}}={\displaystyle \frac{m}{{\mathit{\rho }}_{\text{i}}}}+{\displaystyle \frac{\partial }{\partial z^{\prime }}}\left(\mathit{\kappa }{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z^{\prime }}}\right),\end{array}$$

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+\stackrel{\mathrm{̃}}{w}{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z^{\prime }}}+{\displaystyle \frac{\partial }{\partial z^{\prime }}}\left(qc\right)={\displaystyle \frac{\partial }{\partial z^{\prime }}}\left(\mathit{\kappa }{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z^{\prime }}}\right)-c_{\mathrm{B}}{\displaystyle \frac{\partial w_c}{\partial z^{\prime }}},\end{array}$$

with the ice velocity in the new reference frame $\stackrel{\mathrm{̃}}{w}$ given by

$$\begin{array}{}\text{(28)}& \stackrel{\mathrm{̃}}{w}(z^{\prime },t)=w\left(g\right(t)+z^{\prime })-w\left(g\right(t\left)\right).\end{array}$$

This is the velocity field at time t seen from the (moving) ice material at depth g(t). Figure 3 depicts characteristic curves representing material trajectories on the z−t and $z^{\prime }-t$ plots under the typical compressive flow at a divide, where $-\text{d}w/\text{d}z$ is the vertical strain rate. As a signal descends, layer thinning causes it to narrow, but this is countered by our newly discovered Gibbs–Thomson diffusion and residual diffusion, while anomalous (Rempel) diffusion displaces it from the ice if the temperature gradient is non-zero.

Figure 4 Ice-core background fields used in our signal evolution experiments. GRIP ice core: (a)–(d). (a) Age–depth scale and (b) ice velocity from a Dansgaard–Johnsen model. Dashed line in (a) shows the GICC05 modelext timescale (Seierstad et al., 2014; Rasmussen et al., 2014). (c) Ice temperature from Johnsen et al. (1995). (d) Grain-size data from Thorsteinsson et al. (1997) and spline fit used in our modelling. EPICA Dome C core: (e)–(h). (e) Age–depth scale and (f) ice velocity from the model described in Sect. 3.1. (g) Borehole temperature from Pol et al. (2010). (h) Grain-size data from Durand and Weiss (2004) and spline fit used in our modelling.

To derive a self-contained evolution equation for c_B, we combine Eqs. (27), (26), and (19) by using the same substitutions as those leading to Eq. (20), finding

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& r\left({\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+\stackrel{\mathrm{̃}}{w}{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z^{\prime }}}\right)+w_c{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z^{\prime }}}=\\ & -c_{\mathrm{B}}{\displaystyle \frac{\partial w_c}{\partial z^{\prime }}}+{\displaystyle \frac{\partial }{\partial z^{\prime }}}\left(\mathit{\kappa }{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z^{\prime }}}\right)\\ & -q{\displaystyle \frac{\partial c}{\partial z^{\prime }}}-(\mathrm{1}-r)\left[\mathit{\varphi }\left({\displaystyle \frac{\partial c}{\partial t}}+\stackrel{\mathrm{̃}}{w}{\displaystyle \frac{\partial c}{\partial z^{\prime }}}\right)+c{\displaystyle \frac{\partial }{\partial z^{\prime }}}\left(\mathit{\kappa }{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z^{\prime }}}\right)\right].\end{array}\end{array}$$

The terms on the right-hand side here involving c and ϕ are negligible, as scaling shows that they are $\sim \mathit{\lambda }(\mathrm{1}-r)/H\sim {\mathrm{10}}^{-\mathrm{5}}$ times those terms on the left. This is verified numerically in all of our experiments. Therefore we approximate Eq. (29) as

$$\begin{array}{}\text{(30)}& r{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+(r\stackrel{\mathrm{̃}}{w}+w_c){\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z^{\prime }}}=-c_{\mathrm{B}}{\displaystyle \frac{\partial w_c}{\partial z^{\prime }}}+{\displaystyle \frac{\partial }{\partial z^{\prime }}}\left(\mathit{\kappa }{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial z^{\prime }}}\right).\end{array}$$

In our simulations, Eq. (30) is solved by the explicit finite-difference method, with the anomalous velocity computed from $w_c=-(D/c)\partial c/\partial z$ (Appendix A details the solution for c from Eqs. 4 and 5), the diffusivity κ calculated from Eq. (9), and $\partial c_{\text{B}}/\partial z=\mathrm{0}$ prescribed at the z^′-domain boundaries (this leads to no evolution there), far from the signal of interest. The doped signal is introduced at a depth near the firn–ice transition (≈ 100 m), in ice whose age t corresponds to that depth. The subsidiary variables ϕ, r_v, m, and q are also calculated from c and c_B at each time step. Rempel et al. (2001) argued that $\left|c_{\text{B}}\right(\partial w_c/\partial z\left)\right|\ll \left|w_c\right(\partial c_{\text{B}}$/$\partial z\left)\right|$, but we do not ignore the term $-c_{\text{B}}(\partial w_c$/$\partial z^{\prime })$ on the right-hand side of Eq. (30), because the full flux divergence $\partial \left(w_c{c}_{\text{B}}\right)/\partial z^{\prime }$ is needed for solute conservation, i.e. no leakage.

We experiment with two sets of background profiles (Fig. 4), based on the glaciological conditions at the GRIP ice-core site in central Greenland and the EPICA Dome C core site in Antarctica. In the GRIP runs, we use the depth–age scale t(z) and velocity w(z) from a Dansgaard–Johnsen model with the ice thickness H=3029 m, the kink at 1000 m above the bed, and surface accumulation rate a= 0.23 m yr⁻¹ ice equivalent. In the EPICA runs, we use t(z) and w(z) from the model $w=m_{\text{base}}+(a-m_{\text{base}})\left[\right(H-z)/H){]}^n$ (Ritz, 1992) with H=3275 m, n=1.7, a=0.023 m yr⁻¹, and the basal melt rate m_base=0.0008 m yr⁻¹, which yields a depth–age scale approximating the one published by Parrenin et al. (2007). Smoothed versions of T(z) and d_g(z) measured at the ice-core sites are used (Fig. 4c–d, g–h). The prescribed profiles are exemplative only. In reality, ice at different depths has experienced different glaciological conditions due to changing accumulation, ice-sheet elevation, and climatic temperature over interglacial–glacial timescales. Our interest is not in reconstructing the histories of these conditions.

Figure 5 Modelled evolution of a signal peak in (a)–(c) the GRIP control run and (d)–(f) an otherwise identical run where the Gibbs–Thomson effect is turned off (γ=0). Snapshots are shown in the material reference frame, with displacement z^′ measuring how far the signal has moved from ice of the same age (which lies at $z^{\prime }=\mathrm{0}$). (a, d) Bulk solute concentration c_B; (b, e) vein solute concentration c at one time; (c, f) anomalous velocity w_c, porosity ϕ, and vein curvature r_v at one time. Grey curves in (a) and (d) indicate the initial doped peak. Panel (b) illustrates the Gibbs–Thomson perturbation. See Movie S1 for the full simulations.

3.2 Results: single-peak experiments

Figure 5a presents snapshots from a GRIP run of the evolution of a decimetre-scale signal doped as a Gaussian peak (grey curve: $c_{\text{B}}=\mathrm{1}+\mathrm{5}\exp [-(z^{\prime }/\mathrm{\Delta }{)}^{\mathrm{2}}$], with Δ=0.08 m) in ice 500 years old (z=112.4 m). Initially the peak, centred at $z^{\prime }=\mathrm{0}$, has a “full width at half maximum” (FWHM) of 0.13 m. Its set amplitude, 5 µM, is based on the size of commonly observed peaks in ice-core records (e.g. ∼ 600 µg L⁻¹ for SO${}_{\mathrm{4}}^{\mathrm{2}-}$; ∼ 150 µg L⁻¹ for Cl⁻; ∼ 80 µg L⁻¹ for Na⁺). The peak decays rapidly in the first 20 kyr (upper 2 km at the GRIP site) with negligible migration and migrates into $z^{\prime }>\mathrm{0}$ more noticeably afterwards, as the ice section descends deeper where the temperature gradient increases (Fig. 4c). Movie S1 shows the full evolution of this control run. Strong diffusion of the signal is evident not just from the peak's decay, but also from its broadening, which overcomes the effect of vertical compression. Recall that in the material reference frame, compression shortens the section continually, so ice enters the simulation domain at both ends. Figure 5b and c exemplify the perturbation on c (caused by the c_B peak) and the resulting large wiggle on the velocity w_c, which represents the γ contribution in Eq. (22) and is what causes the Gibbs–Thomson diffusion. As in Rempel's theory, the signals on ϕ and r_v are collocated with the c_B peak throughout the evolution.

To check this diagnosis for the origin of signal diffusion, another run is conducted (Fig. 5d–f; Movie S1) with everything unchanged except that the Gibbs–Thomson term in Eq. (5) is turned off by setting γ=0. As expected, the strong diffusion in the control run disappears, as no perturbation now arises on c and w_c, but there is still residual diffusion from vein motion. The peak narrows under vertical compression without much amplitude reduction until t≈20 kyr. It subsequently decays because strong c_B gradients on its steepening sides amplify the residual diffusion, despite κ being small (∼ 10⁻⁸ m² yr⁻¹; Fig. B2). The peak's migration trajectory in this run is identical to that in the control run because migration is independent of the peak shape and the diffusion mechanisms. By t≈100 kyr (≈ 2800 m depth) it has displaced from the ice by ≈0.6 m. A further experiment with $\mathit{\kappa }=\mathit{\gamma }=\mathrm{0}$ (not shown) reproduces the “Rempel limit” of a migrating peak with no diffusion, as far as its diminishing width can be resolved by our z^′-grid spacing, 0.0025 m. This implies that the simulated signal behaviour in the experiments is not due to numerical diffusion in our finite-difference scheme. Finally, repeating the control run with κ=0 modifies the results in Fig. 5a only slightly, confirming that residual diffusion becomes important only when a signal becomes very narrow.

Figure 6 Modelled evolution of a signal peak in (a)–(c) the EPICA control run and (d)–(f) an otherwise identical run where the Gibbs–Thomson effect is turned off (γ=0). (a, d) Bulk solute concentration c_B; (b, e) vein solute concentration c at one time; (c, f) anomalous velocity w_c, porosity ϕ, and vein curvature r_v at one time. Grey curves in (a) and (d) indicate the initial doped peak. Panel (b) illustrates the Gibbs–Thomson perturbation. See Movie S2 for the full simulations.

Figure 6a–c and Movie S2 present the control run for EPICA, where ice 4 kyr old (z= 89.9 m) is doped with the same peak. The simulated behaviour is similar to that in the GRIP run, but occurs on a much longer timescale due to the low accumulation rate at the EPICA site. The peak migrates from the start because a sizeable temperature gradient spans the ice column (Fig. 4g). Low compressive strain rate, coupled with slow ice submergence and long time for diffusion, yields a wider peak at all depths than in the GRIP run that has a vastly increased “age span” (i.e. the peak's width in the age domain; discussed later in Fig. 8c, f) compared to the doped signal. Again, comparison against a run with γ= 0 (Fig. 6d–f; Movie S2) confirms the Gibbs–Thomson perturbation as the cause of signal damping and broadening and illustrates the weaker residual diffusion.

The rapidity of signal widening versus migration in distorting the peak in both control runs (Figs. 5a and 6a) is anticipated by the non-small dimensionless number χ in Sect. 2.5. According to Eq. (24), near-constant temperature in the top half of the GRIP column (Fig. 4c) preconditions a large χ there. Indeed, signal diffusion dominates that part of the GRIP control run, also confirming its independent operation from migration. Deeper in both cores, migration becomes more significant as χ is reduced by higher T and higher $\text{d}T/\text{d}z$. In theory, larger grain sizes near the bed (Fig. 4d, h) also slow the rate of the simulated signal decay, but the peaks in our control runs have long dissipated before reaching such depths.

These initial runs demonstrate the signal migration of Rempel's theory, but paradoxically highlight that signals may not survive deep into the ice where it predicts their displacement to become so large to be palaeoclimatically important. More precisely, some remnant signals always survive, but with such small amplitudes and such large age spans compared to the original signals that all essential palaeoclimatic information has been lost. There is an apparent problem to resolve, as distinct deep ionic peaks are found in many ice cores (Fig. 1) (e.g. Röthlisberger et al., 2008; Traversi et al., 2009; Svensson et al., 2013; Tison et al., 2015; Schüpbach et al., 2018), although they may be due to impurities outside veins.

Sticking with the vein model for now, can peaks with a different shape survive damping and broadening to reach deep ice? We study this by changing the width of the doped peaks, as this alters their flank gradient, which is a key control of their diffusion rate. Sensitivity experiments are conducted by varying the width parameter Δ of the Gaussian function between 0.02 and 0.32 m (with the control run parameters unchanged), and by tracking the amplitude, FWHM (full width at half maximum), and age span of the peak in each simulation. The age span is found by dividing the FWHM by the local ice velocity w.

Figure 7 Changing morphometry of the signal peak – its amplitude, full width at half maximum (FWHM), and age span – in the GRIP ice core for different model parameters, plotted against depth (a–c) and age of the ice (d–f). Black curves plot the control run of Fig. 5a, and orange curves plot the γ=0 run of Fig. 5d. Grey curves plot the results of altering the width parameter Δ of the doped peak from 0.08 (control) to four other values. Blue curves plot the outcomes of suppressing molecular diffusivity D in the control run by the multiplicative factors 0.1, 0.01, and 0.001 to simulate vein blockage. Parameter labels use the same colours as the curves. Peak width becomes difficult to measure as amplitude diminishes, explaining the jittery appearance of some curves at depth.

Figure 8 Changing morphometry of the signal peak – its amplitude, full width at half maximum (FWHM), and age span – in the EPICA ice core for different model parameters, plotted against depth (a–c) and age of the ice (d–f). Black curves plot the control run of Fig. 6a, and orange curves plot the γ=0 run of Fig. 6d. Grey and blue curves document the same sensitivity tests as conducted for the GRIP core (see Fig. 7 caption for details). The FWHM and age-span axes are scaled to focus more on the blue and orange curves, rather than the deep ends of the grey/black curves, as the corresponding signal amplitudes decay to near zero.

Figures 7 and 8 plot – for GRIP and EPICA, respectively – the evolving peak morphometry in these “Δ experiments” (grey curves). The control runs (black) and the runs where the Gibbs–Thomson effect has been turned off (orange) are included for comparison. As shown by the grey curves, doping a narrower initial peak hampers its survival, as its steep sides cause strong diffusional draw-down of amplitude; broader peaks retain amplitude for longer but meet the same fate as they narrow under vertical compression. We observe an interesting feedback between width and amplitude evolution. Compression steepens the flanks of signals to accelerate their damping, whereas amplitude reduction makes them shallower and less prone to damping and broadening. Thus the compressive strain rate is a key driver of signal diffusion. The balance of compression and broadening causes different peaks to end on similar width trajectories at depth (panels b and e, Figs. 7 and 8). Accordingly, peaks with different initial time durations acquire near-equal age spans increasing down core (panels c and f, Figs. 7 and 8), which define the minimum time resolution for deep climate signals. These interactions are absent from the study of Rempel et al. (2001), who did not simulate signal shape evolution. Their companion papers (Rempel et al., 2002; Rempel and Wettlaufer, 2003) did so but excluded layer thinning and diffusion.

Figure 9 Modelled displacement ${z^{\prime }}_{\text{pk}}$ and age offset Δt_pk of signal peaks (from ice of the same age) at the (a, b) GRIP and (c, d) EPICA core sites for different parameters, plotted against depth and age of the ice. “Control” labels the control runs in Figs. 5a and 6a; 0.1 and 0.01 label those runs in Figs. 7 and 8 where the molecular diffusivity D is suppressed by these factors to simulate vein blockage. The control curve in panel (a) is equivalent to the curve in Fig. 4 of Rempel et al. (2001), except these authors assumed a different age–depth scale from ours.

So, can single peaks survive into deep ice? The Δ experiments show that peaks at the decimetre or centimetre scale struggle to do so. Even for initially wide peaks (e.g. Δ= 0.32 m) near the firn–ice transition, the Gibbs–Thomson diffusion has reduced their amplitude 4-fold by the time they reach z≈ 2300 m at GRIP (where the age is ≈ 25 ka) and 2000 m at EPICA (≈ 175 ka). Setting γ=0 prolongs the signals' survival (Figs. 5–8), but residual diffusion still prevents them from reaching the lowest several hundred metres with a sizeable fraction of their original amplitude, not to mention that ignoring the Gibbs–Thomson effect is unphysical.

In the present theoretical framework, is there any way for signals to reach deep ice without losing integrity (amplitude, narrowness)? One possibility is the suppression of solute transport by partial vein blockage/disconnection, which we simulate here in a crude manner by artificially decreasing the molecular diffusivity D – this cannot capture heterogeneous vein transport at the grain scale. In Figs. 7 and 8, the blue curves plot the results of simulations with D suppressed by different factors. The same doped peak and parameters of the control runs are used otherwise. A suppression factor of 0.001–0.01 postpones signal decay to a similar extent as turning off the Gibbs–Thomson effect. The lesser factor (0.001) allows the peak to reach 2750 m with half its original amplitude. Even with such strong suppression, however, peak survival is hindered in deeper ice because the low strain rate there (Fig. 4b, f) provides ample time for signal diffusion to occur, and because rising temperature near the bed increases κ. Note that in Figs. 7 and 8, a perfectly preserved peak signal that does not diffuse would have constant amplitude and age span, and its FWHM would decrease towards the bed as a result of vertical compression.

The simulated displacement, age offset, and age span of the peaks are of potential palaeoclimatological interest. Figure 9 shows that in the control runs the peaks displace by ∼ 1 m or more in deep ice, causing their apparent age to exceed their true age by hundreds of years in ice ≈ 100 ka at GRIP, and by several thousand years in ice ≈ 500 ka at EPICA. Since the migration rate is independent of the signal shape, the Δ experiments yield the same displacements and offsets as the control experiments. For both sets of experiments, Figs. 7 and 8 (panels c and f) show that peaks arriving in deep ice have age spans of several hundred years at GRIP and several thousand years at EPICA (approaching the precession timescale in ice 700–800 kyr old); however, we caution against using these results to evaluate deep climatic histories retrieved from ice cores, because these sets of experiments predict near-zero signal amplitude at such depths. In contrast, decreasing D suppresses both signal migration and diffusion (see Eq. 23), so the corresponding peaks remain much narrower during their evolution (Figs. 7 and 8, panels b and e) and migrate much less than in the control/Δ runs (Fig. 9, numbered curves). A suppression factor of 0.001–0.01 enables a peak with FWHM < 0.2 m, age span < 200 yr, and potentially detectable amplitude to reach ≈ 2900 m depth, with an age offset of < 50 yr at GRIP and < 300 yr at EPICA. These numerical findings are illustrative, as they depend on the depth–age scale assumed for each site (notably its precise behaviour at depth) and are limited by the fact that we are not solving the inverse problem with time-varying palaeoclimatic forcing. Using them to interpret specific details of the ice-core records is not advisable at this stage because how much matrix or grain-boundary impurities contribute to those records is also unknown (Sect. 4).

For completeness, all of the above experiments have been repeated with doped peaks with twice the amplitude (10 µM), to cater for some especially high (rarer) peaks in the observed records, which may have more chance to survive. Although the corresponding remnant signals retain greater bulk concentrations at all depths than before, their pattern of decay relative to the initial amplitude and the FWHM and age span results are only marginally altered (Figs. S1 and S2 – see Supplement; cf. Figs. 7 and 8). Our single-peak experiments thus confirm the difficulty for palaeoclimate information to be preserved at depth when the veins are connected.

Figure 10 Snapshots (at four times) of the evolution of two neighbouring peaks in a GRIP run that uses the control parameters of the run in Fig. 5a. Diffusional spreading causes the peaks to merge as they approach each other under vertical compression. See Movie S3 for the full simulation.

3.3 Results: multiple-peak experiments

The diffusion of c_B means that neighbouring peaks can merge as they descend the ice column. This process is illustrated in Fig. 10 and Movie S3 by a simulation with two peaks. Their merging begins at ≈ 7 kyr; the deeper peak moves towards $z^{\prime }=\mathrm{0}$ due to vertical ice compression; a bimodal signal ceases at ≈ 12 kyr. Such merging suggests a second explanation for why distinct peaks can feature in deep ice even with strong damping: instead of deriving from a single peak high up in the column, a deep peak might form by the agglomeration of multiple signals or peaks as these merge under compression. This signal-forming mechanism may not be evident from the c_B profile measured from ice cores, which gives an instantaneous record of the signals.

To test this idea, in the next experiments we simulate the evolution of multiple signals doped in shallow ice stretches 20 m long at GRIP and 80 m long at EPICA. Three runs are made for each site, one with the control-case parameters and two with D suppressed by 0.1 and 0.03 (to simulate vein blockage), with both the Gibbs–Thomson effect and residual diffusion included, and using an initial c_B profile formed by adding many Gaussian peaks (numbering 300 at GRIP and 1200 at EPICA) of random amplitudes, widths, and positions onto a 1 µM base level (Fig. 11a). Movies S4 and S5 document these runs.

Figure 11 Simulated evolution of a long sequence of doped signals at the EPICA site, in three experiments assuming the molecular diffusivities D (Table 1), 0.1D, and 0.03D. (a) Initial c_B sequence. Used in all three runs and made from the superposition of 1200 decimetre-scale Gaussian peaks, it is not meant to recreate the actual signals at EPICA. (b)–(e) Snapshots of c_B at later times. Labels near the vertical line indicate the depth and age of the ice at $z^{\prime }=\mathrm{0}$. The dashed boxes trace a group of signals as they evolve into new signals through compression–diffusion merging (Sect. 3.3). As D is reduced, signal persistence into deep ice improves, and signal displacement decreases. The duration of the signal sequence is ≈ 3600 yr in all panels. See Movie S5 for the full simulations.

We focus our analysis on the EPICA runs (Fig. 11; Movie S5), as the GRIP findings are qualitatively similar (although things occur faster there). In the control run (black curves), strong damping and merging smooth the signals rapidly, so c_B(z^′) retains long-scale variations only – and no peaks – at depth. This outcome is consistent with what we learned from the single-peak experiments. When D is reduced (blue and red curves), compressional shortening, with the now slower diffusion, causes a bundle of peaks to merge into new signals that subsume their solute content. This process operates continuously on all signals, with stretches having a high density of peaks turning into peaks and stretches having a low density turning into troughs. The vertical compression is crucial in helping signals maintain their integrity against diffusion.

When D is suppressed by 0.03 (Fig. 11, red curves), we see distinct peaks persisting in deep ice, many traceable back in time to predecessor groups of peaks, rather than a single peak (e.g. dashed boxes). The balance of diffusion and shortening here is such that the deep peaks have similar widths as their shallow counterparts (∼ dm), despite an overall reduction of signal amplitude with depth. The ice in Fig. 11d has shortened by approximately a factor of 10 since the start of the run, so each peak there encapsulates the signals and solute of an original interval some 10 times longer. Signal survival here is aided by the enhanced survival of single peaks due to decreased molecular diffusivity (Sect. 3.2), but also involves the lumping of solute from neighbouring peaks. We find in further experiments (not shown) that when D is reduced even more (suppression factor ≲ 0.001), the peaks continue to narrow into the centimetre range at depth. Figure 11 shows an effect known from the single-peak experiments: a decrease in D reduces signal displacement as well as signal broadening.

The foregoing experiments demonstrate how long-scale averages on c_B at shallow depths – reflecting long-term background levels of impurity input at an ice-core site – evolve to become meaningful variations at depth, as signals are compressed and their fine details filtered out by diffusion. In Fig. 11, the mean level of the 3600-year-long signal sequence is preserved at depth as a bump (of the same duration) about 7 µM above the surrounding ice. Ice-core analyses of the major ions frequently interpret deep features of this kind as reflecting real palaeoclimatic variations on timescales of 10¹–10² kyr (e.g. Mayewski et al., 1997; Wolff et al., 2006; Schüpbach et al., 2018); it is understood that fewer high-frequency palaeoclimatic details are retrievable from deeper ice, due to the finite resolution of ice-core sampling, alongside layer thinning, which causes more time to be encapsulated in a given ice thickness. Our simulations highlight the Gibbs–Thomson effect in vein impurity transport as a further cause of low-pass filtering. Moreover, they show that long-scale signals will migrate under Rempel's anomalous diffusion mechanism unless the vein network is partially disconnected (see the bump's locations in Fig. 11e).

3.4 Grain-size fluctuations and how they create c_B signals

Our experiments have used smooth d_g profiles so far, as the lack of a robust model of grain-size evolution precludes intimate study of how this process intercouples with vein impurity transport. Ice-core records exhibit grain-size fluctuations (e.g. Fig. 4d, h), although continuous measurements of d_g at a high (e.g. decimetre) resolution over long core sections remain rare. How might rapid fluctuations in d_g impact c_B? Here we study this topic with more analysis and a few tentative simulations, discovering a mechanism by which new impurity signals can form post-depositionally as a result of grain-size fluctuations.

In this connection, in the Interactive Discussions of our paper (see interactive comments on The Cryosphere Discussions, RC2 (https://doi.org/10.5194/tc-2020-217-RC2, 2020) and RC3 (https://doi.org/10.5194/tc-2020-217-RC3, 2020) and associated author comments AC2, AC3, and AC4) it was clarified that Rempel et al. (2001) neglected the Gibbs–Thomson effect from the liquidus relation by making the assumption that the vein radii r_v were spatially uniform – the justification for this being an anticorrelation between mean grain size and impurity loading, which has been observed in ice-core records (see Thorsteinsson et al., 1995, 1997; Alley and Woods, 1996, for evidence related to soluble ions). The concept is that grain recrystallisation processes cause d_g to respond to c_B signals at the decimetre scale to prevent variations in r_v at that scale. Thus, whereas Rempel et al. (2002) and Rempel and Wettlaufer (2003) ignored the interfacial effect by invoking smallness of the Gibbs–Thomson term, the alternative assumption addressed here invokes its constancy. Equations (7) and (8) in our model and their general form based on Eq. (5) show that this assumption holds if $d_{\text{g}}^{\mathrm{2}}{c}_{\text{B}}=$ constant. In the Interactive Discussions (see interactive comment on The Cryosphere Discussions, RC2, https://doi.org/10.5194/tc-2020-217-RC2, 2020), it was suggested that the processes in glacier ice might ensure this inverse–square coupling ($d_{\text{g}}\propto c_{\text{B}}^{-\mathrm{1}/\mathrm{2}}$) along the ice column, so that signals migrate by anomalous diffusion, as conceived in Rempel's theory, without the strong (Gibbs–Thomson) diffusion found in our study.

While such a theory is difficult to prove or disprove until firm quantitative modelling has been offered to explain the causal mechanisms in ice linking grain size to impurities, some obstacles for it are apparent: (i) for diffusion to vanish, $d_{\text{g}}^{\mathrm{2}}{c}_{\text{B}}=$ constant needs to be obeyed identically. Anticorrelation between d_g and c_B does not generally imply suppression of diffusion, since the analysis of Sect. 2.5, when repeated with $d_{\text{g}}\propto c_{\text{B}}^{-\mathit{\beta }}$ (β>0, $\ne \mathrm{1}/\mathrm{2}$), shows that a diffusion term still arises in Eq. (23), just with a different form (and with either higher or lower rate) than before. (ii) Whether d_g and c_B are anticorrelated in ice cores is in fact unknown, because the observed anticorrelation used to support the theory concerns their total impurity loading, which includes an unknown contribution of matrix and grain-boundary impurities: it is the vein impurity component, not the total, that must satisfy $d_{\text{g}}^{\mathrm{2}}{c}_{\text{B}}=$ constant. (iii) Existing glaciological theories on how soluble ions affect the mean grain size consider the effect of such impurities at grain boundaries, through drag production to limit grain-boundary mobility (e.g. Alley et al., 1986a, b), not the effect of impurities in the veins. Also, as r_v is uniform in the suggested theory, the solute concentration of the vein water (c) will be uniform even where c_B has a signal (see Eq. 7); hence it is elusive how c_B can influence grain-boundary motion to control d_g. It is also unclear how $d_{\text{g}}^{\mathrm{2}}{c}_{\text{B}}=$ constant is to be maintained despite other grain-size controls, e.g. microparticle abundance and stored strain energy (Faria et al., 2014; Ng and Jacka, 2014). A further complexity that hampers such maintenance in glacier ice is that the melting point depression is regulated by the overall mix of multiple species of dissolved chemicals, while individual species influence the grain size differently. In the light of these considerations, the idea of signal migration with diffusion eliminated by a highly specific coupling between d_g and c_B seems problematic.

The general question about the impact of grain-size fluctuations on c_B is still of interest. Partial insights into this can be gained with our model. In the diffusion analysis of Sect. 2.5, if d_g varies on short distances, then the anomalous velocity found by substituting for c from Eq. (8) into Eq. (14) (using the chain rule when differentiating $d_{\text{g}}^{-\mathrm{1}}{c}_{\text{B}}^{-\mathrm{1}/\mathrm{2}})$ has a longer form

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{\mathit{u}}_c={\displaystyle \frac{D}{(T_{\mathrm{0}}-T)}}& \left[\mathrm{\nabla }T-{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{\mathrm{2}{\mathit{\rho }}_{\text{i}}L{d}_{\text{g}}}}\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\alpha }(T_{\mathrm{0}}-T)}{\mathrm{\Gamma }}}}{c}_{\mathrm{B}}^{-\mathrm{3}/\mathrm{2}}\mathrm{\nabla }{c}_{\mathrm{B}}\right.\\ & \left.-{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{{\mathit{\rho }}_{\text{i}}L{d}_{\text{g}}^{\mathrm{2}}}}\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\alpha }(T_{\mathrm{0}}-T)}{\mathrm{\Gamma }}}}{c}_{\mathrm{B}}^{-\mathrm{1}/\mathrm{2}}\mathrm{\nabla }{d}_{\text{g}}\right]\end{array}\end{array}$$

(cf. Eq. 22), and the signal evolution equation in Eq. (21) becomes

$$\begin{array}{}\text{(32)}& \begin{array}{rl}{\displaystyle \frac{\partial c_{\mathrm{B}}}{\partial t}}+\mathrm{\nabla }\cdot \left({\displaystyle \frac{D\mathrm{\nabla }T}{T_{\mathrm{0}}-T}}{c}_{\mathrm{B}}\right)& =D{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{\mathrm{2}{\mathit{\rho }}_{\text{i}}L}}\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\alpha }}{\mathrm{\Gamma }}}}\mathrm{\nabla }\cdot \left({\displaystyle \frac{c_{\text{B}}^{-\mathrm{1}/\mathrm{2}}\mathrm{\nabla }{c}_{\text{B}}}{d_{\text{g}}\sqrt{T_{\mathrm{0}}-T}}}\right)\\ & +D{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{{\mathit{\rho }}_{\text{i}}L}}\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\alpha }}{\mathrm{\Gamma }}}}\mathrm{\nabla }\cdot \left({\displaystyle \frac{c_{\text{B}}^{\mathrm{1}/\mathrm{2}}\mathrm{\nabla }{d}_{\text{g}}}{d_{\text{g}}^{\mathrm{2}}\sqrt{T_{\mathrm{0}}-T}}}\right).\end{array}\end{array}$$

As in Sect. 2.5, this model equation treats the scenario without ice deformation or compression. The extra term involving ∇d_g on its right-hand side (cf. Eq. 23) indicates that grain-size fluctuations will generate signals even if c_B has no initial signals. This source term balances with the diffusion and advection terms at large t to cause a steady-state signal locked to the d_g fluctuation. The term representing Gibbs–Thomson diffusion still operates, although it does not damp the signal completely when a source term is present. If ∇T=0, steady state is achieved with $d_{\text{g}}^{\mathrm{2}}{c}_{\text{B}}=$ constant (so c_B mirrors the d_g fluctuation). This is the outcome of the evolution, not the external process coupling that is prerequisite for signal migration without diffusion (discussed above), because ∇T=0 implies zero advection. If ∇T≠0, the steady-state signal will not satisfy $d_{\text{g}}^{\mathrm{2}}{c}_{\text{B}}=$ constant.

Figure 12 Four modified GRIP control runs, demonstrating how fluctuations in the mean grain size d_g (a, c, e) at decimetre scale cause signals to form on the bulk concentration c_B (b, d, f). All runs begin at t=500 yr with c_B≡1 µM without an initial impurity signal. Evolution snapshots are shown in the material reference frame. (a, b) Experiment imposing a negative d_g fluctuation of a fixed width. (c, d) Experiment imposing a negative d_g fluctuation that has the same initial form as in (a), but which narrows due to vertical compression of the ice. (e, f) Two experiments with a positive d_g fluctuation, set up as in the last two experiments. One run assumes that the fluctuation does not experience compression (black); the other run assumes that it does (grey). In those runs where the d_g fluctuation is compressed, only results up to 30 kyr are shown, because soon afterwards its width becomes too narrow to be resolved by the numerical grid spacing (0.0025 m). In all runs, the new signal in c_B is localised by the grain-size fluctuation and does not displace into $z^{\prime }>\mathrm{0}$ by anomalous diffusion. See Movie S6 for the full simulations. Figure S3 and Movie S7 present the equivalent experiments for the EPICA core site.

The mechanism creating a new signal can be explained in terms of vein interactions, most easily by considering ice with uniform c_B initially, without compression or temperature gradient. A short stretch with smaller grain size has higher vein length density and lower r_v (see Eqs. 1 and 4). This implies a stronger Gibbs–Thomson effect (in Eq. 6), so that, at a given temperature, the melting-point depression due to solute must be less; i.e. c must be lower in the stretch than in the surrounding ice. The smaller grain size thus sets up a negative perturbation in c – by lowering the curve of Eq. (4) and the intersection point in Fig. 2f – and the resulting concentration gradients drive solute diffusion that raises c_B over the stretch. In this example, the c_B signal will grow until r_v becomes spatially uniform. This happens when $d_{\text{g}}^{\mathrm{2}}{c}_{\text{B}}$ attains the same value everywhere; the vein conditions inside and outside the stretch are then characterised by the same intersection point in Fig. 2f. If the ice deforms or a temperature gradient exists, a signal would still emerge, but its evolution and final shape (if steady) will be different.

To verify these expectations, we perform modified GRIP and EPICA control runs, using the same simulation model and background profiles as in Sect. 3.1, without doping an initial signal on c_B. For d_g, we impose a constant background of 4 mm for simplicity, with a fluctuation $\mathrm{\Delta }{d}_{\text{g}}=\pm \mathrm{3}\exp [-(z^{\prime }/\mathrm{0.1}{)}^{\mathrm{2}}$] superposed on top (Fig. 12a, c, e). Simulations are made with Δd_g negative or positive – representing, respectively, grain fining or coarsening in a decimetre-wide band – and with Δd_g having a fixed width or narrowing under compression (the c_B signal is always vertically compressed). Thus four runs are made for each site. These runs are artificial as our d_g prescriptions are not physically based or necessarily realistic.

Figure 12 and Movie S6 show the GRIP results. In the two runs with negative Δd_g (Fig. 12a–d), a W-shaped signal in c_B emerges in the first few kiloyears and grows more slowly to a height of ∼ 10 µM. Its form evolves over the kilometre depth scale and does not reach steady state due to continued compression and temperature changes at depth (r_v does not become uniform in long time because of these and ∇T≠0; Movie S6). Its peak is narrower than the d_g fluctuation at all times. Compression of the grain-size fluctuation creates a narrower signal and accelerates signal growth (Fig. 12c, d) by causing high spatial gradients in d_g (see the source term in Eq. 32). In the two runs with positive Δd_g, a much smaller S-shaped signal of ∼ 1 µM forms (Fig. 12e, f; Movie S6). The EPICA runs show similar evolution over a longer time (Fig. S3; Movie S7). Note that the Rempel et al. (2001) equations cannot uncover these findings because they do not contain grain size as a variable, and all signals simulated here are localised by the grain-size fluctuation and do not migrate relative to the ice despite non-zero advection. Furthermore, Eq. (32) indicates that vein blockage or disconnection (suppression of D) would slow the growth of signals.

In the mechanism found here, grain fining creates large peaks in c_B, whereas grain coarsening creates small signals, so an anticorrelation between c_B and d_g is expected for signals with this origin. This property – demonstrated for signals that do not migrate by anomalous diffusion – gives another reason why anticorrelation between c_B and d_g does not guarantee signal migration without damping. The c_B signals here can begin to form at any depth, wherever recrystallisation (e.g. strain-induced) processes perturb d_g at a short scale; signals where d_g varies slowly still decay by Gibbs–Thomson diffusion. The associated concern that ice-core records may contain some signals unrelated to palaeoclimatic events is considered in the next section.

4 Discussion and conclusions

For two decades, Rempel's theory has raised concerns that palaeoclimatic signals in the soluble ion records of ice cores may have displaced by anomalous diffusion and suffer large age discrepancies, especially in the older, deeper parts of the records. Objections to signal migration invoke impeded or insignificant solute transport through veins – that the veins are partially disconnected (blocked by solid impurities and bubbles) or that most chemical signals reside outside veins, in the ice matrix or at grain boundaries.

A more fundamental issue with Rempel's theory is explained herein. While signals on c_B (the bulk concentration of an ion species in the veins) unrelated to grain-size fluctuations can migrate in a connected vein system, a strong Gibbs–Thomson diffusion damps them, preventing those signals at decimetre or shorter scales from penetrating into deep ice. Only much longer background variations in c_B can survive the diffusion to exhibit migration, and signals created by grain-size fluctuations do not migrate (not unless these move relative to the ice by extraneous mechanisms). As the physics predicts that no/few migrating short signals survive into deep ice where their displacement matters, the original concerns are no longer valid. In our revised theory, signal damping is aided by a weaker residual diffusion due to stochastic vein motion. Modifying the derivation of this residual diffusion (Appendix B) to include accelerated grain-boundary motion during migration recrystallisation (Duval and Castelnau, 1995) or grain-growth-driven diffusion (Barnes et al., 2003) strengthens our primary conclusion.

The conclusion is unaffected if we consider multiple solute species interacting via the liquidus (Rempel et al., 2002; Rempel and Wettlaufer, 2003). These authors showed that the c_B signals of different species would line up as they evolve, with periodic signals becoming in phase and peaks in each species inducing collocated “sympathetic peaks” in other species; these adjustments occur in a time of ${\mathit{\lambda }}^{\mathrm{2}}/D\sim$ 1–10 yr for short signals ($\mathit{\lambda }\sim {\mathrm{10}}^{-\mathrm{1}}$ m). Extending our theory for multiple species would thus add to the outcome of an initial fast alignment of signals, before they evolve by the mechanisms studied herein.

What do our findings mean for the integrity and interpretation of ion records from ice cores? And what explains the occurrence of well-defined signal peaks deep in those records? To ponder these, it is useful to start with two end-member scenarios:

Scenario 1. Vein-dissolved ionic impurities (c_B) comprise the main contribution to an observed ice-core record, with matrix and grain-boundary impurity contributions negligible. In this scenario, observed deep peaks can only be explained by (i) widespread vein disconnection (which we modelled by reducing the molecular diffusivity D; Sect. 3.2) and/or (ii) signal formation induced by grain-size fluctuations (Sect. 3.4). The GRIP and EPICA simulations show that in (i), the suppression factor on D needs to be ≲ 0.001 for shallow single peaks to survive into depth, but weaker suppression (∼ 0.01) allows some deep peaks to persist via compression–diffusion merging of signals. Figures 7–9 suggest the possibility of determining the suppression factor from the width, amplitude, and position of the deep peaks and then finding their age offset. In practice, the unknown initial peak size and shape will introduce uncertainty to this estimation.

This scenario spells good news for a key aspect of ice-core palaeoclimatic interpretation. Because reduced vein transport limits signal migration (Fig. 9) and because signals created by grain-size fluctuations are “localised” by them (Fig. 12), the deep peaks would not have displaced far – or at all – to accrue large age offsets. However, diffusional merging means that some signals may be a distorted, lumped signature of multiple climate events. The Gibbs–Thomson diffusion can smooth details more than the low-pass filtering caused by finite resolution of the impurity measurements – typically, ∼ 10 cm in traditional ice-core sampling and ≈ 1 cm or less if using continuous-flow analysis (Kaufmann et al., 2008; Bigler et al., 2011). At the model ice-core sites, the highest time resolution of climatic information retrievable from c_B is quantified approximately by the depth-dependent age spans in Figs. 7 and 8.

Our modelling reveals a new issue if Scenario 1 applies: some observed peaks may have formed as a result of grain-size fluctuations (on short length scales) instigated by glaciological processes during or after deposition. This process, which could impact multiple ion species at the same positions down-core, arises from the effect of grain size on vein equilibrium thermodynamics. It may be misleading to read such signals for palaeoclimate information in the usual way for a given ion, although the grain-size fluctuations producing the signals may have a palaeoclimatic origin.

Scenario 2. Impurities in the ice matrix (e.g. in salt micro-inclusions) and/or at grain boundaries dominate the ice-core record, and their relative immobility explains the presence at all depths of prominent peaks, which experience vertical compression but do not diffuse or migrate by our mechanisms (they may be modified by slow diffusion through ice or along grain boundaries). In this scenario, signals on the minor vein component (c_B) still evolve if the veins are connected – where grain size does not fluctuate rapidly – to migrate, decay, broaden, and merge into long-scale variations.

The reality may be a mixture of these scenarios, with vein and matrix and grain-boundary impurities responsible for different signals on an ion record; their contributions may vary down-core and between cores. Thus generally a record may be the sum of an evolved (e.g. diffused, migrated) component and a largely unmodified component. Shallow signals could source from both components, as firn metamorphism apportions impurities to crystal grains and the premelted liquid. Matrix impurities may dominate deeper, as vein signals decay where the grain size varies slowly. If grain growth and recrystallisation relocates some matrix impurities to grain boundaries and hence to the veins, the signals in c_B could be continually refreshed (at any depth). Such impurity transfer has been suggested (Glen et al., 1977; Mulvaney et al., 1988; Rempel et al., 2001), but also debated (Ohno et al., 2005; Eichler et al., 2019), while it is understood that the apportioning and transfer depend on the ion species (Wolff et al., 1988). The creation of new signals by grain-size fluctuations complicates the vein impurity contribution, as part of the signal evolution becomes coupled to recrystallisation processes that determine historical changes of the ice texture. These considerations caution against interpreting all observed ionic signals directly for palaeoclimatic events and variations: some signals may be distorted in form and duration, and some peaks may be caused by local grain fining (induced by recrystallisation processes (Faria et al., 2014) or high levels of dust/microparticles in the ice (Alley et al., 1986a)). Disentangling the vein, matrix, and grain-boundary impurity components of a given ion record and their post-depositional evolution histories, and comparing the record with grain-size profiles to discern peaks with a potential “grain-size-induced” origin, may be necessary for robust reconstruction. Note that our ideas do not oppose the view that most/much ionic impurity occurs in the matrix (Ohno et al., 2005), which does not strictly rule out the presence of any vein impurities.

Based on these considerations, we conclude that distinct deep peaks seen on a record may indicate (i) matrix impurities dominating the record; (ii) relevance of both matrix and vein impurity signals, with the latter damped out at depth or preserved by vein disconnection; (iii) a dominance of vein impurities in a disconnected vein network; or (iv) a dominance of vein impurities in connected veins that receive recent/sporadic impurity input from the matrix and/or where grain-size fluctuations create and maintain impurity peaks. In each case, the observed peaks would not have migrated or migrated far, although some diffusional merging of signals may occur. Unless many peaks owe their origin to grain-size fluctuations, the limited distortion of peaks inferred for all four cases here is consistent with the signal replicability observed between nearby ice cores (Wolff et al., 2005; Gautier et al., 2016) and supports the use of major ion records for synchronising ice-core age scales. We expect the glaciological conditions at different ice-core sites to cause contrasting distortion of signals and different patterns of grain-size induced peaks. Even if two sites receive identical peak signals at the firn–ice transition (there are no differences in atmospheric transport and deposition of the species and in its modification in firn), variations between their records – for the same ion – will result from differences in strain rate, ice temperature, and dust content/bubble density (which affects the level of vein blockage and recrystallisation processes). Some peaks in one record may be absent or more damped in the other; a group of peaks in one record may appear as a single merged peak in the other record. These variations, which are well known in ice-core studies, impact the identification and matching of peaks and peak sequences. The potential corruption of some records by abundant grain-size-induced signals is an emergent problem that should be studied further.

As we have modelled vein impurity transport only, it is beyond our scope in this paper to disentangle the vein, matrix, and grain-boundary components of any specific (observed) ion records. In fact, some ionic impurities in the veins may not be dissolved. An example is NaCl, whose eutectic point at ≈ −23 ^∘C (pressure-dependent) suggests that it could precipitate in cold ice near the surface but re-dissolve in vein water in warmer ice at depth. Also, our analysis is not aimed at resolving where ionic impurities reside in glacier ice. Recent investigations of ice-core samples using Raman spectroscopy (Barletta et al., 2012; Eichler et al., 2019) have yielded a varied picture regarding the distribution of vein versus matrix impurities in ice, with the former authors finding abundant sulfate and nitrate, but the latter inferring a lack of ionic impurities, at triple junctions (i.e. veins). With an extended literature reporting different results on the subject (see Barnes and Wolff, 2004, and the review parts of the two papers cited above), this state of knowledge suggests that one should not generalise any particular impurity distribution to all ice. Striving to understand the range of processes of impurity movement and segregation and their controls, and how they can cause different impurity distributions, is more important.

Hopefully, with better understanding in this direction, research will be able to develop the theory further by coupling our solute and porosity evolution equations for the veins with equations for the formation, transport, and modification of matrix impurity sites – going beyond a static partitioning of the vein and matrix impurities (Rempel et al., 2002; Rempel and Wettlaufer, 2003). Another direction alluded to before is coupling with grain-size evolution. An extended theory could help palaeoclimatic studies more directly, in terms of quantifying post-depositional changes of measured records, revealing their artefacts, and developing refined palaeoclimatic inversions. We know too little at this stage to envision the details, but such a theory will need to address a multi-directional transfer of ionic impurities between matrix, grain boundaries, and veins. Other foreseeable complications include chemical reactions, grain-boundary motion influenced by impurities, and precipitation and dissolution of different impurities at different temperatures. Work that aids the theory development includes (i) a systematic study of solute signals in ice cores that quantifies their depth-varying spectral content, peak density, and peak size and width frequency statistics and compares the signals to high-resolution grain-size measurements and (ii) controlled laboratory experiments on ice samples to recreate evolving signals for testing the theory.

Finally, our model may be used to quantify the degradation of vein ionic signals during ice-core storage. The rates of residual diffusion, Gibbs–Thomson diffusion, and signal migration (if storage imparts non-zero temperature gradient through a sample) are all minimised at low temperatures (Fig. B2 and Eq. 23). Model runs (Movie S8) show that in ice with a mean grain size of 5 mm and containing a 10 µM high, ≈ 10 cm wide signal in c_B, the total diffusion in 100 years would reduce the signal amplitude by only 5 % and 16 % (with negligible broadening) if the sample is stored at −15 and −5 ^∘C, respectively. These upper-bound reductions assume no vein blockage, and the signal diffusion timescale is still several kiloyears (Movie S8).

Appendix A: Solving for vein conditions

From Dani et al. (2012), relevant empirical formulas for the temperature depressions θ_c and θ_p in Eq. (5) (in kelvin) are

$$\begin{array}{}\text{(A1)}& {\mathit{\theta }}_c=k_{\mathrm{1}}c+k_{\mathrm{2}}{c}^{\mathrm{2}}+k_{\mathrm{3}}{c}^{\mathrm{3}}\text{and}{\mathit{\theta }}_p=a_{\mathrm{1}}p+a_{\mathrm{2}}{p}^{\mathrm{2}},\end{array}$$

where c is measured in molarity, pressure p is in pascals, and the constants are k₁=4.7971, $k_{\mathrm{2}}=-\mathrm{1.188}$, k₃=0.685, $a_{\mathrm{1}}=\mathrm{7.61}\times {\mathrm{10}}^{-\mathrm{8}}$, and $a_{\mathrm{2}}=\mathrm{1.32}\times {\mathrm{10}}^{-\mathrm{16}}$. At a depth z, we calculate the overburden pressure as p=ρ_igz without correcting for firn density. Then combining Eqs. (4) and (5) yields

$$\begin{array}{}\text{(A2)}& f=(T_{\mathrm{0}}-T-{\mathit{\theta }}_p)-k_{\mathrm{1}}c-k_{\mathrm{2}}{c}^{\mathrm{2}}-k_{\mathrm{3}}{c}^{\mathrm{3}}-{\displaystyle \frac{\mathit{\gamma }{T}_{\mathrm{0}}}{{\mathit{\rho }}_{\text{i}}L}}\sqrt{{\displaystyle \frac{\mathrm{3}\mathit{\alpha }c}{c_{\text{B}}{d}_{\text{g}}^{\mathrm{2}}}}}=\mathrm{0}.\end{array}$$

This equation has one positive real root for c. We calculate it numerically with Newton's method, choosing $c=(T_{\mathrm{0}}-T-{\mathit{\theta }}_p)/k_{\mathrm{1}}$ as the initial guess.

Appendix B: Residual diffusion of c_B and ϕ

Consider, in polycrystalline ice, a three-dimensional network of veins with random orientations, which, as the grain boundaries migrate, move in random directions (Fig. B1). Here we show that this can cause diffusion of the bulk solute concentration c_B and porosity ϕ. For simplicity, we assume the vein motion to be isotropic.

Figure B1 Velocities v of different vein segments crossing a plane (outlined by bold rectangle) and their distribution in spherical coordinates (upper right). Only crossings in one direction are shown. In Appendix B, a statistical theory is used to calculate the net transport of porosity and vein impurity resulting from this motion.

As each vein segment migrates, its motion takes along pore space and solute (we ignore vein-water flow, ice deformation, and other processes considered in the main model; Sect. 2.2). Transport arises from vein segments moving in myriad directions. Segments with the same size and solute content moving in opposite directions cancel in terms of contribution. A statistical description is needed to calculate the net effect. Suppose their migration velocities $\mathit{v}=(v_x,v_y,v_z)$ across a given plane (Fig. B1) follow the probability density function f, so that the proportion of vein segments with velocity near v (in an incremental box $\text{d}{v}_x-\text{d}{v}_y-\text{d}{v}_z)$ is

$$\begin{array}{}\text{(B1)}& f(v_x,v_y,v_z)\text{d}{v}_x\text{d}{v}_y\text{d}{v}_z=F\left(v\right)v^{\mathrm{2}}\sin {\mathit{\theta }}^{\prime }\text{d}{\mathit{\theta }}^{\prime }\text{d}{\mathit{\varphi }}^{\prime }\text{d}v.\end{array}$$

The right-hand side puts dv in spherical coordinates (we dashed the symbols of the polar angle θ^′ and azimuth ϕ^′ to distinguish them from ϕ and θ). Under isotropic migration, F is a function of speed $v=\left|\mathit{v}\right|$, independent of direction, and one may suppose F decays to 0 as v→∞. This formulation resembles the kinetic theory of gases, where F is the Maxwell–Boltzmann distribution (Chapman and Cowling, 1953). Regardless of the exact form of F for veins, their mean migration speed is

$$\begin{array}{}\text{(B2)}& \stackrel{\mathrm{‾}}{v}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}F\left(v\right)v^{\mathrm{3}}\sin {\mathit{\theta }}^{\prime }\text{d}{\mathit{\theta }}^{\prime }\text{d}{\mathit{\varphi }}^{\prime }\text{d}v=\mathrm{4}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{v}^{\mathrm{3}}F\left(v\right)\text{d}v.\end{array}$$

We expect $\stackrel{\mathrm{‾}}{v}$ to be similar to the speed of grain boundaries. An estimate for the latter can be found from the normal grain growth law (Gow, 1969; Duval, 1985; Cuffey and Paterson, 2010),

$$\begin{array}{}\text{(B3)}& {\displaystyle \frac{\text{d}\left(d_{\text{g}}^{\mathrm{2}}\right)}{\text{d}t}}=K,\end{array}$$

where K has been defined in Sect. 2.2. Specifically, following the results of Hillert (1965) and Ng (2016), we let

$$\begin{array}{}\text{(B4)}& \stackrel{\mathrm{‾}}{v}={\displaystyle \frac{K}{c_{\mathrm{1}}{d}_{\text{g}}}},\end{array}$$

where c₁≈2 to 3. The assumption of normal grain growth gives a low-end estimate for $\stackrel{\mathrm{‾}}{v}$, because strain-induced dynamic recrystallisation can accelerate grain boundary migration (Duval and Castelnau, 1995). On the other hand, impurities and bubbles may reduce grain-boundary mobility (Alley, 1986a, b). In the present formulation, we exclude these complications, as well as anisotropic vein motion resulting from recrystallisation processes. We ignore any influence on K by the solute concentration c, because c refers to impurities dissolved in vein water, rather than impurities at grain boundaries.

Next we calculate the transport fluxes. Imagine a region of uniform porosity and uniform mean grain size where the vein network continually evolves. The incremental flux of porosity (vein space) crossing an area dA in the direction perpendicular to the plane is

$$\begin{array}{}\text{(B5)}& \text{d}J=\mathit{\varphi }v\cos {\mathit{\theta }}^{\prime }F\left(v\right)v^{\mathrm{2}}\sin {\mathit{\theta }}^{\prime }\text{d}{\mathit{\theta }}^{\prime }\text{d}{\mathit{\varphi }}^{\prime }\text{d}v\text{d}A.\end{array}$$

Integrating this over $v\in [\mathrm{0},\mathrm{\infty }]$, ${\mathit{\varphi }}^{\prime }\in [\mathrm{0},\mathrm{2}\mathit{\pi }]$, and ${\mathit{\theta }}^{\prime }\in [\mathrm{0},\mathit{\pi }/\mathrm{2}]$ (for unidirectional flux; Fig. B1) yields the porosity flux density

$$\begin{array}{}\text{(B6)}& j={\displaystyle \frac{\text{d}J}{\text{d}A}}=\mathit{\varphi }\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}F\left(v\right)v^{\mathrm{3}}\text{d}v={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\stackrel{\mathrm{‾}}{v}\mathit{\varphi }.\end{array}$$

The same method applied to the bulk solute content c_B=ϕc gives its flux density as

$$\begin{array}{}\text{(B7)}& j_c={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\stackrel{\mathrm{‾}}{v}{c}_{\mathrm{B}}.\end{array}$$

Given how these macroscopic fluxes originate from microscopic interactions, they are valid on timescales longer than the timescale of vein-crossing events, $\sim d_{\text{g}}/\stackrel{\mathrm{‾}}{v}$.

Figure B2 Diffusivity κ at different temperatures, calculated with Eq. (9) for c₁=2.5.

In a uniform region, the fluxes through dA in opposite directions cancel. But if ϕ (or c_B) varies spatially, a differential flux occurs; diffusion then arises from the vein motion. The net diffusion rate across a plane, say, at elevation z=z₀, is found by subtracting the opposite fluxes at a distance dz on either side, j₊ at z₀−dz and j₋ at z₀−dz, where dz locates the planes for evaluating the fluxes from continuum properties. We determine dz by using an argument similar to that in the kinetic theory. Moving vein segments merge and reconfigure on distances on the order of the mean grain size, so vein segments arriving at the plane typically come from a distance ∼d_g since their last “collision”, which caused them to switch direction and gain or lose solute; d_g is akin to the particle mean free path in the kinetic theory. Accordingly, we choose dz to be the flux-averaged value of the perpendicular distance d_gcos θ^′:

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}\text{d}z& ={\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\int }_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\int }_{\mathrm{0}}^{\mathit{\pi }/\mathrm{2}}{d}_{\text{g}}\cos {\mathit{\theta }}^{\prime }F\left(v\right)v^{\mathrm{3}}\cos {\mathit{\theta }}^{\prime }\sin {\mathit{\theta }}^{\prime }\text{d}{\mathit{\theta }}^{\prime }\text{d}{\mathit{\varphi }}^{\prime }\text{d}v}{{\int }_{\mathrm{0}}^{\mathrm{\infty }}{\int }_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\int }_{\mathrm{0}}^{\mathit{\pi }/\mathrm{2}}F\left(v\right)v^{\mathrm{3}}\cos {\mathit{\theta }}^{\prime }\sin {\mathit{\theta }}^{\prime }\text{d}{\mathit{\theta }}^{\prime }\text{d}{\mathit{\varphi }}^{\prime }\text{d}v}}\\ & =d_{\text{g}}{\displaystyle \frac{{\int }_{\mathrm{0}}^{\mathit{\pi }/\mathrm{2}}{\cos }^{\mathrm{2}}{\mathit{\theta }}^{\prime }\sin {\mathit{\theta }}^{\prime }\text{d}{\mathit{\theta }}^{\prime }}{{\int }_{\mathrm{0}}^{\mathit{\pi }/\mathrm{2}}\cos {\mathit{\theta }}^{\prime }\sin {\mathit{\theta }}^{\prime }\text{d}{\mathit{\theta }}^{\prime }}}={\displaystyle \frac{\mathrm{2}{d}_{\text{g}}}{\mathrm{3}}}.\end{array}\end{array}$$

It follows that the net diffusive transport across the plane is

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}{j}_{\text{net}}& =j_{+}-j_{-}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\stackrel{\mathrm{‾}}{v}\left({\left.\mathit{\varphi }\right|}_{z_{\mathrm{0}}-\mathrm{2}{d}_{\text{g}}/\mathrm{3}}-{\left.\mathit{\varphi }\right|}_{z_{\mathrm{0}}+\mathrm{2}{d}_{\text{g}}/\mathrm{3}}\right)\\ & \approx {\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\stackrel{\mathrm{‾}}{v}\left(-{\displaystyle \frac{\mathrm{4}{d}_{\text{g}}}{\mathrm{3}}}{\left.{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right|}_{z_{\mathrm{0}}}\right)=-{\displaystyle \frac{\stackrel{\mathrm{‾}}{v}{d}_{\text{g}}}{\mathrm{3}}}{\left.{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right|}_{z_{\mathrm{0}}}.\end{array}\end{array}$$

Combining this result with Eq. (B4) gives the diffusivity $\mathit{\kappa }=\stackrel{\mathrm{‾}}{v}{d}_{\text{g}}/\mathrm{3}=K\left(T\right)/\mathrm{3}{c}_{\mathrm{1}}$, as given in Eq. (9) in the text. Figure B2 plots κ against temperature. Interestingly, κ is independent of d_g because smaller grains lead to faster grain boundary migration but proportionally shorter mean free path for the vein motion.