Microstructures and Fabric Transitions of Natural Ice from the Styx Glacier, Northern Victoria Land, Antarctica

We investigated the microstructures of five ice core samples from the Styx Glacier, northern Victoria Land, Antarctica. Evidence of dynamic recrystallization was found in all samples: those at 50 m mainly by polygonization, and those at 170 m, largely by grain boundary migration. Crystallographic preferred orientations of all analyzed samples (view from the surface) typically showed a single cluster of c-axes normal to the surface. A girdle intersecting the single cluster occurs at 140–170 m with a tight cluster of a-axes normal to the girdle. We interpret the change of crystallographic preferred orientations (CPOs) at <140 m as relating to a combination of vertical compression, and shear on a horizontal plane, and the girdle CPOs at depths >140 m, as the result of horizontal extension. Based on the data obtained from the ground penetrating radar, the underlying bedrock topography of a nunatak could have generated the extensional stress regime in the study area. The results imply changeable stress regimes that may occur during burial as a result of external kinematic controls, such as an appearance of a small peak in the bedrock.


Introduction
Understanding the rheological behavior of Antarctic ice is essential for predictions of sea-level change in the future. Ice in the polar regions flow continuously from the continents to the oceans driven by gravity. The flow of ice is usually controlled by two processes: frictional sliding on the bedrock, which often occurs where the base is at pressure melting or pre-melting conditions [1][2][3][4][5][6][7][8], and the internal plastic deformation of the polycrystalline ice [9][10][11][12][13][14][15].
The plastic deformation of Antarctic ice results in the development of crystallographic preferred orientations (CPOs), which may induce a strong mechanical anisotropy that can enhance or impede the flow of ice and can be used to study the kinematics of flow [16]. The CPOs of ice vary with the kinematics (e.g., compression or shear), stress, temperature, and strain [10,[17][18][19]. In nature, usually three types of ice c-axes fabrics form: single cluster, double clusters, and a girdle. Single and double clusters of ice c-axes can form due to simple shear [20]. Single cluster can also form due to uniaxial compression (e.g., [21]). Girdles have been interpreted as due to uniaxial extension [22]. The ice CPOs are often reorganized by kinematic changes, such as the progressive development of folds, shear zones, and an increased depth of burial (confining pressure) [21][22][23][24][25]; however, studies related to the CPO evolution are limited, especially for natural ice.
Single-maxima clusters are observed in uniaxial compression laboratory experiments only at low temperatures and/or fast strain rates [17,26]. The vast majority of laboratory experiments give an open cone (small circle) CPO in uniaxial compression. Laboratory experiments have shown that a transition from an open cone to a single cluster of c-axes [26][27][28] with increasing stress, decreasing temperature, and increasing strain. This has been interpreted as an increase in pyramidal slip systems by Piazolo et al [27] and to a change in balance between lattice rotation and grain boundary migration recrystallization by Qi et al [28]. A series of simple shear experiments have reported a transition of double cluster to single clustered c-axes with increasing strain, and have documented different patterns of evolution with strain dependent upon deformation temperature [29][30][31]. The transition in the CPOs of natural ice is often ascribed to environmental variables, such as the concentration of air bubbles, layers of volcanic tephra, snow accumulation rates, and the structures of bedrock [32][33][34].
Although recent experimental studies have reported stress and temperature as prevailing factors that lead to the formation CPOs of ice [28,29], microstructural evolution of natural ice, which occurs at lower strain rates, and proceeds to higher strain than experiments, probably involves further complexities that are not adequately understood.
In this study, we report a microstructural evolution of ice in cores drilled from the Styx Glacier in northern Victoria Land, Antarctica, which is a renowned site for snow accumulation. The slow ice flow allows for a high probability of the development of diverse stress regimes. We describe the deformation behavior of the ice with increasing depth based on calculations of physical properties and careful investigations of CPOs using electron backscatter diffraction (EBSD) in a scanning electron microscope. We also compare the deformation features of the sample to conventional ice flow models, which is followed by a discussion of glaciological implications.

Geological Outline and Sample Descriptions
The Styx Glacier in northern Victoria Land lies along the western margin of the Ross Sea ( Figure  1) and is a renowned site for snow accumulation that is minimally influenced by katabatic winds. Hence, several ice cores have been retrieved from the Styx Glacier in recent decades [35][36][37][38]. The Styx Glacier consists of a flat surface with a relief of <100 m over a total area of >100 km 2 . The velocity of its ice surface is ca. 2.3 m yr −1 [39], which was measured from multiple satellite (Sentinel-1) synthetic aperture radar interferometry data (InSAR) (Figure 2). As a part of the Korean Antarctic Research Program based at the Jang Bogo Station, a drilling program led by the Korea Polar Research Institute, was conducted at the Styx Glacier to obtain an ice core that developed during the past 2 kyr [40]. An ice core with a total depth of 210.5 m was collected from the Styx Glacier at 73°51.10S, 163°41.22E, at an elevation of 1623 m, during the 2014-2015 austral summer season. The ice base at the sample site was determined to be at 550 m depth by a ground-penetrating radar (GPR) investigation [41]. Bedrock signals in the GPR data confirmed the existence of a nunatak, located ~4 km southeast (SE) of the sample site. The presence of the nunatak is consistent with the reduction in ice flow velocity and divergent change of the ice flow direction, observed from the satellite data (Figures 2 and 3), as ice flows past the site in the south-southeast (SSE) direction.
Due to low densities at <40 m and rich in cracks at >180 m depth, samples of the ice core were selected from depths of about 50, 80, 110, 140, and 170 m, labeled as SG_50, _80, _110, _140, and _170, respectively [40]. The samples were stored at ca. -25 °C for preserving the fabric (e.g., [42]). Characteristics of the cores are summarized in Table 1. The ages were estimated to be 0.26 to 1.83 kyr, based on a combined approach of an empirical firn densification model [43] and the CH4 gas age with the correction for the difference between the gas and ice ages of ~0.3 kyr [44]; densities ranged from 780 to 920 kg/m 3 (Figure 4a,b). The borehole temperature, measured 2 y after the drilling, was -31.7 °C at a depth of 15 m and −28.3 °C at the base of the hole (210m). A model calculation reported a mean annual temperature of -32.5 °C, and a basal temperature of -20.5 °C at 550m [44]. Snow accumulation rate is reported as 203 kg m −2 yr −1 [38].
Vertical longitudinal strain rates can be estimated by assuming that finite longitudinal strain is accommodated by density increase from the surface. The strain rate separating the sample and the surface is calculated by dividing the natural log of the change in ice density (Δ ) by the age of the sample, with shortening being positive. The age is derived from the depth-age scale [44] (Table 1). Initial density of 350 kg/m 3 was assumed at the surface. Finite longitudinal strain increases to 0.97 at 170 m depth. Longitudinal strain rate decreases from 10 −10 s −1 at shallow depths to 2 × 10 −11 s −1 at 170 m depth (Figure 4c). An upper bound horizontal shear strain rate is estimated by assuming a linear decrease in velocity from the surface to the base of the ice: a horizontal shear strain rate of ~10 −10 has been calculated by dividing the ice flow velocity, 2.3 m year −1 , by thickness, 550 m ( Table 1). A lower bound shear strain rate estimate was made by assuming the driving force for horizontal shear strain is at a maximum at the base of the glacier (550 m) and reduces linearly to zero at the surface. Shear strain rate at any depth is proportional to depth to the power of n, where n is the stress exponent in the flow law. Using this relationship and the fact that the average of shear strain rates for all depths must correspond to the average shear strain rate averaged across the entire ice thickness (2.3 m yr −1 /550 m) allows the shear strain rate and shear strain at any given depth to be calculated. A similar approach can be used to account for the warmer temperatures at the base of the glacier using Q = 60 kJmol −1 . Using an n value of 3 (appropriate for low strains), lower bound horizontal shear strain rate increases from 10 −13 at 50 m to 7 × 10 −12 at 170 m depth (Table 1).      ) where Δ is the difference between the density of the surface and the density of the sample. e Upper bound horizontal shear strain rate was calculated by ice flow velocity (2.3 m/year) divided by thickness (550 m). f Lower bound horizontal shear strain rate was estimated with an assumption that the driving force for horizontal shear strain is at a maximum at the base of the glacier and reduces linearly to zero at the surface.

