Tritium signal in Ice Cores

Tritium is the radioactive isotope of hydrogen, with a half-life of 12.32 years (Reference Lucas and UnterwegerLucas and Unterweger, 2000). As its natural concentration is extremely low it is commonly expressed in tritium units (TU). One TU corresponds to a ³H/¹H ratio of 10⁻¹⁸. Tritium is formed naturally in the upper atmosphere, due to cosmic rays interacting with atmospheric gases. Natural tritium is mainly produced by a reaction between nitrogen and thermal neutrons (Reference LibbyLibby, 1946):

The tritium formed is rapidly oxidized to H₂O and enters the hydrological cycle. The natural level of tritium in precipitation is not well known, as measurements started only after the first anthropogenic emissions in the early 1950s (see below). Measurements on vintage wines indicate a natural level of ∼5 TU for mid-latitude areas (Reference Kaufman and LibbyKaufman and Libby, 1954). Ice-core studies indicate a higher natural level of tritium in polar areas, 10–20TU for Greenland and the Canadian Arctic and 25–70TU for Vostok, Antarctica (Reference Kotzer, Kudo, Zheng and WorkmanKotzer and others, 2000; Reference Fourré, Jean-Baptiste, Dapoigny, Baumier, Petit and JouzelFourré and others, 2006). A higher tritium level in polar areas can be explained by an enhanced transport of tritium from the stratosphere to the troposphere at high latitudes (Reference Rozanski, Gonfiantini and Araguás-AraguásRozanski and others, 1991). The spatial variation of tritium in precipitation is also influenced by the distance from the ocean. Water vapour exchange with the ocean (which is low in tritium content) leads to low concentrations for coastal locations.

In the 1950s and early 1960s the level of tritium in the stratosphere increased by a factor of 10³, due to aboveground nuclear bomb tests. Most of the extra tritium is due to the direct injection of ³H into the stratosphere by thermonuclear bombs. In addition, the bombs produce a large number of neutrons, which enhance the occurrence of the reaction given in Equation (1). In the Northern Hemisphere the tritium content of precipitation reached a peak in 1963 of ∼5000TU. Since the Test Ban Treaty of 1963 only a few above-ground tests have been performed and the tritium activity in precipitation has slowly decreased. This decrease is not only the result of the decay of tritium, but also (and mainly) due to the mixing of stratospheric and tropospheric water. Once the water enters the troposphere it precipitates and most of it is stored on Earth as groundwater and in the oceans. Both these reservoirs have a long residence time and have taken up large amounts of thermonuclear tritium.

The mixing between stratosphere and troposphere is not constant throughout the year, but has a maximum in spring. The enhanced mixing leads to a higher tritium content in the precipitation. This was especially so in the 1960s, since most of the bomb tritium was produced in the stratosphere and the stratosphere/troposphere tritium gradient was large. In the tritium records, this seasonal variation with a maximum in May/June is observed. The 1963 tritium peak in ice-core records is often used for dating purposes. Despite the radioactive decay, this peak can still easily be found and is therefore an accurate time marker (Reference PinglotPinglot and others, 2003).

Spitsbergen Ice Cores

In spring 1997 a 121 m long ice core was drilled at the highest point (1250 m a.s.l.) of the Lomonosovfonna plateau (78°51′53″ N,17°25′30″ E) (Reference IsakssonIsaksson and others, 2001). This location (Fig. 1) was chosen based on results of earlier studies on seven ice cores from the Svalbard archipelago drilled by Soviet expeditions. One of the cores was drilled in 1976 on Lomonosovfonna at 1000 m a.s.l. Of the seven cores this was the least disturbed by meltwater percolation. The 1997 core indicates an average accumulation rate of 0.38 m w.e. a⁻¹ for the period 1963–97. Earlier studies of the effect of melt on the stable-isotopic and geochemical signals of this ice core were performed by Reference PohjolaPohjola and others (2002a) and Reference Moore, Grinsted, Kekonen and PohjolaMoore and others (2005). Both concluded that despite high melt ratios and percolation, the seasonal signals in the core were preserved. Reference PohjolaPohjola and others (2002a) found a melt index of 55% for the upper 36 m of the core, which corresponds to the time period AD1920–97. The melt index or melt-ice percentage was introduced by Reference Koerner and FisherKoerner and Fisher (1990) and is defined as the ratio between ice affected by melt (clear ice) and ice not affected by melt (bubbly ice). A 55% melt index in combination with an accumulation rate of 0.38 m w.e. a⁻¹ means that 0.21 m w.e. of the annual layer is snow filled with meltwater. Assuming an average snow density of 350 kg m⁻³ at the surface, this means that 0.37 m of snow has melted and the water has filled the pores of the underlying 0.23 m and frozen to form an ice layer with a density of 900 kg m ^{− 3}.

