Fractionation of O2/N2 and Ar/N2 in the Antarctic ice sheet during bubble formation and bubble-clathrate hydrate transition from precise gas measurements of the Dome Fuji ice core

The variations of δΟ2/Ν2 and δΑr/N2 in the Dome Fuji ice core were measured from 112 m (bubbly ice) to 2001 m (clathrate hydrate ice) at high precision. Our method, combined with the low storage temperature of the samples (-50 °C), 15 successfully excludes post-coring gas-loss fractionation signals from our data. From the bubbly ice to the middle of the bubbleclathrate transition zone (BCTZ) (112 – 800 m) and below the BCTZ (>1200 m), the δΟ2/Ν2 and δΑr/N2 data exhibit orbitalscale variations similar to local summer insolation. The data in the lower BCTZ (800 – 1200 m) have large scatters, which may be caused by mm-scale inhomogeneity of air composition combined with finite sample lengths. The insolation signal originally recorded at the bubble close-off remains through the BCTZ, and the insolation signal may be reconstructed by 20 analyzing long ice samples. In the clathrate hydrate zone, the scatters around the orbital-scale variability decrease with depth, indicating diffusive smoothing of δΟ2/Ν2 and δΑr/N2. A simple gas diffusion model was used to reproduce the smoothing and thus constrain their permeation coefficients. The relationship between δAr/Ν2 and δΟ2/Ν2 is markedly different for the datasets representing bubble close-off (slope ~0.5), bubble-clathrate hydrate transformation (~1), and post-coring gas-loss (~0.2), suggesting that the contribution of the mass-independent and mass-dependent fractionation processes are different for those 25 cases. The method and data presented here may be useful for improving the orbital dating of deep ice cores over the multiple glacial cycles and further studying non-insolation-driven signals (e.g., atmospheric composition) of these gases.

The size-dependent gas fractionation may be related to the size of channels within ice crystals (Huber et al., 2006) or the differences in the dominant diffusion mechanism . Deep ice cores from Antarctic inland show δO2/N2 of -5 to -10 ‰ relative to the atmosphere due to the close-off fractionation (Bender, 2002;Suwa and Bender, 2008a, 35 b; Kawamura et al., 2007;Extier et al., 2018). δAr/N2 has been rarely reported, and it is less depleted than δO2/N2 (Sowers et al., 1989;Bender et al., 1995;Severinghaus et al., 2009).
High correlations were found between local summer insolation and δO2/N2 records of the Vostok and Dome Fuji (DF) cores, which have been used for orbital dating of the cores (Bender, 2002;Suwa and Bender, 2008a;Kawamura et al., 2007). Summer 40 insolation may influence the physical properties of snow, which later control the close-off fractionation (Bender, 2002;Kawamura et al., 2007;Fujita et al., 2009). The phasing between the local summer insolation and δO2/N2 signals and its variability are typically neglected in the dating practice. However, the δO2/N2 data from the EDC, Vostok and DF cores around the last interglacial period show discrepancies and high-frequency variabilities, challenging the assumption of stable phasing (Bazin et al., 2016). 45 The current understanding of the gas fractionation and hence the assessment of δO2/N2 and δAr/N2 as the local insolation proxy are hindered, at least partly, by the scarceness of the data that are unaffected by post-coring gas loss (during ice-core drilling and storage), which also fractionates δO2/N2 and δAr/N2 (Bender et al., 1995;Kawamura et al., 2007;Kobashi et al., 2008;Severinghaus et al., 2009;Sowers et al., 1989). For example, configurational diffusion may occur in microcracks to produce 50 size-dependent and weak mass-dependent fractionations (Bender et al., 1995;Huber et al., 2006). Molecular diffusion through ice matrix may also produce size-dependent fractionation . The DF, Vostok and GISP2 deep cores indeed showed gradual δO2/N2 depletion over several years at -25 --35 °C Suwa and Bender, 2008a, b).