Analytical Methods
The microstructural features of rocks are usually observed in sections cut perpendicular to the foliation and parallel to the lineation (e.g., [47]). The foliation (bedding) of ice, in this setting, can be considered as the plane parallel to the surface, which forms due to the annual snow accumulation and densification during burial. Lineation cannot usually be observed owing to the transparency of ice. Azimuthal orientations are not normally preserved in rotary coring and for this study we cut an arbitrary vertical (normal to foliation) observation plane. The samples were cut by a ring blade (a circular saw) in a cold room at -15 to -20 °C, and mounted on a copper ingot at ~5 °C, thus, producing a bond between the sample and the ingot by melting and re-freezing a thin layer. As soon as a bond was formed, the ingot and the sample were cooled rapidly. To strengthen the bond, water droplets at 0 °C were added into the holes at the back of the ingot which froze onto the bottom of the sample. We used a microtome to cut the samples to thickness of ~1 cm in a freezer at -15 to -20 °C. The final polishing was achieved by using sandpaper (1500 mesh) attached to a polishing plate sit 1-2 cm above liquid nitrogen at the base of a polystyrene box, so that the plate was about -20 to -30 °C (for further details, see [48]).
To study the microstructures of ice, high-quality crystallographic orientation maps were acquired by indexing the Kikuchi bands (diffraction patterns) using an Oxford Instruments EBSD camera on a Zeiss Sigma variable-pressure, field-emission, scanning electron microscope (FE-SEM), housed at the University of Otago, Dunedin, New Zealand. Each polished sample (ca. 30 × 45 mm) was transferred through a nitrogen-filled glove box into the SEM to minimize condensation. After the sample was inserted to the cold stage in the SEM chamber at ≤-80 °C, the chamber was evacuated to a high vacuum and the sample was cooled to ≤-100 °C. A sublimation process was then carried out by bringing the chamber pressure to the atmospheric pressure and warming the sample to approximately -80 °C to remove frost and flatten the surface of the sample, and then the chamber was evacuated again to a high vacuum and the sample was cooled again to less than -100 °C. The orientation maps were generated with an accelerating voltage of 30 kV, a working distance of ~20 mm, a beam current of ~10 nA, and a step size of 30-50 µm with ~10 Pa nitrogen gas pressure to prevent charging (for further details of the process, see [48]). The raw data thus obtained, were processed using the software Channel 5 and the MATLAB-based toolbox MTEX [49]. Contoured pole figures of [0001], <11-20> and <10-10> were generated from all indexed pixels using MTEX. The Jand M-indices were determined [50][51][52] as measure of fabric strength using MTEX. Eigenvalues for calculating shape factor, K, were computed using the MTEX [53].

