Inclusions in ice layers formed by melting and refreezing processes in a Greenland ice core

Abstract

In recent decades, surface melting on the inland Greenland ice sheet has increased, leading to significant meltwater-refreezing in the snow and firn. Increased knowledge of the physical and chemical characteristics of these melt features (i.e., ice layers) is needed to help estimate future global sea-level rise. Here, using a combination of microscopy and spectroscopy, we investigate the size, shape, location, and chemical compositions of inclusions in 12 ice layers of the SIGMA-A ice core in the northwestern Greenland ice sheet (78°03′06″N, 67°37′42″W, 1490 m a.s.l.). In the ice layers, we found inclusions exceeding 30 μm diameter that could only be formed by melting–refreezing, which we classify into columnar-, particle-, and rod-like inclusions. We propose that the smaller columnar- and particle-like inclusions of solid Na₂SO₄⋅10H₂O and CaSO₄⋅2H₂O form first, within the ice grains, followed by the larger rod-like inclusions of brines with mainly Na⁺ and Cl⁻ in grain boundaries. Our results suggest a new proxy that may help identify past warm climates in deeper ice cores in Greenland and for studying future ice sheet melting behavior.

1. Introduction

In recent decades, snow melting has been increasing on the surface of the Greenland Ice Sheet (van den Broeke and others, 2016; Otosaka and others, 2020). In the percolation zones, summer meltwater percolates through the snowpack and refreezes to form horizontally oriented ice slabs, ice layers (0.1–1 m thick), and ice lenses (<0.1 m thick) (Benson, 1962; Koerner, 1977; de la Peña and others, 2015; MacFerrin and others, 2019). Hereafter, we refer to all such features as ‘ice layers’. Ice layers have been studied remotely using radar (e.g., Otosaka and others, 2020) and satellite remote sensing (e.g., Tedesco and others, 2013), as well as climate models (e.g., MacFerrin and others, 2019), due to their existence as an indicator of surface melting.

Direct measurements of ice layers involve analyzing ice cores to determine the melt-feature percentage (MFP), which is then used to reconstruct the history of such surface melting (e.g., Herron and others, 1981; Koerner and Fisher, 1990; Graeter and others, 2018; Fujita and others, 2021). In the deeper part of ice cores with higher densification, the ice layers are difficult to discern, making the MFP there less reliable. Therefore, in deeper ice, past surface melting cannot presently be reconstructed from MFP. Nevertheless, in the firn, the process of redistributing chemical impurities associated with melting and refreezing can be a proxy for ice sheet surface melting.

After examining recent Greenlandic ice cores, Graeter and others (2018) argued that a warming of 1.2°C in the summertime led to a near doubling of the MFP in 1995–2015 compared to that in 1870–1900. Moreover, on 14 August 2021, rain fell at the summit area of Greenland, the area's first rain in recorded history that dates from 1950 (Scambos and others, 2021), implying that even such a dry snow zone may someday be a percolation zone (e.g., McGrath and others, 2013). In such a percolation zone, ice layers can cause meltwater to runoff much more rapidly, increasing the mass loss (Culberg and others, 2021) and associated sea-level rise. Thus, it is important to determine how such layers respond to a warming trend, particularly with the occurrence of Arctic amplification of global warming (e.g., AMAP, 2017). However, determining the ice-layer response requires knowledge of their chemical composition and inclusion properties, both of which are presently lacking.

Many studies (e.g., Eichler and others, 2001; Iizuka and others, 2002; Avak and others, 2019) have examined chemical-impurity redistribution in snow, firn, or refrozen ice during the melting. Most such studies focused on the redistribution before and after melting, not the distribution of inclusions and impurities within the ice layers. A few studies, such as Iizuka and others (2002), examined impurity redistribution in refrozen ice, but the method involved bulk ion concentrations after melting, not in the ice-layer's inclusions. The recent Arctic warming should affect not only the thickness but also the chemical impurity content of the ice layers. For example, if meltwater refreezing is incomplete (e.g., Ashmore and others, 2019), then the resulting runoff flushes out chemical impurities from the refrozen meltwater (e.g., Cragin and others, 1993). Such flushing points to the importance of understanding the chemical characteristics of ice layers, particularly in the inclusions formed by melting-refreezing processes.

Two promising micro-analysis methods for studying impurity chemistry within ice layers are SEM (scanning electron microscope)–EDS (energy dispersive X-ray spectroscopy) and Raman spectroscopy (Stoll and others, 2020). These approaches have been applied to dry ice cores without melting. For example, Mulvaney and others (1988) used SEM–EDS microanalysis on polar ice from Antarctica to show that H₂SO₄ concentrations are several orders of magnitude higher at triple junctions than in grain interiors. A similar finding was later made for H₂SO₄ and HNO₃ using Raman spectroscopy (Fukazawa and others, 1998). In addition, solid Na₂SO₄ and CaSO₄ salts in ice grains have been found in polar ice from Antarctica using Raman spectroscopy (Ohno and others, 2005, 2006; Sakurai and others, 2011). The use of SEM–EDS also showed the impurities were mainly in triple junctions and grain boundaries, with Na and Cl primarily in the Greenland ice core (GISP2), whereas Mg and S were primarily in an Antarctic ice core (Byrd) (Cullen and Baker, 2001; Obbard and others, 2003). Using a similar method, Stoll and others (2020) reported that insoluble inclusions are often found in the grain interior, whereas soluble inclusions occur at both the grain boundaries and grain interior. For a dry ice core, the inclusions were typically a micrometer in diameter. However, for an ice core with melting, these methods have yet to be applied.

With melt-refreezing, the impurities in meltwater become concentrated in the liquid phase during freezing, collecting mainly in the triple junctions in the last region of freezing (e.g., Takenaka and others, 1996; Bartels-Rausch and others, 2014). Therefore, in these ice layers, this approach should be particularly useful. Moreover, we have little understanding of the chemical characteristics of inclusions distributed by the refreezing. Therefore, it is essential to investigate how inclusions are distributed and preserved in the ice layer.