55
Another obstacle to understanding the close-off fractionation lies in the transformation of bubbles to clathrate hydrates. In the bubble-clathrate hydrate transition zone (BCTZ: e.g., 450 -1200 m at Dome Fuji, (Narita et al., 1999), extreme gas fractionations on the order of several 100 ‰ can occur between individual bubbles and clathrate hydrates due to lower dissociation pressure and larger permeation coefficient of O2 relative to N2 (e.g., Chazallon et al., 1998;Ikeda-Fukazawa et al., 2001). The microscopic gas fractionation, combined with the post-coring gas-loss fractionation, creates high variabilities of 60 δO2/N2 and δAr/N2 from gas measurements in the BCTZ (Bender, 2002;Kawamura et al., 2007;Kobashi et al., 2008;Lüthi et al., 2010;Shackleton et al., 2019). Below the BCTZ, the variabilities of δO2/N2 and δAr/N2 decrease with depth over several hundred meters, possibly due to diffusive smoothing (Lüthi et al., 2010;Bereiter et al., 2014;Shackleton et al., 2019). The understanding of the smoothing process has been insufficient because of the lack of high-resolution gas data for a proper numerical simulation (Bereiter et al., 2014). Moreover, the fact that reliable δO2/N2 records have been available only from the clathrate hydrate ice (i.e., older than ~100 kyr) has hindered the discussion of how the original close-off fractionation signals are preserved through the BCTZ.
A better understanding of the gas fractionation associated with bubble close-off and clathrate hydrate formation is important for the ice-core dating and eventual reconstruction of the past atmospheric ratios. Bereiter et al. (2009) predicted that the 70 original gas composition is preserved in the inner part of the cores for a few decades at low temperature (e.g., -50 ˚C), based on a numerical simulation of gas diffusion. The prediction was verified by analyzing the DF core, which has been stored at -50 ˚C for ~20 years, after sufficiently removing its fractionated outer part (Oyabu et al., 2020). In this work, we employ the method by Oyabu et al. (2020) and measure δO2/N2 and δAr/N2 from the bubbly ice to clathrate hydrate ice in the DF core. In addition to the discrete sampling, we obtain high-resolution, continuous gas and chemical data at selected depth intervals to 75 understand the smoothing process of δO2/N2 and δAr/N2 below the BCTZ. Moreover, the outer ice is also measured for discussing the artifactual fractionation during storage. We analyze the data mainly based on the relationships between δO2/N2, δAr/N2 and δ 18 O of O2 to distinguish the gas fractionation processes associated with the bubble close-off, transformation of bubbles to clathrate hydrates, and post-coring gas loss. A simple gas diffusion model is also employed to constrain the permeation coefficients of N2, O2 and Ar, and discuss the preservation of insolation signals in the gas records through the 80 BCTZ.

Ice core
We use the first Dome Fuji ice core (DF1 core) drilled from 1993 to 1997 (Watanabe et al., 2003). After transporting to Japan, the samples had been stored at −50 °C at the Institute for Low Temperature Science, Hokkaido University (until 2009) and the 85 National Institute of Polar Research (NIPR, 2009 -present). Within the DF ice core, 105 -450 m is bubbly ice, 450 -1200 m is bubble-clathrate transition zone (BCTZ), and below 1200 m is bubble-free ice (Narita et al., 1999;Ohno et al., 2004).

Air extraction and measurements
A typical ice sample is ~60 g and ~110-mm long (corresponding to ~10 years) (Fig. 1a), which is cut out from a bulk section (~500-mm long). Because δO2/N2 and δAr/N2 near the sample surface are depleted by gas loss during storage, we removed the 90 surface by more than 8 mm for clathrate ice (Oyabu et al., 2020) and 5 mm for bubbly ice (see below). For some samples in the BCTZ (1000 -1200 m), we removed the surface by 20 mm (Fig. 1b) to minimize the signal of gas loss from high-pressure bubbles (Bender, 2002).
The thickness of surface removal (5 mm) for bubbly ice was determined as follows. Within a pair of samples with 3-mm 95 removal, the sample with larger outer surface ("a2" in Fig. 1a) had generally lower δO2/N2 and δAr/N2, and slightly higher δ 18 O than the other sample ("a1" in Fig. 1a) (Fig. B1). No such tendency is observed with >5 mm surface removal (Fig. B1).

100
The experimental procedures at NIPR are described elsewhere (Oyabu et al., 2020). Briefly, ice samples were evacuated for 2 hours and melted, and the released air was cryogenically collected into a sample tube. The air sample was split into two aliquots and measured with a mass spectrometer (Thermo DELTA V, for δO2/N2, δAr/N2, δ 15 N and δ 18 O) and gas chromatographs (two Agilent 7890A, for CO2, CH4 and N2O concentrations; not shown in this study). Corrections were made for non-linearity of the elemental and isotopic ratios upon sample pressure, and for the dependency of δ 15 N, δ 18 O and δAr/N2 upon δO2/N2 (Bender 105 et al., 1994b). Then, the values were normalized to the modern atmosphere (Oyabu et al., 2020). Gravitational corrections were applied to δO2/N2, δAr/N2 and δ 18 O using δ 15 N of the same sample (Sowers et al., 1989). We measured 522 and 162 depths with single and replicates, respectively, between 112.88 to 2001.12 m (3 to 173 kyr BP).
We tested the possibility of gas loss during the sample evacuation (Craig et al., 1988). For a bubbly ice (446.4 m) and a 110 clathrate hydrate ice (2001.1 m), we split each sample into five thin pieces of 13 -18 g (Fig. 1d), and evacuated them for different durations (2 pieces for 20 minutes, 1 for 2 hours, and 2 for 5 hours). There are no dependence of δO2/N2, δAr/N2, δ 15 N and δ 18 O on the pumping time (Fig. B2). The bubbly ice data show somewhat low δO2/N2 and δAr/N2 for the 20-min evacuation, but the result is opposite to the expectation from gas loss. Also, larger variabilities of δO2/N2 and δAr/N2 are expected in the bubbly ice than in the clathrate hydrate ice (section 4.1). Thus, our method and core quality do not significantly 115 fractionate the gases.
Eleven samples from 179.69 to 2496.55 m, stored for 11 years at -25 ˚C, were measured for δO2/N2 and δ 15 N at Tohoku University. The analytical method was slightly modified from the previous method by Kawamura et al. (2007) as follows. The typical sample size was 250 g, and the measurements were made with a Delta plus XP (Thermo Fisher Scientific) with 64 120 changeover cycles. Analytical precisions (1σ) were estimated to be 0.4 and 0.011 ‰ for δO2/N2 and δ 15 N, respectively.

High-resolution analyses
The air composition and ion concentrations were measured on five 50-cm segments at the resolutions of 1 -2 cm. Three segments were measured for the air (1258.51 -1258.99 m, 1389.79 -1390.32 m and 1893.51 -1893.99 m), and two segments 125 were measured for the air and ions (1292.29 -1292.82 m and 1399.03 -1399. 48 m). In Fig. 1c, the X and Y sections were divided into 25-and 12.5-mm pieces for air and ion measurements, respectively. For the ion measurements, the ice was decontaminated by shaving-off by ~3 mm, and then chopped off with a ceramic knife and sampled in particle-free plastic bags. They were melted in a clean laboratory (class 10,000) and analyzed for the 130 concentrations of Cl -, SO4 2-, NO3 -, F -, CH3SO3 -, Na + , NH4 + , K + , Mg 2+ and Ca 2+ using ion chromatography (Thermo Fisher Scientific Dionex ICS-5000+) (Goto- Azuma et al., 2019).

Diffusion model
We simulate diffusive smoothing of δO2/N2 and δAr/N2 in the clathrate hydrate ice with a one-dimensional diffusion model 135 (Ikeda-Fukazawa et al., 2005;Bereiter et al., 2009;Bereiter et al., 2014). The model assumes that the molecular diffusion through ice lattice is driven by the concentration gradient of gas molecules dissolved in ice, which are in equilibrium with clathrate hydrates in respective layers of the model. Permeabilities of O2 and N2 are poorly known, thus we used three sets of There are no published permeability of Ar in ice, thus we used two formulations proposed by Kobashi et al. (2015) (Fig. 2).
The first permeability !"($) uses diffusion coefficients of N2, O2 and Ar at 270 K from the molecular dynamics simulations (1) 145 The second permeability !"($$) is based on the observations that δAr/N2 is often depleted about half of δO2/N2 (e.g., Severinghaus et al., 2009), and is given by: (2) To take into account the actual ice sheet condition, the temperature and grid size were revised every ~100 model-years according to the temperature and thinning function. Initial depth profiles of δO2/N2 and δAr/N2 were constructed by repeating 150 the 50-cm segment of the 2.5-cm continuous data at 1258 m (81.7 kyr BP) (Fig. 6). We run the model from 1258 to 2000 m (for ~100 kyr) to cover the depth range of our data. Further details of the model are described in Appendix A.
Variations in δO2/N2 and δΑr/N2 for the depths shallower than ~800 m and deeper than ~1200 m have similarity with local 160 summer insolation curve, while little similarity is found for 800 -1200 m with extremely large scatters (Fig. 4). We assess the scatters in the data by taking residuals of δO2/N2 and δAr/N2 from their low-pass filtered curves ( Fig. 3d and 3e). The average residuals of δAr/N2 also show similar pattern (0.5 ‰ for 112 -450 m, 0.5 ‰ for 450 -800m, 2.8 ‰ for 800 -165 1200 m, 1.3 ‰ for 1200 -1480 m, and 0.3 ‰ for ~1480 -2000 m). Below, we divide our dataset into four depth ranges guided by the scatters as well as the classic boundaries between bubbly ice zone, BCTZ and clathrate ice zone (Narita et al., 1999;Ohno et al., 2004).
Figs. 2f and 2g show differences of δO2/N2 and δAr/N2 between ice samples cut from the same depth ( Fig. 1a) (hereafter referred to as pair difference or ΔδO2/N2 and ΔδAr/N2). Ranges of ΔδO2/N2 and ΔδAr/N2 are ~0.02 to 0.97 ‰ and 0.01 to 0.96 ‰, respectively. Pair differences are usually used to calculate the pooled standard deviation, a metric of analytical precision. Instead, we use it to evaluate spatial variability in air composition for a given depth in the ice sheet (pair differences 175 of δ 15 N and δ 18 O are plotted in Fig. B3 and pooled standard deviations are summarized in Table C1).
Our new δO2/N2 data are compared with previously obtained data at Tohoku University using the samples stored at -25 °C for various durations (Fig. 5). As expected, the samples stored at -25 °C generally show lower δO2/N2 than those stored at -50 °C, with dependency on the storage duration. Also, while the data from the -50 °C samples slightly increase with depth in the 180 bubbly ice zone as expected from the insolation variation, the data from the -25 °C samples tend to decrease with depth possibly due to the increase in bubble pressure (leading to larger gas-loss) (Ikeda-Fukazawa et al., 2005;Kawamura et al., 2007).

450 -800 m (upper BCTZ)
The number and size of clathrate hydrates gradually increase with depth, but air bubbles are still the dominant form of air 185 inclusion (Fig. 3i, 3j) . In the lowermost part of this range (below 720 m), clathrate hydrates with extremely enriched δO2/N2 (> ~1000 ‰) are found by Raman spectroscopy (Fig. 3h)  δΑr/N2 values from the gas analyses smoothly connect with those in the bubbly ice zone, with a slightly increasing trend with depth ( Fig. 3b, 3c). The pair differences of δO2/N2 and δΑr/N2 are smaller than 1 ‰ until ~650 m, below which they show high values (> 3 ‰, Fig. 3f, 3g). 190

800 -1200 m (lower BCTZ)
Major transition from air bubbles to clathrate hydrates occur in this depth range (Fig. 3i, 3j) . The δO2/N2 values of clathrate hydrates decrease from extremely high values towards the atmospheric value, and those of remaining air bubbles decrease to extremely negative values (~ -500 ‰) ( Fig. 3h) (Ikeda-Fukazawa et al., 2001). 195 In contrast to the upper BCTZ, the bulk δO2/N2 and δΑr/N2 in the lower BCTZ show elevated scatters ( Fig. 3b-e). Very high δO2/N2 scatters have also been reported from the Vostok (Suwa and Bender, 2008a), GISP2 (Kobashi et al., 2008) and WAIS Divide (Shackleton et al., 2019) ice cores, with numerous positive values (higher than the atmosphere). The previous data were affected by the post-coring gas loss, thus they were interpreted only as artifacts created by the preferential gas loss from 200 extremely N2-rich air bubbles (Bender, 2002;Ikeda-Fukazawa et al., 2001;Kobashi et al., 2008;Shackleton et al., 2019). In our data, relatively few samples (six and ten samples for δO2/N2 and δΑr/N2, respectively) show positive values, and even the samples with twice the surface removal (20 mm) show large scatters. Thus, the large δO2/N2 and δΑr/N2 variabilities naturally occur in this zone. The pair differences of δO2/N2 and δΑr/N2 are also large (up to ~10 ‰, Fig. 3g, 3h; Table C1), indicating large inhomogeneities also in the horizontal direction. 205
Below 1480 m, the δO2/N2 and δΑr/N2 variations show remarkable similarities with the local summer insolation with small 215 scatters (Fig. 4). The pair differences are also small ( Fig. 3f, 3g, Table C1). The previous δO2/N2 data corrected for the gasloss fractionation (measured at Tohoku University; Kawamura et al., 2007) agree well with our new data (Fig. 5), suggesting that the inner part of the ice core stored at -50 °C has kept the original δO2/N2 for two decades.

Outer ice 220
In almost all outer ice samples, δO2/N2 and δAr/N2 are significantly lower and δ 18 O is significantly higher than those in the inner ice, as expected for gas loss (Ikeda-Fukazawa et al., 2005;Severinghaus et al., 2009). The outer ice stored at -30 ˚C for ~1 year shows larger gas-loss signals than those stored at -50 ˚C (Fig. B4), indicating that the gas-loss fractionation strongly depends on the storage temperature and duration. In the lower BCTZ, some δO2/N2 and δAr/N2 values from the outer ice are positive ( Fig. B4a, B4b), possibly due to the preferential gas loss from N2-rich bubbles (Bender et al., 1995). We do not find 225 significant differences between the outer and inner pieces for δ 15 N, CH4 concentration and N2O concentration ( Fig. B4d -f).
The results are consistent with the earlier notion that only the molecules smaller than 3.6 Å in collision diameter (O2 and Ar among the above species) significantly fractionate (Huber et al., 2006;Severinghaus and Battle, 2006).

235
The data are averaged over 4 to 5 consecutive samples and compared with the data from the normal (11 cm) samples ( Fig. 6).
Two normal samples at ~1293 and 1390 m are taken from the same depths as the high-resolution measurements, and they agree with the averages of the high-resolution data. The standard deviation of the averaged data (blue curves in Fig. 6) for the upper four depths are 1.4 to 2.8 ‰ for δO2/N2 and 0.8 to 1.5 ‰ for δΑr/N2, which are significantly higher than that for the deepest sample (0.1 and 0.2 ‰). High variabilities at 10-cm scales have also been found in the continuous measurements of 240 the GRIP core (in the BCTZ, Huber and Leuenberger, 2004) and EDML core (just below the BCTZ, Lüthi et al., 2010). Also, in the GRIP data, the variabilities of δO2/N2 well below the BCTZ (~2500-m depth) are much smaller than those in the lower BCTZ (~1100 -1400 m) (Huber and Leuenberger, 2004).
Major ion concentrations (Na + , Mg 2+ , Ca 2+ , Cl -, SO4 2and NO3 -) at 12.5-mm resolution from the two depths ( smooth variations possibly due to their higher mobility in firn and ice. The correlation coefficients between the ion 250 concentrations and gas ratios are summarized in Table 1. For calculating the correlation coefficients, the ion data are resampled to the depths of the gas data. Significant correlations with the gas data are found for the Na + , Mg 2+ and Ca 2+ concentrations.

Fractionation in firn
The δO2/N2 and δΑr/N2 data for the bubbly ice, upper BCTZ and below BCTZ show variations similar to the local summer 255 insolation (Fig. 4). For the first time, we find the insolation signals in the top part of the ice sheet, supporting the link between the local summer insolation and close-off fractionation through the effects on the snow metamorphism (Bender, 2002;Fujita et al., 2009).
In the bubbly ice zone, the scatters of the data are significantly larger than those of the pair differences, probably reflecting 260 high variabilities of natural close-off fractionation in the vertical direction. In the Dome Fuji region, firn layerings of several 10 cm are found in δ 18 O of ice, ion concentrations and density (Hoshina et al., 2014;Fujita et al., 2016). The layers closed-off deeper may be less depleted in Ο2 and Ar because (1) air bubbles spend a shorter time in the close-off zone, and (2) δΟ2/Ν2 and δAr/Ν2 are enriched deeper in the open pores (e.g., Severinghaus and Battle, 2006). Also, microbubbles may form near the surface (Lipenkov, 2000) and become extremely depleted in Ο2 and Ar (Ohno et al., 2021), also causing the δO2/N2 and 265 δAr/N2 variations (Kobashi et al., 2015). Hereafter, we collectively call the natural fractionation during the firn densification and bubble formation as "close-off fractionation".
As discussed above, the relationship between δAr/N2 and δO2/N2 (after gravitational correction using δ 15 N) in our data is expected only to reflect natural fractionations. For bubbly ice, we find a high correlation between δAr/N2 and δO2/N2 with the 270 slope of 0.50±0.01 (Fig. 8, Table C2), which agrees with that of pair differences (0.53±0.04) (Fig. 9, Table C2). The smaller fractionation for δAr/N2 is consistent with the size-dependent fractionation. We find similar slopes in the upper BCTZ, lower BCTZ and below the BCTZ (0.45±0.01, 0.61±0.01, 0.42±0.02, respectively, Table C2), despite large scatters in the lower BCTZ and just below the BCTZ (~800 -1480 m). We interpret the similarity of the slopes as preservation of the close-off fractionation signals through the BCTZ. The slight differences between the slopes in the different zones might arise from 275 different temperature or accumulation rate in the past (affecting the close-off fractionation), or natural variations of O2 in the atmosphere.
A firn-air study at WAIS Divide (Battle et al., 2011) provided the evidence of mass-dependent fractionation associated with the close-off process with the slope of δ 18 O vs. δO2/N2 of -0.0090 (‰ / ‰). For the DF ice-core data, the pair differences in 280 the bubbly ice zone (Fig. B5), as well as the data for the last 2000 years (Oyabu et al., 2020)  between δ 18 O and δO2/N2. Thus, our data cannot verify the mass-dependent fractionation, and much higher precision is required to detect the δ 18 O change of ~ 0.01 ‰ expected for the range of DF δO2/N2.
From the comparison of three Antarctic inland cores (EDC, Vostok, DF) that showed large short-term δO2/N2 variabilities and 285 discrepancies between the cores around the last interglacial period, Bazin et al. (2016) discussed the possibility of non-orbitalscale variabilities of close-off fractionation. In our data, the residuals of δO2/N2 from the filtered curve are small (1σ = ~0.5 ‰) below ~1480 m (Fig. 3d), whereas the previous data (after gas-loss correction) show the variability of 1.3 ‰ for the similar depth range   (Fig. 5). Thus, the large short-term variabilities in the previous datasets may be mostly attributed to the poor ice core quality, and the actual short-term variability of close-off fractionation may be rather small, at 290 least for the DF core.

Bubble-clathrate transformation fractionation
In the upper BCTZ, the residuals of δO2/N2 and δAr/N2 from the low-pass filtered curves exhibit scatters of ~3 ‰ and ~2 ‰ (peak-to-peak), respectively, which are much smaller than those in the lower BCTZ (Fig. 4). Pair differences are similar to 295 those at the bottom of the bubbly ice zone, suggesting insignificant effects of the clathrate hydrate formation on the measured δO2/N2 and δAr/N2 until ~40 % of air bubbles are transformed to clathrate hydrates. This may be reasonably explained by the dominance of direct conversion of air bubbles to clathrate hydrates (thus little displacements of molecules) in the early stages of bubble-clathrate transformation (Lipenkov, 2000). It may also be the case that the distance of gas diffusion from the bubbles to clathrate hydrates, even if it occurs, is too short to create scatters in the bulk δO2/N2 and δAr/N2. 300 In the lower BCTZ, the scatters around the orbital scale variations dramatically increase for δO2/N2 and δAr/N2 despite the sufficient removal of the outer ice (Fig. 3, 4). Microscopic observations by Ohno et al. (2004) found layered distributions of air bubbles and clathrate hydrates in the lower BCTZ, as well as high spatial variability of total number of air inclusions on a few mm-scales, possibly due to diminishing bubbles by transferring their molecules to nearby clathrates. Using their number 305 concentrations of air bubbles and clathrate hydrates, the average distance between air inclusions is estimated to be 1.1 mm.
Thus, the distance of air migration associated with the bubble-clathrate transformation should be on this order. Ikeda-Fukazawa et al. (2001) observed extremely fractionated bubbles and clathrate hydrates in the lower BCTZ, up to about +1000 ‰ for clathrate hydrates and -740 ‰ for bubbles for δO2/N2. From these observations, we suggest that the highly fractionated bubbles and clathrate hydrates may be stratified in mm-scale layers, and that the scatters in our dataset may be produced by random 310 inclusion of such fractionated layers at the top and/or bottom of the ice samples.
The slopes of δAr/N2 vs. δO2/N2 from the pair differences are close to 1 for both the upper and lower BCTZ ( Fig. 9c and 9d, Table C2). They are strikingly different from the that in the bubbly ice (Fig. 9b, Table C2) and should be created by the fractionations associated with bubble-clathrate transition, because the amplitudes of ΔδO2/N2 and ΔδAr/N2 are larger than 315 those in the bubbly ice by an order of magnitude. We also investigate the pair differences of δ 18 O vs. δO2/N2, which may be an indicator of mass-dependent fractionation. We find no correlation between Δδ 18 O and ΔδO2/N2 with relatively wide ranges of ΔδO2/N2 (up to 10 ‰, Fig. B5), indicating that the fractionation by the bubble-clathrate transformation is not massdependent. This is clearly different from the known correlations between Δδ 18 O and ΔδO2/N2 for bubble close-off (Battle et al., 2011) and post-coring gas loss (Severinghaus et al., 2009). Thus, we suggest that the net fractionation associated with the 320 bubble-clathrate transformation and molecular diffusion between air inclusions within the ice sheet is mass-independent, and that the bubble close-off and post-coring gas-loss fractionations include mass-dependent processes for molecules smaller than 3.6 Å such as He, O2 and Ar (Craig and Scarsi, 1997;Severinghaus et al., 2009). For the BCTZ, the slower permeation of Ar than that of O2 may be counteracted by the lower dissociation pressure of Ar than O2, which should produce a steeper Ar partial pressure gradient from the bubbles to the growing clathrate. This hypothesis may explain the observed similarity of δAr/N2 325 and δO2/N2 enrichment in clathrates. For the bubble close-off, such a cancellation cannot occur, explaining the observed δAr/N2 that is only half as enriched as δO2/N2. Instead, the bubble close-off fractionation of small molecules (<3.6 Å molecular diameter) appears to be caused by permeation through the thin ice wall of freshly-formed bubbles via both the bond-breaking mechanism (mass-independent) and the interstitial mechanism (mass-dependent; Ikeda-Fukazawa et al., 2004;Severinghaus and Battle, 2006;Battle et al., 2011). 330 We find similar slopes (~1) of ΔδAr/N2 vs. ΔδO2/N2 in the BCTZ of the WAIS Divide (Seltzer et al., 2017) and South Pole (Severinghaus, 2019) ice cores, although the samples have possibly experienced post-coring gas loss to some extent (Table   C2). The fact that the different cores from a wide range of temperature show similar relationships between ΔδAr/N2 and ΔδO2/N2 may suggest that the permeation coefficients of O2 and Ar in ice have similar dependence on temperature. 335

Preservation of insolation signal through BCTZ
Below the BCTZ, the low-pass filtered δO2/N2 and δAr/N2 show close resemblance to the variations of local summer insolation, although the residuals are rather scattered until 1480 m (Fig. 3, 4). The large scatters probably originate in the extreme layered fractionation in the lower BCTZ. Until the conversion of bubbles to clathrate hydrates complete (at ~1200 m), gas flux from 340 the remaining (extremely fractionated) bubbles to clathrate hydrates should continue and enhance the inhomogeneity of gas composition. Once all bubbles disappear, the gas flux occurs only between clathrate hydrates; thus, their compositions become gradually homogenized. The major homogenization finishes within 300 m (25,000 years) below the BCTZ, and further homogenization continues to ~500 m (~50,000 years) below the BCTZ where the scatters become small and stable (±0.2 ‰) (Fig. 3d, 3e and 12). 345 In contrast to the scatters from the 11-cm samples discussed above, the pair differences decrease sharply at the bottom of BCTZ (Fig. 3f, 3g), suggesting smaller spatial variability of gas composition in the horizontal direction. The 2.5-cm continuous data from the four 50-cm samples (1258 to 1399 m) also show much higher variabilities (5 -10 ‰ difference between neighboring samples) than the pair differences (Fig. 3f, 3g and 6). Thus, our data consistently show the large anisotropy of gas 350 fractionation, as expected from the layered clathrate hydrate formation in the BCTZ . The positive correlations between the gas ratios and Na + , Mg 2+ and Ca 2+ concentrations are consistent with the suggestion by Ohno et al.
(2010) that nucleation of clathrate hydrates tends to occur in the layers with abundant micro-inclusions, because the elements Na, Mg and Ca likely exist as water-soluble salts (Oyabu et al., 2014 and references therein). Thus, the high-micro-inclusion layers may attract O2 and Ar from air bubbles in the adjacent layers with fewer micro-inclusions. The same conclusion was 355 reached by Huber and Leuenberger (2004) (albeit with ~50 times larger measurement uncertainties than ours), based on the quasi-annual cycles in the δO2/N2 and δΑr/N2 in the GRIP ice core, Greenland, and their similarity with Ca 2+ concentrations.
The molecular diffusivities should depend on the direction in the ice sheet and also contribute to the horizontal homogeneity.
The c-axes orientation gradually becomes clustered around the vertical, and the crystal grains become horizontally elongated, 360 by the vertical compression within the ice sheet (Azuma et al., 2000). In terms of molecular diffusion through ice matrix, the diffusivity in the direction perpendicular to the c-axis is greater than that in the direction parallel to the c-axis . For the diffusion along grain boundaries, the pathway should be shorter in the horizontal direction due to the crystal elongation. In addition, lattice defect and dislocation within ice crystals tend to develop along the basal plane ), possibly providing another pathway for gas diffusion. All these mechanisms should contribute to preferential 365 homogenization of gases in the horizontal direction.
The results of the diffusion model for 2.9, 13.7, 14.5 and 67.5 kyr (at 1292, 1390, 1399 and 1894 m, respectively) after the initial state (81.7 kyr BP at 1258 m, Fig. 6) are resampled at 2.5-cm intervals and compared with the data (Fig. 10). The IkFk05 permeation parameters give the smoothest δO2/N2 profiles with a poor agreement with the data (Fig. 10c). The IkFk01 370 parameters produce the profiles similar to the data at 1292 m, but the model results are too smooth at 1390 and 1399 m (Fig.   10b). The results with the Salm01 parameters agree well with the data, including rapid changes (>10 ‰) within a few consecutive samples and the standard deviation (~3 ‰) at 1399 m (Fig. 10a). We note that the data at 1390 and 1399 m show rather different standard deviations (5.8 and 3.1 ‰) probably reflecting the original fractionations in the BCTZ. Thus, the model-data comparison is inadequate with the small number of cases. For Ar, the model results with the scaling function Ar(II) 375 are closer to the data than those with Ar(I) (Fig. 11).
The standard deviations of the model results resampled at 11-cm intervals are compared with the residual δO2/N2 and δAr/N2 data from the low-pass filtered curves (Fig. 12) (following Bereiter et al., 2014). Exponential fitting curves through the residual data (black line) are in close agreements with the model results with the Salm01 and Ar(II) permeation parameters. On the 380 other hand, the model results with the other parameters (IkFk01, IkFk05 and Ar(I)) show too rapid decrease in comparison with the data. Therefore, our datasets (high-resolution and normal datasets) consistently support the Salm01 and Ar(II) permeation parameters at around 240 K (temperature at DF for the simulated depths).
From the Salm01 parameters, the rate of diffusive migration is on the order of 0.1 mm per 10 kyr (10 -10 m s -1 ). Therefore, we 385 favor the interpretation that the extreme scatter of δO2/N2 and δAr/N2 in the BCTZ in our datasets are caused by mm-scale inhomogeneity of the compositions of air inclusions combined with the finite sample length, rather than by cm-scale migration of gas molecules. We also suggest that the original insolation signal on δO2/N2 and δΑr/N2 in the BCTZ may be reconstructed by analyzing long ice samples (>50 cm) to average out the inhomogeneity.

Gas-loss fractionation
In the outer ice, the δO2/N2 and δΑr/N2 are significantly depleted and δ 18 O is enriched due to post-coring gas loss, with higher storage temperature leading to larger fractionation (Fig. B4). The slope of ΔδAr/N2 vs. ΔδO2/N2 is 0.22 (Fig. 9f, Table C2), which is significantly smaller than those from the pairs of inner ice, suggesting that the process for gas loss is different from natural fractionation processes within the ice sheet. The slope of Δδ 18 O vs. ΔδO2/N2 is -0.0083 (Fig. B5), which is similar to 395 those previously reported for artifactual gas loss from the Siple Dome (Severinghaus et al., 2009), WAIS Divide (Seltzer et al., 2017) and EDC (Extier et al., 2018) ice cores. Severinghaus et al. (2009) extensively measured the Siple Dome ice core and found a tight relationship between ΔδAr/N2 and ΔδO2/N2 with the slope of ~0.5, which was also typical for other ice cores that experienced large gas loss (Sowers et al., 1989;400 Bender et al., 1995). The slope of ~0.5 may be the combination of size-dependent fractionation (e.g., 1:1 for δO2/N2 and δAr/N2) and mass-dependent fractionation (1:3 for δO2/N2 and δAr/N2) (Severinghaus et al., 2009). The slope for the DF gasloss fractionation (~0.2) is significantly smaller than most of the previous values, implying that mass-dependent fractionation may be more important for the storage condition of the DF core. Significant mass-dependent fractionation is expected for gas loss through micro-cracks (Bender et al., 1995). Thus, there may have been micro-cracks in the outer ice, although we did not 405 observe visible cracks. The insensitivity of δ 15 N to gas loss, consistent with earlier findings, may be explained by viscous flow for N2, or molecular-size-limited diffusion through ice matrix (Severinghaus et al., 2009).

Conclusions
The variations of δΟ2/Ν2 and δΑr/N2 within the ice sheet from bubbly ice to clathrate hydrate ice at Dome Fuji are reconstructed without post-coring gas-loss signals. Their variations in the bubbly ice zone, upper BCTZ and below BCTZ show close 410 similarity to the local summer insolation. Our δΟ2/Ν2 data from the clathrate hydrate ice zone agree with the previous data after the gas-loss correction , with much less scatters, demonstrating that the original air composition is preserved in the ice core stored at -50 ˚C and may be measured with our technique for precise orbital tuning. In the future, the high-precision δΟ2/Ν2 and δAr/N2 data may also provide non-insolation signals on the gases and eventually be useful for the reconstructions of past atmospheric oxygen and argon concentrations. 415 The large scatters in the lower BCTZ may be created by mm-scale vertical inhomogeneity of air composition combined with finite sample length. The insolation signal originally recorded at the bubble close-off remains through the BCTZ, and the insolation signal may be reconstructed by analyzing long ice samples (>50 cm).

420
Below the BCTZ, the scatters around the orbital-scale variability decrease with depth. The high-resolution analyses of five 50cm segments show decreasing cm-scale variability of δΟ2/Ν2 and δΑr/N2 with depth, suggesting diffusive smoothing. A onedimensional gas diffusion model reproduces the smoothing in this zone with the permeation coefficients of Salamatin et al. (2001).

425
The slope of the pair differences of δAr/Ν2 to δΟ2/Ν2 is about 1 in the BCTZ, and no correlation is found between those of δ 18 O and δO2/N2, suggesting that the O2 and Ar fractionations associated with the bubble-clathrate transition are massindependent. On the other hand, the slope for the bubble close-off process is around 0.5, suggesting a combination of massindependent and mass-dependent fractionation for O2 and Ar. For the BCTZ, the slower molecular diffusion of Ar than O2 may be canceled by the lower dissociation pressure of Ar than O2, which should produce a steeper Ar partial pressure gradient 430 from the bubbles to the growing clathrate. For the bubble close-off of small molecules, such a cancellation cannot occur, and both the bond-breaking mechanism and the interstitial mechanism may play roles for the molecular diffusion through the thin ice wall of fresh bubbles. The slope is small (0.2) in ice that experienced large gas loss, suggesting that mass-dependent fractionation may be important for the storage condition of the DF core.

A1 Model construction
We simulate diffusive smoothing of δO2/N2 and δAr/N2 in the clathrate hydrate ice zone with a simple one-dimensional gas diffusion model (Ikeda-Fukazawa et al., 2005;Bereiter et al., 2009;Bereiter et al., 2014). The model assumes that the molecular diffusion through the ice lattice is driven by the concentration gradient of dissolved gas in the ice, which is in equilibrium with 440 clathrate hydrates. The governing equation is (see Table A1 for a full list of symbol meanings) /0 ( given by (Miller, 1969;Kuhs et al., 2000) (Fig. A1) where 3 and 3 are constant and shown in Table 2. Because there are no published values for !" and !" , we found them 450 by fitting the dissociation pressures of argon hydrate vs. temperature measured by Nagashima et al. (2018) with eq. (A3) (Fig.   A1). The !" and !" are 3.63 and 739.5, respectively. The molar fraction of the m-molecule in the clathrate hydrate is given where 3 is the concentration of m-molecule in total air content. 914:"; ) is the concentration of minor gases, which is assumed 455 to be constant and given by 914:"; where ) is concentration of total air content in ice, 3 ) is concentration of m-molecule with the atmospheric ratio in the total air content, is total air content, ? ! * is molar mass of ice (H2O), @7A is molar volume of a gas at standard temperature and pressure, and 3 is the atmospheric ratio of m-molecule. 460 The model uses an initial box height (∆ ) of 5×10 -3 m. The downward diffusive flux ( 3 ) of m-molecule per unit area at the top boundary of i-th box is the product of the diffusivity and concentration gradient: " is relative thinning function (thinning function divided by that at the top of the model domain). By substituting eq. A2 into 465 eq. A6, 3 is expressed as The net flux of m-molecule is and the concentration change of m-molecule in total air content becomes 470 where ∆ is time step (~11.6 days).