Fig. 1. Map of Svalbard illustrating the location of the two drill sites on Spitsbergen: Lomonosovfonna and Holtedahlfonna. The tritium ice-core records are compared with precipitation data from the coastal station Isfjord. Longyearbyen and Ny-Ålesund are denoted L and N, respectively.

The second ice core investigated in this study was drilled in April 2005 at a saddle point on Holtedahlfonna (79°8′15″ N, 13°16′20″ E) at 1150 m a.s.l. Due to its location, closer to the west coast of Spitsbergen, the annual accumulation rate at this site is higher than at Lomonosovfonna. Based on the depth of the 1963 bomb peak and density measurements we find a mean accumulation rate for the period 1963–2005 of 0.50 m w.e. With the higher accumulation rate this ice core is less sensitive to melt as the annual layers are thicker. However, the melt at this location is expected to be higher, due to the slightly lower altitude and location closer to the west coast, which lead to higher temperatures in this region.

Measurements

The tritium activity of the ice-core samples was measured using the gas proportional counting system of the Centre for Isotope Research, Groningen, The Netherlands. First, the water sample (5 mL, corresponding to a layer thickness of 5 cm) is reduced to hydrogen gas in a magnesium furnace at 600°C:

Using palladium as a catalyst, the hydrogen reacts with tritium-free ethylene to form the counting gas ethane:

The activity of the gas samples is then counted for 1 day. For samples with low activities, the counting time is extended to 2 or 3 days to obtain an uncertainty in the measured value of the samples of, at most, 2.5TU. If the activity of the sample is very low, as is the case for the deeper ice-core samples from 1950 and earlier, it is necessary to enrich the sample before measurement. This is done by electrolysis of the water using tritium-free NaOH as the electrolyte (Reference GröningGröning and others, 2009). The typical enrichment achieved is a factor 9, by electrolysing typically 90% of the water. After the enrichment, the remaining water is distilled to remove the electrolyte and measured in the same way as described above. This results in a typical measurement uncertainty of 0.2TU in the original sample.

Seven samples of the Lomonosovfonna ice core dated to be from the period ad1949–52 were enriched before measurement, both as a quality control for our measurements and in an attempt to determine the natural background level at the time of deposition. After correcting for the decay (a factor of ∼16), we find an average value of 48TU (ranging from 27 to 74TU) for the precipitation water from this period, with a measurement uncertainty of 4 TU.

This value can also be used as a quality control as it gives an upper limit for the level of contamination. The decay-corrected value of 48TU corresponds to a measured value of 3 TU. This means that before enrichment the tritium concentration of the sample was 0.3TU. This level of contamination is possible as the samples (not originally intended for tritium analysis) were stored in plastic containers through which tritium diffusion is possible. However, thanks to the high levels of tritium still in the water in the years of interest (note that these samples did not need to be enriched), even the worst possible case of 0.3 TU modern contamination is fully acceptable for our goals.

In total, 189 samples from the Lomonosovfonna core and 196 from the Holtedahlfonna core were analysed. Every sample covers a 5 cm thick layer. As the annual accumulation is 0.38 and 0.50 m w.e. for Lomonosovfonna and Holtedahlfonna, respectively, seasonal variations should be visible in the record even when the snow is transformed to ice as a result of compression and refreezing of meltwater.

