1. Introduction
The global mean sea level has risen faster with acceleration in recent decades. It has been reported that the sea-level changes occur due to ice loss particularly in the Greenland ice sheet. However, it has also been reported that the acceleration in mass loss between 2006 and 2015 was at an average rate of 278 ± 11 Gt yr−1 from Greenland and 155 ± 19 Gt yr−1 from Antarctica (IPCC, Reference Pörtner, Roberts, Masson-Delmotte, Zhai, Tignor, Poloczanska, Mintenbeck, Alegríam, Nicolai, Okem, Petzold, Rama and Weyer2019: Summary for Policymakers). One of the reasons for mass loss is believed to be related to the ice sheet flows from Greenland and Antarctica. As mentioned above the ice sheet mass loss is particularly significant in Greenland. Such flows can occur due to basal slip on the soil surface and deformation of ice (Cuffey and Paterson, Reference Cuffey and Paterson2010). The latter is related to the creep phenomenon with the assistance of dislocation motions under applied stress and temperature. The concept of creep of ice is well established for natural ice sheets and experimental artificial ice, although it is known that a universal constitutive law for ice does not exist. For example well-known Glen’s flow law explains how ice plastically deforms (Glen, Reference Glen1955). During creep deformation particularly in a tertiary creep regime dynamic recrystallisation usually takes place (Duval, Reference Duval1981; Duval and others, Reference Duval, Ashby and Anderman1983) with fabric development. Budd and Jacka (Budd and Jacka, Reference Budd and Jacka1989) summarised different creep behaviours between ice sheets and artificial ice due to different stress fields and crystal orientations. Although anisotropic ice deformation was common in the ice sheet, it was ill-considered in the flow-law model. Azuma investigated the deformation behaviours of anisotropic polycrystalline ice and then constructed a power-law creep equation with a steady-state creep rate, including the mean Schmid factor in the pre-exponential factor (Azuma, Reference Azuma1995), and then extended it to the Cartesian coordinate system (Azuma, Reference Azuma1994; Azuma and Goto-Azuma, Reference Azuma and Goto-Azuma1996).
One of the important factors affecting deformation behaviour in ice is the existence of microparticles. The reasons for the presence of such impurities in ice could be related to minerals of continental origins, sea salts, biological origins and so on. In an earlier study ultrafine amorphous silica was artificially dispersed in ice (Nayar and others, Reference Nayar, Lenel and Ansell1971). To model the base of an ice sheet which usually includes a layer of dirty ice several meters in thickness, sand was included in artificial ice (Hooke and others, Reference Hooke, Dahlin and Kauper1972). These creep experiments show strengthening effects due to the existence of sands in ice on creep deformation. On the other hand the presence of particles in ice also affects dynamic recrystallisation or grain boundary migration during creep (Humphreys and Hathery, Reference Hori, Hondoh, Oguro and Lipenkov2004), depending on the stress level (Song and others, Reference Song, Cole and Baker2004, Reference Song, Baker and Cole2008), resulting in softening of the ice. This discrepancy between the effects of particles on hardening or softening of creep deformation of ice could be related to the average grain size of the ice and particles. We then tested making a fine-grained ice with a few tens of micrometres in diameter, containing micro- particles with 300 nm in diameter and found that the addition of microparticle softens the ice during creep, and the corresponding strain rate-strain curves do not show clear appearance of a steady state, irrespective of the microparticle additions. Thus, it is concluded that the recovery of dislocation pile-ups due to lattice diffusion based on a dislocation creep mechanism is the rate-controlling process, and the microparticles insignificantly influence the deformation of the fine-grained ice, except for their effect on grain boundary pinning (Saruya and others, Reference Saruya, Nakajima, Takata, Homma, Azuma and Goto-Azuma2019). The details of the creep conditions mentioned above are given in Table 1.
Experimental conditions of creep deformation using artificial ice
When considering a dislocation glide or climb during creep deformation, the creep rate 
$\dot \varepsilon $ in both cases can be described by the well-known Orowan equation,
\begin{equation}\dot \varepsilon = \rho bv\end{equation} where 
$\rho $ is dislocation density, 
$b$ the Burgers vector and 
$v$ the average velocity of dislocations (Poirier, Reference Poirier1976; Cuffey and Paterson, Reference Cuffey and Paterson2010). In the case of steady-state creep under a constant-stress condition the steady-state creep rate is expressed as 
$\dot \varepsilon \propto {\sigma ^m}$, according to Weertman (Weertman, Reference Weertman, Whalley, Jones and Gold1973). Here 
$m$ is a stress exponent. More details will be discussed later. Based on these equations we understand that the dislocation density is one of the important factors in determining the steady-state creep rate. The parameters used in equations are shown in Table 2.
Symbols and definitions in equations
First systematic deformation experiments of water-frozen ‘single-crystal’ ice were conducted by Glen (Glen, Reference Glen1952). He observed that slips appear after a tensile test. Nakaya was the first to report plastic deformation behaviours in single crystalline ice from the Mendenhall Glacier in Alaska with the aid of dislocations, showing slip lines during a bending test, which were observed by Foucault’s method of shadow photography (Nakaya, Reference Nakaya1956). He then found that gliding planes appeared after bending the single crystalline ice. He attributed this to the occurrence of slips related to the dislocations. In the early stage of the research the existence of dislocations was indirectly confirmed using an etch-pit method and electron microscopy for lake ice (Muguruma, Reference Muguruma1961). However the observations of the dislocations in ice were successfully made using X-ray photographs (X-ray diffraction topography) as reported by Hayer and Webb (Hayes and Webb, Reference Hayes and Webb1965) and were then extensively analysed by Higashi et al. using the X-ray diffraction topography (Fukuda and Higashi, Reference Fukuda and Higashi1969; Higashi, Reference Higashi1988) based on the Lang method (Lang, Reference Lang1959). Nevertheless, while the X-ray diffraction topography was a powerful tool to directly reveal dislocations in the single crystals, the measurable dislocation densities were limited to around 102 cm−2 (106 m−2). Instead of measuring dislocation densities by X-ray diffraction topography Hori et al. applied an X-ray diffraction method to obtain dislocation densities in the Vostok ice core from the rocking-curve measurements with Cu-K
$\alpha $ radiations (Hori and others, Reference Hori, Hondoh, Oguro and Lipenkov2004). The instrument was installed in a cold room kept at −20ºC, and a 
$\left( {10\bar 12} \right)$ reflection in each crystal was analysed. It was found that the measured dislocation densities ranged from 1.9 × 108 to 2.4 × 1010 m−2, and large numbers of dislocations were detected in the glacier ice extracted from a depth of around 3286 m.
The measurements of dislocation densities were significantly developed by Krivoglaz and Wilkens and Ungár (Krivoglaz, Reference Krivoglaz1969; Wilkens, Reference Wilkens, Simmons, de Wit and Bullough1969; Ungár and others, Reference Ungár, Ott, Sanders, Borbély and Weertman1998). It is well known that both contributions from grain sizes of the matrix of polycrystalline materials and strains are involved in the diffracted-peak intensities. As described by Hori and others (Hori and others, Reference Hori, Hondoh, Oguro and Lipenkov2004), they only chose the 
$\left( {10\bar 12} \right)$ reflections for the diffracted intensities of the ice cores. In their appendix, they considered the effect of the dot product of 
${\boldsymbol{g}}$ and 
${\boldsymbol{b}}$, where 
${\boldsymbol{g}}$ and 
${\boldsymbol{b}}$ are reciprocal lattice and Burgers vectors, respectively, in order to have contrast. Krivoglaz and Wilkens included the dislocation contrast factor as an effect of the strains, which will be discussed in more detail in the experimental methods section of this paper. This method, known as the modified Warren-Averbach and modified Williamson-Hall plots, can obtain more precise dislocation densities in strained materials, even in deformed ice. Nevertheless, such an approach has never been applied to ice.
In order to evaluate the dislocation densities in deformed ice samples using the modified Warren-Averbach and modified Williamson-Hall plots of X-ray line broadening we developed an apparatus that to be located inside a commercial XRD instrument with Cu-K
$\alpha $ radiations to maintain ice at low temperatures. By using this method we succeeded in accurately measuring the dislocation densities in ice as a function of the creep strain at −20ºC under an applied stress of 2 MPa. The effect of microparticles on the dislocation densities is also discussed.
2. Experimental methods
2.1. Sample preparations
Ultra-pure water (18.2 M
$\Omega $-cm) was sprayed into liquid nitrogen in order to produce ice powders. Additionally SiO2 (silica) powders, simulating impurities in ice cores, with a diameter of 300 nm and an amorphous structure were mixed with the ice powders. Two different artificial ice samples were prepared; (1) pure ice powders and (2) pure ice powders mixed with microparticles containing ice. The latter contained 0.1 mass% silica in the ice. These powders were sieved with a mesh size of 710 
$\mu $m so as to obtain a uniform powder size. The powders were loaded into a metallic die and then sintered at a cold room maintained at −10ºC for 1 h and pressed at 70 MPa using a hydraulic jack. During the sintering the inside of the die was evacuated by a rotary pump to remove air bubbles as much as possible (Azuma and others, Reference Azuma, Miyakoshi, Yokoyama and Takata2012). The average dimensions of the sintered ice with a diameter of 33 mm and a height of 65 mm were cut. The obtained artificial ice had polycrystalline microstructures.
Creep specimens were cut from a sintered sample into four pieces with equivalent volumes, and then each piece was turned on a lathe to shape the specimen into a cylinder with a diameter of 15 mm and a height of 30 mm.
2.2. Uniaxial compressive tests of the bulk specimens
Creep specimens were immersed in a container filled with silicone oil to prevent sublimation during the test and to suppress friction between the specimen and jig. A load was applied to the cylindrical creep specimen parallel to its long axis and controlled to maintain the applied stress at 2 MPa. The creep test machine was placed in a constant temperature bath inside a freezer maintained at −20ºC. The strain was measured using the digital displacement meter (Mitutoyo ID-F125) attached to the creep machine.
2.3. Microstructural observations
The microstructures before and after the creep tests were observed using an optical microscope (Olympus BX51). The samples were cut from the cylindrical specimens used for the creep tests, and the observation directions were fixed perpendicular to the long axis of the cylindrical specimens. The upper and bottom parts of the cylinder were cut about 10 mm from the edges, and a remaining 10 mm part located near the centre of the cylinder was used. The average grain diameter was calculated based on the equivalent circle diameter after measuring the area of each grain using Image J. Additionally, subgrain boundary densities were also calculated by measuring the length of the subgrain boundaries using Image J within the area of a field of view.
2.4. Sample holder for X-ray diffraction analyses
The X-ray diffraction analyses were performed using an commercial XRD machine (Rigaku Ultima IV) installed in a room maintained at ambient temperature. A sample holder shown in Figure 1a was made of a Cu plate (thermal conductivity: 400 Wm−1K−1), with blocks of thermal insulators made of Teflon covering the top surface, except for the space parallel to the beam line. The ice was maintained between the Teflon blocks, which were 8 mm thick, on the Cu plate and cooled at −20ºC. After fixing the ice sample to the Cu plate in a cold room kept at −10ºC, the space was sealed using a Kapton tape with a width of 20 mm so that it did not disturb the X-ray emissions to the surface of the ice at large Bragg angles. The Cu plate was cooled using a container which was made of an A5052P aluminium alloy (thermal conductivity: 200 Wm−1K−1) and filled with cooling nitrogen gas (purity of 99.9%) flowing into a helical and 750-mm-long Cu tube kept in liquid nitrogen at −196ºC as shown in Figures 1c and d. The sample holder was easily removable and placed on the cooling container just before starting the X-ray analysis. The container could be adjusted to the central position by an XYZ linier translation stage, and the optimum height was adjusted using a stereomicroscope (Figure 1b). The surface of the Kapton tape was gently covered with dry air to prevent frosting due to moisture absorption from the air in the experimental room maintained at ambient temperature. The X-ray radiation used was Cu-k
$\alpha $ with a wavelength of 0.15406 nm. The tube voltage and current used were 40 kV and 20 mA, respectively. The step size and scan speed used were 0.01º and 0.2º min−1, respectively. The schematics of the developed cooling system for the X-ray measurements of ice samples are summarised in Figure 2.
(a) design of removable sample holder, (b) arrangement of x-ray generator, sample holder and x-ray detector in ultima IV, (c) at the top of the dewar vessel poured liquid nitrogen and (d) a helical and 750-mm-long cu tube which was immersed in liquid nitrogen.
Schematic of the developed cooling system for x-ray measurements of ice samples using N2 gas and liquid N2. Symbols (a) to (d) correspond to Figures 1a to 1d, respectively.
2.5. Dislocation density measurements
At least two independent parameters of grain size and strain in polycrystalline materials are included in diffracted intensities of the X-ray (Warren and Averbach, Reference Warren and Averbach1950). The shape of a diffracted intensity can be represented by a Fourier series, and it is expressed as
\begin{equation}{A_n}\!\left( l \right) = A_n^PA_n^D\left( l \right)\end{equation} where 
$n$ is an index ranging from 
$ - \infty $ to 
$\infty $ including zero, 
$l$ is the Miller index given e.g. by 
$\left( {00l} \right)$ and 
$A_n^P$ and 
$A_n^D\left( l \right)$ are particle size and distortion coefficients, respectively (Warren and Averbach, Reference Warren and Averbach1952). Williamson and Hall originally reported the composite broadening of diffracted intensities in X-ray diffraction (XRD) produced by particle size and strain effects (Williamson and Hall, Reference Williamson and Hall1953), expressed by
\begin{equation}\delta = \frac{\lambda }{{t{\text{cos}}\theta }} + 2\xi {\text{tan}}\theta \end{equation} where 
$\delta $ is a spectral integral line breadth, 
$\lambda $ the X-ray wavelength, 
$t$ a mean linear dimension of the grain, 
$\theta $ the Bragg angle and 
$\xi $ the strain distribution. If the well-known Eq. (3) is rearranged, e.g., to
\begin{equation}\frac{{\delta {\text{cos}}\theta }}{\lambda } = \frac{1}{t} + \frac{{2\xi {\text{sin}}\theta }}{\lambda }\end{equation} then from the intercept and slope of the 
$\delta {\text{cos}}\theta /\lambda - 2{\text{sin}}\theta /\lambda $ diagram, the particle size 
$t$ and strain 
$\xi $ can readily be extracted. However, the scattering of the X-rays is significantly affected by the incident beam direction, dislocation direction, the Burgers vector and elastic anisotropy of the crystals. To improve the accuracy of the measurement the dislocation contrast factor was included in the theory of the broadening of the X-ray intensities by Krivoglaz and Wilkens (Krivoglaz, Reference Krivoglaz1969; Wilkens, Reference Wilkens, Simmons, de Wit and Bullough1969). This was then extended to be analysed from the actual XRD patterns of cubic crystals (Groma and others, Reference Groma, Ungár and Wilkens1988; Ungár and others, Reference Ungár, Ott, Sanders, Borbély and Weertman1998), and was applied for hexagonal crystals (Klimanek and Kužel, Reference Klimanek and Kužel1988; Ungár and others, Reference Ungár, Gubicza, Ribárik and Borbély2001; Dragomir and Ungár, Reference Dragomir and Ungár2002).
The intuitive expression of the diffracted intensity 
${I_g}\!\left( S \right)$ including the effect of elastic strains can be explained by (Wilkens, Reference Wilkens, Simmons, de Wit and Bullough1969)
\begin{equation}{I_g} = \mathop {\smallint} {A_g}\left( n \right){\text{exp}}\left( {2\pi inS} \right)dn\end{equation} where 
${A_g}\!\left( n \right)$ is the Fourier transform of 
${I_g}\!\left( S \right)$ which is normalised to unity, 
$S = 2\left( {{\text{sin}}\theta - {\text{sin}}{\theta _0}} \right)/\lambda $, 
$\theta $ diffraction angle, 
${\theta _0}$ the Bragg angle of the undistorted lattice and 
$\lambda $ the X-ray wavelength. The effect of strains 
${\varepsilon _n}$ is inserted in 
${A_g}\!\left( n \right)$ as
\begin{equation}{A_g}\!\left( n \right) = \langle {\text{exp}}\!\left( {2\pi ign \cdot {\varepsilon _n}} \right)\rangle\end{equation} The brackets 
$\langle \,\rangle $ indicate a mean value and 
$g$ the reciprocal lattice vector. According to the theory of dislocations the elastic strain energy of dislocations per unit length can be expressed as (Wilkens, Reference Wilkens1967; Hirth and Lothe, Reference Hirth and Lothe1968)
\begin{equation}W = C\frac{{\mu {b^2}}}{{4\pi }}\rho {\text{log}}\left( {\frac{{{R_e}}}{{{r_0}}}} \right)\end{equation} where 
$\mu $ is the shear modulus. Additionally, Wilkens derived the following equation for the mean square of the differential strain using Eq. (7) and a rough approximation of Young’s modulus 
$E = 2\mu $ by considering Poisson’s ratio as 0.325 (Petrenko and Whitworth, Reference Petrenko and Whitworth1999),
\begin{equation}\langle \varepsilon _0^2\rangle = \frac{{{b^2}}}{{4\pi }}\rho C{\text{log}}\left( {\frac{{{R_e}}}{{{r_0}}}} \right)\end{equation} here 
$C$ is the dislocation contrast factor, 
${R_e}$ and 
${r_0}$ are the outer and inner cutoff radii of the dislocation, respectively. Thus we realise from Eqs. (5), (6) and (8) that the intensity profiles of XRD patterns can be broadened due to the presence of strains related to the dislocations whose effect is included in the parameter of 
$C$.
The obtained X-ray profiles were fitted using the Lorentz function. The modified Warren-Averbach (WA) and modified Williamson-Hall (WH) plots were then used to evaluate the dislocation densities as a function of the creep strain. Here we assign 
$K = 2{\text{sin}}\theta /\lambda $ and 
$\Delta K = 2{\text{cos}}\theta \left( {\Delta \theta } \right)/\lambda $. The full width at half-maximum (FWHM) of line profiles can be given as
\begin{equation}\Delta K = \frac{{0.9}}{D} + {\left( {\frac{{\pi {H^2}{b^2}}}{2}} \right)^{1/2}}{\rho ^{1/2}}K{\bar C^{1/2}} + \left( {\frac{{\pi H'{b^2}}}{2}} \right){U^{1/2}}{K^2}\bar C\end{equation} where 
$D$ is the volume averaged grain size, 
$H$ a constant, 
$\bar C$ the average dislocation contrast factor, 
$H{^{\prime}}$ a constant and 
$U$ a correlation factor (Ungár and others, Reference Ungár, Ott, Sanders, Borbély and Weertman1998). The last term in Eq. (9) is usually small and hence it can be neglected. Eq. (9) is known as the modified WH plot. In the case of hexagonal materials the average dislocation contrast factors are expressed as (Dragomir and Ungár, Reference Dragomir and Ungár2002).
\begin{equation}{\bar C_{hk.l}} = {\bar C_{hk.0}}\!\left( {1 + {q_1}x + {q_2}{x^2}} \right)\end{equation} where 
${q_1}$ and 
${q_2}$ are fitting parameters and 
$x = \left( {2/3} \right){\left( {c/a} \right)^2}$: 
$a$ and 
$c$ are lattice parameters of the hexagonal materials, and 
$c$ is measured parallel to the 
$c$ axis. From Eqs. (9) and (10)
\begin{equation}\frac{{{{\left( {\Delta K} \right)}^2} - \alpha }}{{{K^2}}} = \kappa {\bar C_{hk.0}}\left( {1 + {q_1}x + {q_2}{x^2}} \right)\end{equation} is derived, where 
$\alpha $ is 
$0.9/D$ and 
$\kappa = {\left( {{H^2}{b^2}/2} \right)^{1/2}}{\rho ^{1/2}}$.
The modified WA plot is given by
\begin{align}lnA\left( L \right) &\approx ln{A^P}\left( L \right) - \left( {\frac{{\pi \rho {b^2}{L^2}}}{2}} \right){\text{ln}}\left( {\frac{{{R_e}}}{L}} \right){K^2}\bar C + \left( {\frac{{U{\pi ^2}{b^4}{L^4}}}{4}} \right)\nonumber\\ &\quad{}\times{\text{ln}}\left( {\frac{{{R_1}}}{L}} \right){\text{ln}}\left( {\frac{{{R_2}}}{L}} \right){K^4}{\bar C^2}\end{align} where 
$L = n{a_3}$ is the Fourier length, 
${a_3} = \lambda /2\left( {{\text{sin}}{\theta _2} - {\text{sin}}{\theta _1}} \right)$ is in the direction of the diffraction vector 
${\boldsymbol{g}}$ and the line profile is measured from 
${\theta _1}$ to 
${\theta _2}$ of an angular range, 
${A^P}\!\left( L \right)$ is the size contribution to the Fourier coefficient for the intensity distribution and 
${R_1}$ and 
${R_2}$ are auxiliary constants which are not interpreted physically (Ungár and others, Reference Ungár, Ott, Sanders, Borbély and Weertman1998). Eq. (12) can be simplified such as
\begin{equation}{\text{ln}}A\left( L \right) \approx \tau - \varphi {K^2}{b^2}\bar C + o{K^4}{b^4}{\bar C^2}\end{equation} and the last term is omitted and then the parameter 
$\varphi $ becomes a fitting parameter when a 
$lnA\left( L \right) - {K^2}{b^2}\bar C$ diagram is drawn by selecting several Bragg peaks in the line profile
\begin{equation}\varphi = \rho \frac{{\pi {L^2}}}{2}{\text{ln}}\left( {\frac{{{R_e}}}{L}} \right)\end{equation} From Eq. (14) we can create a diagram of 
$\varphi /{L^2} - {\text{ln}}L$, and it reveals the dislocation density from the slope of 
$ - \pi \rho /2$. The dislocation density was indeed measured as the sum of mobile and immobile dislocations using this technique.
3. Results and discussion
3.1. Creep curves
The measured strain rate and strain diagrams associated with the creep experiments are given in Figure 3, which show samples crept to 20% strain. Both samples show concave upward behaviours as a function of the strain: as the strain increases, the strain rate decreases. The initial trends of the strain rates are very similar to those reported earlier in our paper (Saruya and others, Reference Saruya, Nakajima, Takata, Homma, Azuma and Goto-Azuma2019); that is, the pure ice shows slightly concave downward behaviour up to ≈ 0.36% strain, while the silica-containing ice shows a rapid decrease in the strain rate until around 0.5% stain. This means that silica may promote creep deformation particularly in the early stages of the creep deformation in the primary creep region. However pure ice shows continuous decrease in the strain rate within our experimental conditions. As a result the silica-containing ice is much softer than the pure ice above 0.56% strain. This result is consistent with that reported by Saruya and others (Reference Saruya, Nakajima, Takata, Homma, Azuma and Goto-Azuma2019).
Strain rate and strain diagrams of pure and silica-containing ice measured at-20ºC and 2 MPa after the creep tests.
3.2. Microstructural observations and quantifications of grain size and subgrain boundary density
The creep tests were conducted up to 1, 10 and 20% strains in both the pure and silica-containing ice based on Figure 3. The observed microstructures are shown in Figure 4. The quantitative results associated with the average gran size and subgrain boundary density as a function of the strain are given in Figure 5. In the initial state before the creep deformations grain boundaries (GBs) indicated by black arrows are well distinguished. The shape of GBs is equiaxed and appears straight. Additionally, some spherical and weak contrasts are seen in the silica-containing ice, as indicated by white arrowheads. They are the silica particles, but they appear aggregated when observed with the optical microscope (OM).
Optical microscope images before and after creep deformations were obtained from pure and silica-containing ice. ‘Initial’ indicates the initial state before the creep deformations, and 1%, 10% and 20% indicate interrupted creep strains as shown in Figure 3.
Average grain size or subgrain boundary density as a function of strain measured at − 20ºC and 2 MPa.
Once stress is applied to the ice (1% strain), subgrain boundaries (sGBs) shown by white arrows appear within the grains particularly in the pure ice. They divide the crystal into several subgrains. Since there are no silica particles, dislocations can move easily. Large amounts of dislocations are generated as the 1% strain is within the primary creep regime. Thus both grain refinements and increases in sGB density occur in the pure ice as indicated in Figure 5. However the sGBs are rarely observed in the silica-containing ice at the 1% strain. The grain size also hardly changes, due to the pinning effect of the silica particles on the GBs (Figure 5). The drastic changes in the microstructures occur in the 10% strain. In both the samples the appearances of dynamic recrystallisation are obvious, and average grain sizes decrease. SGB densities also increase; it is readily estimated that large amounts of dislocations are introduced in the samples. Both microstructures are very similar, as shown in Figure 4, and therefore the quantitative results in Figure 5 also reveal similar trends. At 20% strain apparent grain growth can be confirmed in the pure ice, as shown in Figure 4. This is evident even from the grain size measurement shown in Figure 5. Nevertheless, while the grain size continuously decreases in the silica-containing ice, the size increases from 10 to 20% strains in the pure ice. Surprisingly, the sGB densities in both types of ice are almost saturated.
3.3. Dislocation densities in the crept ice
An example of the XRD profile obtained from the as-sintered ice is shown in Figure 6. We assume that the crystal is hexagonal ice (Ih) with the space group symbol of 
$P{6_3}/mmc$ (No. 194), with lattice parameters of a = 0.4501 nm and c = 0.7348 nm: only oxygen sites of 
$1/3,\,2/3,\,{z_O}$; 
$1/3,\,2/3,1/2 - \,{z_O}$; 
$2/3,\,1/3,\,1/2 + \,{z_O}$ and 
$2/3,\,1/3,\,1 - \,{z_O}$ are considered (Hobbs, Reference Hobbs1974; Petrenko and Whitworth, Reference Petrenko and Whitworth1999). The parameter 
${z_O}$ represents a small shift of the hexagonal rings in the plane normal to the c-axis. 
${z_O}$ is reported to be 0.0622. Based on this crystal structure the Miller indices are assigned to all the diffracted peaks shown in Figure 6, although they do not usually represent the four indices for hexagonal materials. Some small deviations from the ideal Bragg angles for ice Ih indeed appear. We neglect them as they are trivial and measure the modified WH and WA plots from some independent Bragg peaks appearing at more than 20 º.
XRD profile obtained from as-sintered pure ice.
During the calculations we used the hypothesis that the dislocation slips almost occur exclusively based on the basal slip system 
$ \lt 11\overline20 \gt \left(0001\right)$ (Hondoh, Reference Hondoh and Hondoh2000), even during the creep deformations under our experimental conditions. Then 
${\bar C_{hk.0}}$ parameter in Eq. (10) was calculated using a computer program called ANIZC (Borbély and others, Reference Borbély, Dragomir-Cernatescu, Ribárik and Ungár2003) which provided us with the average dislocation contrast factor by inputting the second order elastic constants 
${C_{ij}}$ or elastic stiffness constants 
${S_{ij}}$. The values used for 
${C_{11}}$, 
${C_{12}}$, 
${C_{13}}$, 
${C_{33}}$ and 
${C_{44}}$ were 13.20, 6.69, 5.84, 14.42 and 2.89 GPa, respectively, and c/a was 1.63 (Hobbs, Reference Hobbs1974). Here, c/a indicates the ratio of the c lattice parameter to the a lattice parameter.
The measured dislocation densities as a function of the stain are summarised in Figure 7. The minimum and maximum dislocation densities obtained in this study range from 1011 to 1016 m−2. The absolute values are quite high, and some reported dislocation densities for metallic alloys are on the order of 1014 m−2 using the modified WH and WA methods (Yin and others, Reference Yin, Hanamura, Umezawa and Nagai2003; Shintani and Murata, Reference Shintani and Murata2011). We must note that the dislocation density reaches 1.4 × 1014 m−2 even in the as-sintered pure ice in our case. This may indicate that dislocation densities significantly increase during the pressure-sintering process. Indeed metallic alloys subjected to severe plastic deformation, such as cold-working, show dislocation density greater than 1015 m−2 (Saymour and others, Reference Saymour2017). Though the same process was used for the fabrication of the silica-containing ice, the density is 3.8 × 1011 m−2. The decrease in the dislocation densities in the silica-containing ice just after sintering could be explained by the behaviour of creep deformations of the silica-containing ice at 20% strain, as mentioned below. At least we can discuss the relative dependence of the silica addition on the changes of dislocation densities by comparing between pure and silica-containing ice.
Dislocation densities measured using modified WH and WA plots for pure and silica-containing ice as a function of strain.
At 1% strain dislocation densities in both types of ice rapidly increase. In the silica-containing ice it shows a two orders of magnitude increase compared to the initial state. Nevertheless, although the dislocation densities increase, the grain sizes do not alter. The sGBs rapidly increase in the pure ice, but they do not increase in the silica-containing ice (Figure 5). The strain rate-strain diagrams shown in Figure 3 reveal a raid decrease in the strain rate in both samples, and the pure ice is harder. Hence, it can be concluded that the rapid decrease in the strain rate is related to work hardening, particularly in the pure ice in the primary creep regime. Similar conclusions are drawn in our previous paper (Saruya and others, Reference Saruya, Nakajima, Takata, Homma, Azuma and Goto-Azuma2019); the hardening behaviours of the strain rate-strain diagrams are related to dislocation pile-ups. However the silica-containing ice shows a rapid decrease in the strain rate until the 0.5% strain. Hence dispersion hardening (Ashby and Jones, Reference Ashby and Jones1980) may occur only in the initial stage of the deformation, but it does not continue in the later stage of the deformation. The dislocation density at 10% strain slightly decreases or remains almost consistent with that in the 1% strained sample in the pure ice and work hardening continues, while it exhibits the maximum density in the silica-containing ice. Simultaneously the grain sizes decrease and the sGB densities increase. This means recrystallisation occurs in addition to a sufficient increase of dislocations (Humphreys and Hathery, Reference Hori, Hondoh, Oguro and Lipenkov2004): it is confirmed in Figure 4. Some dislocations are used for the formation of sGBs (Figure 5) as well as strain hardening; therefore the strain rates decrease compared with those in the initial stage of the creep.
In both types of ice the stain rates continuously decrease as shown in Figure 3. This is also consistent with what we reported previously (Saruya and others, Reference Saruya, Nakajima, Takata, Homma, Azuma and Goto-Azuma2019). However, the dislocation densities show different behaviour. Although the measured dislocation densities are almost constant in the pure ice after the 1% strain, the density rapidly decreases in the silica-containing ice at 20% strain. The decrease is more than two orders of magnitude. Figure 5 suggests that the grain size of the silica-containing ice continuously decreases from 10% strain, while in pure ice it remains almost constant within the error bars. The sGB densities do not significantly change with or without the silica addition. The grain refinement must be related to the occurrence of dynamic recrystallisation, which consumes dislocations for the nucleation of new crystalline grains (Humphreys and Hathery, Reference Hori, Hondoh, Oguro and Lipenkov2004), but at the same time relatively large grains also remain (Figure 4); we can see zig-zag shapes of the GBs, i.e., occurrence of strain-induced boundary migrations (SIBM, Faria and others, Reference Faria, Weikusat and Azuma2014). The occurrence of SIBM is also confirmed in the pure ice. Again the hardening effect in the pure ice at 20% strain compared with that of the silica-containing ice (Figure 3) can be explained by work hardening, the occurrences of which is directly confirmed in Figure 7. Thus the silica addition promotes dislocation annihilations near the grain boundaries of the fine grains (Gifkens, Reference Gifkens1976) due to the pinning effect of the silica additions, leading to a decrease in the dislocation density at 20% strain. Hence the silica-containing ice is softer than the pure ice (Figure 3).
3.4. Rough calculations of dislocation velocities in the crept ice
Lastly the dislocation velocities in pure and silica-containing ice were roughly calculated based on Eq. (1) by considering only the basal glide during the creep deformation. The velocities were evaluated at the 10% strains. The dislocation velocities in the pure and silica-containing ice are 3.10 × 10−12 and 4.17 × 10−14 m s−1, respectively, assuming that the measured dislocations include both geometrically necessary dislocations and statistically stored dislocations (Ashby, Reference Ashby1970). It should be emphasised that the dislocation velocity in the silica-containing ice is two orders of magnitude lower than that in the pure ice. The delay in the dislocation velocity could be related to the presence of nanoscale silica particles, which could also be dispersed on the basal planes. At the creep temperature dislocation glide and climb should occur simultaneously. It is difficult to consider uniform distributions of the silica in the ice crystals, but some segregations can easily be expected, as such uniformly dispersed particles are usually unrealistic in the case of composite materials (for example, Homma and others, Reference Homma, Nagai, Katou, Arai, Suganuma and Kamado2009). This means that dislocations might be pinned at the densely distributed silica regions. Such regions can be seen in the 10% and 20% strained samples shown in Figure 3. The orders of magnitude of the dislocation velocities are smaller than those reported by other researchers, although they are calculated or measured (Perez and others, Reference Perez, Maï and Vassoille1978; Higashi, Reference Higashi1988). The reported minimum velocity was around 10−6 to 10−9 m s−1 (Forouhi and Bloomer, Reference Forouhi and Bloomer1978; Whitworth, Reference Whitworth1978) using single crystals and X-ray diffraction topography methods. This indicates that since dislocation densities in the single crystal ice are lower, dislocation tangles in the polycrystalline artificial ice sufficiently reduce the velocities.
In the case of the steady-state creep with a constant-stress condition the steady-state creep rate is expressed as
\begin{equation}\dot \varepsilon = B{\sigma ^m}{\text{exp}}\left( { - \frac{Q}{{kT}}} \right)\end{equation} according to Weertman (Weertman, Reference Weertman, Whalley, Jones and Gold1973). Here 
$B$ is a constant, 
$\sigma $ the applied stress, 
$Q\,$an activation energy, 
$k$ Boltzmann’s constant and 
$T$ the absolute temperature. Weertman then estimated:
\begin{equation}\rho = {\left( {\frac{{\beta \sigma }}{{\mu b}}} \right)^2}\end{equation} if the average internal stress generated by dislocation tangles is equivalent to dislocation densities, where 
$\beta $ is a constant. The velocity of dislocations is given by
\begin{equation}v = \gamma \sigma {\text{exp}}\left( { - \frac{Q}{{kT}}} \right)\end{equation} where 
$\gamma $ is a constant. Glen identified the stress exponent as 3.7 (Reference Glen1955), which is slightly larger than that reported by Weertman (Reference Weertman, Whalley, Jones and Gold1973). Based on the equations we understand that the dislocation density is one of the important factors in determining the steady-state creep rate. Using Eq. (17) the constant 
$\gamma $ for the pure and silica-containing ice can roughly be determined using activation energies reported by Saruya and others (Reference Saruya, Nakajima, Takata, Homma, Azuma and Goto-Azuma2019). These are 60 and 66 kJmol−1 for the pure ice and 0.1% silica-containing ice, respectively. By replacing Boltzmann’s constant with the gas constant we obtain the 
$\gamma $ values in Eq. (17) as 3.8 × 10−6 and 8.8 × 10−7 m Pa−1 s−1 for the pure and silica-containing ice, respectively. These results suggest that the dislocation velocities with microparticles should almost always be lower than those in pure ice at the stead-state creep region in our experimental conditions.
4. Conclusions
The cold stage and sample holder for X-ray diffraction measurements of ice were developed using a nitrogen gas combined with liquid nitrogen. The obtained X-ray diffraction patterns were analysed based on the modified Warren-Averbach (WA) and modified Williamson-Hall (WH) plots of X-ray line broadening; this is the first report to accomplish dislocation measurements with the modified WH and WA methods for ice samples.
Dislocation densities as a function of strains were measured from the modified WA and WH plots using artificial ice with and without silica. The silica models the existence of impurities in ice sheets. The X-ray diffraction system was set in a room maintained at ambient temperature. The creep tests were conducted at −20ºC and 2 MPa.
In the initial ice large amounts of dislocations are stored due to the sintering process, which uses 70 MPa of the sintering pressure. In the primary creep region dislocation densities increase with the increase in the grain size, particularly in pure ice. Since the aggregates of silica act as pinning sites for grain growth, the grain sizes in the silica-containing ice are relatively smaller than those in the pure ice. The maximum dislocation densities are achieved at 10% strain, with high subgrain boundary densities, while at 20% strain, the density decreases. However the dislocation densities in the pure ice remain almost constant values. The silica dispersions act as pinning sites for the suppression of grain growth. Such small grain sizes may introduce sinks for dislocation annihilations.
The small grain size due to the silica dispersions increases dislocation densities because of the pile-up at grain boundaries at 10% strain. The dislocation velocity was roughly calculated based on the Orowan equation. As can easily be estimated, the silica addition retards the dislocation velocity compared with that in pure ice due to the grain refinement and large number of dislocations.
Acknowledgements
We acknowledge all the members of the Arctic Challenge for Sustainability (ArCS). This research was supported by ArCS (program grant number JPMXD1300000000), JSPS KAKENHI (grant numbers JP15K13567, JP17H02957 and JP18H04140), Japan Arctic Research Network Center (J-ARC Net) in 2017 and NIPR through General Collaboration Project, Japan no. 27-16 and 30-13.
We also acknowledge Dr. Morimasa Takata of Nagaoka University of Technology, Japan for his enthusiastic mentorship regarding the research and developments of the instruments for X-ray diffraction measurements and cooling systems at low temperature.
Author contributions
This study was conceptualized by Nobuhiko Azuma and Tomoyuki Homma with guidance from Kumiko Goto-Azuma. Tomoyuki Homma led the experimental and theoretical work, directed experimental plans and co-authored the paper with Kumiko Goto-Azuma and Nobuhiko Azuma. Tomoyuki Homma and Kazuteru Hirai developed the systems for X-ray diffraction measurements, including the cooling systems and conducted the experimental work. Particularly Kazuteru Hirai was responsible for the data analysis. Kumiko Goto-Azuma and Nobuhiko Azuma were supervised, reviewed and edited. All authors contributed to the discussions of the results and drafting of the manuscript.
Open access
