Mechanical properties of bi- and poly-crystalline ice

A sound knowledge of fundamental mechanical properties of water ice is of crucial importance to address a wide range of applications in earth science, engineering, as well as ice sculpture and winter sports, such as ice skating, ice ﬁshing, ice climb-ing, bobsleighs, and so on. Here, we report large-scale molecular dynamics (MD) simulations of mechanical properties of bi- and poly-crystalline hexagonal ice (I h ) under mechanical loads. Results show that bicrystals, upon tension, exhibit either brittle or ductile fracture, depending on the microstructure of grain boundaries (GBs), whereas they show ductile fracture by amorphization and crystallographic slips emit-ted from GBs under compression. Under shearing, the strength of bicrystals exhibits a characteristic plateau or sawtooth behavior drawn out the initial elastic strains. Nanograined polycrystals are destabilized by strain-induced amorphization and col-lective GB sliding. Their mechanical responses depend on the grain size. Both tensile and compressive strengths decrease as grain size decreases, showing inverse Hall-Petch weakening behavior. Large fraction of amorphous water structure in polycrystals with small grain size is mainly responsible for the inverse Hall-Petch softening. Dislocation nucleation and propagation are also identiﬁed in nanograined ice, which is in good agreement with experimental measurements. Beyond the elastic strain, a combination of GB sliding, grain rotation, amorphization and recrystallization, phase transformation, and dislocation nucleation dominate the plastic deformation in both bicrystals and polycrystals. © 2018 Author(s). All article content, except where oth-erwise noted,


I. INTRODUCTION
Water is one of the most essential substances to life on Earth. A water molecule is composed of one oxygen atom and two hydrogen atoms, and the oxygen atom is covalently bonded to the two hydrogen atoms. The two hydrogen atoms of each water molecule can be weakly bonded to other two water molecules in hydrogen bond form, respectively. Upon cooling, water tends to freeze into ice. So far, over 15 different solid-state crystal structures of water, including ice XVI, XVII from clathrates and ice XVIII, have been identified with the aid of experimental and computational at 263.15 K and a strain rate of 5.5 × 10 -4 s -1 . High values of peak stress were resulted for sample size/crystal size ratios less than 12. 62,63 Furthermore, deformation behaviors and rheology of ice have been investigated in the past few decades. 64 By uniaxial compressive experiments on polycrystalline ice, Jacka et al. demonstrated that the minimum flow rate of isotropic polycrystalline ice was of little or no crystal size dependence. 65 Visible cracks have not been found in laboratory-prepared randomly-oriented ice at low strain rates. 66 Instead, Duval et al. 67 reported that the diffusional flow and the transient behavior in ice show grain size-dependence. More than four independent slip systems govern the compatible deformation of ice polycrystals. Both deformation of the grains and slips at the GBs were found to contribute to the deformation process of ice at temperature ranging from 233.15 to 270.15 K. 51 Muguruma 68 reported that GBs were a prominent source of dislocations in the early stage of deformation of columnar-grained ice, and the maximum stress of columnar-grained ice was grain size-dependent. Under compressive stresses ranging from around 0.1 to 1 MPa at a temperature from 263.15 K to the melting-point, a transient creep component and a continuing or quasi-viscous component have been found by Glen 69 in ice creep. In addition, ice rheology at low stress was also demonstrated to exhibit a strong dependence on the interactions between ice crystals. 70 Moreover, ice properties have also been studied in situ for the analysis of flow of glaciers and ice sheets. 71,72 Creep experiments on fine-grained ice 73 demonstrated that superplastic flow of ice at stresses less than 0.1 MPa governs the rate-limiting creep mechanism under a wide range of temperatures and grain sizes.
Despite many efforts made, our understanding on mechanical properties of ice still remains very limited, particularly, on the deformation mechanisms, strain-induced structural transition and dislocations at the molecular level. As the Hall-Petch relationship is established experimentally in polycrystalline ice with millimeter-sized grains, it is interesting to understand whether this Hall-Petch strengthening behavior is expected to continue or switch to an inverse Hall-Petch softening in nanograined ice (grain size ranging from 1 to 100 nm). However, experimental investigation is challenging because of the difficulties in preparing high-quality nanograined samples and experimental nano-visualizing of structure of GBs and deformation responses. Classic molecular dynamics (MD) simulations as an important tool have been successfully applied to understand formation mechanisms, thermodynamic and mechanical properties of a variety of solid-state materials at the molecular level, 9,16,74,75 This study provides critical knowledge of mechanical properties of bi-and poly-crystalline I h under mechanical loadings. The failures and plastic deformation mechanisms, such as GB sliding, GB accommodation mechanisms and intragranular deformation, are elucidated by classic MD simulations at nanoscale (1 -100 nm). These results are of key importance in glacier studies, ice engineering, and frozen ground mechanics and operations in polar area. Moreover, this work provides the molecular insight to understand the differences in mechanical mechanisms of ice, clathrate hydrate and ice-like water-dominant materials that often coexist with ice in the permafrost regions.