The time between the first and last measurements of the ice-core samples was several years. To ensure that the effect of decay is not larger for the samples measured later, all measured values are corrected to represent the activity of the sample on 1 May 1997 (the time of the drilling of the Lomonosovfonna ice core, the older of the two cores). Unless otherwise stated, all tritium activities reported in this paper refer to the activity on that date. We explicitly avoid computing the tritium activity at the time of deposition, as the uncertainty in the depth/age relation would add an additional error. To facilitate a comparison between the measurements and precipitation record, the latter is also decay-corrected to May 1997.

Figure 2 shows the measured tritium profiles for the Lomonosovfonna and Holtedahlfonna ice cores. The timescale of the cores is determined by counting peaks in the stable-isotope and ion records (Reference IsakssonIsaksson and others, 2001), the tritium data presented here and other radioactive measurements (Reference PinglotPinglot and others, 2003). The 1963 bomb peaks are clearly visible in both records as the highest peak with a tritium ratio up to 450TU. Correcting this value for decay gives a value of ∼3000TU at the time of deposition. The actual value of the precipitation at that time will have been higher, as the tritium signal is smoothed due to diffusion and meltwater percolation. For Holtedahlfonna the 1963 peak is found at greater depth, which is caused by a higher accumulation rate in this area and by the 8 year time difference between the drilling of the two cores. Below the 1963 bomb peak other sharp peaks are visible in the profile. These are caused by bomb tests in the 1950s and early 1960s. For the period after 1963 no sharp peaks are found in the tritium profiles. Also, there is no clear seasonal variation in the record for this period, which is contrary to expectations based on the precipitation record, as shown in the next section.

Fig. 2. The measured tritium concentrations for (a) Lomonosovfonna and (b) Holtedahlfonna as a function of depth. The highest peak in both profiles corresponds to the year 1963. The 13–22 m core section from Lomonosovfonna spans the period ad1975–55. For Holtedahlfonna the 23–33.5 m section corresponds to ad1973–54. Correcting the 1963 peak values for decay yields a tritium concentration of ∼3000 TU at the time of deposition. However, as the profile is altered as a result of diffusion and percolation of meltwater it is likely that the tritium concentration in the precipitation was higher.

The Virtual-Ice-Core Model

To investigate the influence of melt on the ice-core record, we compare the measured tritium record with precipitation data. A model is constructed in which the main processes that alter the tritium record after deposition of the precipitation are included. In the model, layer after layer of precipitation, following a prescribed monthly precipitation record, are stacked upon each other. For each time-step, the tritium activity of every layer in the core is calculated using existing theories of diffusion, decay and densification. This produces a record which is valid for a situation where there is no melt. The effect of melt on the ice-core record is then investigated by including melt in the model, with some basic assumptions about how meltwater is redistributed over the underlying firn layers. The output of the model is a ‘virtual ice core’, which is used for comparison with measurements on the two Spitsbergen ice cores.

The main input of the virtual ice core is precipitation, which is added to the top of the ice core every month. The monthly amount of precipitation and the average ³H activity of the precipitation is retrieved from the global network for isotopes in precipitation (GNIP) database (International Atomic Energy Agency/World Meteorological Organization). GNIP consists of several stations where precipitation water is collected on a monthly basis. The total amount of precipitation water is recorded and samples are taken for isotope analysis. The first tritium analyses on these samples were done in the mid-1950s on samples from Ottawa, Canada, and Valentia, Ireland. However, as the tritium content in precipitation water varies spatially over the Earth, these samples only give an indication of the total tritium content in the stratosphere. Since the early 1960s, precipitation water has also been collected at Isfjord radio station, a coastal station on Spitsbergen. Tritium measurements of precipitation from this station started in July 1961 and continued until May 1975. The monthly amount of precipitation is recorded for a slightly longer period (January 1960 to December 1976). The tritium precipitation records from this station and from Ottawa and Valentia are given in Figure 3. They show that, although the tritium precipitation signal is qualitatively similar for the different locations, the magnitude of the signal varies with latitude and distance from the ocean, as discussed above.

Fig. 3. Tritium content in precipitation water for three different GNIP stations. The longest record (top) is from Ottawa, Canada. Comparing this record with those of Isfjord radio station (Spitsbergen, middle) and Valentia, Ireland (bottom) illustrates the large spatial variability in tritium content in precipitation water. The tritium concentrations given here are the actual values measured at the time (not corrected for decay).