In this study, we aim to establish a new proxy for ice layers formed by refreezing processes. Using an ice core from northwestern Greenland, we examine the shape and chemical forms of inclusions in the ice layers by applying liquid analyses, as well as the micro-inclusion measurements of SEM–EDS and Raman spectroscopy. Here we describe how the inclusions redistribute in the thick ice layers during the refreezing process.

2. Methods

2.1. SIGMA-A ice core and stratigraphic analyses

Between 23 and 28 May 2017, we extracted a 60.19 m ice core near the SIGMA-A observation site (78°03′06″N, 67°37′42″W, 1490 m above sea level (a.s.l.)) in northwestern Greenland (Matoba and others, 2018). The site lies on the ridge of the Hayes Peninsula, 70 km northeast of the seaside village of Qaanaaq, in a percolation zone with an estimated annual average accumulation rate of 0.25–0.40 m yr⁻¹ between 1903 and 2005 (Kurosaki and others, 2020). The temperature at 20 m depth is −18.9°C. The SIGMA-A ice core covers from the period from 1903 to 2017 and was dated by examining the measured δ ¹⁸O, high resolution density, electric conductivities, and tritium content (Kurosaki and others, 2020).

Each SIGMA-A core section was cut to ~30 mm thickness with a band saw (Ryowa, BSW-200), and its surface was then polished with a microtome (Leica, 2400) to obtain parallel planes for physical observation in a cold room of the Institute of Low Temperature Science at Hokkaido University, Japan. The sections were photographed under transmitted light (Elpa, ALT-1090IR(D)) using a Nikon D3300 camera. In the photograph insets of Figure 1b, the firn is darker than the white ice layers due to the lower transmittance of light in firn. Ice layers were also identified by high (1 mm) resolution density by X-ray transmission method (Fig. 1a; Kurosaki and others, 2020).

Fig. 1. Depth profiles of physical properties and major ions in the ice. Depth profiles of (a) high resolution density, (b) melt feature percentage (%), and (c) ion concentrations (SO₄²⁻, Cl⁻, Na⁺, Ca²⁺, NH₄⁺). The inset photographs in (b) show ice layers formed in 2006 (6.056–6.211 m, 155 mm thick) and 1981 (20.581–20.596 m, 15 mm thick). The whiter (greater light transmission) parts are ice layers.

The SIGMA-A ice core has many ice layers. We chose 10 thin ice layers at random depths (Table 1) plus the two thickest ice layers, from 4.330–4.480 to 6.056–6.211 m (Table 1). With the following methods, we examined the shape and chemical forms of inclusions in the ice layers. In addition, we investigated the depth distribution of inclusions in the two thick ice layers at 4.330–4.480 m and 6.056–6.211 m depths.

Table 1. The methods applied to ice layers in different depths and years

2.2. Microscope observations and Raman spectroscopy analyses of inclusions

Using a microscope (Olympus STU-UM) in a cold room of −22°C, we observed micro inclusions in five cuboid samples that we chose as three thin ice layers at random from 2.780–2.805 m, 20.581–20.596 m, 34.610–34.661 m depth, as well as the two thick ice layers from 4.330–4.480 to 6.056–6.211 m (Table S1). (The objective lens is an Olympus ULWD Neo SPlan 50, with settings of numerical aperture (NA) 0.55, magnification 50× and 500×). We scanned the sample (field of view is 1.6 mm by 2.4 mm) while adjusting the z-axis by 5 mm to scan the depth. We examined any inclusion greater than 1 μm across, noting its location, shape, and size.

We then used micro-Raman spectroscopy (T64000, Horiba Jobin-Yvon) with a laser (Quantum Torus 532, Horiba) to determine the chemical form of selected inclusions. Our method is that described by Ohno and others (2005, 2006) and Sakurai and others (2010a). Briefly, the ice layer samples were put in a cold chamber at −30°C. Laser light (532 nm wavelength, software NGS LabSpec ver.5.64.15) of power 140 mW was focused on an inclusion using a long-working-distance objective lens with 6 mm focal length (M Plan Apo 100×; Mitutoyo). To identify chemical compositions in the inclusion, we compared the spectra to reference data (Ohno and others, 2005, 2014, 2016; Sakurai and others, 2010b, 2011).

2.3. Observation of nonvolatile inclusions by the sublimation–EDS method

To extract the nonvolatile inclusions in the ice layers, we used two ice sublimation methods. In the first method, we cut seven thin ice layer samples collected at random and two thick samples from the depths given in Table S2. Then, we extracted the nonvolatile inclusions at −22°C by the method described in Supplementary Information. In the second method, we pulverized four samples collected at random from the surface snow as well as ice layers in the 2.780–2.805, 20.581–20.596, and 22.611–22.631 m depths. We then extracted the nonvolatile inclusions from each 1 g sample at −50°C by following the sublimation method described in Iizuka and others (2009).

Constituent elements, inclusion diameters, and areas of each nonvolatile residue were measured using a JSM-6360LV (Japan Electron Optics Laboratory) SEM and a JED2201 (JEOL) EDS system. The accelerating voltage was 20 keV with a working distance of 20 mm and a collecting time of 30–60 s. We determined the levels of Si, Al, S, Cl, Na, Mg, and Ca in the inclusions, which were the major detected elements. In this analysis, O, C, and N were also detected but their levels were considered unreliable due to their light masses.

Particles containing Si are assumed to have silicate material, and almost all S and Cl in the nonvolatile particles are from soluble material of sulfate and chloride compounds (Iizuka and others, 2009, 2012b). The dry-air ventilation removes gas-phase acidic compounds such as HCl and HNO₃, and the remnant liquid acid droplets (H₂SO₄) are unstable under the measurement environment in a vacuum at +25°C of the SEM chamber. The coefficient of variation (CV) of the atomicity ratio from the SEM-EDS is based on the average values and standard deviations of these 20 measurements, which gives a value of 0.40 (Iizuka and others, 2012a). The sublimation-EDS method cannot be used to investigate the preservation and location of inclusions in the ice matrix, though it can be used to investigate the diameter, surface shape, and elemental compositions.

2.4. Ion concentration and water isotope ratio measurements