A. Atomic models
The initial atomic structure of monocrystalline I h was obtained by X-ray diffraction results. 76 One unit cell of monocrystalline I h is composed of 16 water molecules. Figures 1a-1c present three geometrical configurations of I h with different crystal planes highlighted by green. Monocrystalline I h shows hexagonal plates and columns with top and bottom faces those are basal planes {0001}, and 6-equivalent side faces which are called the prism faces a. Additionally, planes in I h formed by the sides of the chair-like structure are the secondary prism faces {1120}. Because of the complexity of GB structures in polycrystals, three distinct bicrystals formed by two different single crystals were specifically created to reveal the mechanical characteristics of GBs under tensile, compressive and shearing loads. Figure 1f presents a typical side-viewed relaxed bicrystalline structure with a flat and cohesive GB that is formed by joining two monocrystals with {1120} and {0001} planes. Five-and seven-membered rings identified at GB were green-and yellow-highlighted, respectively. To reveal the effect of grain size on mechanical properties of polycrystalline I h , five samples containing a number of monocrystalline grains with average grain sizes ranging from 4.54 to 15.87 nm were constructed based on a Voronoi construction. 77,78 All polycrystals had identical dimensions of 40 × 40 × 40 nm 3 . The number of grains determined by grain size in polycrystals ranged from 16 to 686. The average grain size in polycrystals is defined as follows: where L and N are the length of samples (nm) and the number of grains in polycrystals, respectively. Figure 1d shows a representative three-dimensional (3D) molecular configuration of polycrystalline I h with an average grain size of 15.87 nm, and each grain is colored for clarification. The orientations of grains in polycrystals are randomly distributed. Figure 1e displays one typical grain with a geometry of the body-centered cubic Wigner-Seitz cell from the polycrystal. To avoid artificial particle overlaps in the polycrystals, molecules that are closer than 0.1 nm to any adjacent one in the neighboring grains were removed. To facilitate further statistical analysis, four polycrystalline samples with identical average grain size but different GBs were generated. Each polycrystal was composed of approximately 1,930,000 water molecules.

B. Forcefields
A mono-atomic model was adopted to describe water molecules. The tetrahedral short-ranged interactions among mono-atomic water are modeled by the available Stillinger-Weber (SW) model. 79 This coarse-grained water model 80,81 is over 2 orders of magnitude more computationally efficient than a fully atomistic model in reproducing a range of properties of liquid and solid phases of water. 80 Moreover, this model was successfully applied to investigate the mechanical properties of methane