The precipitation record used in the model is mostly based on the measured data from Isfjord. However, the model starts just before the first major bomb tests in 1953. Therefore, for the first period (1953–61) there are no data available for Isfjord and we create an approximate precipitation record based on measurements of Ottawa precipitation. For the period 1965–71 the tritium activity in Isfjord’s precipitation water is, on average, 61% of the tritium activity of Ottawa precipitation. Therefore, the tritium content of Ottawa precipitation for the period 1953–61 is scaled by a factor of 0.61 to obtain an approximate record for Isfjord for this time period.

The monthly amount of precipitation in this period is obtained from the average annual amount of precipitation at Isfjord in the 1960s, which is 0.48 m. As there is no clear seasonal variation present, we assume a fixed precipitation rate of 40 mm month⁻¹ for those months where no data are available. Using these assumptions we have created a continuous artificial precipitation record for Isjford for the period of interest (1953–75).

We assume that the tritium content in the precipitation is the same for the different locations in Spitsbergen. For the amount of precipitation, we scale the constructed Isfjord precipitation record based on the annual-layer thickness measured in the ice-core records. The average annual-layer thickness in the Lomonosovfonna ice core is 0.38 m w.e. As this is only 79% of the annual amount of precipitation at Isfjord, the constructed Isfjord precipitation record is scaled by 0.79 to obtain a record for Lomonosovfonna. Similarly, for the Holtedahlfonna virtual ice core a correction of 1.04 is used to compensate for the difference between the annual accumulation at this site (0.50 m w.e.) and the annual accumulation at Isfjord.

As precipitation is added on top of the virtual ice core, all layers below sink to greater depth. With increasing depth, the density of the firn also increases. Therefore, after each precipitation event a new density, and thus also a new thickness, is assigned to every layer. The density/depth profile used is based on measurements on the two ice cores (Reference PohjolaPohjola and others, 2002a; Reference SjögrenSjögren and others, 2007) (Fig. 4). In the measured density profiles we can easily identify the layers at low depth that have been affected by melt. Layers in which meltwater has refrozen are identified by a higher density. In the model we prescribe the density as a function of depth in the absence of melt. To obtain this function, a fit of the measurements is made excluding the higher-density layers. The data are fitted to a theoretical expression obtained from the Herron–Langway densification model (Reference Herron and LangwayHerron and Langway, 1980). Although the Herron–Langway model is developed for dry snow conditions and not for sites that experience melt on a regular basis, we find that the model fits the measured values well.

Fig. 4. Density/depth profiles for (a) Lomonosovfonna and (b) Holtedahlfonna. These data are used to obtain an approximate density/depth relation in the absence of melt. This relation is given by the solid curve which is a fit to Equation (5). In this fit the high-density layers at shallow depth (the gray dots in the graph) are excluded.

In the Herron–Langway model, the change in firn density, ρ, with depth, z, for a constant accumulation rate and temperature is given as

where ρ _ice is the density of ice (917 kg m⁻³) and K a constant. Integrating Equation (4) from the surface to some depth, z, gives the density as a function of depth:

where ρ ₀ is the density at the surface (z = 0). In the Herron–Langway model, two values for K are used (K ₁ and K ₂ for ρ < 550 kg m⁻³ and 550 < ρ < 800 kg m⁻³, respectively) to reflect a change in densification rate. As we observe no clear difference in the densification rate at different depths, K is kept constant over the entire depth. The function is fitted to the measured density profile using K and ρ ₀ as the fit parameters (Fig. 4). For Lomonosovfonna the best fit is obtained with K = 1.16 × 10⁻⁴ m² kg⁻¹ and ρ ₀ = 317.9 kg m⁻³, whereas for Holtedahlfonna the values for the fit parameters are K = 1.43 × 10⁻⁴ m² kg⁻¹ and ρ ₀ = 318.0 kg m⁻³. The fit parameters are then used in Equation (5) to assign a density to every layer in the model. This is, however, the density profile for the situation in which no melt occurs. The increase in density caused by refreezing meltwater is calculated separately in the melt module of the model.

Adding precipitation and recalculating the density of the firn is done every month, but the effects of diffusion and decay are calculated for a smaller time interval. The effects of diffusion on the isotope ratio can be described by the following differential equation (Reference JohnsenJohnsen, 1977):