To investigate the depth distribution of ion concentration and water isotope ratio in the thick ice layers, we used the two thickest ice layers in the core at 4.330–4.480 m and 6.056–6.211 m depth, which are 150 and 155 mm thick. To measure the ion concentrations and the water isotope ratio of these ice layers, we first cut the ice layers from the core sections in the cold room (−22°C). Using a clean ceramic knife in the cold clean room (class 10 000), the samples were further divided into 15 samples 10 mm thick for the 150 mm ice layer and into 31 samples 5 mm thick for the 155 mm ice, then put into a clean polyethylene bag. The bagged ice was then melted at room temperature in a clean room (+25°C). We measured the ion concentrations of CH₃SO₃⁻, Cl⁻, SO₄²⁻, NO₃⁻, Na⁺, Ca²⁺, NH₄⁺, Mg²⁺, and K⁺ using ion chromatography (Thermo Scientific ICS-2100). For the cations, the column was Dionex CS-12A, and the eluent was 20 mM MSA. For the anions, the column was Dionex AS-17A column with 1–18 mM KOH gradient eluent. The injection volumes were 1000 μL for anion and 500 μL for cation. The analytical precision of the ion concentration was 10%. We also analyzed the stable isotope ratios of the samples using a water isotope analyzer (Picarro, L2120-i) with an evaporating device (Picarro, A0212 vaporizer). The analytical precisions of δ ¹⁸O and δD were 0.08 and 0.8‰, respectively.

For the ion concentrations and isotope measurements, the entire SIGMA-A ice core was cut into 894 samples with depth intervals 50–160 mm (average thickness of 67 mm) in a cold room (−22°C) and analyzed (Kurosaki and others, 2020).

3. Results and discussion

3.1 Physical and chemical characteristics of the SIGMA-A ice core

First, we considered the vertical profiles of density. Figure 1a shows that the x-ray-determined density increases steadily with increasing depth. Above 6.5 m, which dates back to 2004 (Kurosaki and others, 2020), a few particularly large peaks exceed 830 kg m⁻³. Such high values exist sporadically through the core and indicate thick ice layers corresponding with ice layers observed in the cold room.

There are 243 observed ice layers thicker than 1 mm since 1903, with an average thickness of 13.8 mm. These visual observations and density characteristics show that the frequency and the thickness of the ice layers has increased since 2004 (81 layers with an average thickness of 24 mm). The two thickest layers from 2012 and 2006 may be the result of the warming trend after 2000 in Greenland (Furukawa and others, 2017: MacFerrin and others, 2019). The MFP profile (Fig. 1b) shows the same trends, all being consistent with an increase of meltwater refreezing to ice layers as found by de la Peña and others (2015) in some regions of the Greenland ice sheet.

Most MFP values are less than 20%. The low values likely occurred during years when the ice layers formed within the year that the snow fell. The SIGMA-A site is in northwestern Greenland where temperatures in the winter are typically below −20°C (Aoki and others, 2014). Hence, the meltwater tends to refreeze within the year. On the other hand, the high MPF values in 2012 and 2006 likely are a result of meltwater percolating through several years of snowpack before freezing. These two years also contain the two thickest ice layers: 150 mm in 2012 and 155 mm in 2006. Figure 1b shows an image of a thick ice layer from 2006 and a thin layer from 1981.

Ion concentration profiles in Figure 1c are averaged over 10-year periods in Table S3. In the table, SO₄²⁻ has high concentrations during 1950–1980s (average of 3.35 μeq L⁻¹), a period when the emission of anthropogenic SO₄²⁻ was at a maximum (Fischer and others, 1998). In contrast, the core-average concentrations of NH₄⁺ and Ca²⁺ are much less, at 1.08 and 1.05 μeq L⁻¹. Their values increase after 2000, a trend also observed in other ice cores obtained from Greenland (e.g., Iizuka and others, 2018; Nagatsuka and others, 2021). The highest core-averaged concentrations come from Na⁺ and Cl⁻, at 4.51 and 5.29 μeq L⁻¹, generally over twice the averages of the other four ion species, suggesting a high impurity contribution from sea salt (Curtis and others, 2018), which is likely due to the site being in a coastal region. These Na⁺ and Cl⁻ concentrations vary by roughly a factor of two between decades, though the average decadal values hardly change over the past 100 years.

Snowmelt and refreezing processes affect ion concentrations (e.g., Iizuka and others, 2002), so we expect them to differ between firn and ice layer. The average concentrations (μeq L⁻¹) of Na⁺ and Cl⁻ are 5.25 and 5.65 in the upper firn, which are just above of the ice layers, 4.07 and 4.97 in the ice layers, and 3.99 and 4.90 in the lower firn, which are just below of the ice layers. In Figure 2, we compare the values in a given depth region to that in the region above; in particular, the ratios of the ice layers to the upper firn values, and of the lower firn to the ice layers values. For Na⁺ and Cl⁻ ratios for the ice layers to the upper firn (magenta), there are a relatively large number of values for ratios over 1.0, compared to the same ratios for the lower firn to the ice layers. This shift suggests ion enrichment from the upper firn to the ice layers, presumably due to ion migration in meltwater. Figure 2c also shows the general trend in Mg²⁺/Na⁺ for the upper firn, ice layers, and lower firn. If Mg²⁺/Na⁺ for the precipitation is assumed to equal that of the sea salt ratio (0.11), then a high ratio for the ice layers can only be attributed to the inflow of meltwater (e.g., Iizuka and others, 2002). From Figure 2c, the Mg²⁺/Na⁺ in the upper firn tends to be less than 0.11, whereas that in the ice layers tends to exceed 0.11. These results suggest that the chemical compositions tend to concentrate in the ice layers from the upper firn.