C. Mechanical MD tests
Prior to MD calculations of uniaxial straining, all molecular structures were quasi-statically relaxed to a local minimum energy configuration by the conjugate gradient method with an energy tolerance of 1.0× 0 -4 eV and a force tolerance of 1.0×10 -4 eV/Å. 83 MD relaxations were conducted with 1.0×10 6 timesteps at 223.15 K and atmospheric pressure in the NPT (constant number of particles, constant pressure, and constant temperature) ensemble, and then with another 1.0×10 6 timesteps at 223.15 K under NVT (constant number of particles, constant volume, and constant temperature) ensemble. The Nosé-Hoover thermostat and barostat method with damping times of 0.1 and 0.5 ps were used to control the temperature and pressure, respectively. Our MD calculations allowed ice structures to relax in order to obtain the stable structures. The velocity-Verlet integration algorithm was employed to integrate the equations with timestep of 1 fs. 3D periodic boundary conditions were applied in three directions to obtain data. Initial velocities of water particles in the systems were assigned following a uniform distribution according to the given temperature. In this study, the tensile (compressive) strain is defined as where l 0 and l are the original length and the final length of ice sample along the deformational direction, respectively. For mechanical tests, the uniaxial loadings were conducted by the deformation control technique in the NpT ensemble. This procedure was carried out on the relaxed samples with a constant strain rate of 1.0×10 8 s −1 by uniformly rescaling the z-coordinates of all molecules every 1000 timesteps. Deformational MD simulations were performed in a modified NpT ensemble, specifically, NVT in the loading direction but NpT in the lateral directions. The deformational MD simulations were performed in NL z p x p y T ensemble, where L z is the length of simulation box along the loading direction. The Nosé-Hoover anisotropic barostat and thermostat was employed in NL z p x p y T ensemble to control pressure only in the lateral directions independently and temperature. This allowed ice samples to experience expansion/contraction in transverse directions due to the Poisson's effect. The stress of water particle was calculated based on the virial definition of stress using forces on the particles collected during MD simulations. Both stress and potential energy of molecules in systems were averaged over 5000 timesteps to eliminate thermal oscillations. All MD calculations were performed by using the Large-scale Atomic-Molecular Massively Parallel Simulator (LAMMPS) software package. 84 To identify and visualize dislocations and stacking faults of hexagonal structures during MD simulations, an extended common neighbor analysis (CNA) method 85 was employed. The extended CNA is able to classify pairs of particles based on their local environments. Furthermore, the topology of water rings in GBs was also identified by using the "shortest path ring" algorithm developed by Franzblau. 86 The development of dislocation activities during deformation was also observed. The dislocation results were obtained by means of the Dislocation Extraction Algorithm (DXA) method. 87

Mechanical properties of bicrystals under both tension and compression
Mechanical properties of bicrystalline I h with three distinct GB structures under both tension and compression are examined. Figure 2a shows the simulated stress-strain curves of bicrystals with typical {1120}|{00 0 1}, {1120}|{0110} and {00 0 1}|{0110} GB structures subjected to tension along the directions perpendicular to GBs. Bicrystals show initially a linear response, followed by a nonlinear elastic response. However, differences in their peak tensile strengths are detected. The peak tensile strength from the highest to the lowest are sorted as {1120}|{0110} > {1120}|{00 0 1} > {00 0 1}|{0110}, which indicates different cohesive energy of those GB structures. For bicrystals with {1120}|{00 0 1} and {00 0 1}|{0110} GB structures, tensile stresses almost drop to zero after the peak strength, suggesting a brittle fracture mechanism. For the one with {1120}|{0110} GB structure, however, it shows a ductile fracture pattern as indicated by Figure  at around 340 MPa after the sudden drop of stress. The fracture strains are also found sensitive to bicrystalline structures. Upon compression, as shown in Figure 2b, the mechanical response is similar to the case of bicrystal with {1120}|{0110} GB structure under tension. The order of peak compressive strength is found identical to that of peak tensile strength. To understand the differences in ultimate mechanical strength, the energy of selected crystal surface that forms the GB structures and GB binding energy are accordingly calculated, as presented in Figures 2c and d. Apparently, the {1120} plane has the lowest energy, while the {0110} plane shows the highest energy. For binding energy of GBs, their order is found to be {00 0 1}|{0110} > {1120}|{00 0 1} > {1120}|{0110}, which is consistent with the calculated ultimate mechanical strengths. In addition, the Young's moduli of the bicrystal determined from initial slopes of stress-strain curves are around 7.6-9.1 GPa.
To unravel tensile and compressive deformation mechanisms, a number of typically localized snapshots of bicrytalline I h are captured and shown in Figure 3. By examining the microstructures of bicrystals at equilibrium state (zero-strain), it is observed that all water molecules in GBs are identified and GB structures have comparable potential energy to the bulk counterparts, indicating highly cohesive GB structures. Moreover, structural defects, 88 including five-and seven-membered rings that commonly exist in graphene and diamond, are observed in GBs. These defects also explain the superior stiffness and mechanical strength of bicrystals. When the tensile strain exceeds the critical values, molecular structures of bicrystals with {1120}|{00 0 1} and {00 0 1}|{0110} GBs destabilize, as demonstrated by Figures 3a-3, a-7, c-3 and c-7. This corresponds to a sharp drop of stress in stress-strain curves of Figure 2a. Further strain causes the complete decohesion of GBs, suggesting brittle fracture mechanism. For a bicrystal with {1120}|{0110} GB, however, a ductile fracture mechanism is observed, as illustrated by Figure 3b. This plastic deformation mechanism differs from those of monocrystal. 82 The bicrystal destabilizes by crystallographic slip along the (1120) crystal plane that emits from one site of the pre-existed GB. This results in the formation of a new GB that is less cohesive than pre-existed GB. As a result of its small cohesive force, the following plastic deformation is controlled by sliding of this new GB accompanied by recrystallization and dislocation AIP Advances 8, 125108 (2018) at GBs. The two distinct fracture patterns are mainly attributed to the difference in crystal face binding energy. Upon compression, all bicrystals show ductile fracture behaviors. As shown in Figure 3d, the plastic deformation mechanism is similar to that of bicrystal with {1120}|{0110} GB structures under tension.