where δ is the isotope ratio, Ω_fi the firn diffusivity and ε_Z the vertical strain. The second term on the right-hand side accounts for the thinning of the layers due to vertical strain, ε_z . Close to the surface the thinning of the layers is mainly caused by densification, which is discussed above. The thinning continues at larger depth where it is caused by gravitational spreading. To account for this effect we use a simple linear Nye thinning model (Reference NyeNye, 1963), which relates the layer thickness at a height h from the bed, L_h , to the original layer thickness at the surface, L_s :

where H is the ice thickness and all length scales are mw.e. This correction to the depth scale, as well as the correction for densification, is made every month in the model. Therefore, the second term on the right-hand side of Equation (6) is accounted for. What remains is Fick’s second law, which is solved numerically using a forward difference method. Typical values for the time and spatial step used in this method are 1 day and 5 mm, respectively.

The diffusion coefficient or diffusivity, Ω_fi, in the diffusion equation is a function of several parameters. A thorough description of the diffusivity and its dependencies on the temperature and structure of the firn (density, ρ, tortuosity, τ) was given by Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others (2000). Here we summarize their final result:

where m is the molar mass of the water molecule, R the gas constant and p _sat the water vapour saturation pressure which can be expressed as a function of temperature, T (Reference Murphy and KoopMurphy and Koop, 2005):

Following Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others (2000) we parameterize the tortuosity, τ, as a function of firn density, ρ:

Subscript i is added to those parameters that are different for different isotopes. In the diffusivity these are the ice-vapour fractionation factor, α_i , and the diffusivity of water vapour in air, Ω _ai . As neither of them is well known for tritiated water (¹H¹⁶O³H), we use estimates based on their values for deuterated water (¹H¹⁶O²H). The fractionation effects arise as a consequence of the mass difference between the rare isotopic molecule and the abundant isotopic molecule. Therefore, in general, the fractionation for the heaviest isotope is taken to be twice as large as that for the less heavy isotope (Reference MookMook, 2000); for example, for the carbon isotopes:

(Here α is the fractionation factor, while the fractionation is given by ε = α−1. As ϵ ≪ 1 this implies that when α ₁₄ = α ₁₃ ², ϵ ₁₄ ∼ 2ϵ ₁₃.) Also here, we assume that the fractionation for the heaviest isotope (tritium) is twice as large as for the less heavy isotope (deuterium). In reality the fractionation ratio may slightly deviate from 2, as was seen in the oxygen isotopes of water. Measurements of ¹⁷O and ¹⁸O of several natural waters by Reference Meijer and LiMeijer and Li (1998) showed a fractionation ratio of 1.8935 ± 0.005.

For the diffusivity in air of the most abundant water molecule (¹H₂ ¹⁶O) we have (Reference Hall and PruppacherHall and Pruppacher, 1976)

where T₀ and p ₀ are a reference temperature and pressure, equal to 273.15K and 1 atm, respectively. The diffusivity decreases when one of the atoms is replaced by a heavier isotope. For deuterium Reference MerlivatMerlivat (1978) found Ω_a2 = Ω_a/1.0251. This difference between normal and deuterated water is mainly caused by the mass difference between the deuterium atom and the hydrogen atom. As the mass difference between ³H and ¹H is twice the mass difference between ²H and ¹H, we assume that the reduction in diffusivity for tritium is twice the reduction for deuterium. Therefore we adopt for the diffusivity of ¹H¹⁶O³H in air:

Also for the fractionation for the phase transition from ice to vapour we assume that fractionation effects are twice as high for tritium compared to deuterium. Therefore we use:

where the parameterization of the deuterium fractionation factor as a function of temperature is taken from Reference Merlivat and NiefMerlivat and Nief (1967).

