The Cryosphere, 19, 5693–5717, 2025 https://doi.org/10.5194/tc-19-5693-2025 © Author(s) 2025. This work is distributed under the Creative Commons Attribution 4.0 License. Doomed descent? How fast sulphate signals diffuse in the EPICA Dome C ice column Felix S. L. Ng1, Rachael H. Rhodes2, Tyler J. Fudge3, and Eric W. Wolff2 1School of Geography and Planning, University of Sheffield, Sheffield, UK 2Department of Earth Sciences, University of Cambridge, Cambridge, UK 3Department of Earth and Space Sciences, University of Washington, Seattle, WA, USA Correspondence: Felix S. L. Ng (f.ng@sheffield.ac.uk) Received: 2 April 2025 – Discussion started: 19 May 2025 Revised: 28 August 2025 – Accepted: 27 September 2025 – Published: 14 November 2025 Abstract. The loss of climate information due to smooth- ing of ionic impurity signals in ice provides a strong mo- tivation for understanding their diffusion rates at ice-core sites. By analysing sulphate signals in the EPICA Dome C (EDC) core, recent studies estimated the vertical profile of effective diffusivity Deff at that site. However, Deff crudely approximates the local diffusivity D in the ice, it being a nonuniform-weighted average of D over large intervals. We formulate the mathematical inversion for retrieving the D profile from observed signals, which reconciles the findings of the earlier studies as well as elucidating the averaging ap- proximation. Inversion for EDC sulphate reveals a rapid de- crease in D through the firn layer – from ≈ 10−6 m2 yr−1 at the surface to ≈ 1.7× 10−8 m2 yr−1 at the firn-ice tran- sition (≈ 100 m depth, ≈ 2.5 ka), followed by a gradual de- cline to ≈ 10−10 m2 yr−1 through 100–2700 m (2.5–390 ka). This profile enables new interpretation of sulphate transport in the EDC column. We propose vapour diffusion of H2SO4 through interconnecting air pores as the cause of the high firn diffusivity. By evaluating the mechanisms controlling D be- low the firn (diffusion through ice crystals, liquid veins and grain boundaries and diffusion arising from interfacial mo- tion), we infer a dominant partitioning of signals immedi- ately below the firn to a connected vein system, and progres- sive smoothing of vein signals by Gibbs–Thomson diffusion down to ≈ 2000 m depth, which leaves more and more of the remaining signals to grain boundaries. We conclude that those sulphate signals that survive the initial fast diffusion in the firn to “punch through” to its base might survive into deep ice, and that EDC sulphate preserves a strongly filtered history of volcanic and climatic forcing that underrepresents changes and events shorter than a few years. For the Beyond EPICA – Oldest Ice and Million Year Ice Core drilling sites on Little Dome C, calculations assuming a diffusivity pro- file like our EDC profile and not exceeding 10−10 m2 yr−1 in ice older than 450 ka constrain the sulphate diffusion length in ice 1–2 Ma old to 2 cm at most, and probably as low as ≈ 1 cm, for atmospheric-sourced signals that experienced only diffusion and mechanical shortening in the column. 1 Introduction Ionic impurities in ice cores provide valuable records of climate and environmental change (e.g. Legrand and Mayewski, 1997). The realisation that impurity signals in ice may be altered – not necessarily carrying climatic in- formation “written in stone” – motivates study of the post- depositional processes threatening their integrity. Diffusion attenuates and broadens signals as they descend the ice col- umn, potentially causing severe signal loss at depth, where the diffusion rate may be enhanced by higher temperature. The vertical pattern of the diffusion rate is of interest to ques- tions about the reliability of ice-core ion records, the amount of climatic information retrievable from their signals, and the methods of reconstructing past forcings at the ice-sheet sur- face – questions that matter the more as ice-coring campaigns seek older and older records, such as in the Beyond EPICA – Oldest Ice project (BE-OI, 2017) and the Million Year Ice Core project (MYIC, 2020). Recently, Fudge et al. (2024) and Rhodes et al. (2024) quantified diffusion on the high-resolution sulphate record Published by Copernicus Publications on behalf of the European Geosciences Union. 5694 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice of the EPICA (European Project for Ice Drilling in Antarc- tica) Dome C or “EDC” ice core (EPICA community mem- bers, 2004) from Antarctica. This record of sulphate concen- tration, measured by fast-ion chromatography (FIC) on bulk ice samples at ≈ 4 cm spacing down to 770 m depth and 1– 2 cm spacing at greater depths (Traversi et al., 2002, 2009), is shown in Fig. 1a. By analysing how its signals vary along the core, together with a signal-evolution model that accounts for diffusion and vertical mechanical shortening of the ice, Fudge et al. (2024) and Rhodes et al. (2024) estimated the “effective diffusivity” Deff of sulphate at EDC. Sulphate is relevant in the diffusion context because vol- canic events, which occur as sharp peaks on such records, provide data for synchronising ice-core timescales (e.g., Sev- eri et al., 2012; Svensson et al., 2020) and inferring the his- tory of volcanism – the record in Fig. 1a has been used to study eruption frequency back as far as 200 ka (Castellano et al., 2004; Lin et al., 2022; Wolff et al., 2023). Sulphate may also experience rapid transport in the liquid veins of polycrystalline ice, given the low eutectic temperature of sul- phuric acid (−73 °C) implies its likely dissolution in vein wa- ter located at grain triple junctions (Mulvaney et al., 1988; Wolff et al., 1988; Nye, 1989; Mader, 1992), and given the- oretical modelling which shows that ionic signals residing in a network of connected veins diffuse rapidly due to the Gibbs–Thomson effect (Ng, 2021). However, when studying impurity transport in ice, it is difficult to know how the bulk concentration of an ion partitions into contributions from dif- ferent impurity sites – the ice-crystal lattice, grain bound- aries, veins, and micro-inclusions; the mechanisms of impu- rity transfer between these sites also remain elusive (Barnes et al., 2003; Ng, 2021; Stoll et al., 2021). Thus, our under- standing of how signals on the bulk concentration evolve is incomplete. Because the model used by Fudge et al. (2024) and Rhodes et al. (2024) in their diffusivity inversions tracks sulphate bulk concentration without resolving the partition- ing, their effective diffusivity (Deff) estimates for the EDC site reflect the overall outcome of different grain-scale trans- port processes. Yet, for this reason, their estimates provide global constraints on how these processes operate. In this paper, we formulate a theory of diffusivity inversion that extends the methods of Fudge et al. (2024) and Rhodes et al. (2024), and which may be applied to other ions and to other ice cores besides EDC. Their studies referred to the “effective diffusivity” in part because of the caveat about im- purity partitioning, but more specifically because their inver- sions assumed constant diffusivity acting on each signal as it evolves. Accordingly, they recognisedDeff as some weighted average of the true diffusivity. The averaging process has not been made clear though. We show mathematically that their Deff estimates, owing to the averaging approximation, de- viate significantly from the true diffusivity D (unless noted otherwise, all diffusivities in this paper pertain to sulphate). We improve upon their results to obtain the vertical profile ofD in the EDC ice column, deriving new information about ionic impurity transport there. Notably, we discover high D values localised to the firn layer, whose cause is discussed towards the end. We also briefly consider what the findings mean for signal survivability at the sites of the BE-OI and MYIC projects. For convenience, we abbreviate Fudge et al. (2024) and Rhodes et al. (2024) as “F2024” and “R2024”, respectively, given how often they are referenced below. Figure 1 illustrates their diffusivity inversions. R2024’s approach utilised the decay of signal peak amplitude, whereas F2024 employed two approaches, one based on sig- nal peak widening and the other on the decay of signal vari- ability down core (Fig. 1b). In R2024, 537 sulphate peaks were identified in the record down to ≈ 2800 m depth (0– 450 ka). For each peak, the height of the corresponding orig- inal peak at deposition on the ice-sheet surface is recon- structed, by assuming that it held the same amount of sul- phate as the observed peak (after removing local background concentration due to non-volcanic sources of sulphate such as marine biogenic emissions) and that it was Gaussian- shaped, with a duration of 3 years at “full width at tenth maximum” (FWTM), as is typically found for the width of volcanic sulphate peaks in Antarctic snow; see R2024 for de- tailed justification. Then, using their model, which we give in Eq. (1) below, R2024 numerically simulated the evolution of the reconstructed peak forward in time, tuning the diffusivity in multiple model runs to match the observed peak’s height at its recorded age, to find Deff for the peak. We denote by DR their amplitude-based Deff estimate. In contrast, F2024 studied only signals in interglacial and glacial maximum periods (red and blue shading in Fig. 1a) and made separate inversions for these period types, to cater for the possibility of interglacial ice and glacial ice having different diffusivities. This is motivated by the idea that the different ice-column conditions (e.g. strain rate, mean crystal size) in these periods might affect impurity transport differ- ently. Their width-based inversion, which gauges each peak’s width by its “full width at half maximum” (FWHM), per- forms best-fit numerical simulations as R2024 did, but uses two peaks below the surface (Fig. 1b) rather than one peak and its reconstructed surface counterpart. Specifically, for in- terglacials and glacial maxima separately, they ran simula- tions to evolve a Gaussian signal with an initial width equal to the median width of observed peaks in the most recent period (either the Holocene or LGM) to match the median width of observed peaks in the earlier interglacials or glacial maxima, thus backing out Deff for the intervening intervals. We denote by DF1 their width-based Deff estimate. The other approach of F2024 uses a method pioneered by Barnes et al. (2003) for quantifying signal variations in terms of “mean absolute gradient” (explained in Sect. 2.3) to estimate Deff from the decrease of signal variability down core. Using the method, Barnes et al. (2003) had estimated Deff = 3.9± 0.8× 10−8 m2 yr−1 for the Holocene part (top 350 m) of the sulphate record in Fig. 1a. F2024 essentially applied the method to older parts of the core, focussing on The Cryosphere, 19, 5693–5717, 2025 https://doi.org/10.5194/tc-19-5693-2025 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice 5695 Figure 1. Approaches and results of the inversions for effective diffusivity Deff by Fudge et al. (2024) and Rhodes et al. (2024), for sulphate at the EPICA Dome C ice-core site. (a) Depth profile of sulphate concentration from fast ion chromatography (Traversi et al., 2009), showing abundant peaks, many of them recording volcanic eruptions. Interglacial and glacial maximum periods are highlighted by red and blue shading, respectively; for their age and depth ranges, see Table A1 of Fudge et al. (2024). Grey shading marks the record > 2800 m, which is not studied herein. (b) Schematic of the approaches of Rhodes et al. (2024) and Fudge et al. (2024) for finding their Deff estimates – DR, DF1, and DF2, which are based on peak-amplitude decay, peak widening and signal-variability reduction, respectively. (c) Plot of their DR, DF1, andDF2 results versus age back to 450 ka. The depth scale is indicated on the top axis. Horizontal bar shows the age range of eachDeff estimate, and vertical bar its uncertainty. Green point plots the Deff estimate of Barnes et al. (2003) for the Holocene part of the record. the sequence of interglacials and glacial maxima. We denote by DF2 their gradient-based Deff estimate (Fig. 1b). The effective diffusivities of R2024 and F2024 (Fig. 1c) show striking differences. Although DR, DF1 and DF2 in the deeper record≈ 200 to 450 ka (≈ 2100–2800 m) have similar magnitudes, ∼ 10−9–10−8 m2 yr−1, DR is much higher (up to 10−6 m2 yr−1) than DF1 and DF2 in ice. 50 ka, where it decays with age and depth. As R2024 reported, their median DR value for Holocene ice (0–10 ka), 2.4× 10−7 m2 yr−1, is nearly ten times the Deff estimate of Barnes et al. (2003) (green data point in Fig. 1c). Beyond its initial decay, DR averages at ≈ 10−8 m2 yr−1 in 50–200 ka, still about twice of DF1 and DF2. On seeing that DF1 and DF2 (≈ 5× 10−9 m2 yr−1) are not much higher than the self- diffusivity of ice (≈ 3× 10−10 to 3× 10−9 m2 yr−1 at −50 to −35 °C; Ramseier, 1967), F2024 inferred that the fast sig- nal diffusion in liquid veins modelled by Ng (2021) occurs only to a limited extent for sulphate in the upper ≈ 90 % of the ice column, and hence most sulphate there resides within ice crystals and at grain boundaries – not in the veins. On the other hand, R2024 interpreted the initial high (falling) DR values for significant (diminishing) diffusion of sulphate in interconnected veins in the top quarter of the ice column. Resolving these differences is imperative because the dif- fusivity profile is key to understanding how the crystal-scale diffusion mechanisms vary with depth and the factors in- volved, such as impurity partitioning. Besides adopting dif- ferent inversion approaches, R2024 and F2024 processed the FIC data differently. R2024 only analysed sulphate peaks that are certainly volcanic by omitting others coincident with dust peaks, whereas F2024 applied the scaling procedure of Barnes et al. (2003) to the sulphate record to reduce the influ- ence of background climate variations before extracting sig- nals for analysis. These methodological differences can only explain minor discrepancies, not the overall incompatibility, between DR and DF1,2. The results in Fig. 1c also raise in- triguing questions, notably the cause of the near-surface de- cay in DR in ≈ 0–50 ka, which seems to continue through ≈ 100–450 ka at lower rate, and why (as both their studies pointed out)Deff does not increase with depth, against the ex- pectation that molecular diffusivity increases with tempera- ture. The ice temperature at the EDC site increases monoton- ically from ≈−53 °C at the surface to ≈−12 °C at 2800 m (Fig. S1 in the Supplement). Herein, our theory not only allows estimating the true dif- fusivity D, which is a more fundamental quantity than Deff for probing impurity transport mechanisms; it also shows how the DR, DF1 and DF2 estimates may be reconciled on account of their underlying averaging and two needed cor- rections in theDF2 inversion. A key insight is that the signal- evolution model of F2024 and R2024 can be solved analyt- ically, so the inversions can be done without numerical sim- https://doi.org/10.5194/tc-19-5693-2025 The Cryosphere, 19, 5693–5717, 2025 5696 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice ulation. While our inversion results draw interest to the firn diffusivity, their signal-evolution model ignores firn densifi- cation; we therefore also examine its validity when used for inversions within the firn. We focus on the EDC record in 0–2800 m (Fig. 1) by using the data collected by R2024 and F2024 without reprocess- ing the FIC sulphate concentrations. The record at depths > 2800 m (which features in part of F2024’s study) is ex- cluded for the reason given by R2024: there, some sulphate peaks may be non-volcanic and shaped by post-depositional processes other than diffusion and vertical mechanical short- ening. This is shown by the presence of (i) anomalous peaks below 2800 m depth that have been chemically modified, as evidenced by ion association (Traversi et al., 2009), and (ii) other anomalous peaks starting from ≈ 2700 m (perhaps as shallow as 2500 m) that exhibit side troughs, indicating sulphate being “sucked” from neighbouring background lev- els towards zones with high cation concentration to form the peaks (Wolff et al., 2023). These artefacts reflect added com- plexity in the evolution of signals in deep ice at EDC that makes their origin uncertain. Our theory and analyses strictly concern signals without such artefacts, which give the ideal input data for diffusivity inversion. While R2024’s data mit- igate the issue by excluding potential artefact peaks during data collection, the deepest data of F2024 used by us may contain artefact signals, especially anomalous peaks of type (ii); but, for reasons explained later, this should not affect our conclusions. 2 Mathematical theory 2.1 Signal evolution We begin with the advection–diffusion equation for signal evolution down the ice column, used by F2024 and R2024. In a coordinate frame moving with the ice, where z denotes dis- tance below a material horizon descending towards the bed, signals in the bulk impurity concentration C(z, t) (measured in µg kg−1, or µg L−1 of meltwater) evolve according to ∂C ∂t =D(t) ∂2C ∂z2 − ε̇z(t)z ∂C ∂z . (1) Here, t is the age of the horizon,D is the impurity diffusivity, and ε̇z(< 0) is local vertical strain rate. Equation (1) encap- sulates the effects of mechanical shortening and diffusional spreading. Table A1 lists other mathematical symbols used in the paper. Following F2024 and R2024, we use Eq. (1) to model sul- phate signals, assuming an invariant strain-rate profile and constant surface accumulation rate at the core site – thus, a steady-state column with constant thickness and vertical ve- locity profile. In this system, signals travel through fields that are functions of depth in the column only, not time, so the age–depth scale allows translation between D(t) and its ver- tical profile. Material at age t has shortened from its original thickness at the surface by the thinning factor S, given by S(t)= exp t∫ 0 ε̇z(η) dη, (2) where η denotes the variable of integration. Differentiating Eq. (2) gives dS/dt = ε̇zS. The thinning function S decays with age t from its value at the surface, S0 = S(t = 0)= 1. The inversion methods of R2024 and F2024 (elaborated in Sect. 2.2 and 2.3) use Eq. (1) as the basis, but as noted earlier, assume a constant D for each signal as it evolves down col- umn. The resulting effective diffusivities DR, DF1 and DF2 do not strictly represent the true (local) diffusivityD, instead averages measuring its cumulative effect over finite age and depth intervals; as we shall see, these intervals are large. By solving Eq. (1) analytically below, we develop exact inver- sions for D(t) that circumvent this assumption, at the same time deriving equations linking D(t) to DR, DF1 and DF2. How Eq. (1) is affected by firn densification will be exam- ined in Sect. 3.5, after we glimpse high firn diffusivity from our inversions. 2.2 Theory: peak-based inversions 2.2.1 The inversion possibility A key property we exploit is that a Gaussian signal stays Gaussian under the combined mechanical shortening and dif- fusional spreading described by Eq. (1). F2024 and R2024 both initialised their simulations with Gaussian peaks, but did not harness this property. As alluded to by R2024, the sul- phate flux from eruptions reaching the ice sheet often varies asymmetrically in time, but the deposited peaks rapidly relax to near-Gaussian. This motivates a Gaussian approximation to their shape. To see the property, define the transformed depth ζ = z S(t) (3) and define the variable τ(t)= t∫ 0 D(η) S2(η) dη+ τ0, (4) where τ0 is the value of τ at zero age. Here, ζ is the de- strained or unthinned thickness, and τ , an indirect proxy of age or time, accounts for the histories of diffusion and layer thinning. On letting C(z, t)= f (ζ,τ ), these changes of vari- able convert Eq. (1) to the classical heat equation ∂f ∂τ = ∂2f ∂ζ 2 , (5) which has the well-known (Gaussian) similarity solution The Cryosphere, 19, 5693–5717, 2025 https://doi.org/10.5194/tc-19-5693-2025 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice 5697 f (ζ,τ )∝ 1 √ τ e−ζ 2/4τ . (6) In the ζ -direction, this Gaussian’s width expressed as a stan- dard deviation is σ = (2τ)1/2. Equations (5) and (6) mean that in τ–ζ space, signals experience uniform diffusion at unit rate, and a Gaussian peak decays in amplitude following the factor 1/ √ τ and widens following √ τ . Consequently, observations of peak widening or amplitude reduction down core, which provide data on τ(t), can be used to recoverD(t) via Eq. (4). This idea forms the basis of the peak-based in- versions. For example, consider an amplitude-based inversion, where the “relative peak amplitude” (the ratio of a peak’s observed amplitude to its original amplitude on deposition at the surface at t = 0) has been compiled for different peaks along the core, as done by R2024. Suppose the relative am- plitudes vary with age to trace out the function α(t). Then we have α(t)= √ τ0/τ according to Eq. (6), and differentiating Eq. (4) with respect to t gives the inversion D(t)= S2(t) dτ dt = S2 d dt ( τ0 α2 ) =− 2τ0S 2α′ α3 (7) (the ′ denotes derivative). This inversion requires the value τ0 = τ(t = 0). For each observed peak, R2024 reconstructed the amplitude of the original peak by assuming it to be Gaussian, with a 3-year duration at FWTM and carrying the same total impurity load as the observed peak. They ig- nored firn densification effects and used the ice density in the reconstruction. We therefore set τ0 = σ 2/2, with σ = 3 years× a/4.2919, where σ (m) is the standard deviation mentioned above, a is the ice-equivalent accumulation rate (m yr−1), and the factor 4.2919 converts the FWTM of the Gaussian to σ . Positive τ0 ensures a finite amplitude for the initial peak in Eq. (6). The differentiation in Eq. (7) assumes α to be smoothly varying; in practice, one fits a curve to the α-data prior to inversion. Similarly, in a width-based inversion, where data on “rel- ative peak width” (the ratio of observed width to origi- nal width) trace out the function β(t), such that β(t)=√ τ(t)S2(t)/τ0S 2 0 = S(t) √ τ(t)/τ0, we derive D(t)= S2(t) dτ dt = S2 d dt ( τ0β 2 S2 ) = 2τ0β(β ′ −βε̇z). (8) Applying this inversion necessitates an assumption for the original peak width (for compiling β and τ0). However, for comparing against the width-based inversion of F2024, Eq. (8) first needs to be adapted for use on two peaks be- low the surface, rather than one at the surface and one below (Fig. 1b). We attain the relevant result via a different route below. 2.2.2 Full-fledged theory Before applying the above theory to data, we expand the mathematical analysis to unravel how peak-based inversions work and establish the relationship between D(t) and the ef- fective diffusivities of F2024 and R2024. The ratio D/S2 recurring in Eqs. (4), (7) and (8) relates a physical effect. As a signal shortens mechanically, its vari- ations steepen, so it diffuses faster than if shortening were absent. With S < 1 below the surface, D/S2 represents the amplified diffusivity. Another way of picturing this effect is to imagine the signal experiencing the diffusivityD, but over a longer time – longer by 1/S2 times. This motivates us to in- troduce another age variable, ψ , defined by dψ = dt/S2(t). Specifically, we set ψ = 0 at t = 0, so that ψ(t)= t∫ 0 S−2(η)dη. (9) This function has unit slope at t = 0 (since S(t = 0)= 1) and curves upward (e.g. Fig. 3c). We call ψ the dilated age be- cause it accounts for thinning but excludes diffusion, unlike the proxy variable τ , which accounts for both. On moving from t–z toψ–ζ space (Fig. 2), the transforma- tion z to ζ geometrically destrains the signal to track material horizons, whereas the transformation t to ψ stretches time to capture the mechanically-induced enhanced diffusion on the signal. With coordinate stretching absorbing both effects, the transformed signal obeys ∂C/∂ψ =D(ψ)∂2C/∂ζ 2 with- out a shortening term (Fig. 2b). Crucially, under the move, Eq. (4) is converted to τ(ψ)= ψ∫ 0 D(ψ) dψ + τ0, (10) which shows that inversion for D fundamentally involves D(ψ)= dτ dψ . (11) In other words, as a Gaussian peak evolves in ψ–ζ space, its unthinned width squared and its inverse squared amplitude (recall that unthinned width ∝ √ τ and amplitude ∝ 1/ √ τ) increase at a rate with respect to ψ that equals the instanta- neous or local diffusivity. Equivalently, the local diffusivity is given by the rate of change of these peak-form parameters with respect to dilated age ψ (Fig. 2b). Not surprisingly, the age-domain inversions in Eqs. (7) and (8) also involve rates of change. Given these insights, we can calculate the effective diffu- sivitiesDR andDF1 of R2024 and F2024 analytically, which obviates need to integrate Eq. (1) numerically and perform multiple simulations to fit data. As noted before, their inver- sions assumed constant D during each signal’s descent. If https://doi.org/10.5194/tc-19-5693-2025 The Cryosphere, 19, 5693–5717, 2025 5698 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice Figure 2. Evolution of a Gaussian signal in (a) t–z space and (b)ψ– ζ space. Solid black curves signify material trajectories. Dashed curve in (b) marks the unthinned signal width, whose square in- creases at a rate with respect to ψ that reflects the instantaneous diffusivity. Eq. (10) is used to reproduce their inversions, then we set D ≡Deff (constant) in its integral, which gives τ(ψ)=Deffψ + τ0, (12) or Deff = τ − τ0 ψ = τ0 ψ ( τ τ0 − 1 ) , (13) which describes the inversion based on a subsurface peak and its reconstructed original (surface) peak. Where the inversion uses a pair of subsurface points, say, τ1 at dilated age ψ1 and τ2 at dilated age ψ2, differencing the application of Eq. (12) to these data yields Deff = τ2− τ1 ψ2−ψ1 . (14) From these results, it follows that the R2024 inversion is equivalent to DR(t)= τ0 ψ(t) ( 1 α2(t) − 1 ) , (15) with ψ given by Eq. (9) and the data for α and τ0 gathered as before (Sect. 2.2.1), whereas the width-based inversion of F2024 for two peaks of age t1 and t2 has the analytical coun- terpart DF1 = τ(t2)− τ(t1) ψ(t2)−ψ(t1) , (16) in which the τ values derive from observed peak widths. F2024 measured the unthinned FWHM of each peak, so τ = σ ∗2/2, where σ ∗ = FWHM/2.3548 is the destrained standard deviation of the Gaussian. The effective diffusivity from Eq. (15) or (16) is valid for the specific interval bracketed by the paired data, as in R2024 and F2024’s simulation-based inversions. The inter- val in R2024’s inversion spans each peak’s entire history. In F2024’s inversion, which uses paired data between the Holocene and earlier interglacials or between the LGM and earlier glacial maxima, the intervals exceed ∼ 100 kyr. Thus, DR of R2024 andDF1 of F2024 are effective diffusivity esti- mates for different periods – this is a key reason behind their discrepancy, which we will point out again when analysing results in Sect. 3. Next we relate the effective diffusivities to the true diffu- sivity D(t). Applying Eq. (10) to paired data (ψ1, τ1) and (ψ2, τ2), eliminating τ0, and using Eq. (14), yields Deff = 1 ψ2−ψ1 ψ2∫ ψ1 D(ψ) dψ, (17) or Deff(ψ)= 1 ψ ψ∫ 0 D(ψ) dψ (18) if the upper data point lies at the surface. These results show that Deff is the interval average of D, not over t but over the dilated age ψ . Since ψ(t) curves upward (Fig. 3c), Deff is biased towards D in the older part of the averaging interval; but it is influenced byD in the younger part. The larger is the interval, the more crudely Deff approximates D at the lower (deeper) data point. Equation (18) leads to further insights on the profileDR(t) retrieved by the R2024 inversion. By evaluating its integral, working with time rather than ψ as the integration variable, we derive DR(t)= 1 ψ(t) t∫ 0 D(η)ψ ′(η) dη =D(t)− 1 ψ(t) t∫ 0 D′(η)ψ(η) dη. (19) According to this expression, DR found from an observed peak not only reflects the local diffusivity D at its depth, but also inherits a signal from the variations in D throughout its earlier shallower history: we call this the “memory effect”. Notably, DR(t) overestimates (underestimates) D(t) if D(t) is a decreasing (increasing) function. Differentiating the first equation in Eq. (19) gives the op- posite conversion from DR to D, D(t)= 1 ψ ′ d dt [ψDR(t)] =DR+D ′ RψS 2, (20) which shows that D is less (greater) than DR wherever DR decreases (increases) down core. R2024 recognised DR as a “time-weighted diffusivity” and took care when interpret- ing their DR(t) profile; but without the analytical result in The Cryosphere, 19, 5693–5717, 2025 https://doi.org/10.5194/tc-19-5693-2025 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice 5699 Figure 3. Functions used in our diffusivity inversions for the EPICA Dome C ice-core site: (a) the AICC2012 age–depth scale (Bazin et al., 2013; Veres et al., 2013) and the corresponding (b) depth profile of thinning factor S, (c) dilated age ψ versus age t , and (d) ψS2 versus age. Black curves derive directly from the AICC2012 scale. Magenta curves, used in our inversions, are smooth approximations based on an ice-flow model assuming the submergence velocity wi = as(h/H) 1.2, where h is height above the bed, H = 3165 m (mean ice thickness chosen by Rhodes et al., 2024), and surface accumulation rate as = 0.0195 m yr−1. Grey triangle in (c) illustrates the misamplification factor (Sect. 2.3). Eq. (20), inferring D from DR is challenging. The present analysis also reveals the weighting to be highly nonlinear. In Sect. 3, we use Eq. (20) to estimate D(t) from DR(t) and Eq. (19) to predict DR(t) from D(t), discovering a marked difference between these curves. 2.3 Theory: gradient-based inversion We turn to F2024’s inversion for the effective diffusivityDF2, which calculates the “mean absolute gradient” m̄ of signals with the Barnes et al. (2003) method, which in turn is based on the diffusion-length theory of Johnsen (1977). We extend this framework to derive an exact inversion for D from m̄, exposing the averaging approximation behind DF2. We find that the Barnes et al. (2003) method – and thus the DF2 esti- mates – require two corrections. In the Barnes et al. (2003) method, the concentration record C is first destrained and processed to suppress un- wanted signals from background climate variations. Signal peaks on the processed record, Cp, are thought to reflect the sulphate input from volcanic events more reliably (with less bias) than C. To quantify the signal variability on Cp, they studied different 10 m long sections down core by calculat- ing their signal mean absolute gradient, m̄= 1 n1ζ n∑ i=1 ∣∣Cp,i+1−Cp,i ∣∣ , (21) where Cp,i denotes individual processed concentration mea- surements, 1ζ is the destrained interval between measure- ments, and n is the number of intervals in each 10 m. Equa- tion (21) is the same as their Eq. (1), despite written with different symbols. Their method quantifies the rate of diffusive smoothing by using the observed decrease in m̄ down core (Fig. 1b) and retrieves Deff from the rate. Notably, they regard the prin- cipal signals on Cp as periodic, with a wavenumber k∗ that does not vary with depth on the destrained record (we use ∗ to signify destraining). Accordingly, the ratio of m̄ of a core section at depth to the mean absolute gradient m̄0 of a ref- erence section higher in the column measures the amplitude decay of the signals, and they equate this ratio to the signal attenuation predicted for Eq. (1) by Johnsen (1977) – thus, m̄ m̄0 = exp(−k∗2σ ∗2/2), (22) in which σ ∗ is the destrained value of the diffusion length σ ; that is, σ ∗ = σ/S(t). In Johnsen’s theory, the diffusion length σ evolves accord- ing to the ordinary differential equation dσ 2 dt = 2D(t)+ 2ε̇z(t)σ 2, (23) and transforming this to the destrained coordinate system yields dσ ∗2/dt = 2D(t)/S2(t). (24) However, Barnes et al. (2003) took dσ ∗2/dt = 2D without the final 1/S2, assuming Eq. (23) with ε̇z set to zero to be a valid diffusion-length equation for unthinned records. On taking a constant (effective) diffusivity, they then found σ ∗2 = 2Defft , which, together with Eq. (22), led them to the inversion formula Deff =− 1 k∗2t ln ( m̄ m̄0 ) . (25) When rewritten for a pair of subsurface data points, this gives the F2024 inversion formula: DF2 =− 1 k∗2(t2− t1) ln ( m̄2 m̄1 ) . (26) https://doi.org/10.5194/tc-19-5693-2025 The Cryosphere, 19, 5693–5717, 2025 5700 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice These formulas are approximate because of the missing 1/S2 in the underlying diffusion-length model: strictly, Eq. (24) should be used instead. In particular, for the EDC core site, the approximation is reasonable for signals in t . 104 years (because S ≈ 1 up to that age; Fig. 3a, b) but not beyond. It follows that the Deff estimate of Barnes et al. (2003) for the Holocene ice (Fig. 1c) is approximately valid, but the DF2 estimates of F2024 for older sections of ice suffer large inac- curacies. Having explained the Barnes et al. method, we modify it to derive an exact inversion for D. In the ψ–ζ coordinate system, Johnsen’s diffusion-length equation (Eqs. 23 or 24) takes the form1 dσ ∗2 dψ = 2D(ψ), (27) and substituting for σ ∗2 from Eq. (22) gives D(ψ)=− 1 k∗2 d(lnm̄) dψ . (28) This inversion formula involves a rate of change, as in the peak-based inversions, and shows that D can be estimated from the slope of the logarithmic plot of m̄ versus dilated age ψ (we will explore this with EDC data in Sect. 3.3). One can again relate the effective diffusivity to D. Sup- pose D in Eq. (27) equals a constant, DE; this is the effec- tive diffusivity that is found by the inversion without the ap- proximation of Barnes et al. (2003) described above. Then σ ∗2 = 2DEψ . Using this together with Eq. (22) for paired data leads to the inversion formula DE =− 1 k∗2(ψ2−ψ1) ln ( m̄2 m̄1 ) ( ≡ 1 ψ2−ψ1 ψ2∫ ψ1 D(ψ)dψ ) . (29) We see that DE is the average of D over the dilated age ψ , as for DR and DF1. On comparing DE against the effective diffusivities in Eqs. (25) and (26), we find DF2 = ψ2−ψ1 t2− t1 DE, (30) which means that DF2 of F2024 (also Deff of Barnes et al. 2003) is misamplified by (ψ2−ψ1)/(t2− t1) and strongly overestimates the effective diffusivity in deep intervals (see 1Eq. (27) can also be derived from the τ–ζ formulation (Sect. 2.2.1). It is well known that given the Gaussian solution of Eq. (5), its general solution can be written as the convolution inte- gral f (ζ,τ )= 1 σ ∗ √ 2π ∫ ∞ −∞ F(η)e−(ζ−η) 2/2σ ∗2 dη, where F is the initial condition (e.g. Johnsen, 1977). Substituting this into Eq. (5) yields dσ∗2/dτ = 2, which, after a change of variable from τ to ψ , gives Eq. (27). triangle in Fig. 3c). The issue stems from the missing 1/S2. The misamplification ratio allows the effective diffusivities of Barnes et al. (2003) and F2024 to be corrected to give the desired value, DE. The other correction in the mean absolute gradient ap- proach concerns the signal wavenumber k∗. Barnes et al. (2003) and F2024 estimated it via k∗ = 2π/w̄, where the mean destrained wavelength w̄ of the signals is found by cal- culating w̄ = 4 m̄L n∑ i=1 ∣∣Cp,i − C̄p ∣∣1ζ (31) for a long record (F2024 used the Holocene or LGM part of the EDC record for this); L is the length of the record and C̄p its mean impurity concentration. Barnes et al. (2003) idealised the signals as triangular-shaped when deriving Eq. (31), but our repeat derivation in Appendix B shows that its right-hand side should be doubled, or π /2 times larger if one assumes sinusoidal signals. Consequently, their method overestimates k∗ by ≈ 1.6–2 times, and the DF2 estimates of F2024 and the Deff estimate of Barnes et al. (2003) for Holocene ice (3.9± 0.8× 10−8 m2 yr−1) are too small by a factor of k∗2 ≈ 2.5–4 times. In our DF2 inversions below, we remedy both issues by correcting the results with this factor and the misamplification ratio. 3 Diffusivity inversions: results and analysis We proceed to estimate the true diffusivity profileD(t) at the EDC site, using the theory of Sect. 2 and data from F2024 and R2024 as input. The work is done in stages. In Sect. 3.1, 3.2, and 3.3, we undertake inversions from peak amplitude, peak width, and signal gradient in turn, exploring avenues including the conversion of Deff to D and direct inversion of data, as well as finding the effective diffusivities DR, DF1 and DF2 with analytical formulas. While these sections al- low gleaning information about D(t), we further constrain its form by forward modelling in Sect. 3.4. In Sect. 3.5, after inferring high sulphate diffusivity confined to the firn layer, we examine how firn densification impacts the diffusivity in- version. F2024 and R2024 used the AICC2012 chronology of the EDC site (Fig. 3a; Bazin et al., 2013; Veres et al., 2013) throughout data compilation and analyses. To max- imise compatibility of our results with theirs, we employ the same chronology, rather than the newer AICC2023 chronol- ogy (Bouchet et al., 2023). In particular, our inversions use what we call “AICC2012-based” functions – smoothed forms of the thinning factor S and dilated age ψ(t) (Fig. 3, ma- genta curves), which we derive from a power-law model of the ice submergence velocity in the EDC column fitted to the AICC2012 age–depth scale; see the caption of Fig. 3 for the details. Although the smoothing injects minor differences be- tween our Deff estimates and those of R2024 and F2024, it The Cryosphere, 19, 5693–5717, 2025 https://doi.org/10.5194/tc-19-5693-2025 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice 5701 is desirable because the thinning function provided with the AICC2012 dataset is non-monotonic (black curve, Fig. 3b), with small bumps that imply negative strain rate at various depths. 3.1 Inversions from peak-amplitude decay Figure 4a shows the relative amplitudes α of the peaks stud- ied by R2024, obtained from their Supplementary data by dividing observed peak heights by original peak heights. The values show considerable scatter but generally decay with age. To compute the effective diffusivity DR for each peak, we apply Eq. (15) to its α-value, setting τ0 as described in Sect. 2.2.1. When calculating the strain rate ε̇z for simulating Eq. (1), R2024 adopted an ice submergence velocity profile derived not from the AICC2012 scale, instead from a Nye model with an ice thickness of 3165 m and a surface accu- mulation rate that puts the peak at its observed depth, so that its age and depth agree with the AICC2012 scale. Thus, their inversion ofDR envisages a slightly different steady-state ice column for each peak. Their use of the Nye model, which does not resolve the details of firn compaction near the top of the column, seems consistent with their choice of work- ing with ice-equivalent depths when compiling original peak amplitudes (Sect. 2.2.1). We shall say more about the effect of the firn processes in Sect. 3.5. To show that our analytic approach can reproduce the DR estimates of R2024, in our first use of Eq. (15) we adopt their ice-flow approximation by using ψ(t) based on their peak-specific Nye model instead of our AICC2012- based model. The correspondingDR results (Fig. 4b, circles) agree closely with their estimates (Fig. 1c), decaying from ∼ 10−6 to ∼ 10−9 m2 yr−1, rapidly in ≈ 0–50 ka and slowly beyond. We find four values exceeding 10−6 m2 yr−1 near t = 0 and four values below 10−9 m2 yr−1 not reported by R2024 (Fig. 4c and b). These and other minor discrepan- cies between our results arise because Eq. (15) is an exact formula for DR, whereas their DR estimates are constrained to 50 graded values (visible from the banding in Fig. 4c) on the log scale between 10−9 and 10−6 m2 yr−1 (values outside these bounds are clipped to them). Performing the same inversion with our AICC2012-based function ψ(t) (Fig. 3c) yields lowerDR estimates, especially for deeper peaks (Fig. 4b and d, magenta points). This is be- cause their Nye model tends to overestimate S (underesti- mate the amount of thinning) at depth; less of the thinning- induced enhancement in signal diffusion (Sect. 2.2.2) is cap- tured, making their DR estimates larger. The lowering helps explain some of the difference betweenDR and F2024’sDeff estimates. Next, we attempt to estimate the true diffusivity profile D(t) by two approaches. The first applies the time-domain inversion in Eq. (7) (same as Eq. (11) in theψ–ζ domain) to a smoothed version of α, which we derive by fitting the α-data with the sum of two exponentials (solid curve, Fig. 4a). This function is preferred to a single exponential (dashed curve, Fig. 4a) because it captures the high α-values near t = 0 bet- ter. In Eq. (7) we use the AICC2012-based thinning function S (Fig. 3b). The second approach converts DR to D with Eq. (20), assuming the AICC2012-based function ψS2 (Fig. 3d). To derive a smooth input for Eq. (20), in which the derivative D′R appears, we spline-fit our AICC2012-basedDR estimates from Fig. 4b on log-10 scale. These estimates show pro- nounced fluctuations and scatter on time scales shorter than ≈ 20 kyr that indicate uncertainty and noise on the relative amplitudes α, so we choose a level of spline smoothing to suppress these fluctuations; see Spline 2 in Fig. 5b, e, h. How- ever, the exact time scales on which fluctuations inDR reflect true changes in diffusivity is unknown, so we experiment also with less and more smoothing by using Spline 1 and Spline 3 (the left-hand and right-hand columns of panels in Fig. 5). Spline 3 strongly suppresses fluctuations in DR shorter than about 50 kyr. The curves of D(t) computed by this second approach (Fig. 5, solid black curves) indicate high, steeply-decaying diffusivity in the first ≈ 2 to 8 ka – from ∼ 10−6 m2 yr−1 to well below 10−7 m2 yr−1, followed by generally low diffu- sivity beyond (∼ 10−8 m2 yr−1) and even negative diffusivity in some age ranges (see comments below). Although stronger spline smoothing lengthens the initial fast decay, all three curves portrayD as greatly diminished from its surface value by several ka (Fig. 5d–i). Because the firn-ice transition at EDC lies at≈ 100 m depth (e.g. Landais et al., 2006; Calonne et al., 2022), where t ≈ 2.5 ka, much of the initial steep drop inD apparently occurs in the firn layer; we explore the cause of this later in Sect. 4.1. In contrast,D(t) retrieved by the first approach (Fig. 5, dashed black curve in all panels) shows a much more subdued decay over the first 30 ka, starting from a lower surface diffusivity≈ 6× 10−8 m2 yr−1. We think that this is because the two-term exponential in Fig. 4a does not adequately capture the high negative slope of the α-data near t = 0, which is necessary for Eq. (7) to reconstruct the details of D there. In the second approach, the curves of D(t) lie below the DR estimates in many places. As expected from our theory in Sect. 2.2.2, D(t) lies above (below) the spline curve of DR where this curve rises (drops). Where DR increases with age, high D (>DR) is retrieved because peaks with much lower amplitude than overlying peaks in the column imply high diffusivity in the intervening depth interval. Where DR decreases with age, lowD ( tc. The inclined floor is assumed to have a slope equal to the mean slope of the fully-corrected median DF2 estimates (Fig. 7b), but its level is fixed by the corner location, not by the es- timates. For either floor type, we find the combination of Dc and tc that best-fits the predicted DR(t) profile to our AICC2012-based DR estimates (Fig. 4b, magenta points). The DF1 results (Fig. 6c) are excluded from the exercise, given their large uncertainties. Figure 8 shows the best-fit profiles and maps of mis- fit over the tc–Dc parameter space from the forward mod- elling. In the flat-floor experiment (Fig. 8a, c), the pre- dicted DR(t) profile fits the DR estimates moderately well, and D decreases from its surface value to the corner dif- fusivity Dc ≈ 2.1× 10−9 m2 yr−1 in 9.1 ka. In the inclined- floor experiment (Fig. 8b, d), DR(t) fits the DR esti- mates better, capturing their gentle decay trend at large t . Here, D(t) has a shorter initial drop (2.8 ka), Dc is higher (≈ 1.74× 10−8 m2 yr−1), and the floor shoots through the fully-corrected DF2 values even though their level is not a fitting target; D(t) also lies slightly below their trend, as an- ticipated. Thus, this D(t) profile (Fig. 8b) explains the mean absolute gradient data as well as the peak-amplitude data and yields the better reconstruction of the two experiments. It also gives a more plausible estimate of the true diffusiv- ity than the opposite conversion from DR to D (Sect. 3.1), which reconstructed negative D intervals. Its steep initial drop is mainly constrained by the DR decay in 0–50 ka, and its inclined-floor level by the deeper DR values. Importantly, the corner age (2.8 ka) confirms our finding from Sect. 3.1 that the high D decaying through the upper column does not extend far below the firn-ice transition. Note that D(t) in Fig. 8b is reliable to a maximum age of only ≈ 390 ka https://doi.org/10.5194/tc-19-5693-2025 The Cryosphere, 19, 5693–5717, 2025 5706 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice Figure 8. Forward modelling of the effective diffusivity profile DR(t) from the true diffusivity profile D(t), and least-squares fitting to constrain D(t). As described in Sect. 3.4, panels (a) and (c) report an experiment assuming a flat floor for D; and panels (b) and (d), an inclined floor. (a, b) Plot of log diffusivity versus age, showing D(t) (dashed), the predicted DR(t) profile (solid curve), and DR data from peak-amplitude inversion (points). Panel (b) includes the fully-corrected DF2 results of Fig. 7b for comparison. Grey shading about D(t) shows its maximal variation as found from the confidence intervals of the best-fit parameters, tc and Dc. (c, d) Root mean square (RMS) mismatch in log-10 scale between predicted and estimatedDR values, as a function of the corner age tc and corner diffusivityDc of theD(t) profile. White crosses locate the optimal tc and Dc values in (a) and (b), which yield RMS mismatches of 7.16 and 6.70, respectively. (≈ 2700 m) because the deepest DR and fully-corrected DF2 data constraining the fit are interval-based results. In these experiments, the long tails on DR caused by the high initialD values confirm the memory effect, and theD(t) profiles are broadly consistent with the Holocene Deff es- timate of Barnes et al. (2003), ≈ 1.3± 0.6× 10−7 m2 yr−1 after applying the k∗-correction (Sect. 2.3). That D(t) in Fig. 8b lies below the DR, DF1 and DF2 estimates of R2024 and F2024 by up to 1–2 orders of magnitude confirms that these effective diffusivities approximate D crudely, and that the discrepancy between them stems from the underlying av- eraging, the different intervals used to evaluate them, and er- rors in the mean absolute gradient method. Finally, the DR estimates being fitted depend on the assumed 3-year FWTM duration for the surface peaks (Sect. 2.1). To gauge the impact of this assumption, we conducted an ensemble of 105 best-fit forward model runs, where each of the 537 DR estimates serving as fitting target in each run was picked randomly from its three values based on FWTMs of 1, 3, and 5 years. The maximal ranges found for tc and Dc are 2.3–2.9 ka and 1.70–1.78× 10−8 m2 yr−1, respectively, narrowly bracketing the results in Fig. 8b. This is not surprising as DR is weakly sensitive to the FWTM, as found by R2024. 3.5 Diffusivity inversion in the firn The preceding inversions highlight the firn diffusivity as a key interest, but their methods ignore firn densification and assume an EDC column consisting of ice only. Might the high D found for the firn be an artefact of this neglect? Should the methods be adjusted to account for firn density change? The DR inversion is especially relevant in this re- gard, as it involves signals descending through the firn layer; its results are also used in the estimation of D(t) (Sect. 3.1 and 3.4). Here we show that because of the way the bulk con- The Cryosphere, 19, 5693–5717, 2025 https://doi.org/10.5194/tc-19-5693-2025 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice 5707 centration C is defined and used, Eq. (1) correctly describes signal evolution in the firn as well as in the ice, so that the DR inversion and our findings for D are valid. The concerns are two-fold. R2024’s reconstruction of the original surface peaks, which provides input data for the DR inversion and takes each peak’s width (FWTM) to be 3a, uses the ice-equivalent thickness and assumes surface ma- terial with the ice density ρi (917 kg m−3) throughout cal- culation (Sect. 2.2.1). If the firn surface density ρ0 (<ρi) is used in the reconstruction, each original peak would be wider (3aρi/ρ0), its height proportionally less, so α may have been underestimated and DR overestimated. A second concern is that Eq. (1) might not conserve the amount of impurity in densifying firn. With D defined as the diffusivity of the bulk material (ice–air composite in the case of firn), D∂2C/∂z2 in Eq. (1) describes impurity flux divergence only if C is the impurity amount per unit volume of bulk material, not if C is impurity amount per unit mass – as used by us and R2024 in theDR inversions. Consequently, one fears that Eq. (1) might be the wrong model, and it is unclear whetherD presently re- trieved for the firn describes its bulk diffusivity or some other quantity. We dispel these concerns in the following by deriving Eq. (1) from first principles. Consider the firn layer in the Cartesian coordinates (X, Y , Z), with depth Z measuring down from the surface. Let U = (U , V ,W) be the firn veloc- ity, and ρ = ρ(Z) be the firn density profile, assumed time- invariant. We define the bulk concentration cB as the impurity amount per volume, i.e., cB = ρC. Then, impurity conserva- tion obeys ∂cB ∂t +∇.(UcB)=∇.(D∇cB), (32) in which D =D(Z) describes the vertical diffusivity profile, and the equation for water mass conservation is ∂ρ ∂t +∇.(ρU)= 0. (33) At ice-sheet divide or summit locations, ρ, cB and W have negligible horizontal variations (they are functions ofZ only) so the above equations become ∂ρ ∂t +W ∂ρ ∂Z + ρ∇.U = 0, (34) ∂cB ∂t +W ∂cB ∂Z + cB∇.U = ∂ ∂Z ( D ∂cB ∂Z ) , (35) where W =W(Z) is the submergence velocity profile in the column. Now suppose the depth–age scale Z = g(t), with the func- tion g given by t = g−1(Z)= Z∫ 0 dη W(η) . (36) We define z= Z− g(t) in order to use the reference frame of Eq. (1), which follows the material as it descends. The variable change fromZ to z gives ∂/∂Z→ ∂/∂z and ∂/∂t→ ∂/∂t −W(g(t))∂/∂z, and Eqs. (34) and (35) become ∂ρ ∂t + w̃ ∂ρ ∂z + ρ ( ∂U ∂X + ∂V ∂Y + ∂W ∂z ) = 0, (37) ∂cB ∂t + w̃ ∂cB ∂z + cB ( ∂U ∂X + ∂V ∂Y + ∂W ∂z ) = ∂ ∂z ( D ∂cB ∂z ) . (38) In this reference frame, material seen from the horizon with age t has the velocity w̃(z, t)=W(g(t)+ z)−W(g(t)). (39) On the short length scale of the signals (∼ 10−1–100 m), the vertical gradient of velocity W can be approximated by the strain rate, so w̃ ≈ ε̇z(t)z. Density and diffusivity varia- tions across individual signals can be assumed to be small on this scale, so we take ∂ρ/∂z≈ 0 and ∂(D∂cB/∂z)/∂z≈ D∂2cB/∂z 2. Applying the first of these approximations in Eq. (37) yields the compaction relation ∂U/∂X+ ∂V/∂Y + ∂W/∂z=−(∂ρ/∂t)/ρ, which, when used in Eq. (38), con- verts it to ∂cB ∂t − cB ρ ∂ρ ∂t =D ∂2cB ∂z2 − ε̇z(t)z ∂cB ∂z . (40) This is an approximate general evolution model for signals on cB in firn or ice. In the ice, where ρ ≡ ρi (constant), Eq. (40) loses the com- paction term and reduces to the same form as Eq. (1). This means that Eq. (1) is valid in the ice whether the bulk impu- rity concentration is defined in per mass or volume terms. In the firn, Eq. (1) is missing the compaction term of Eq. (40) so it cannot be used to track cB. But by substitut- ing cB = ρC into Eq. (40), we recover Eq. (1) exactly after some algebra. This means that Eq. (1) is valid in the firn and is the right model for formulating the diffusivity inversion, provided that C measures the bulk impurity concentration in per mass terms, as done by R2024 and us here. Inversions with the impurity concentration in per volume terms must use Eq. (40) instead. It follows that R2024’s reconstruction of the original peaks gives the right inputs, and the DR inversion is valid for both peaks in the firn (which experienced diffusion in a densifying material) and peaks in the ice (which experienced diffusion in a densifying material and then diffusion under constant density). A further realisation unknown to the earlier studies is that D retrieved for the firn by inversions based on Eq. (1) automatically quantifies its bulk diffusivity. It thus turns out fortunate that R2024 used the chemical measurements ex- pressed as sulphate concentration in per mass terms directly as C in Eq. (1). https://doi.org/10.5194/tc-19-5693-2025 The Cryosphere, 19, 5693–5717, 2025 5708 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice 4 Physical controls on sulphate diffusion at EDC Armed with D(t) from Fig. 8b, we discuss the mechanisms of sulphate transport in the EDC column, going beyond the interpretations made by R2024 and F2024 from their ef- fective diffusivities. D drops steeply from its surface value ≈ 10−6 to ≈ 1.7× 10−8 m2 yr−1 at 2.8 ka, an age coincid- ing roughly with the firn base (≈ 2.5 ka at 100 m depth); this drop is much shorter in duration than the initial decay in DR (Figs. 1c and 4). A slower decay in D to ∼ 10−10 m2 yr−1 follows from 2.8–390 ka, although its real form may not be exactly log-linear as posited in our forward model;D is sim- ilar at ≈ 125 ka to the DF1 and DF2 estimates of F2024 (1.6– 6.1× 10−9 m2 yr−1) but much lower in deeper ice. We con- sider these intervals in turn. 4.1 Vapour diffusion in the firn Recall that D retrieved for the firn reflects its bulk diffusiv- ity (Sect. 3.5). We interpret the high D on the steep drop as being due to diffusion of H2SO4 vapour through intercon- necting air pores in the firn. This mechanism is plausible be- cause H2SO4, though often viewed as nearly non-volatile, does have a vapour pressure (Tsagkogeorgas et al., 2017). A back-of-the-envelope calculation of the diffusion rate involving the H2SO4 vapour pressure pv and H2SO4 dif- fusion coefficient �a in firn air supports the interpreta- tion. We assume H2SO4 transport by vapour diffusion to be much faster than solid-state diffusion of sulphate through ice grains, but slow compared to sulphate exchange between ice and air, such that vapour diffusion is rate limiting. This as- sumption, which is justified by the time scales found below, features also in the Whillans and Grootes (1985) model for water isotope diffusion in firn, except evaporation replaces fractionation of the species here. To estimate pv, we as- sume sulphuric acid to be available on firn grain surfaces to exchange with vapour; then, Eq. (11) of Tsagkogeorgas et al. (2017) gives pv ∼ 5× 10−9 Pa at−53 °C. To estimate�a, we extrapolate the H2SO4 diffusion coefficients measured by Brus et al. (2017) at 278–298 K under laminar conditions to ≈−50 °C, finding �a ∼ 0.015 to 0.05 cm2 s−1 or ∼ 45 to 150 m2 yr−1 when bearing in mind the uncertainty in its power-law temperature dependence noted by these authors (see their Fig. 5 and Table 1). We adopt �a ∼ 150 m2 yr−1 from the top of the range, as it is based on a weaker tempera- ture dependence (a power of 1.5) that is more consistent with the one (1.75) found across the literature on gas diffusion (e.g., Tang et al., 2014). Now, if we focus on the top few metres of the firn and account for the relative density ρ0/ρi ∼ 0.4 there, then the diffusion coefficient for the bulk firn would be less, ≈�a(1− ρ0/ρi), but this correction is offset by a strong en- hancement of diffusion by firn ventilation and wind pumping (e.g., Colbeck, 1989; Waddington et al., 1996). We there- fore proceed by using the free-air diffusivity �a without correction. A ballpark estimate of D from vapour diffusion is found by scaling �a by the abundance ratio of SO2− 4 in the air to ice. The vapour pressure pv converts via the ideal gas law to 2.7× 10−12 mol m−3, whereas for a vol- canic signal in the firn with bulk concentration peaking at 200 ppb (≈ 2 nmoles g−1), the peak sulphate abundance is ≈ 2 mmol m−3 in the ice grains or ≈ 0.8 mmol m−3 in the bulk firn. The resulting abundance ratio, ≈ 3× 10−9, leads to the bulk diffusivityD∼ 4.5× 10−7 m2 yr−1, similar to the shallow high values on our D(t) profile. The assumption regarding time scales may be checked. For signals on the decimetre scale l ∼ 0.1 m, the vapour- diffusion time scale, l2/D ∼ 104 year, is much longer than the solid-diffusion time scale, d2 g/Ds ∼ 400 years. These val- ues are based on the mean grain diameter dg in the upper firn at Dome C (≈ 0.1–0.2 mm; Gay et al., 2002) and the as- sumption that the H2SO4 diffusivity within ice grains, Ds, is similar to the H2O self-diffusivity of monocrystalline ice (∼ 10−10 m2 yr−1 at −53 °C). The vapour diffusion model also explains the steep drop ofD in the range≈ 0–2.8 ka, because firn metamorphism re- duces the porosity and seals off interconnecting airways to lower the bulk diffusivity (Calonne et al., 2022) on descent through the firn layer, and because as grain growth occurs in the firn, larger grains slow the H2SO4 diffusion from their interior to their surfaces, progressively limiting the bulk dif- fusion rate. At pore close-off (ρ ≈ 845 kg m−3 at Dome C; Calonne et al., 2022) vapour diffusion terminates, and other processes must control D thereafter. 4.2 Sulphate transport below the firn-ice transition Below pore close-off at EDC, signal smoothing may result from (i) solid-state diffusion within ice grains, (ii) “Gibbs– Thomson diffusion” of vein signals (Ng, 2021; i.e., dif- fusion of the part of sulphate bulk concentration in liq- uid veins due to thermodynamic interactions including the Gibbs–Thomson effect and melting-point depression by dis- solved ions), (iii) diffusion through the grain-boundary net- work, (iv) “residual diffusion” caused by the stochastic three- dimensional motion of veins and grain boundaries carrying impurities (Ng, 2021), and any combination of these pro- cesses. We expect suppression of Gibbs–Thomson diffusion where the veins are disconnected or blocked by microparti- cles or dust (Ng, 2021), and suppression of residual diffusion where such particles impede grain-boundary motion (Durand et al., 2006). In terms of how sulphate is partitioned in EDC ice across crystal lattice, grain boundaries and liquid veins, direct ob- servations are limited to a small number of shallow ice sam- ples, but they indicate its presence at grain boundaries and triple junctions. Barnes and Wolff (2004) analysed 6 sam- ples in the 140–501 m depth range with scanning electron microscopy, finding sulphur at grain-boundary sites, and at triple junctions only in samples with high sulphate concen- The Cryosphere, 19, 5693–5717, 2025 https://doi.org/10.5194/tc-19-5693-2025 F. S. L. Ng et al.: Diffusivity of sulphate signals in Antarctic ice 5709 tration; they suggested that the veins could carry much im- purity only when the grain boundaries were saturated. Re- cently, Bohleber et al. (2025) used laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) to map the abundances of S, Cl and Na at tens-of-microns resolution in samples from 281.6, 585.2 and 1000 m depth (9, 27.3, 64 ka), in places away from volcanic spikes. Their elemen- tal maps show strong localisation of S at grain boundaries, with little of it in grain interiors, and no apparent trend in the partitioning with depth. They argued diffusion through grain boundaries (as well as veins) as potentially playing a key role in signal evolution. Measurements targeting triple junctions with a laser spot size of 1 µm also revealed S there, but the ablated material volumes were too small for deter- mining whether S was concentrated at the junctions – as was found for Na and Cl – and its abundance ratio between triple junctions and grain boundaries. If Ds, Dvn, Dbn and Dres symbolise the respective “com- ponent diffusivities” contributing to D from processes (i) to (iv) above, then one might regard D as their signal- partitioning weighted sum. In the following, we assess how they conspire with signal partitioning to different impurity sites to govern D(t) in 2.5≤ t ≤ 390 ka, by estimating their profiles down column. For reasons given below, we evaluate Dbn only qualitatively. We calculate Ds, Dvn and Dres by using published equa- tions (Appendix C) and temperature and grain-size data from EDC as input (Fig. S1 in the Supplement). These diffusivities increase with temperature. The Gibbs–Thomson diffusivity Dvn decreases with the mean grain size and sulphate bulk concentration cB (Eq. C2). We calculate it for cB from 1 µM (the typical background concentration at EDC) to 10 µM (or- der of magnitude of large volcanic peaks), noting that this range constrains a lower-bound Dvn because not all of cB may reside in veins. For finding Ds and Dvn, empirical esti- mates of the molecular diffusivities of sulphate in ice single crystals and water are desired but lacking; we approximate them by the H2O self-diffusivities in those materials. Turning to the grain-boundary network diffusivity Dbn, we note that grain boundaries, each bounded by triple junc- tions, must form a discontinuous transport network that is in- terrupted repeatedly by triple junctions, where we envisage liquid veins to exist. Sulphate can diffuse within each grain boundary according to its molecular diffusivity, Dgb; this is probably several orders of magnitude larger than the solid- state/lattice diffusivity, Ds (Lu et al., 2009; Ng, 2024). But on the centimetre or longer scale of signals of interest (over multiple grains), the veins intersect and strongly short-circuit the grain-boundary transport; the vein-water molecular dif- fusivity is several orders of magnitude above Dgb (Lu et al., 2009; Ng, 2024). Consequently, grain-boundary signal diffu- sion is inherently coupled to vein-signal diffusion. With Dbn representing diffusion of only the part of the bulk signal in grain boundaries, we expect Dbn