Fig. 2. Ratios of ion concentrations in the ice layers and firn. (a) For Na⁺ ratios. Not shown are the ratios exceeding 5.0, which comprise 16% of the cases for [Na⁺] _{Ice layer / Upper firn} and 10% of the cases for [Na⁺] _{Lower firn / Ice layer}. For [Na⁺] _{Ice layer / Upper firn} and [Na⁺] _{Lower firn / Ice layer}, the median and Std dev. are 1.26 ± 11.86 and 1.14 ± 13.47. (b) For Cl⁻ ratios. Not shown are the ratios exceeding 5.0, which comprise 13% of the cases for [Cl⁻] _{Ice layer / Upper firn} and 10% of the cases for [Cl⁻] _{Lower firn / Ice layer}. For [Cl⁻] _{Ice layer / Upper firn} and [Cl⁻] _{Lower firn / Ice layer,} the median and Std dev. are 1.41 ± 7.85 and 1.18 ± 5.90. (c) Mg²⁺/Na⁺ ratios. Not shown are Mg²⁺/Na⁺ ratios over 0.4. These comprise 0.9, 7.3, and 4.6% of the cases for the upper firn, the ice layer, and the lower firn, respectively.

Considering now the Na⁺ and Cl⁻ ratios of the lower firn to that in the ice layers, the values seem to show two peaks in Figures 2a and 2b, one of which is less than 1.0 and the other is ~ 1.5. The peak below 1.0 suggests that the ice layers can collect ions that would otherwise flush to the lower firn; on the other hand, the peak near 1.5 suggests further ion migration occurred to the lower firn from the meltwater, which will become ice layers. Finally, the Mg²⁺/Na⁺ values in the ice layers and the lower firn tend to exceed 0.11, suggesting that these regions are partly formed by meltwater from the upper firn.

3.2. Shape and chemical components of inclusions in the ice layers

3.2.1 Microscope observation and Raman analyses

We chose five ice layers, dated 1955, 1981, 2006, 2012, and 2014 in which to examine the inclusions using microscope observations and Raman analyses (Table S1). In these layers, we first examine the shape and size of 513 inclusions. Of the 513, the most common shapes are particle-like (n = 253) and rod-like (n = 230), with the remaining 30 classified as columnar-like.

Concerning this classification, we first check if the inclusion is angular, and if so, it is classified as columnar. If not, then we check if the noncolumnar inclusion has an aspect ratio (major-to-minor axes) exceeding 2:1, and if so, the inclusion is classified as rod-like. What remained were particle-like, an example of which is shown in Figure 3a. The particle-like inclusions occurred both in ice grains and on grain boundaries. The rod-like inclusions were found in grain boundaries and triple junctions (Fig. 3b) and had diameters of about 10–100 μm. This type generally had the longest lengths of all inclusions. Almost all rod-like inclusions were detected in thick ice layers such as those from 2012 and 2006 (e.g., Fig. 1b). The columnar-like inclusions were found in ice grains (Fig. 3c). Some of the particle- and columnar-like inclusions are greater than 30 μm in diameter.

Fig. 3. Microscope images and Raman spectroscopy analyses of the three types of inclusions. (a) Particle-like inclusion in an ice grain from the 2014 ice layer, with peaks at 990 and 1008 cm⁻¹. At right is the distribution of chemical compositions of 85 particle-like inclusions from Raman analyses from all ice layers. Colors, defined at bottom, are based on peaks at 982, 984, 990, 1008 cm⁻¹. ‘Others’ means peaks found at other wavenumbers. (b) Rod-like inclusion in a grain boundary from the 2006 ice layer, with a peak of 984 cm⁻¹. Right side is the same as that in (a) except for 97 rod-like inclusions. (c) Columnar-like inclusion in a grain interior from the 2006 ice layer, with a peak of 1008 cm⁻¹. Right side is the same as that in (a) except for 30 columnar-like inclusions.

The chemical forms of the inclusions are measured by Raman spectroscopy. In the spectra, we look for peaks at 982, 984, 990, and 1008 cm⁻¹. The 982 and 984 cm⁻¹ peaks have a broad half-width, suggesting liquid inclusions. On the other hand, the 990 and 1008 cm⁻¹ peaks are sharp, suggesting solid inclusions. Following Ohno and others (2005), we assume the 990 cm⁻¹ peak is Na₂SO₄⋅10H₂O, and the 1008 cm⁻¹ peak is CaSO₄⋅2H₂O. The broad peaks at 984 cm⁻¹ indicate SO₄²⁻ in a liquid phase (Sakurai and others, 2017). However, we found no prior reference for the broad 982 cm⁻¹ peak, so we ran a few tests that identified the peak as SO₄²⁻ and NH₄⁺ in a liquid phase (Supplementary information).

Among our 525 measurements, 221 showed Raman peaks, whereas 304 measurements did not (Table S4). Of the 221 with peaks, 96% (n = 213) show at least one chemical compound at liquid (NH₄)₂SO₄, liquid H₂SO₄, solid Na₂SO₄⋅10H₂O, and solid CaSO₄⋅2H₂O. The breakdown of the 213 peaks is such that 31% have the liquid (NH₄)₂SO₄, 31% have the liquid H₂SO₄, 10% have the solid Na₂SO₄⋅10H₂O, and 28% have the solid CaSO₄⋅2H₂O (Table S4).

The specific peaks depend on the inclusion's shape-type as shown by the pie charts on right side of Figure 3. Among the particle-like type, individual inclusions had one or two compounds, distributed between all four compounds (Fig. 3a, right). Nearly all of the rod-like inclusions are either liquid (NH₄)₂SO₄ or liquid H₂SO₄ (35 and 46‰ of cases in Fig. 3b). On the other hand, all columnar-like inclusions are solid CaSO₄⋅2H₂O (gypsum) (Fig. 3c). There are many inclusions with no Raman peaks, with 56% of them being the particle-like inclusions and the remaining 44% being the rod-like inclusions. The chemical composition of these inclusions is discussed next.

3.2.2 Shape and elemental components of nonvolatile inclusions in ice layers

The microscope observation method used above has the advantage of determining where the inclusions are located in an ice layer (e.g., within ice grains, along grain boundaries, in triple junctions). However, the method cannot detect all the inclusions because bubbles in an ice layer make it difficult to locate every inclusion. On the other hand, the sublimation–EDS method can measure the shapes and sizes of almost all inclusions formed in the ice layers.