From the equations above it can be seen that the diffusivity depends strongly on the temperature of the firn (mainly through the saturation pressure of water vapour). Higher temperatures are associated with a higher diffusion coefficient, leading to a stronger smoothing of the isotopic signal. The temperature at the top of the firn column is influenced by the surface air temperature and by the latent heat released by refreezing meltwater. At greater depth the firn is more isolated from the surface and changes at the surface will not affect the temperature. The temperature profile used in this model is based on surface temperature measurements from Isfjord and on measured temperatures in the borehole of the core. Below 10 m depth the temperature in the borehole is nearly constant at −2.5°C (Reference Van de WalVan de Wal and others, 2002). In the model, we use this temperature for the firn below 10 m depth. At the surface the monthly temperature record of Isfjord radio station is used, corrected with a lapse rate of 0.44°C (100 m)⁻¹ (Reference PohjolaPohjola and others, 2002a). The surface and 10 m depth temperatures are linearly interpolated to obtain a first-order approximation for the temperature of the layers in between. If the temperature obtained in this way exceeds 0°C at any depth, the temperature is set to zero. In reality the temperature profile will be affected by the release of latent heat during refreezing of meltwater (e.g. Reference Pfeffer and HumphreyPfeffer and Humphrey, 1996), but for our purposes this simplified method is sufficient.

At every time-step the effect of diffusion is calculated and the tritium activity is corrected for decay. The activity changes according to the formula

where Δt and T_1/2 are the time-step and the half-life of tritium, respectively. The half-life of tritium is taken as 12.32 years (Reference Lucas and UnterwegerLucas and Unterweger, 2000).

The effects of melt and the percolation and refreezing of meltwater on the tritium profile can be included in the model by making some basic assumptions about the redistribution of water. Here the aim was not to include a full melt model, but only to quantify the effects of varying melt parameters. In the model it is assumed that all meltwater is taken up in the pores of the snow below the melting zone. This means that there is no run-off of meltwater. Melt is accounted for in the model by removing a layer at the surface (the melt layer) and redistributing the water and the corresponding tritium concentration over the underlying layers of firn (the percolation layer). The water can be redistributed over the percolation zone evenly or in such a way that most water stays at the top. This can be varied as a free parameter in the model, as well as the melt rate and the percolation depth. In the next section we explore the sensitivity of the resulting tritium profile to these model parameters.

Model Results

Before we include melt in the model and investigate the different parameters in the melt module, we first compare the measured data with a model run in which no melt takes place. Any differences observed in this comparison are then, in principle, due to melt effects. In Figure 5 the model run without melt (middle panel) together with the measured tritium profiles (left panel) for both ice cores are shown. The right panels of this figure show the precipitation records that are used as inputfor the model. The precipitation records shown here are corrected for decay effects in order to facilitate comparison with the measured and modelled data. For the same reason the depth scale of all the plots is given in mw.e. A slight shift upward of the measured and modelled profiles, compared to the precipitation record, is clearly visible. The thickness of a certain layer decreases with increasing depth due to the weight of the overlying ice. This is not accounted for in the stacked precipitation record, which explains the shift to lower depth of the precipitation record. A smaller shift can be observed in the modelled profile with respect to the measurements. This is the result of uncertainties in the assumptions for the strain rate and the monthly amount of precipitation in the model.

Fig. 5. Comparison of the measured profile (left panel) with a model run in which no melt is present (middle panel) for (a) Lomonosovfonna and (b) Holtedahlfonna. The right panels show the precipitation record that is used as input for the model. The difference in depth scale between the stacked precipitation and both the model results and measured data is the result of strain due to ice flow, which is not present in the precipitation.

We also observe differences in the variability of the signal between the measured data and the model output. For the period after 1965 (depths lower than 12.5 m w.e. for Lomonosovfonna (Fig. 5a) and lower than 20 m w.e. for Holtedahlfonna (Fig. 5b)) there is a seasonal signal visible in the modelled tritium concentration which is notvisible in the measured results. This is most obvious in the model run for Holtedahlfonna, where the layers are thicker and therefore the signal is better preserved in the model. In addition, for the deeper parts (14–16 m w.e. depth for Lomonosovfonna and 21–24 m w.e. depth for Holtedahlfonna) the sharp peaks visible in the measured data are much less pronounced in the model runs.

The melt module in the virtual-ice-core model has several parameters that can be varied. In the following we investigate the influence of these parameters in model runs for Holtedahlfonna. Model runs for Lomonosovfonna are qualitatively very similar to those for Holtedahlfonna. The main difference is the lower amplitude in the seasonal variation for Lomonosovfonna, due to the lower accumulation rate at this location.