A2 Model input
The diffusivity 3 , solubility 3 , or their product ( where 3 ) is a constant, 3 is activation energy of permeation for m-molecule, R is the gas constant. The permeability (m 2 s -1 ) at temperature T (K) and 1MPa of Salamatin et al. (2001) (hereafter Salm01) is given by 490 where 3 '') is dissociation pressure of m-molecule at 220 K. The diffusivity 3 or permeability 3 (m 2 s -1 ) at temperature T

(K) of Ikeda-Fukazawa et al. (2005) (hereafter IkFk05) is given by
The solubility at 1MPa of Ikeda-Fukazawa et al. (2005) is given by 495 where 3 ) is a constant for m-molecule, 3 @ is activation energy of solubility for m-molecule. We used those permeation parameters for our model (parameters are summarized in Table A2 and each permeability is shown in Fig. 2 and Table A3).
There are no published values of !" , thus we estimated it from &' and *' in Salamatin et al. (2001) with two formulations 500 by Kobashi et al. (2015). The first one !"($) uses diffusion coefficients of N2, O2 and Ar at 270 K from the molecular dynamic The second permeability !"($$) is based on the observations that δAr/N2 is often depleted about half of δO2/N2 (e.g., Severinghaus et al., 2009). The permeability of Ar(II) is expressed as 505 Initial values of δO2/N2 and δAr/N2 are made by repeating the 2.5-cm continuous data at 1258 m (81.7 kyr BP) (Fig. 6), and initial total air content is constant over the model domain. Borehole temperature profile (T) (Buizert et al., 2021) and thinning function (Nakano et al., 2016) were used to simulate the evolution of layer thickness and temperature of the model domain (Fig. A2). We run the model for 100 kyr to cover the depth range covered by our data (until about 2000 m). 510 Figure A2: Thinning function (red, Nakano et al., 2016), temperature (blue, Buizert et al., 2021) and age (black, Kawamura et al., 2007) in the ice sheet at Dome Fuji used for the diffusion model. Grey shading indicates the depth range of the model run.