The SEM measurements revealed 1230 nonvolatile inclusions in the three shape types through the 12 ice layers. Of these types, the particle-like inclusions occur in almost all ice layers and account for 94.4% (n = 1161) of the total inclusions (light and dark blue bars in Fig. 4b). Almost all rod-like inclusions occur in just the two thickest ice layers, though they account for 3.1% (n = 38) of the total inclusions (pink and red bars in Fig. 4b). The columnar-like inclusions exist in both thin and thick ice layers, but account for only 2.5% (n = 31) of the total inclusions (light and dark green bars in Fig. 4b).

Fig. 4. Depth profiles of composition and inclusion type, as well as the elemental compositions in the nonvolatile inclusions. Analyses cover the surface snow and selected ice layers (Raman: 5 layers; SEM & EDS: 12 layers). (a) The number of times the chemical composition is detected by Raman analyses. (b) The number of each type of nonvolatile inclusion detected by SEM. (c) The number of times the element combinations are detected by EDS analyses in nonvolatile inclusions.

The cross-sectional area S (μm²) of an inclusion is measured by tracing the inclusion's perimeter on a photomicrograph. For the particle-like inclusions, we then calculate the radius r (μm) as

(1)$$\matrix{ {r = \sqrt {\displaystyle{{\rm S} \over {\rm \pi }}} \;} \cr } $$

We then define the grain diameter as twice the radius. For rod-like and columnar-like inclusions, the grain length is defined as the major axis. We refer to the grain diameter and length as D_inclusion (μm).

The two methods (microscope and SEM) similarly find the particle type to be relatively common and the rod type to be more common in thick ice layers. Moreover, both methods show columnar to be the least common type. The proportion of particle-like is much larger in the SEM method than in the microscope method. The difference is likely due to the particle-like inclusions being generally smaller (see below) and more easily detected with the SEM method. In addition, the rod-like inclusions with small D_inclusion (e.g., Fig. 8a, see ‘(X)’) are difficult to distinguish from particle-like inclusions in the SEM method, so some rods may have been counted as particle type.

Almost all nonvolatile inclusions in the surface snow sample at the site are particle-like (Fig. 4b), with an average grain diameter of 9.0 ± 8.7 μm. Also, nearly all inclusions in the surface snow are smaller than 30 μm in diameter, with 66% being smaller than 10 μm. These particle types and sizes are consistent with nonvolatile inclusions in dry snow and ice cores in Antarctica and Greenland analyzed previously with the sublimation method (Iizuka and others, 2012a, 2012b, 2013; Oyabu and others, 2014, 2015, 2020). Given that aerosols are generally less than 30 μm in diameter (e.g., Whitby, 1978), we conclude that the inclusions greater than 30 μm in diameter formed by refreezing meltwater. So, we divide the inclusions further into two size categories: smaller and larger than 30 μm in diameter. For particle-like inclusions, 22% are greater than 30 μm in diameter (dark blue bar in Fig. 4b). For rod-like inclusions, this percentage is 87% (red bar in Fig. 4b), whereas for columnar-like inclusions, the percentage is 61% (dark green bar in Fig. 4b). Thus, most rod-like and columnar-like inclusions are over 30 μm in diameter, and therefore likely formed by refreezing.

The elemental analyses of the nonvolatile inclusions in the surface snow sample and 12 ice-layer samples are shown in Figure 4c. We do not discuss Mg because its abundance is very small (5%). Additionally, we do not discuss Al because it is mainly in dust (insoluble terrestrial materials), which has Si as the main component. As a result, we classified the content of the inclusions into 12 categories based on their constituent elements of Na, Ca, S, and Cl, which are associated with Na⁺, Ca²⁺, SO₄²⁻, and Cl⁻ species. For example, inclusions containing all of Na, Ca, S, and Cl are classified as ‘Na∩Ca∩S∩Cl’. Those with Na and S but without Ca nor Cl are classified as ‘Na∩S’. Those with only S or only Cl are classified as ‘only S or Cl’. Finally, ‘Others’ have neither S nor Cl. (Others mainly contain Si, suggesting silicate minerals.) In the surface snow sample, 36% of the inclusions are classified as ‘Na∩S’, and nearly all of the remaining are ‘Others’ (Fig. S1a). Compared to the surface snow, inclusions in the ice layers contain much more ‘Na∩Ca∩S∩Cl’ (dark blue; 11%), ‘Na∩S∩Cl’ (dark red; 11%), and ‘Na∩Cl’ (dark green; 10%), and much less ‘Na∩S’ (light pink; 5%) (Fig. S1b). These differences are more pronounced in the larger inclusions (Fig. S1c). So, the combination of ‘Na∩Cl’ with or without S and Ca is common in the inclusions, particularly those formed by refreezing. The meltwater contains Cl⁻ (4.97 μeq L⁻¹) and Na⁺ (4.07 μeq L⁻¹) with a few SO₄²⁻ (2.85 μeq L⁻¹) and Ca²⁺ (1.28 μeq L⁻¹), and these ions are considered to aggregate as brine or salt during the refreezing.

The most common signal in the particle type is ‘Others’ (grey; 42% in Fig. 5a). Most of the ‘Others’ inclusions contain Si without Na, Ca, S, and Cl, suggesting that the particle-like inclusions mainly consist of mineral dust. The ‘only S’ inclusions may be NH₄⁺ and SO₄²⁻ as liquid brine because Raman measurements showed the liquid 982 cm⁻¹ peak, not the solid 975 cm⁻¹ peak. In this case, liquid brine may exist at the ice temperature of −18.9°C due to the eutectic temperature of (NH₄)₂SO₄ being −19.0°C, or due to an ion balance of SO₄²⁻ ≫ NH₄⁺ (Fig. S2c).

Fig. 5. SEM-EDS analyses. (a) Particle-like inclusion from the 2012 ice layer and composition of 1161 particle-like inclusions from EDS analyses. (b) Rod-like inclusion from the 2012 ice layer and composition of 38 rod-like inclusions. (c) Columnar-like inclusion from the 2006 ice layer and composition of 31 columnar-like inclusions.