One of the parameters that can be varied in the melt module is the redistribution scheme. This determines how the water from the melt layer is distributed over the underlying layers in the percolation zone. In the model the percolation zone is divided into four layers of equal thickness. A fraction of the total meltwater is assigned to each of these layers according to the redistribution scheme. Model runs with different distribution schemes are depicted in Figure 6a. The solid curve represents the model run in which the water is equally distributed over the percolation layer. For the other curves we assume a more realistic scenario, in which the amount of meltwater assigned to a certain layer decreases with depth. From this figure it is clear that changing the way the water is distributed over the four layers has only limited influence. A more equal distribution of the water over the four layers leads to a slightly higher amplitude and a small shift to larger depths, but the results are very similar.

Fig. 6. Model runs with different schemes for (a) redistributing the percolating meltwater over the underlying layers and (b) varying the distribution of the total melt over the different summer months. The numbers in the redistribution scheme indicate the fraction of the meltwater assigned to layers 1–4 in the percolation layer. Here the first layer is at the top and the fourth at the bottom of the percolation layer. The melt fraction is the monthly melt expressed as a percentage of the annual melt. These model runs show that both parameters have only limited influence on the resulting tritium profile.

The second parameter that can be varied within the melt module is the period over which melt takes place. Melt effects are calculated on a monthly basis (on the first day of the month) for certain months of the year. The amount of melt for each month as a percentage of the annual melt is varied between the different model runs depicted in Figure 6b. This shows that varying the melt fraction per month has negligible influence on the tritium profile.

Figure 7a shows model runs in which the depth to which the meltwater percolates is varied between 0.5 and 4 m with an annual melt of 40 cm. From these runs we can see that the annual oscillation is better preserved in the case of a shallow percolation depth. This can be explained by the density profiles for these runs (Fig. 7b). A shallow percolation depth leads to high-density layers, which effectively block the diffusion. This effect is incorporated in the model through Equations (8) and (10).

Fig. 7. (a) Different runs of the virtual-ice-core model in which the maximum depth to which the meltwater percolates is varied, and (c) runs in which the annual melt is varied. (b, d) A small section of the density profiles of the runs in (a) and (c), respectively.

Another parameter that can be varied within the melt module is the thickness of the melt layer. From Figure 7c it is obvious that a change in melt depth leads to a shift of the profile towards lower depths for increasing melt. Note that the depth scale in this figure is mw.e.; this shift can therefore not be explained by a change in density. One of the reasons for this shift to occur is the meltwater redistribution. It is the summer precipitation layer that is subject to melt, and the water in this layer is distributed over the underlying layers. As a result, the spring snow (which contains the largest tritium concentration) is found at shallower depth.

The density profiles of the different model runs are shown in Figure 7b and d. The main pattern of this profile is determined by the minimum density as a function of depth, which is calculated using the Herron–Langway model. Added to this minimum density, most layers have an additional amount of meltwater which leads to the scattering in the density profiles. In the model this extra meltwater is not allowed to increase the density of the firn above 917 kg m⁻³ (the density of ice). As soon as this density is reached, adding extra meltwater to the layer will instead increase the thickness. The density profiles in Figure 7b and d illustrate how strongly the density is affected by melt. Increasing the melt depth leads, on average, to higher density values, whereas increasing the percolation depth leads to a smoother density profile.

The density profiles obtained are very useful in our interpretation of the measured data. They show that the high-density layers measured for the cores (e.g. in the Lomonosovfonna core in the depth range 5–15 m; Fig. 4) can only occur with low percolation depth (Fig. 7b). Thus most of the meltwater refreezes in the layers just below the surface.

Discussion and Conclusions

In the previous section we explored the sensitivity of the different parameters in the melt module of the virtual ice core on the tritium record. In the following we compare the model results with the two Spitsbergen ice-core records, in as much detail as needed to pursue our goals: (1) to determine how far melt has destroyed the ability of the ice core to preserve the tritium signal present in precipitation and (2) to check in a semi-quantitative way the validity of our virtual-ice-core concept, especially the melt module.

