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 103 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.
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.
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, 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 is used in our calculations). We let cB denote its bulk concentration, using the unit mol L−1 or M to refer to the amount of impurity in a unit volume of ice. Note that cB 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.
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 cB, 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 cB encapsulating vein impurities, the relation cB=cϕ holds1, 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 increases this flux with distance to draw down cB, whereas on the leading (lower) edge of the peak, reduces the flux with distance to bump up cB. Thus the peak signal in cB 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 cB. However, we find that when this seemingly reasonable approximation is not made, the cB 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 cB 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 and m2 yr−1, 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 cB 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 cB 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 cB. This mechanism does not require a temperature gradient or grain-size variations to operate and is unaffected by vein disconnection.
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 dg (grain diameter) has the vein length density
we adopt equality in this expression herein. If the vein faces have the radius of curvature rv (Fig. 2e; Nye, 1989), then the ice porosity ϕ is
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,
and using Eqs. (1) and (2), leads to
For ice at a temperature T below the reference temperature T0, thermodynamic equilibrium between ice and vein water means that the liquidus relation is satisfied:
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−1 for SO. The third term is typically small and may be absorbed into T0 by accounting for glaciostatic overburden. Accordingly, the version of Eq. (5) that we use for the present analysis is
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).
Given dg, T, and cB, Eqs. (4) and (6) can be solved for rv 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 rv makes the Gibbs–Thomson term negligible so that ; 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 cB – albeit weakly – for fixed T and dg. This dependence lies at the heart of our “Gibbs–Thomson diffusion”. Specifically, when , a first-order approximate solution of Eqs. (4) and (6) is
where , or
This differs from Rempel's approximation by the γ (Gibbs–Thomson) term, which causes and a perturbation to appear on c(z) where cB(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 dg and T could perturb c(z) to drive solute transport and influence cB. 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 ([ ]) with [ ] as given in Eq. (8) except that replaces ( 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, , the incompressibility condition 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
where is the temperature-dependent grain growth rate, R is the gas constant, Q is activation energy, and c1≈2 to 3 (we use 2.5 in our simulations).
Accordingly, porosity conservation is described by
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
and we deduce from this together with Eq. (10) that
where is the ratio of ice density to water density.
The local solute transport flux is . 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
Following Rempel's theory, we define the “anomalous velocity”
and use cB=ϕc to rewrite the D contribution in Eq. (13) to derive
If both water flux q and the κ-term (“residual diffusion” due to vein motion) are ignored, this simplifies, after using the incompressibility condition , to the Rempel et al. (2001) equation:
As is well rehearsed in their theory, the advection term here predicts migration of chemical signals at the velocity uc relative to the ice, with uc 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 () and horizontal variations in ϕ and cB are negligible, Eqs. (10) and (15) become (on using again)
where is the anomalous velocity downward. Equation (18) tracks the evolution of the bulk impurity profile cB(z,t) and its signals, but the water flux q needs to be found via
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 cB, T, and dg. There are five equations (Eqs. 4, 5, 17–19) for the five unknowns cB, 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 wc, because of water percolation induced by the migration. To see this, we need to put all advective parts of Eq. (18) (notably ) in terms of cB. Let c(z) be time-invariant as in Rempel's theory, and ignore the κ-term (vein motion), to write Eq. (18) as
Substituting for from Eq. (19), using Eq. (17) for m, converts this to
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−1 m). In Eq. (20), the prefactor r of the time derivative means that the migration speed is , which exceeds wc by ≈ 9 %.
What causes this correction? As introduced in Sect. 1, a migrating cB 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 cB will diffuse, as we emphasised at the outset. In Rempel's theory, cB 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:
For the anomalous velocity, substituting for c from Eq. (8) into Eq. (14) yields, after some algebra,
Note that this result assumes that dg varies slowly. If grain size fluctuates strongly over short distances (e.g. decimetres), the evaluation of ∇c would introduce another term involving ∇dg (Sect. 3.4). With the present result in Eq. (22), Eq. (21) becomes
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 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):
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 K, |∇T|∼20 K km−1 = 0.02 K m−1, dg∼5 mm, cB∼1 µM, and |∇cB|∼1 µM/10 cm = 10−5 M m−1 (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 m2 yr−1, 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.
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 dg 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 dg(z) in most of our experiments. In our final experiments in Sect. 3.4, we impose short-scale fluctuations on dg 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−1 m) is much shorter than the core (ice thickness H∼103 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 , where g(t) is the core's depth–age scale, defined by
The change of variable from z to z′ involves and , with , so Eqs. (17) and (18) become
with the ice velocity in the new reference frame given by
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 plots under the typical compressive flow at a divide, where 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.
To derive a self-contained evolution equation for cB, we combine Eqs. (27), (26), and (19) by using the same substitutions as those leading to Eq. (20), finding
The terms on the right-hand side here involving c and ϕ are negligible, as scaling shows that they are times those terms on the left. This is verified numerically in all of our experiments. Therefore we approximate Eq. (29) as
In our simulations, Eq. (30) is solved by the explicit finite-difference method, with the anomalous velocity computed from (Appendix A details the solution for c from Eqs. 4 and 5), the diffusivity κ calculated from Eq. (9), and 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 ϕ, rv, m, and q are also calculated from c and cB at each time step. Rempel et al. (2001) argued that /, but we do not ignore the term / on the right-hand side of Eq. (30), because the full flux divergence 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−1 ice equivalent. In the EPICA runs, we use t(z) and w(z) from the model (Ritz, 1992) with H=3275 m, n=1.7, a=0.023 m yr−1, and the basal melt rate mbase=0.0008 m yr−1, which yields a depth–age scale approximating the one published by Parrenin et al. (2007). Smoothed versions of T(z) and dg(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.
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: ], with Δ=0.08 m) in ice 500 years old (z=112.4 m). Initially the peak, centred at , 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−1 for SO; ∼ 150 µg L−1 for Cl−; ∼ 80 µg L−1 for Na+). The peak decays rapidly in the first 20 kyr (upper 2 km at the GRIP site) with negligible migration and migrates into 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 cB peak) and the resulting large wiggle on the velocity wc, which represents the γ contribution in Eq. (22) and is what causes the Gibbs–Thomson diffusion. As in Rempel's theory, the signals on ϕ and rv are collocated with the cB 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 wc, 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 cB gradients on its steepening sides amplify the residual diffusion, despite κ being small (∼ 10−8 m2 yr−1; 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 (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 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 . 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.
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.
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.
3.3 Results: multiple-peak experiments
The diffusion of cB 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 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 cB 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 cB 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.
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 cB(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 cB 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 101–102 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 cB signals
Our experiments have used smooth dg 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 dg at a high (e.g. decimetre) resolution over long core sections remain rare. How might rapid fluctuations in dg impact cB? 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 rv 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 dg to respond to cB signals at the decimetre scale to prevent variations in rv 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 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 () 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, constant needs to be obeyed identically. Anticorrelation between dg and cB does not generally imply suppression of diffusion, since the analysis of Sect. 2.5, when repeated with (β>0, ), 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 dg and cB 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 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 rv is uniform in the suggested theory, the solute concentration of the vein water (c) will be uniform even where cB has a signal (see Eq. 7); hence it is elusive how cB can influence grain-boundary motion to control dg. It is also unclear how 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 dg and cB seems problematic.
The general question about the impact of grain-size fluctuations on cB is still of interest. Partial insights into this can be gained with our model. In the diffusion analysis of Sect. 2.5, if dg 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 has a longer form
(cf. Eq. 22), and the signal evolution equation in Eq. (21) becomes
As in Sect. 2.5, this model equation treats the scenario without ice deformation or compression. The extra term involving ∇dg on its right-hand side (cf. Eq. 23) indicates that grain-size fluctuations will generate signals even if cB 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 dg 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 constant (so cB mirrors the dg 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 constant.
The mechanism creating a new signal can be explained in terms of vein interactions, most easily by considering ice with uniform cB initially, without compression or temperature gradient. A short stretch with smaller grain size has higher vein length density and lower rv (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 cB over the stretch. In this example, the cB signal will grow until rv becomes spatially uniform. This happens when 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 cB. For dg, we impose a constant background of 4 mm for simplicity, with a fluctuation ] superposed on top (Fig. 12a, c, e). Simulations are made with Δdg negative or positive – representing, respectively, grain fining or coarsening in a decimetre-wide band – and with Δdg having a fixed width or narrowing under compression (the cB signal is always vertically compressed). Thus four runs are made for each site. These runs are artificial as our dg prescriptions are not physically based or necessarily realistic.
Figure 12 and Movie S6 show the GRIP results. In the two runs with negative Δdg (Fig. 12a–d), a W-shaped signal in cB 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 (rv does not become uniform in long time because of these and ∇T≠0; Movie S6). Its peak is narrower than the dg 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 dg (see the source term in Eq. 32). In the two runs with positive Δdg, 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 cB, whereas grain coarsening creates small signals, so an anticorrelation between cB and dg 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 cB and dg does not guarantee signal migration without damping. The cB signals here can begin to form at any depth, wherever recrystallisation (e.g. strain-induced) processes perturb dg at a short scale; signals where dg 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.
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 cB (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 cB 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 cB 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 1–10 yr for short signals ( 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 (cB) 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 cB 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 (cB) 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 cB 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 cB, 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).
From Dani et al. (2012), relevant empirical formulas for the temperature depressions θc and θp in Eq. (5) (in kelvin) are
where c is measured in molarity, pressure p is in pascals, and the constants are k1=4.7971, , k3=0.685, , and . At a depth z, we calculate the overburden pressure as p=ρigz without correcting for firn density. Then combining Eqs. (4) and (5) yields
This equation has one positive real root for c. We calculate it numerically with Newton's method, choosing as the initial guess.
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 cB and porosity ϕ. For simplicity, we assume the vein motion to be isotropic.
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 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 is
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 , 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
We expect 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),
where K has been defined in Sect. 2.2. Specifically, following the results of Hillert (1965) and Ng (2016), we let
where c1≈2 to 3. The assumption of normal grain growth gives a low-end estimate for , 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
Integrating this over , , and (for unidirectional flux; Fig. B1) yields the porosity flux density
The same method applied to the bulk solute content cB=ϕc gives its flux density as
Given how these macroscopic fluxes originate from microscopic interactions, they are valid on timescales longer than the timescale of vein-crossing events, .
In a uniform region, the fluxes through dA in opposite directions cancel. But if ϕ (or cB) varies spatially, a differential flux occurs; diffusion then arises from the vein motion. The net diffusion rate across a plane, say, at elevation z=z0, is found by subtracting the opposite fluxes at a distance dz on either side, j+ at z0−dz and j− at z0−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 ∼dg since their last “collision”, which caused them to switch direction and gain or lose solute; dg 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 dgcos θ′:
It follows that the net diffusive transport across the plane is
Combining this result with Eq. (B4) gives the diffusivity , as given in Eq. (9) in the text. Figure B2 plots κ against temperature. Interestingly, κ is independent of dg because smaller grains lead to faster grain boundary migration but proportionally shorter mean free path for the vein motion.
Our MATLAB code and the simulated data of our control runs are archived at https://doi.org/10.15131/shef.data.12735191 (Ng, 2021).
Movies S1–S8 are available at: https://doi.org/10.15131/shef.data.12739169. The supplement related to this article is available online at: https://doi.org/10.5194/tc-15-1787-2021-supplement.
The author declares that there is no conflict of interest.
I thank Rita Traversi for providing the data in Fig. 1a and b and Andrew Sole and Adam Hepburn for comments on the manuscript.
This paper was edited by Nanna Bjørnholt Karlsson and reviewed by Alan Rempel, Eric Wolff, and one anonymous referee.
Alley, R. B., Perepezko, J. H., and Bentley, C. R.: Grain growth in polar ice: I. Theory, J. Glaciol., 32, 415–424, 1986a.
Alley, R. B., Perepezko, J. H., and Bentley, C. R.: Grain growth in polar ice: II. Application, J. Glaciol., 32, 425–433, 1986b.
Alley, R. B. and Woods, G. A.: Impurity influence on normal grain growth in the GISP2 ice core, Greenland, J. Glaciol., 42, 255–260, 1996.
Barletta, R. E., Priscu, J. C., Mader, H. M., Jones, W. L., and Roe, C. W.: Chemical analysis of ice vein microenvironments: II. Analysis of glacial samples from Greenland and the Antarctic, J. Glaciol., 58, 1109–1118, https://doi.org/10.3189/2012JoG12J112, 2012.
Barnes, P. R. F. and Wolff, E. W.: Distribution of soluble impurities in cold glacial ice, J. Glaciol., 170, 311–324, 2004.
Barnes, P. R. F., Mader, H. M., Röthlisberger, R., Udisti, R., and Wolff, E. W.: The evolution of chemical peak shapes in the Dome C ice core, Antarctica, J. Geophys. Res., 108, 4126, https://doi.org/10.1029/2002JD002538, 2003.
Bigler, M., Svensson, A., Kettner, E., Vallelonga, P., Nielsen, M. E., and Steffensen, J. P.: Optimization of high-resolution continuous flow analysis for transient climate signals in ice cores, Environ. Sci. Technol., 45, 4483–4489, https://doi.org/10.1021/es200118j, 2011.
Chapman, S. and Cowling, T. G.: The Mathematical Theory of Non-Uniform Gases, Cambridge University Press, London, 1953.
Cuffey, K. M. and Paterson, W. S. B.: The Physics of Glaciers, 4th edn., Butterworth-Heinemann, Oxford, 2010.
Dani, K. G., Mader, H. M., Wolff, E. W., and Wadham, J. L.: Modelling the liquid-water vein system within polar ice sheets as a potential microbial habitat, Earth Planet. Sc. Lett., 333–334, 238–249, https://doi.org/10.1016/j.epsl.2012.04.009, 2012.
Durand, G. and Weiss, J.: EPICA Dome C Ice Cores Grain Radius Data, IGBP PAGES/World Data Center for Paleoclimatology, Data Contribution Series No. 2004-039, NOAA/NGDC Paleoclimatology Program, Boulder CO, USA, 2004.
Duval, P.: Grain growth and mechanical behaviour of polar ice, Ann. Glaciol., 6, 79–82, 1985.
Duval, P. and Castelnau, O.: Dynamic recrystallization of ice in polar ice sheets, J. Phys. IV [Paris], 5, 197–205, https://doi.org/10.1051/jp4:1995317, 1995.
Eichler, J., Weikusat, C., Wegner, A., Twarloh, B., Behrens, M., Fischer, H., Hörhold, M., Jansen, D., Kipfstuhl, S., Ruth, U., Wilhelms, F., and Weikusat, I.: Impurity Analysis and microstructure along the climatic transition from MIS 6 into 5e in the EDML ice core using cryo-Raman microscopy, Front. Earth Sci., 7, 20, https://doi.org/10.3389/feart.2019.00020, 2019.
EPICA community members: Eight glacial cycles from an Antarctic ice core, Nature, 429, 623–628, https://doi.org/10.1038/nature02599, 2004.
Faria, S. H., Weikusat, I., and Azuma, N.: The microstructure of polar ice. Part II: state of the art, J. Struct. Geol., 61, 21–49, https://doi.org/10.1016/j.jsg.2013.11.003, 2014.
Fujita, S., Parrenin, F., Severi, M., Motoyama, H., and Wolff, E. W.: Volcanic synchronization of Dome Fuji and Dome C Antarctic deep ice cores over the past 216 kyr, Clim. Past, 11, 1395–1416, https://doi.org/10.5194/cp-11-1395-2015, 2015.
Gautier, E., Savarino, J., Erbland, J., Lanciki, A., and Possenti, P.: Variability of sulfate signal in ice core records based on five replicate cores, Clim. Past, 12, 103–113, https://doi.org/10.5194/cp-12-103-2016, 2016.
Glen, J. W., Homer, D. R., and Paren, J. G.: Water at grain boundaries: its role in the purification of temperate glacier ice, Int. Assoc. Hydrogeol., 118, 263–271, 1977.
Gow, A. J.: On the rates of growth of grains and crystals in South Polar firn, J. Glaciol., 8, 241–252, 1969.
Hillert, M.: On the theory of normal and abnormal grain growth, Acta Metall., 13, 227–238, 1965.
Iizuka, Y., Takata, M., Hondoh, T., and Fujii, Y.: High-time-resolution profiles of soluble ions in the last glacial period of a Dome Fuji (Antarctica) deep ice core, Ann. Glaciol., 39, 452–456, 2004.
Johnsen, S. J., Dahl-Jensen, D., Dansgaard, W., and Gundestrup, N.: Greenland palaeotemperatures derived from GRIP bore hole temperature and ice core isotope profiles, Tellus, 47B, 624–629, 1995.
Kaufmann, P. R., Federer, U., Hutterli, M. A., Bigler, M., Schüpbach, S., Ruth, U., Schmitt, J., and Stocker, T. F.: An improved continuous flow analysis system for high-resolution field measurements on ice cores, Environ. Sci. Technol., 42, 8044–8050, https://doi.org/10.1021/es8007722, 2008.
Legrand, M. and Mayewski, P.: Glaciochemistry of polar ice cores: a review, Rev. Geophys., 35, 219–243, 1997.
Mader, H. M.: Observations of the water-vein system in polycrystalline ice, J. Glaciol., 38, 333–347, 1992a.
Mader, H. M.: The thermal behaviour of the water-vein system in polycrystalline ice, J. Glaciol., 38, 359–374, 1992b.
Mayewski, P. A., Meeker, L. D., Twickler, M. S., Whitlow, S. I., Yang, Q., Lyons, W. B., and Prentice, M.: Major features and forcing of high-latitude northern hemisphere atmospheric circulation using a 110,000-year-long glaciochemical series, J. Geophys. Res., 102, 26345–26366, 1997.
Mulvaney, R., Wolff, E. W., and Oates, K.: Sulphuric acid at grain boundaries in Antarctic ice, Nature, 331, 247–249, 1988.
Ng, F. S. L.: Statistical mechanics of normal grain growth in one dimension: A partial integro-differential equation model, Acta Mater., 120, 453–462, https://doi.org/10.1016/j.actamat.2016.08.033, 2016.
Ng, F.: Computer code and simulated data of the paper “Pervasive diffusion of climate signals recorded in ice-vein ionic impurities” [code], https://doi.org/10.15131/shef.data.12735191, 2021.
Ng, F. and Jacka, T. H.: A model of crystal-size evolution in polar ice masses, J. Glaciol., 60, 463–477, https://doi.org/10.3189/2014JoG13J173, 2014.
Nye, J. F.: The geometry of water veins and nodes in polycrystalline ice, J. Glaciol., 35, 17–22, 1989.
Nye, J. F.: Thermal behaviour of glacier and laboratory ice, J. Glaciol., 37, 401–413, 1991.
Ohno, H., Igarashi, M., and Hondoh, T.: Salt inclusions in polar ice core: location and chemical form of water-soluble impurities, Earth Planet. Sc. Lett., 232, 171–178, 2005.
Osman, M., Das, S. B., Marchal, O., and Evans, M. J.: Methanesulfonic acid (MSA) migration in polar ice: data synthesis and theory, The Cryosphere, 11, 2439–2462, https://doi.org/10.5194/tc-11-2439-2017, 2017.
Parrenin, F., Barnola, J.-M., Beer, J., Blunier, T., Castellano, E., Chappellaz, J., Dreyfus, G., Fischer, H., Fujita, S., Jouzel, J., Kawamura, K., Lemieux-Dudon, B., Loulergue, L., Masson-Delmotte, V., Narcisi, B., Petit, J.-R., Raisbeck, G., Raynaud, D., Ruth, U., Schwander, J., Severi, M., Spahni, R., Steffensen, J. P., Svensson, A., Udisti, R., Waelbroeck, C., and Wolff, E.: The EDC3 chronology for the EPICA Dome C ice core, Clim. Past, 3, 485–497, https://doi.org/10.5194/cp-3-485-2007, 2007.
Pol, K., Masson-Delmotte, V., Johnsen, S., Bigler, M., Cattani, O., Durand, G., Falourd, S., Jouzel, J., Minster, B., Parrenin, F., Ritz, C., Steen-Larsen, H. C., and Stenni, B.: New MIS 19 EPICA Dome C high resolution deuterium data: Hints for a problematic preservation of climate variability at sub-millennial scale in the “oldest ice”, Earth Planet. Sc. Lett., 298, 95–103, https://doi.org/10.1016/j.epsl.2010.07.030, 2010.
Rasmussen, S. O., Bigler, M., Blockley, S. P., Blunier, T., Buchardt, S. L., Clausen, H. B., Cvijanovic, I., Dahl-Jensen, D., Johnsen, S. J., Fischer, H., Gkinis, V., Guillevic, M., Hoek, W. Z., Lowe, J. J., Pedro, J. B., Popp, T., Seierstad, I. K., Steffensen, J. P., Svensson, A. M., Vallelonga, P., Vinther, B. M., Walker, M. J. C., Wheatley, J. J., and Winstrup, M.: A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy, Quat. Sci. Rev., 106, 14–28, https://doi.org/10.1016/j.quascirev.2014.09.007, 2014.
Rempel, A. W. and Wettlaufer, J. S.: Segregation, transport, and interaction of climate proxies in polycrystalline ice, Can. J. Phys., 81, 89–97, 2003.
Rempel, A. W., Waddington, E. D., Wettlaufer, J. S., and Worster, M. G.: Possible displacement of the climate signal in ancient ice by premelting and anomalous diffusion, Nature, 411, 568–571, https://doi.org/10.1038/35079043, 2001.
Rempel, A. W., Wettlaufer, J. S., and Waddington, E. D.: Anomalous diffusion of multiple impurity species: Predicted implications for the ice core climate records, J. Geophys. Res., 107, 2330, https://doi.org/10.1029/2002JB001857, 2002.
Ritz, C.: Un modèle thermo-mécanique d'évolution pour le basin glaciaire Antarctique Vostok-Glacier Byrd: Sensibilité aux valeurs des paramètres mal connus, PhD thesis, Laboratoire de Glaciologie et Géophysique de l'Environnement, Université Joseph Fourier, Grenoble, 377 pp., 1992.
Röthlisberger, R., Mudelsee, M., Bigler, M., de Angelis, M., Fischer, H., Hansson, M., Lambert, F., Masson-Delmotte, V., Sime, L., Udisti, R., and Wolff, E. W.: The Southern Hemisphere at glacial terminations: insights from the Dome C ice core, Clim. Past, 4, 345–356, https://doi.org/10.5194/cp-4-345-2008, 2008.
Schüpbach, S. Fischer, H., Bigler, M., Erhardt, T., Gfeller, G., Leuenberger, D., Mini, O., Mulvaney, R., Abram, N. J., Fleet, L., Frey, M. M., Thomas, E., Svensson, A., Dahl-Jensen, D., Kettner, E., Kjaer, H., Seierstad, I., Steffensen, J. P., Rasmussen, S. O., Vallelonga, P., Winstrup, M., Wegner, A., Twarloh, B., Wolff, K., Schmidt, K., Goto-Azuma, K., Kuramoto, T., Hirabayashi, M., Uetake, J., Zheng, J., Bourgeois, J., Fisher, D., Zhiheng, D., Xiao, C., Legrand, M., Spolaor, A., Gabrieli, J., Barbante, C., Kang, J.-H., Hur, S. D., Hong, S. B., Hwang, H. J., Hong, S., Hansson, M., Iizuka, Y., Oyabu, I., Muscheler, R., Adolphi, F., Maselli, O., McConnell, J., and Wolff, E. W.: Greenland records of aerosol source and atmospheric lifetime changes from the Eemian to the Holocene, Nat. Commun., 9, 1476, https://doi.org/10.1038/s41467-018-03924-3, 2018.
Seierstad, I. K., Abbott, P. M., Bigler, M., Blunier, T., Bourne, A. J., Brook, E., Buchardt, S. L., Buizert, C., Clausen, H. B., Cook, E., Dahl-Jensen, D., Davies, S. M., Guillevic, M., Johnsen, S. J., Pedersen, D. S., Popp, T. J., Rasmussen, S. O., Severinghaus, J. P., Svensson, A., and Vinther, B. M.: Consistently dated records from the Greenland GRIP, GISP2 and NGRIP ice cores for the past 104 ka reveal regional millennial-scale δ18O gradients with possible Heinrich event imprint, Quat. Sci. Rev., 106, 29–46, https://doi.org/10.1016/j.quascirev.2014.10.032, 2014.
Severi, M., Becagli, S., Castellano, E., Morganti, A., Traversi, R., Udisti, R., Ruth, U., Fischer, H., Huybrechts, P., Wolff, E., Parrenin, F., Kaufmann, P., Lambert, F., and Steffensen, J. P.: Synchronisation of the EDML and EDC ice cores for the last 52 kyr by volcanic signature matching, Clim. Past, 3, 367–374, https://doi.org/10.5194/cp-3-367-2007, 2007.
Svensson, A., Bigler, M., Blunier, T., Clausen, H. B., Dahl-Jensen, D., Fischer, H., Fujita, S., Goto-Azuma, K., Johnsen, S. J., Kawamura, K., Kipfstuhl, S., Kohno, M., Parrenin, F., Popp, T., Rasmussen, S. O., Schwander, J., Seierstad, I., Severi, M., Steffensen, J. P., Udisti, R., Uemura, R., Vallelonga, P., Vinther, B. M., Wegner, A., Wilhelms, F., and Winstrup, M.: Direct linking of Greenland and Antarctic ice cores at the Toba eruption (74 ka BP), Clim. Past, 9, 749–766, https://doi.org/10.5194/cp-9-749-2013, 2013.
Thorsteinsson, T., Kipfstuhl, J., Eicken, H., Johnsen, S. J., and Fuhrer, K.: Crystal size variations in Eemian-age ice from the GRIP ice core, Central Greenland, Earth Planet. Sc. Lett., 131, 381–394, 1995.
Thorsteinsson, T., Kipfstuhl, J., and Miller, H.: Textures and fabrics in the GRIP ice core, J. Geophys. Res., 102, 26583–26599, 1997.
Tison, J.-L., de Angelis, M., Littot, G., Wolff, E., Fischer, H., Hansson, M., Bigler, M., Udisti, R., Wegner, A., Jouzel, J., Stenni, B., Johnsen, S., Masson-Delmotte, V., Landais, A., Lipenkov, V., Loulergue, L., Barnola, J.-M., Petit, J.-R., Delmonte, B., Dreyfus, G., Dahl-Jensen, D., Durand, G., Bereiter, B., Schilt, A., Spahni, R., Pol, K., Lorrain, R., Souchez, R., and Samyn, D.: Retrieving the paleoclimatic signal from the deeper part of the EPICA Dome C ice core, The Cryosphere, 9, 1633–1648, https://doi.org/10.5194/tc-9-1633-2015, 2015.
Traversi, R., Becagli, S., Castellano, E., Marino, F., Rugi, F., Severi, M., de Angelis, M., Fischer, H., Hansson, M., Stauffer, B., Steffensen, J. P., Bigler, M., and Udisti, R.: Sulfate spikes in the deep layers of EPICA-Dome C ice core: Evidence of glaciological artifacts, Environ. Sci. Technol., 43, 8737–8743, https://doi.org/10.1021/es901426y, 2009.
Wolff, E. W., Mulvaney, R., and Oates, K.: The location of impurities in Antarctic ice, Ann. Glaciol., 11, 194–197, 1988.
Wolff, E. W., Cook, E., Barnes, P. R. F., and Mulvaney, R.: Signal variability in replicate ice cores, J., Glaciol., 51, 462–468, 2005.
Wolff, E. W., Fisher, H., Fundel, F., Ruth, U., Twarloh, B., Littot, G. C., Mulvaney, R., Röthlisberger, R., de Angelis, M., Boutron, C. F., Hansson, M., Jonsell, U., Hutterli, M. A., Lambert, F., Kaufmann, P., Stauffer, B., Stocker, T. F., Steffensen, J. P., Bigler, M., Siggaard-Andersen, M. L., Udisti, R., Becagli, S., Castellano, E., Severi, M., Wagenbach, D., Barbante, C., Gabrielli, P., and Gaspari, V.: Southern Ocean sea-ice extent, productivity and iron flux over the past eight glacial cycles, Nature, 440, 491–496, https://doi.org/10.1038/nature04614, 2006.
Rempel et al. (2001) formulated their model with ϕ as a mass fraction (instead of volume fraction) and used a correspondingly different unit for cB, but the interactions are the same.