In contrast, the ‘Na∩Cl’ with or without S and Ca inclusions are considered to be salt mixtures of NaCl, Na₂SO₄, and CaSO₄. The Na₂SO₄ and CaSO₄, having eutectic temperatures of −1.3 and −0.05°C, are likely to exist as solid salts, whereas the NaCl, having a eutectic temperature of −21.3°C, can exist as a liquid brine of Na⁺ and Cl⁻ at −18.9°C. Thus, by combining the Raman and EDS analyses, we argue that the particle-like inclusions contain various solid salts of Na₂SO₄⋅10H₂O and CaSO₄⋅2H₂O, as well as brines of Na⁺, Cl⁻, NH₄⁺, and SO₄²⁻.

The primary category for rod-like inclusions is ‘Na∩Ca∩S∩Cl’ (dark blue; 45% in Fig. 5b). On the other hand, the Raman measurements show that the rod-like inclusions have mostly broad peaks of 982 and 984 cm⁻¹ (35 and 46% in Fig. 3b), as well as a high fraction of nonpeak-detected inclusions. The broad 982 and 984 cm⁻¹ peaks suggested that inclusions have NH₄⁺ and SO₄²⁻ (Fig. S2c), whereas the high fraction of nonpeak-detected inclusions suggests an ionic compound of Na⁺ and Cl⁻. Sakurai (2010) showed that mixtures of Na₂SO₄ and HCl have an extremely low eutectic temperature of −42.7°C (Na₂SO₄: 4.3 wt%, HCl: 2.1 wt%), which is colder than that of the Raman measurements (−30°C). Thus, a possible chemical form of the rod-like inclusions is a brine of SO₄²⁻, Cl⁻, Na⁺, and NH₄⁺.

For the columnar-like inclusions, the main elemental component is ‘Ca∩S’ (green; 55% in Fig. 5c). Thus, the columnar-like inclusions probably consist of CaSO₄. Also, their shape is consistent with the crystal structure of gypsum (CaSO₄⋅2H₂O), which is the prismatic class of the monoclinic system. These results suggest that the columnar-like inclusions contain impurities with the components of CaSO₄, which is consistent with the results from Raman analyses that indicated that the columnar-like inclusions contain solid gypsum salts (CaSO₄⋅2H₂O).

3.3 Chemical form and shape of inclusions in thick ice layers formed in 2012 and 2006

The two thickest ice layers of 2012 and 2006 have many large rod-like inclusions (brine of SO₄²⁻, Cl⁻, Na⁺, and NH₄⁺). In this section, we examine the chemical forms, shapes, and distributions of inclusions in the thick ice layers.

3.3.1 The refreezing processes of the two thick ice layers