For the period after 1961 the tritium content in this estimate is reliable as it is based on Spitsbergen data (albeit at lower altitude), but for the period before 1961 we have to rely on data from stations much further away. In contrast to the tritium content, the monthly amount of precipitation may also vary over short distances. Thus, the scaling of the amount of precipitation at Isfjord to represent the amount at the drilling locations introduces extra uncertainties. Therefore, we are limited to studying only the main features present in the record, such as the distinct bomb peaks (mainly the 1963 peak) and the seasonal cycle in the precipitation after 1965. Apart from the tritium profiles, we can also use the density profile for our comparison between the measured profile and the model results. As can be seen from Figure 7b and d, the density profile is strongly affected by both the amount of melt and the percolation depth. Therefore, the distribution of the layers with high density due to melt, percolation and subsequent refreezing is a valuable additional tool for the measurements/model comparison.

The direct effect of melt is a broadening of the tritium peaks to higher depths. However, the model results show that even with high annual melt it is still possible to have very sharp peaks. This is caused by the increased density in the layers, which results in a lower firn diffusion rate. In the measured tritium profiles of the two ice-core records, very distinct peaks are present for the period before 1963. The model runs show that these peaks would be much less pronounced or completely absent if no melt had occurred (Fig. 5). Thus, for this period the main effect of melt is that firn diffusion effects are reduced as a consequence of the increased density. We conclude that for this period most meltwater is stored within the annual strata, confirming the suggestions of Reference PohjolaPohjola and others (2002a) and Reference Moore, Grinsted, Kekonen and PohjolaMoore and others (2005).

A clear seasonal signal is not visible in the measured tritium profiles for Lomonosovfonna and Holtedahlfonna for the period 1964–75. This suggests that there was more melt in this period and meltwater may have percolated deeper into the firn pack. The large isolated peak in the Isfjord precipitation record in September 1972 is not visible in the ice cores. As this difference cannot be explained by diffusion and melt, we conclude that this precipitation event did not occur at the drilling locations.

The comparison between the different density profiles produced by the model and the measured profiles points to a relatively low percolation depth. Density profiles of model runs with high percolation depth yield only limited scattering of the density values. Only a high percolation depth with a high value for the annual melt can produce the same scattering as in the measured profiles (Fig. 4). However, if the annual melt is high, the density reaches the density of ice even at very shallow depths. As this is not visible in the measured profiles, we conclude that annual melting in both cores did not exceed 40 cm a⁻¹ of firn. This means that most of the meltwater is stored within the annual layers, supporting our conclusion based on the tritium profiles.

Our main conclusion is thus that the occurrence of melt does not necessarily lead to a strong smoothing of the tritium profile. Depending on the annual accumulation rate and the amount of melt, the occurrence of melt can lead to better preservation of the isotope signals. For the two ice cores discussed in this paper, the refreezing of meltwater seems to stop the effects of firn diffusion for the pre-1963 period and therefore lead to a more structured signal than in the case of no melt. Refreezing of meltwater takes place mainly in the layers directly below the surface. Thus, sub-annual signals in the core may disappear, but in general the annual signal in the core is hardly disturbed. The same process applies to the other isotopes of water and we can therefore conclude that the stable-isotope signal of these Spitsbergen ice cores can be used as a reliable estimate for annual temperatures for the period 1955–63 and for periods with similar or lower temperatures. It would be interesting to apply the model to ice cores from different locations to investigate how well it can simulate the real melting process in two cases: (1) when melt is much higher and meltwater percolates further down and (2) when melt occurs less often.

The whole virtual-ice-core approach can then be applied to the stable isotopes ²H and ¹⁸O, for which the annual signal is near constant and well known. For the cores discussed here the model outcome can be compared with the measurements (Reference PohjolaPohjola and others, 2002b). Especially, the prevention of diffusion by high-density layers caused by melt can be shown in that case. For chemical tracers (Reference Moore, Grinsted, Kekonen and PohjolaMoore and others, 2005), the diffusion process is of an entirely different nature. Their characteristics (diffusion through air channels and in the water and ice phases, co-transport with percolating meltwater) need to be brought into the model in order to interpret the measurements and the effect of percolation.