Results
All results from our microstructural analyses are summarized in Table 2. As illustrated in Figure  4, the age and density of the ice increases, and the longitudinal strain rate decreases with increasing depth. In this section, we describe the microstructures, CPOs, and related data.
The high-quality EBSD maps, figures for misorientation between each pixel and the mean orientation of the parent grain (mis2mean), and the grain-size distributions are displayed in Figure5a-c, respectively, for intracrystalline and intercrystalline deformation features. Analyzed sections have different arbitrary azimuthal orientations as azimuthal information was lost in coring. Sample SG_50 contains a large number of air bubbles. The ice crystals in this sample are characterized by relatively small and uniform grain size, equigranular shapes, and straight grain boundaries (Figure 5a,b). Samples from greater depths have; a lower fraction of bubbles, larger grain size, larger variations in grain size, and more angular and irregular grain boundaries. The figures for mis2mean indicate that a few areas with high misorientations are distributed along the grain boundaries, particularly near the triple junctions (Figure 5b). Mean grain size increases with depth and the peak distribution ranged between 1000-2000 µm (Figure 5c). Grain size at the peak frequency of the distribution is about 1.5 mm at 50 m and decreases to ~1 mm at 80 and 110 m. The crystallographic orientations are shown in a reference frame where the vertical is in the center of the stereonets and the data are rotated around the vertical axis so that the a-axis maxima lie in the east-west (E-W) direction. This approach facilitated a better comparison of data from different depths, but the reader is reminded that the azimuthal orientation is arbitrary (the direction of the aaxis maxima is unknown). The c-axes become more concentrated with increasing depths (Figure 6a). The CPOs are characterized by a vertical cluster of c-axes with weak secondary maxima at ≤110 m. The strong cluster of c-axes lies within a broad vertical girdle at depths ≥140 m (Figure 6b). The girdle is partial at 140 m and almost complete at 170 m. The vertical girdle of c-axes lies perpendicular to the a-axis maxima. The a-axes and poles to m-planes lie sub-horizontally. Weaker a-axis sub-maxima lie at ±60 degrees to the main maximum (normal to c-axis girdle). Maxima in the poles to m-planes lie at ±30 degrees and at 90 degrees to the a-axis maxima: those at ±30 degrees are strongest.  Table 2). Mean grain size and standard deviation increase slightly, corresponding to increasing density, with increasing depth (Figure 7b and Table 2). The shape factor, K, computed from eigenvalues of the c-axis distribution [53] provides quantitative identification of c-axes being more clustered at 50-110 m depth and more girdle-like at 140-170 m. Figure 7. Plots of (a) J-index, (b) mean grain size, and (c) shape factor, K, as a factor of increasing depth. The range presented as black arrows and a blue arrow and in (b) denote the standard deviation and peak grain size, respectively. The shape factor is calculated from eigenvalues [53].
Rotation axes (misorientation axes) of low-angle boundaries (2°-10°) in our samples were analyzed in inverse pole figures (Figure 8). The rotation axes of ice were found to be generally distributed along the basal plane, and/or parallel to c-axes. Strong intensity parallel to the c-axes is activated strongly at 140 m, and more weakly at 50 and 170 m, while intensity along the basal plane is clear in most samples, except the sample from 140 m.  a Success rate represents the initial number of indexed grains divided by all the analyzed samples before noise reduction, and zero solution includes the space represented by the bubbles and cracks; therefore, the number is underestimated.