Consider first the ice layer formed in 2012. Figure 6 shows its structure analyzed to 10 mm resolution. In (a), the values of δ ¹⁸O and δD vary by 2.1 and 17‰, respectively. According to previous studies (e.g., O'Neil, 1968), these values are consistent with meltwater refreezing. δ ¹⁸O and δD have a correlation (Craig, 1961), and Dansgaard (1964) defined the deuterium excess (d−excess) as the following:

(2)$$d-excess = {\rm \delta D}-{\rm \;}8{\rm \;}{\rm \delta }^{18}{\rm O\;}$$

Fig. 6. Depth profiles of the thick ice layer formed in 2012. (a) δ ¹⁸O, δD, d-excess, and major ion species (SO₄²⁻, Cl⁻, NH₄⁺, Na⁺, Ca²⁺). (b) The number of times the chemical compositions are detected by Raman analyses. (c) The number of inclusions of various types. (d) The number of times elements are detected by SEM analyses. Yellow bands are levels that likely froze last.

Lacelle (2011) suggests that the last ice to freeze is likely to have a higher d-excess value, such as that at depths ~ 4.370–4.390 m and 4.430–4.450 m, which we mark light yellow in the depth profiles. Hence, we believe that the 2012 ice layer is formed by at least three refreezing cycles.

We show similar analyses for the 2006 ice layer in Figure 7. (Inclusions above 6.071 m depth and below 6.191 m depth are not analyzed here due to the lack of water isotope data in those regions.) In (a), the values of δ ¹⁸O and δD fluctuate by 1.2 and 6.9‰, respectively. The slope of δ ¹⁸O and δD is 5.6, suggesting that the ice was formed by the refreezing of meltwater (e.g., Jouzel and Souchez, 1982; Souchez and Jouzel, 1984; Souchez and de Groote, 1985). We consider the 6.101–6.126 m depths as the last region to freeze, due to the low δ ¹⁸O and high d-excess there, which are marked in light yellow in Figure 7. We believe that at least two refreezing cycles formed the 2006 ice layer.

Fig. 7. Same as Figure 6, except from the ice layer formed in 2006.

Ion concentrations of all species are relatively high at the depths to freeze last (Figs 6a, 7a). For the 2012 ice layer (~ 4.370–4.390 and 4.430–4.450 m depths), each ion concentration is 1.1–1.8 times higher than that of other depths (Table 2). For the 2006 ice layer, each ion concentration except NH₄⁺ (μeq L⁻¹) in the last region of refreezing (~ 6.101–6.126 m) is also 1.1–1.4 higher than other depths (Table 2). As ice rejects impurities upon freezing (e.g., Halde, 1980; Takenaka and others, 1996), ionic materials in the thick ice layers should be concentrated in the last ice to freeze.

Table 2. The average of major ion concentrations (μeq L⁻¹) in the thick ice layers of 2012 and 2006

^a Last regions to freeze.

3.3.2 Formation mechanism of inclusions during the refreezing processes of the two thick ice layers

The number and size of inclusions in the last region to freeze gives clues to their formation mechanism. To show the spatial distribution of inclusions for the example in the last region of the 2012 ice layer (4.370–4.390 m), we present Fig. S3. For all inclusions in the entire 2012 ice layer, the SEM measurements give the number n and average diameter φ of 574 and 36 ± 99 μm. Considering just the last two regions to freeze (4.370–4.390 m and 4.430–4.450 m), the n and φ are 311 and 34 ± 115 μm. For the entire 2006 ice layer, these values are 224 and 37 ± 78 μm. At the last region of refreezing in the 2006 layer (6.101–6.126 m), the values are 143 and 32 ± 45 μm. Table 3 shows the average diameter (μm) and areas (μm²) of three types of inclusions in the last region of refreezing in the 2012 and 2006 ice layers. Because of long length, the average diameters and areas of the rod-like inclusions are much larger than those of the particle-like and columnar-like types. For the last refreezing of 2012 and 2006 ice layers, the large rod-like inclusions are abundant (red bar in Figs 6c, 7c; red frame on box in Fig. S3). In the 2012 ice layer, 64% of the total rod-like inclusions over 100 μm in diameter occur (n = 7) in the last ice to freeze (4.370–4.390 m and 4.430–4.450 m). In the 2006 ice layer, 45% (only 5) of total rod-like inclusions over 100 μm in diameter occur in the last ice to freeze (6.101–6.126 m).

Table 3. Average of diameters (μm) and areas (μm²) of inclusions in the last stage of refreezing in the 2012 and 2006 ice layers, according to the SEM analysis

Figure 8 shows the photomicrographs of the last region of refreezing in the 2012 and 2006 ice layers. The arrows in the micrographs (Figs 8a, 8b) show inclusions observed in the ice grains or grain boundaries. The rod-like inclusions are found here along grain boundaries (e.g., Fig. 8a X and Y, 8b Z), mostly in the last regions of refreezing. In the 2012 ice layer, the largest ones are about 150 μm long and 7.8 μm wide (Fig. 8a Y), and for the 2006 ice layer, they are about 500 μm long and 15.8 μm wide (Fig. 8b Z). As such, they are much wider than the other filaments in ice cores that never experience melting (e.g., Baker and others, 2003). The columnar-like inclusions instead are found within the ice grain (e.g., Fig. 8a arrow 4 and 8b arrow 6), whereas the particle-like inclusions (e.g., Fig. 8a arrow 1 and 8b arrow 7) are found both on the grain boundary and within the ice grain.

Fig. 8. Photomicrographs and Raman spectra of inclusions. (a, b) Photomicrographs of regions within 4.370–4.380 m (last region to freeze) in the 2012 ice layer (a) and 6.116–6.121 m (last region to freeze) in the 2006 ice layer (b). The left edge is shallower depth, right edge deeper. ‘TJ’ indicates a triple junction. Green arrow is a particle-like inclusion of double salts of Na₂SO₄⋅10H₂O and CaSO₄⋅2H₂O within an ice grain. Orange arrows are liquid particle- and rod-like inclusions of (NH₄)₂SO₄ in grain boundaries. Red arrows: similar to orange arrows except H₂SO₄. Blue arrows are columnar-like inclusions of CaSO₄⋅2H₂O in the ice grain. The purple arrow is a particle-like inclusion of H₂SO₄ and CaSO₄⋅2H₂O in the triple junction. White arrows mark particle- and rod-like inclusions at grain boundaries or ice grains that have no Raman peaks detected. The black lines show a smaller rod-like inclusion (X) and larger rod-like inclusions (Y and Z). (c) Raman spectra of inclusions marked in (a) and (b).

With regard to the inclusion compositions, the Raman measurements (see section 3.2.1) indicate that particle-like types have a brine of NH₄⁺ and SO₄²⁻ (orange and red arrows in Figs 8a, 8b) as well as salts of Na₂SO₄⋅10H₂O (s) and CaSO₄⋅2H₂O (s) (green arrow in Fig. 8a). The columnar-like types have salts of CaSO₄⋅2H₂O (s) (blue arrows in Figs 8a, 8b). Additionally, the rod-like types have a brine of NH₄⁺ and SO₄²⁻ (orange and red arrows in Figs 8a, 8b). However, in the 2012 ice layer, 65% of the Raman measurements of rod-like inclusions in 4.370–4.390 m, one of the last regions to freeze, show no significant peaks (60% for section 4.430–4.450 m). Similarly, in the 2006 ice layer, 55% of the rod-like inclusions in the last region of refreezing show no significant peaks. These rod-like inclusions are marked by white arrows in Figures 8a and b. The lack of peaks is consistent with the rod-like inclusions containing Cl⁻ and Na⁺ as described in section 3.2.

In considering how an inclusion's composition relates to its size, we examined the relationship between the atomic fractions of Na, S, Cl, and Ca on the inclusion's cross-sectional area. The sum of four fractions equals 1. For the inclusions in the last stage of refreezing, the fractions of Ca in Figure 9 decrease significantly as the inclusion areas increase. Similarly, the fractions of S decrease gradually in the 2012 ice layer (4.370–4.390 m) and 2006 ice layer. The fraction of Na also decreases gradually, reaching about 0.4 in the largest areas. For the 2012 ice layer (4.430–4.450 m), the fraction of Cl instead increases to about 0.5 as the areas increase (Fig. 9b). In contrast, in the 2012 ice layer (4.370–4.390 m) and 2006 ice layer, the fraction of Cl has no significant change. However, the fraction of Cl is probably the same as Na, since the fractions of S and Ca decrease to 0 as the inclusion areas increase. Given that the area of the rod-like inclusions tends to be large, these results suggest that large rod-like inclusions in the last region of refreezing consist mainly of Na and Cl elements, with relatively little S and Ca. As described above, the rod-like inclusions show no significant Raman peaks, consistent with a large amount of Na⁺ and Cl⁻ (brine) in the rod-like inclusions. On the other hand, S and Ca increase as the areas decrease.

Fig. 9. The fractions of Na (red), S (green), Cl (blue), and Ca (black) of a given inclusion and the inclusion area (μm²) for inclusions formed in the last region of refreezing. Also shown are the least-square-error straight-line fits for each element. Each of correlation coefficient (r) and p-value are in the box. Inclusions without Na, S, Ca, and Cl have mainly Si. We consider them as insoluble inclusions and do not plot them here.

Generally, the impurities in water concentrate in the liquid phase during freezing, with a tendency to collect in the triple junctions in the last region of freezing (e.g., Takenaka and others, 1996; Bartels-Rausch and others, 2014). The waterfront pressed the brine further down until the porosity or cold front prevents further progress. Then the brine is more concentrated and preserved as rod-like inclusions along the resulting grain boundary or triple junction. Thus, the redistribution of impurities can be explained by fractionation between solid and liquid impurities during the freezing. We summarize the proposed process in Figure 10: (a)–(b) The smaller columnar- and particle-like inclusions of solid Na₂SO₄⋅10H₂O and CaSO₄⋅2H₂O form in the ice grains. (c) Particle-like inclusions of brines with Na⁺, Cl⁻, NH₄⁺, and SO₄²⁻ form in grain boundaries. (d) Finally, if sufficient impurity exists in the meltwater, large rod-like inclusions of brines with Na⁺, Cl⁻, NH₄⁺, and SO₄²⁻ form at grain boundaries. In this way, when the ice layer is particularly thick, as is the case with the 2012 and 2006 ice layers, the grain boundaries can have particularly large rod-like brine inclusions.

Fig. 10. Proposed formation mechanism of inclusions at the last region of refreezing in the two thick ice layers. The sequence of inclusion formation is summarized in the last (green-framed) drawing. White lines are grain boundaries or triple junctions.

3.3.3 Ion amount needed to form a large rod-like brine inclusion

To form a large rod-like brine inclusion, a certain amount of ion flux is needed. The major ion species in a large rod-like brine inclusion are Na⁺ and Cl⁻ (Fig. 9). To evaluate the ion amount, three parameters should be considered: (i) the concentration of impurities in the precipitation, (ii) how ion levels change as the meltwater percolates to the ice layer, and (iii) enrichment in the liquid phase during refreezing. As an example, we consider the large rod-like inclusion in 2012 (Fig. S3a-(1)) with the inclusion size as a cuboid of about 1600 μm in diameter and about 50 μm in width. The elemental mass ratios of Na and Cl compared to all seven elements (Si, Al, S, Cl, Na, Mg, and Ca) are 15.9 and 19.0%, respectively. With these numbers, we estimate that the inclusion contains 1.6*10⁻⁶ g Na and 1.9*10⁻⁶ g Cl.

For (i), the concentration of impurities in the precipitation, Na⁺ and Cl⁻ have average concentrations of 4.51 and 5.29 μeq L⁻¹, respectively, averaged over the entire core (Fig. 1). For (ii), the percolation process, the median values of the ion flushing ratio to ice layers from the upper firn are 1.26 for Na⁺ and 1.41 for Cl⁻, respectively (Fig. 2). Multiplying the ratio by the original concentrations, any ice layer in this ice core should have a potential of about 5.66 and 7.47 μeq L⁻¹ for Na⁺ and Cl⁻ before refreezing.

Finally, we evaluate (iii), enrichment in liquid phase during refreezing. To aggregate the estimated 1.6*10⁻⁶ g Na in the inclusion from the initial water concentration of 5.66 μeq L⁻¹, 12.31 g of water is required. For Cl, 7.17 g of water is required. Assuming uniform freezing in one direction as sketched in Fig. S4, the thickness of ice needed for the Cl case is about 7 cm.

This threshold thickness should depend on characteristics of the core, but in this case, the value is roughly consistent with nearly all rod types found in the 15.0 and 15.5 cm thick layers of 2012 and 2006 except those in thin ice layers. For example, in the earlier warm period of the 1920s, the thickest ice layer was 0.7 cm, and here we chose a 0.5 cm ice layer at random and detected in it no rod-like inclusions.

3.4 Implications of the large rod-like inclusion

Large rod-like inclusions such as those found here have not been reported in ice cores that had never experienced melting. A possible reason for their not appearing in such dry ice cores is that they might only appear in the last regions to freeze of thick ice layers.

Concerning their implications, a high concentration of such inclusions in grain boundaries would decrease the pressure melting point and lead to melting in veins along grain boundaries (Wolff and Paren, 1984). In this way, they could change the electrical properties and thermodynamics of the ice. In addition, Fisher and Koerner (1986) and Paterson (1991) suggest that a high concentration of inclusions produces a higher deformation rate of the ice. Thus, large rod-like inclusions in grain boundaries are a place of special vulnerability in the ice layer, which could affect the motion of the entire glacier or ice sheet.

As a proxy, the rod-like inclusions may be an indicator that a large amount of meltwater has refrozen. Such a proxy may be useful for identifying past warm climates (e.g., MIS5e) in deeper ice cores in the Greenland ice sheet. Furthermore, the millimeter-scale rod-like brine inclusions might become a weak point in the refrozen ice against the Arctic warming. The brine can then change into water veins in the refrozen ice. In this way, the millimeter-scale brine may accelerate the melting of the ice sheet. Further studies should help clarify these possible implications.

4. Conclusion

We measured the size, shape, location, and elemental composition of inclusions preserved in ice layers formed by melt-refreezing in the Greenland ice sheet. In the ice layers, we identified three types of inclusions: (1) Particle-like inclusions in ice grains and grain boundaries that contain solid salts of Na₂SO₄⋅10H₂O and CaSO₄⋅2H₂O, as well as brine with NH₄⁺ and SO₄²⁻. (2) Rod-like inclusions in grain boundaries that contain brine with SO₄²⁻, Cl⁻, Na⁺, and NH₄⁺. (3) Columnar-like inclusions within ice grains that consist of CaSO₄⋅2H₂O (gypsum). We argued that the inclusions over 30 μm in diameter had formed by melt-refreezing.

In thick ice layers that formed in 2012 (4.330–4.480 m depths) and 2006 (6.056–6.211 m depths), we discovered particularly large rod-like inclusions up to about 1 mm in diameter at depths corresponding to the last region of refreezing. In this region, we argued that the process followed the following sequence. (1) First, the smaller columnar-and particle-like inclusions of solid Na₂SO₄⋅10H₂O and CaSO₄⋅2H₂O form within the ice grains. (2) Then, particle-like inclusions of brines with Na⁺, Cl⁻, NH₄⁺, and SO₄²⁻ form in grain boundaries. And finally, in the case that the meltwater has sufficient impurity, then (3), large rod-like inclusions of brines with Na⁺, Cl⁻, NH₄⁺, and SO₄²⁻ form in grain boundaries as the final product. We suggested that the large rod-like inclusions could be used to indicate that a large amount of meltwater refroze. Such a proxy may help identify past warm climates in deeper ice cores in Greenland and for studying future ice sheet melting behavior.