Mechanical properties of bicrystals under shear strain
Shearing responses of bicrystals are also key mechanical properties of interest in solid-state materials. It is instructive to examine mechanical properties of bicrystalline I h under pure shear deformation in GBs planes. Pressures in the directions that are perpendicular to shear directions are controlled to be ambient pressure during the shear deformations. Figure 4 shows the simulated shear mechanical properties of bicrystals with different GBs along two nonequivalent shear directions. Obviously, the features of all shearing curves differ from those of tensile and compressive ones. By comparing MD shearing results of bicrystals, bicrystals with {1120}|{00 0 1} and {00 0 1}|{0110} GBs show very limited elasticity and small first peak shear strengths, indicating their least cohesive forces of GBs. In sharp contrast, for bicrystal with {1120}|{0110} GB, the similarity between elastic shear, tension and compression strains and its large first peak shear strengths indicate the large cohesive forces of {1120}|{0110} GB. Those agree with the results of tensile and compressive MD simulations. Beyond yield points, three distinct categories of shear stress-strain behaviors are classified, depending on bicrystalline structures and shearing directions. Upon To elucidate the shear deformation mechanism of bicrystalline I h , the associated molecular deformation modes are analyzed by examining the microstructural evolution during shearing process.  (Figures 5c-4 and -6) are responsible for the nearly constant large shear strength. Formation and destruction of five-and seven-membered defect pair was also identified as the controlling step in homogeneous ice melting, 88,97,98 similar to other solid materials. 99 It is also found in the nucleation and growth of ice, 100  directions easily break metastable hydrogen bonds, resulting in early GB sliding. Instead, dislocation-free region is found during the entire deformation process (Figure 5e-f). In short, both elasticity and plasticity in bicrystals strongly depend on microstructures of GBs and shearing directions. Mechanical plastic deformation of bicrystals under two nonequivalent shear directions involves GB sliding, the competition mechanism among decomposition, recrystallization, and dislocations.