Stress Regimes and CPO Development
Deformation features can be identified from intracrystalline microstructures, such as the CPOs (e.g., [54][55][56][57]). (Table 1). In terms of density, sample SG_50 is classified as firn-ice transition, and the other samples are defined as ice [58,59]. Deformation temperatures between -29.1 to -31.2 °C estimated from the measurements of the bore hole, are colder than the vast majority of published deformation experiments for which CPO and microstructural data have been presented, but are at a similar temperature to −30 °C experiments presented by [17,26,29].
The CPOs of natural ice could be attributed to both vertical shortening and horizontal shear. The upper bounds of the estimated rate of horizontal shear is equivalent or greater than the calculated rate of finite longitudinal shortening at 80-170 m ( Table 1). Lower bound estimations of horizontal shear strain at 140-170 m (>0.05) could induce a strong CPO [60]. The weak central cluster with or without multiple maxima of c-axes and weak fabric strengths (Figures 6 and 7a) at depths of 110 m, can be ascribed to low strains [20,[60][61][62][63]. Higher longitudinal strain rate than horizontal shear strain rate at 50 m suggests that the single cluster of c-axes is related with the vertical compression under firn-ice transition condition (Figure 4c) [18]. The single cluster with multiple maxima of ice c-axes from 80-110 m can be attributed to a dominant activation of the horizontal shear, as estimated upper bound horizontal strain rates are higher than estimated vertical strain rates in all samples below the 50 m sample (Figure 4c). All samples have estimated shear strains that are significant, given that strong CPOs are well developed by shear strains of ~0.2 [29]. Similar CPOs have been reported as products of deformation due to recrystallization under a simple shear regime by numerical models [29,64]. Therefore, we infer that the CPOs formed at 110 m are influenced by both the finite longitudinal and horizontal shear regimes.
In contrast, the c-axes of deeper samples (≥140 m) show a strong cluster with weak maxima aligned vertically (a girdle) with an abrupt increase in fabric indices, followed by more girdle-like (K < 1) shape factors, and longitudinal strain rates that are about one to two orders of magnitude lower than the horizontal strain rates (Figures 4c, 7a,c and Table 1). The girdle for c-axes is often ascribed to the development of a horizontal extension, based on the orientation of the vertical girdle that intersects the single cluster [19,22,23,[65][66][67]. We suggest that horizontal extension may become important, probably in addition to vertical compression and horizontal shear. The strong alignment a-axis at 140-170 m depths may lend support for the horizontal extension origin of the girdle CPOs. Slip on basal planes to facilitate extension would lead to the slip direction aligning with extension direction. The lower J-index at 170 m rather than at 140 m, may be attributed to the 'countervailing effect' of the stronger girdle against the presence of a strong cluster or a small number of analyzed crystals at 140 m. The J-and M-indices are scalars and they become less useful for comparing CPO strengths when CPO symmetry changes. The occurrence of tephra at depths of 97.01, 99.18, and 165.37 m indicates that the change in CPO may be irrelevant with impurities in the ice [40].

Recrystallization Processes and Slip Systems
The flow of ice can be influenced significantly by microstructural changes, particularly those attributed to recrystallization by grain boundary migration (GBM) and/or by lattice rotation including associated recovery, polygonization, and subgrain rotation (e.g., [28,29,68]). Recrystallization by GBM is driven by reduction in the mean intragranular distortion of the polycrystals. GBM alone results in an increase in grain size as no new orientation (nuclei) are created. Recrystallization by lattice rotations is accompanied by a substantial reduction as defects and associated distortion become localized in subgrain boundaries greating cells that may evolve to become new grains (nuclei) [22,69].
Sample SG_50, representing the firn-ice transition [58,59], shows straight grain boundaries, a relatively uniform grain size, triple junctions of ~120°, weak fabrics, and a very weak single cluster of c-axes (Table 1, Figures 5 and 6), which is typically regarded as a product of polygonization in shallower ice with compression [18]. Scarcity of small domains with a high misorientation suggests that the numerous air bubbles, in or near grain boundaries, may play an important role in generating fabrics at 50 m (Figure 5b). Between 80-110 m, the analyzed ice samples displayed a steady increase in grain size and fabric strength, and a minor change in grain shape (Figures 5a,b and 7b). At depths between 110 m and 140 m, the ice samples exhibited an abrupt increase in fabric strength. This change corresponds to a more diverse grain size distribution (high standard deviation), and the appearance of lobate grain boundaries ( Figure 5), thus indicating recrystallization involves GBM. Small domains with high misorientation, 120 triple junctions, and straight grain boundaries are distributed at all depths, suggesting the formation of new grains and the occurrence of recovery and subgrain rotation and recovery among the analyzed samples ( Figure 5). The evidence suggests, therefore, that recrystallization at depths ≥140 m is dominated by GBM accompanied by recovery and subgrain rotation.
Analysis of low angle boundaries is often useful for delineating the dominant slip systems in ice (e.g., [64,70]). The rotation axes of the analyzed samples, shown in Figure 8, can be classified into two types: one parallel to the c-axes that is activated at the depths of 50, 140, and 170 m, and another along the basal plane that has been highlighted at all depths except at 140 m. Assuming that only the basal slip system is active, rotation axes parallel to the c-axis can be related to twist walls (screw dislocations) of the basal slip system, while rotation axes along the basal plane can be related to tilt walls (edge dislocations) of the basal slip system. Statistically significant data sets from experiments identify both of these rotation axes [26,71] and subgrain boundary trace analyses, and weighted burgers vector analyses [72] of these same samples (Sheng Fan, personal communication) indicate that the rotation axes in basal plane are clearly related to basal plane dislocations. Other studies have constrained the operation of non-basal dislocations in both experimental [73,74] and natural [70,75] ice and further work is needed to determine the contribution of non-basal dislocations in these samples.
There is clearly a change in misorientation axes between 110 and 140 m, corresponding to the change in CPO. Thus we infer that the changes in deformation kinematics that cause the change in CPOs is probably also responsible for the change in misorientation axes. Whether this in turn relates to a redistribution of slip from tilt walls to twist walls or a changing involvement of non-basal dislocations needs further analysis.
In the finite longitudinal/horizontal shear stress regimes, the operation of the basal slip could explain the formation of a single cluster with weak maxima of c-axes at depths 110 m and with a girdle at 140-170 m. Under a compressional (and/or a simple shear) stress regime at 50-110 m, slip occurs along the basal plane and the c-axes align around grains in high Schmid factor orientations at the expense of grains in the low Schmid factor orientations ( Figure 9a). Under an extension regime, at 140-170 m, slip occurs along the basal plane and the c-axes disperse to rotate to form a girdle, with fixed a-axes along the tensile direction (Figure 9b). These processes can explain the formation of distinctive CPOs under different stress regimes.

Origin of Extensional Stress Regime and Implications for Ice Sheet Dynamics
Antarctic ice forms by recurrent patterns of snow accumulation, densification, and transformation through firn to ice [38,58,[76][77][78]. The Styx Glacier is a renowned site of snow accumulation; hence, the c-axes of ice crystals can be expected to be aligned perpendicular to the surface. The results from the microstructural investigations in this study, however, demonstrate the CPO evolution with depth from a single cluster of c-axes to a single cluster with a girdle. The ice in the study area flows and meets a nunatak within ~4 km to the SE. The topography may generate an extensional stress regime that may be responsible for the weak girdle of c-axes at depths ≥140 m (Figure 10a). Our GPR data suggest that the bedrock becomes shallower from the sample location towards the nunatak (Figure 3). Additionally, a small bedrock peak corresponds to the core site and may contribute to an extensional stress regime at the sample location (Figure 10b,c).
We have compared our results to conventional ice-flow models that regard ice as a viscoelastic material. Conventional numerical models of ice flow [76,79] are characterized by a low velocity at the base due to a high friction with the bedrock, and a logarithmic increase in speed toward the surface. In the study area, the evolution of the CPO from a single cluster of c-axes to a single cluster with a girdle, can be attributed to the kinematic control of the underlying bedrock. Although we have concluded that the structure of the bedrock is the most likely cause for the generation of the CPO evolution, other environmental factors, such as the concentration of air bubbles, gas or volcanic tephra due to variations in annual snow accumulation, need to be tested [32][33][34]. Therefore, we suggest that the microstructural analyses of various ice cores should be explored extensively by sampling from kinematically distinct locations within ice sheets. This process is crucial for a thorough review of the conventional ice flow models. Moreover, new techniques must be devised to facilitate preservation of the azimuthal orientations of cored ice samples.

Conclusions
We conducted a microstructural study of five natural ice samples taken at intervals of ~30 m from depths of 50-170 m by a core drilled in the Styx Glacier, located in northern Victoria Land, Antarctica. Our interpretations of the microstructures, based on high-quality EBSD maps combined with GPR, can be summarized as follows: 1. With increasing depth, the analyzed samples show a systematic change from a weak single cluster of c-axes at 50 m, to a single cluster with weak multi-maxima between 80-110 m, and to a single cluster with a girdle between 140-170 m. An a-axis cluster occurs normal to the caxis girdle. The formation of the girdle that accompanied the strong single cluster is characterized by an abrupt increase in fabric indices, and a decrease in the shape factor; thus, suggesting the extension began to occur at the depth of <140 m. 2. Straight grain boundaries, relatively uniform grain size, triple junction of ~120°, and a weak single cluster of c-axes at 50 m (Figures 4 and 6) are typically regarded as products of polygonization under compression, mainly due to bubble extraction. With an increased depth of burial, grain size increased, fabric strengthened, and the irregularity increased in grain shape ( Figures 5 and 6). At the depth of 140 m, highly curved grain shapes, abrupt increase of fabric indices, and the occurrence of a girdled c-axes (K <1) suggest dynamic recrystallization dominated by GBM. 3. The ice generally flows at ca. 2.3 m yr −1 from the northwest (NW) to the SE, dominated under the finite longitudinal and horizontal shear regimes. The presence of a nunatak located ~4 km far from the sample location in the SE direction, and a low surface velocity of ice, allowed the appearance of a girdle of c-axes between 140-170 m. The girdle was the result of an extensional stress regime generated at that depths. These results suggest the occurrence of changeable stress regimes as a response to diverse environments, such as an abrupt appearance of a small peak in the bedrock.
Author Contributions: D.K., D.J.P., and Y.H. collaboratively designed the project and wrote the manuscript. D.K., D.J.P., and C.Q. conducted EBSD analysis, and H.H. and H.T.J provided satellite and GPR data, respectively. All authors contributed to correct the initial version of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by a KOPRI project PE20030.