1.55×10 -11
Diffusivity of O2 (m 2 s -1 ) Permeability of N2 (m 2 s -1 MPa -1 )    Figure 1: Typical cross-sectional cutting plan. The outer black line in each panel is the original surface that was exposed to the atmosphere for ~20 years at -50 °C. Dotted lines are the cutting lines for this study.   , and (c) δO2/N2 and δAr/N2 with low-pass filtered curves (cut-off at ~17 kyr, the filter was designed by Kawamura et al., 2007) (red lines).

Figure 6:
High-resolution data of δO2/N2 and δAr/N2 in 5 depth intervals. Black, blue and red lines are the raw data (2.5-cm resolution), 10cm averages and 25-cm averages, respectively. Standard deviations (SD) of the raw 2.5-cm data (black) and 10-cm averages (blue) are shown in each panel. Thick horizontal black bars with gray shadings in some panels are the independent data from the 11-cm samples with 1σ error (measured by the standard protocol).

Figure 12:
Diffusive smoothing of (a) δO2/N2 and (b) δAr/N2 below the BCTZ. Open markers with dashed lines (gray) are the residuals from the low-pass filtered data. Black lines are exponential fits through the residuals for > 5600 years (corresponding to the model's initial age). The model results (standard deviations of 11-cm resampled outputs) with the best permeation parameters are shown by solid red lines. Blue and green lines are the other model results. Bold numbers indicate significant correlations (p < 0.05). Numbers in brackets are calculated using the data from the lower 20 cm (because of the experimental problem in the Ca 2+ data). 780