E. Mechanical properties of polycrystals
For conventional polycrystalline solid materials such as metals, carbon materials and methane hydrates, grain size plays a key role in their mechanical behaviors. In natural and laboratory settings, I h was commonly observed to be grown well with polycrystalline texture. MD simulations of polycrystalline I h with grain size varying from 4.54-15.87 nm under both tension and compression are performed to investigate mechanical properties. Figure 6a shows the initial and MD relaxed molecular structures of polycrystals with two different grain sizes. It is found from snapshots of #1 and 2 that, after MD relaxations, GBs are able to maintain their flat plane. However, besides hexagonal and amorphous structures, a cubic phase of ice structures is observed at GBs. As shown in snapshots of #3 and 4 of Figure 6a, GBs exhibit higher potential energy than that of the grain interior, indicating their metastability. By analyzing water topological rings at GBs, various structural defects containing four-, five-and seven-membered rings are found, as presented in Figure 6b.
A series of uniaxial tensile and compressive loadings on polycrystalline I h at 223.15 K and ambient pressure were carried out. Figures 7a and b show the MD simulated mechanical stress-strain curves of polycrystals with fine grain size ranging from 4.54-15.87 nm under both tension and compression, respectively. Apparently, all polycrystals show unique mechanical responses that strikingly differ from those of bicrystalline counterparts. Particularly, characteristic features of nonlinear stress-strain curves are remarkably sensitive to the grain size under both tension and compression. Three deformational stages can be roughly observed. At the first stage, all polycrystals in all tests show linear loading responses, corresponding to the initially limited elastic deformations. Young's modulus obtained by fitting the stress-strain curves in this stage varies from 4.9 to 7.7 GPa with grain size ranging from 4.54 to 15.87 nm. The values of elastic moduli are lower than those of bicrystals in our simulations, as well as those of experimental measurements. 101,102 This is mainly attributed to differences in microstructure of samples. For example, bicrystals show infinite long GBs, while polycrystals show finite long GBs and complex triple GB junctions. Crystalline grains of ice in laboratory, engineering or nature are typically millimeter-and micron-sized, 31,47,62,65,103,104 which are orders of magnitude larger than these of polycrystalline ice in the present study. Similar reduction characteristics have been also identified in metals, diamond and methane hydrates. 82,105,106 This indicates that the quality of ice can be dictated by the grain size dependent mechanical stiffness, which indicates its useful applications in winter sports such as ice skating, ice fishing, ice climbing, and bobsleighs. Moreover, in experimental measurements, the used strain rates ranging from 10 -7 to 10 -2 s -1 31,47,103,104 were  lower than that of our MD studies. The strain rates in our MD simulations may be comparable to hail stone impacts or winter sports. The second deformation stage corresponds to nonlinear mechanical responses prior to the yield points. In this stage, mechanical strength strongly depends on grain size, although all polycrystals show monotonic increasing strength, which is indicative of different deformation elastoplasticity in polycrystal with different grain size beyond the early small elastic deformations. In addition, both yield tensile and compressive strengths are also greatly grain sizedependent, as presented in Figure 7c. Intriguingly, under tension and compression, polycrystalline I h exhibits inverse Hall-Petch behavior. The yield stress monotonically increases as the fine grain size increases from 4.54 to 15.87 nm. Experimentally, strong dependence of mechanical strength on grain size in macroscopically grain-textured polycrystalline ice was also observed at various measurement conditions. 31,47,62,63,103,104 For example, it was reported that polycrystalline ice with millimeter-scale grains of 1.5-5.8 mm are able to follow either Hall-Petch or inverse Hall-Petch behavior, depending on the measured strain rates. 47 Currier et al. 31 showed that, under tensile experiments with strain rate of 1.0×10 -6 s -1 at 263.15 K, the fracture strength of millimeter-sized grained ice reduces as the grain size is increased. However, Cole et al. 47 reported that the peak stress increases with increasing grain sizes when the applied strain rate is less than a critical value of around 4.0×10 -6 s -1 . The conflicting findings from different experimental measurements suggest that deformation mechanisms of polycrystals with large grain size are sensitive to strain rate. This is indeed generally accepted in glaciology, 107 where dislocation creep is assumed to be the dominant mechanism in glaciers and ice sheets, but grain-size sensitive flow can occur at small grain sizes or at very slow strain rates. 73,107 Although the experimentally observed inverse Hall-Petch relation is in accordance with our findings, the applied strain rates are different between experiments and our MD simulations. Moreover, opposite observations between experiments (Hall-Petch) and our calculations (inverse Hall-Petch) can be mainly attributed to disparity in grain size. Previous reports via experiments or computer simulations showed very small critical grain sizes at the transition between Hall-Petch and inverse Hall-Petch relationships in metals, metallic alloys and methane hydrates. 82,108,109 It is therefore expected that the mechanical strength of extremely fine nanograined ice follows inverse Hall-Petch relationship as predicted by our simulations. For the last stage that corresponds to the deformation responses after the yield points, it is observed that the mechanical characteristics of nanograined ice are also strongly grain size dependent. Under tension, all polycrystals show strain-induced softening behavior, which becomes more pronounced for ice with larger grain size. Such strain softening is commonly observed in glassy organic materials. [110][111][112] The deformation responses are associated with changes in microstructures. Posterior to strain-softening, ice with grain size of 4.54 and 6.35 nm exhibits evidently steady flow stress, whereas that with grain size of 7.94 nm exhibits monotonic strain-hardening drawn out the first at strain of 0.15. Particularly, tensile strengths in polycrystals with two large grains increase, and then decline with increasing strain, reflecting a transition from strain-hardening to strain-softening, consistent with the previous experimental data. 67 Upon compression, loading responses of polycrystal with the smallest grain are similar to the case under tension. Polycrystals with grain size of 6.35 and 7.94 nm exhibit monotonic strain hardening that continues over the predefined strain of 0.25. Such strain hardening persists to a large degree relative to most metals and ceramics at high temperatures. In contrast, polycrystals with large grains display ultimate compressive strengths at critical strains followed by apparent strain-softening behavior. Figure 7d shows the average flow stress as a function of grain size. Similar to yield stress, the average flow stress increases as the grain size increases, which indicates inverse Hall-Petch behavior. Beyond ice, natural gas hydrates occurring in deep-sea sediments or permafrost regions are also icy crystalline substances. Figure 10 shows a comparison between ice and methane hydrates, 31,47,82,[113][114][115][116][117] and both solids exhibit grain size and temperature dependent mechanical properties. However, differences can be also identified, attributed to the difference in the intrinsic water-dominated framework structure.
For some solid materials (metals, black phosphorus and methane hydrates), a large fraction of disorder molecules at GBs of polycrystals with extremely fine grains is responsible for the inverse Hall-Petch softening behavior. 82,105,118,119 Accordingly, variations of the fraction of amorphous water in polycrystals with imposed tensile and compressive strains are calculated and plotted in Figures 7e and f, respectively. At equilibrated state (zero strain), polycrystal with small grain size shows the largest fraction of disorder water particles. The existence of amorphous water molecules in GBs weakens polycrystals as indicated by the difference in mechanical properties between bicrystals and polycrystals. Similarly, other solid materials exhibit decreasing mechanical strength due to the high defect density in GBs. 99 Upon loadings, all polycrystals show a nonlinear behavior for the fraction of disorder water molecules as a function of the deformation strain. The change tendency in the fraction of disordered water molecules is closely connected with deformation stages. In the first stage, no change in the fraction of disorder water molecules is observed, reflecting limited elastic responses. In the second stage, the fraction of disordered water molecules increases as the strain is increased, implying amorphization of crystalline water in polycrystals to accommodate the imposed strain. Such strain-induced amorphization does not lead to the decrease in strength but degradation in stiffness. Polycrystals are found to amorphize at the maximum damage rate as the stress reaches the yield point. Amorphization continues at a decreasing damage rate to critical points, which corresponds to the first part of strain-softening periods. In the remaining part of strain softening, for polycrystals with small grain size under tension, the fraction of disordered water tend to become a constant. Finally, slight reductions in disordered water are observed in the late deformation stage, implying crystallization occurring in polycrystals. In contrast, polycrystals with large grain size show an apparent decrease of disordered water in the following strain-softening stage, indicating large-scale recrystallization in this phase. Such significant recrystallization strengthens deformed polycrystals, resulting in strain hardening as shown in Figure 7a. However, polycrystals re-amorphize in the late deformation stage, as evidenced by the increasing fraction of disordered water in Figure 7e. Upon pressurization, polycrystals exhibit a similar trend for the fraction of disordered water in the late part of strain softening, except for polycrystal with a grain size of 10.58 nm that display an increase for the fraction of disordered water. In contrast to the case of tension, a linear enhancement in the fraction of disordered water is observed for all polycrstals in the late straining, while polycrystals with large grain sizes demonstrate strain-softening behavior. This indicates that different deformation mechanisms govern the subsequent over-straining in polycrystals. The issue of an amorphous phase and failure strength in nanograined polycrystalline ice under extremely high strain rate is of importance for winter sports (ice-skating, bobsleighs, etc.), as certainly engineering and natural hazards if one thinks of hail stone impacts.
To further elucidate grain size-dependent mechanical responses from microscopic structures of polycrystals, a number of localized snapshots of polycrystals with two grain sizes upon different strains are captured and shown in Figure 8. Apparently, local crystalline textures in polycrystals change with the increase of both tensile and compressive strains. Figures 8a and 8b display the evolution of cross-sectional morphology of polycrystals with grain size of 7.94 and 15.87 nm under tension. Prior to yield strains, there is no GB sliding as evidenced by the highlighted solid segment which crosses GB. Crystalline grains are amorphous at the external surface to react to mechanical straining as indicated by the increasing fraction of disordered water (Figure 7e). Furthermore, minor structural transition from hexagonal to cubic phase locally occurs at high potential energy GBs. When the imposed strains exceed the yield values, collective GB sliding takes place to relax the strain energy as visually confirmed by the breaking of highlighted solid segment. This explains the observed strain-softening behavior as shown in Figure 7a. For polycrsytals with large grains, nanovoids are clearly nucleated at GB-surface multi-junctions, as a consequence of GB sliding that is non-accommodated. However, because of short-distance diffusion from GB to junctions, GB sliding is effectively accommodated by GB diffusion process, which leads to no apparent formation of nanovoids until large strains are imposed. Such nanovoid nucleation has been also observed in other solid materials. [120][121][122] At large strain of 0.24, minor cubic structure of ice is observed in the vicinity of multi-junctions in polycrystals with small grain size (Figure 8a-3), whereas for that with the largest grain size, a locally hexagonal-to-cubic structural transition occurring from multi-junction to interior grain is identified (Figure 8b-5). Upon compression, similar to the case of tension, no GB sliding occurs when the strains are below yield strains. However, changes in microstructure of GB are observed as shown by Figure 8c-1 and -2. Likewise, cooperative GB sliding happens to release the excess strain energy when the strains are beyond the yield points. In contrast to the case of tension, such large-scale GB sliding does not cause nucleation of nanovoids. With increase of compressive strain, structure geometry of some grains in polycrystals is greatly changed. As a consequence of high stress-concentration located at GBs that are perpendicular to loading direction, grains in polycrystals are greatly destabilized by amorphization, recrystallization, and phase transformation at GBs accompanied by roughening of the smooth GBs interfaces and rotation of grains (Figures 8c-3 and -6, d-3 and -6). Metals, upon heating and mechanical loading, are able to undergo such GB interface roughening transition. [123][124][125] Particularly, polycrystals with large grains show dislocation nucleation at the interior grains during AIP Advances 8, 125108 (2018)  Figure 7b. This is the main source of the well-known classic Hall-Petch softening behavior in metals. 126 Experimentally, Schulson et al. 127 also showed that similar transgranular microcrack occurs in polycrystals composed of equiaxed and randomly oriented aggregates of granular ice, but those with small grains do not exhibit transgranular fracture behavior. Moreover, Duval et al. 67 reported that crystallographic slips govern mechanical deformation of polycrystalline ice and distribution of internal stress.
Some studies have been conducted to understand ice plasticity related to dislocation behaviors, 126,[128][129][130][131][132][133][134][135][136] and dislocations have been richly studied in conventional solid materials. [137][138][139][140][141] To provide more information on dislocation in nanograined ice, snapshots of dislocation structures are captured. Figure 9 presents the simulated dislocation structures of polycrystalline ice with grain size of 15.87 nm. The results show rich dislocation behaviors in polycrystalline ice, consistent with the experimental observations. 126,[128][129][130][131][132] A competition between dislocation-nucleation and dislocationannihilation is found to occur during the entire deformation process. From Figure 9a, it is found that plenty of dislocations are uniformly distributed at GBs in the early deformation stage. However, in the strain-softening stage, dislocation nucleation activity is greatly weakened, but major 1/3 1 2 1 0 and minor 1/3 1 1 0 0 dislocation develop to form complex network topology in the entire polycrystal, as shown in Figure 9a-2. Using synchrotron X-ray topography, Liu et al. 130 also revealed that dislocations mainly nucleate at GBs of polycrystal with milimeter-sized grains, resulting in stiffness-softening. The change in dislocation type and location of dislocation nucleation implies the differences in plastic deformation modes between the two deformation stages. Figure 9b shows the detailed development of dislocation occurring in one grain. Initially, dislocations nucleate at the top and bottom damaged GBs that are perpendicular to loading direction. The top dislocations propagate and intersect to form a network motif, accompanied by new dislocations. Once the interaction between dislocation lines is turned on, dislocations come into being in the interior of grain. Meanwhile, the generation of new dislocations is also initiated at the bottom GB, and those dislocations propagate towards the interior of grain. Finally, dislocations coalesce to form complex dislocation network in the interior of grain to accommodate crystal deformation. However, among a variety of dislocations occurring in the interior of ice grain, screw dislocations with Burgers vector 1/3 1 2 1 0 are the most frequently observed as shown in Figure 9b. The similar dislocation behaviors have also been observed in hexagonal close-packed (HCP) solid materials. [137][138][139][140][141] Shōji et al. 142 reported that dislocation glide region in polycrystalline ice with grain size of 2 millimeters can be divided into dislocation glide with crack formation, dislocation creep, and fracture at high stresses. Similar to our findings, plastic deformation in twining-free polycrystalline metals with extremely fine grains was mostly governed by GBs. 105,[143][144][145] Interplays between dislocations and GBs are utilized in GB engineering. 146 However, for nano-twinned polycrystalline metals, transition in plastic deformation mechanisms from dislocation pile-up and cutting through twin planes to a dislocation-nucleation-controlled softening mechanism with twin-boundary migration was revealed. 147 Dislocation behavior in polycrsytalline ice appears very complex, attributed to different mechanisms of generation and annihilation of dislocation. Importantly, plastic responses greatly depend on dislocation behaviors, which contain AIP Advances 8, 125108 (2018)  31,47,82,[113][114][115][116][117] mutual interactions between dislocations, interactions between dislocations and other defects such as GBs.

III. CONCLUSIONS
In summary, large-scale classical MD simulations with the coarse-grained mW water model were performed to reveal failure and deformation mechanism of bi-, and poly-crystalline I h under mechanical loads. For bicrystals, a series of tension, compression and shear MD tests reveal that mechanical characteristics, such as elastic moduli, ultimate strength, failure modes and dislocation behaviors, strongly depend on GB structures that are determined by the plane of jointed monocrystals. MD simulations of polycrystals with extremely fine grains subjected to both tensile and compressive loads demonstrate pronounced grain size effect on mechanical behaviors. Both the yield stress and flow stress decrease with grain size under the simulation conditions, which is similar to the inverse Hall-Petch behavior observed in other polycrystalline textured materials. The high ratio of amorphous crystalline water molecules in polycrystal with small grain size is responsible for the inverse Hall-Petch softening. Mechanical plastic deformation of ice polycrystal is found to be cooperatively governed by strain-induced amorphization and recrystallization, GB sliding, dislocation nucleations, phase transformation, GB roughening and grain rotation. This study provides key insights into the structural, failure and plasticity of solid water ice, shedding light on the distinct deformation mechanism of water ice at the nanoscale.