Carbon dioxide enhances fragility of ice crystals

Ice caps and glaciers cover 7% of our planet, greater than the land area of Europe and North America combined [1]. They reflect 80–90% of the solar radiation and trap a large amount of carbon dioxide. Specifically, the Arctic accounts for 10–15% of the earth's carbon sink [2]. The glaciers' dynamical behaviour and stability play important roles in controlling the global climate [3, 4]. Thereby, the facture mechanism of ice is of paramount importance for the understanding of glacial dynamics [5], and at a fundamental level is controlled by how a single crack propagates in ice crystals via the breaking of chemical bonds [6]. Such growth of cracks eventually mediates the breakdown of ice, as exemplified in recent incidents of large-scale fracture of glaciers [7, 8]. Very large-scale ice fractures occurred recently close to the Pine Island Glacier in the Antarctic region, which generated an iceberg with an area of around 880 square kilometres, the size of the city of Berlin [9]. While the macroscopic mechanical properties of pure ice are well understood by laboratory tests and its behaviour has been characterized by existing fracture mechanics models [10], the effect of environmental conditions such as the concentration of carbon dioxide, have not been incorporated into existing models. Lack of such knowledge prevents us from understanding how changes in the chemical environment affect the stability and movement of glaciers, which is important given the rising levels of carbon dioxide in the atmosphere.

It is known that the chemical environment can have significant effects on the mechanical properties of ice. The ions of salts contained in sea water cause the ice to form internal brine pockets, enabling ice to deform and rupture more easily because emerging vacancies lead to local stress concentrations [11]. It has also been argued that air bubbles may lead to a similar effect, but the mechanism remains elusive since there has been a shortage of measurements of this effect in experimental work, especially for the difficulty to obtain the density, orientation and size distributions of the bubbles [12]. A recent study based on world-wide ice core records suggested that the rising concentration of carbon dioxide in the atmosphere may have been closely associated with the last deglaciation period [13, 14]. However, it is not clear if increased temperatures caused by the presence of carbon dioxide is the only reason to cause deglaciation. The oceans are huge reservoirs for carbon dioxide (CO₂) as it is soluble [15, 16]. The average molar concentration of carbon dioxide in polar ice caps is about 0.02% [17, 18], and the value reaches 0.1% for deep sea ice [18]. The two oxygen atoms in CO₂ have a partial negative charge, which can interact with the positive charged hydrogen atoms in water molecules and thus facilitate the formation of surrounding water clusters [19]. This evidence suggests that the effects of the presence of carbon dioxide in ice may not only be to induce stress concentrations by forming bubbles [20], but may also have more fundamental implications by affecting its fracture property from the atomistic level upwards.

Here we investigate the role of carbon dioxide in ice fracture based on molecular dynamics simulations carried out using the first-principles based ReaxFF reactive force field [21] implemented in the LAMMPS simulation package [22]. We examine the fracture toughness and critical strain of ice as the number of carbon dioxide molecules at the crack tip is increased.

We use the ReaxFF reactive force field to model the water and carbon dioxide [21, 23]. ReaxFF force field is a bond-order potential that uses a general relationship between bond distance and bond order and between bond order and bond energy. This force field thereby has the advantage to model the chemical system where there might be reactions involving the breaking and reforming of the bonds while it is significantly faster than quantum chemical calculations. The total energy of the system is given by

$$\begin{eqnarray} &&\fl E=\sum {E_{{\rm bond}}} +\sum {E_{{\rm over}}} +\sum {E_{{\rm under}}} +\sum {E_{{\rm val}}} +\sum {E_{{\rm pen}}}\nonumber\\ &&+\,\sum {E_{{\rm tor}}} +\sum {E_{{\rm conj}}} +\sum {E_{{\rm vdW}}} +E_{{\rm Coul}}, \end{eqnarray} \tag{ 1 }$$

where each term is defined as the bond energy, atom under-coordination and over-coordination penalty term, valence angle term, penalty energy term for two double bonds with a sharing atom, torsion angle, conjugation effect term, van der Waals interaction and Coulomb interaction. The detailed expression of each energy term and the parameters of the model are included in the original paper [23]. Together with the parameters of oxygen [21], we can model all interactions in our system.

We use an atomistic model of hexagonal ice (ice Ih) as shown in figure 1. The initial overall geometry of the model is $44\,{\mathring{\rm A}}\times 44\,{\mathring{\rm A}}\times 44\,{\mathring{\rm A}}$, the distance between neighbouring oxygen atoms is 2.7 Å (which expend to 2.75 Å in equilibrium), the H–O distance in a water molecule is 0.96 Å (which shrinks to 0.95 Å in equilibrium) and the H–O–H angle is 104.5° (which approaches 106.7 ± 2.6° in equilibrium). By taking this initial geometry of the bulk ice, the system contains 2880 water molecules with a density of ∼1 g cm⁻³ (which decrease to ∼0.94 g cm⁻³ in equilibrium). A wedge-shaped crack with sharp ends is made by removing the water molecules at the middle of the bulk ice through the z direction as depicted in the figure. The carbon dioxide molecules are randomly seeded within the crack to mimic the natural structure of an air bubble. The C–O distance is initially chosen as 1 Å and expands to 1.25 Å in equilibrium, and the O = C = O angle is 180°.

Zoom In Zoom Out Reset image size

The concentration of the carbon dioxide is given by the mole fraction of the carbon dioxide as $C={n_{{\rm CO}_2}}/{n_{{\rm tot}}}$, where $n_{{\rm CO}_2}$ is the number of carbon dioxide molecules and n_tot is the total number of molecules in the entire system. Periodic boundary conditions are applied in all directions of the simulation box. The temperature and pressure is controlled by Nose-Hoover thermostat and barostat, respectively. The ice and carbon dioxide system are fully equilibrated in a constant temperature (250 K) and isotropic pressure of one atmosphere (NPT ensemble for 100 000 steps) before loading. We increase the loading in a quasi-static way as we minimize the energy and equilibrate the material after each strain increment at constant temperature. The boundary condition is defined by the strain–time relationship and is composed of several increments as shown in figure 1. For each strain increment, we firstly apply 5% additional strain in 4000 steps by uniformly deforming the periodic length in the y direction with a constant strain rate. During this stage the temperature is controlled to be constant at 250 K and the pressures in x and z directions are controlled independently to be one atmosphere. Energy minimizations are carried out using the conjugate gradient algorithm until the energy or the force reaches convergence. During equilibration we apply a constant temperature and fixed boundary condition (NVT ensemble) and equilibrate the system for 100 000 steps in 250 K, which is found to be sufficient for convergence of the stress in the system. We compute the average stress during this stage and record the stress–strain curve. The temperature and pressure parameters are chosen to match estimated climate conditions in the polar area. The time step through the simulation is set to be 0.25 fs.

Metadynamics simulations are performed using LAMMPS and the PLUMED package [22, 24]. Energy deposition events occur every 100 steps in the form of a Gaussian with height of 0.01 kcal mol⁻¹ and a width of 0.35 Å. The resulting free energy landscape is obtained by summing up all Gaussians deposited. We run the simulation for 8000 000 steps until there are no further changes in the shape and depth of the free energy landscape.

A hexagonal ice (ice Ih) model is used for this study as it is the most common form of natural snow and ice on Earth [1]. The basal plane is composed of hexagonal plates, each of which is formed by six water molecules at the vertexes connecting by hydrogen bonds. For natural growth of ice on top of water, the basal plane is parallel to the surface and accounts for the radius growth of ice and its growth in the thickness is in the direction along the axis of basal plane [25]. The initial crack in our model penetrates the ice plate along the axis of the basal plane as shown in figure 1. We apply uniaxial tensile strain on the (1 1 0) surface which is perpendicular to the basal plane and the long axis of the initial crack. This computational set up allows us to investigate the dynamics behaviour of ice under tensile (mode I) loading. The ultimate stress leading to crack propagation in the ice can be caused by thermal expansions or bending by crushing in contact with other materials.

We observe an interesting behaviour of carbon dioxide molecules, which are seen to move towards the boundary and the tips of the crack. As depicted in the series of snapshots focused on one representative carbon dioxide molecule in figures 2(a)–(d), the molecule attaches to the crack surface due to attractive Van der Waals interactions. As it gets closer to the boundary, the molecule forms hydrogen bonds with nearby water molecules of the ice crystal. The hydrogen bonding energy is measured to be 5.5 kcal mol⁻¹ with an O–H distance of 2.0–2.2 Å, while it is a little smaller than the hydrogen bonding energy between two water molecules (measured as 5.8 kcal mol⁻¹). Once the loading is applied, the molecule adjusts its position and moves towards the crack tip, a high-stress region, through a flipping motion. It is found that most of the carbon dioxide molecules have the tendency to move towards the crack tip as is shown by the Markov Chain of each carbon dioxide location shown in figure 2(e). It is shown by the chains that the movements of the carbon dioxide molecules are consistently along the crack surface. We also calculate the relative density distribution along the long axis of the crack as shown in figure 2(f), and it is seen that the local carbon dioxide density around the crack tip can reach twice as much as that of the centre region under stress, supporting the generality of the above described mechanism.

Zoom In Zoom Out Reset image size

We now investigate how the interactions between carbon dioxide and water molecules at the crack tip affect the mechanical properties of the ice crystal. As shown in figure 3(a), the presence of the carbon dioxide molecule significantly alters the strain–stress relationship of cracked ice crystals. In particular, we find that it makes ice more brittle by decreasing the material strength as obtained by peak forces in figure 3(a). Using linear elastic fracture mechanics, the fracture toughness of the ice can be estimated as [26]

$$\begin{equation} K_{{\rm Ic}} =\sigma _c \sqrt {W\tan \left( {\frac{\pi a}{W}} \right)}, \end{equation} \tag{ 2 }$$

where σ_c is critical stress state at which the crack begins to grow, a = 10 Å is the half crack length and W = 44 Å is the system length in x direction. We find that higher concentrations of carbon dioxide lower the toughness of the material, as directly shown in figure 3(b). The fracture toughness of pure ice is 19.4 kPa m^1/2 without adding carbon dioxide (C = 0), while it decreases to 12.0 kPa m^1/2 with C = 2% carbon dioxide. This pure ice result agrees with the surface energy obtained as 0.065 ± 0.003 J m⁻² for ice crystals as measure by experiments [27] (the fracture toughness is 24.5 ± 5.3 kPa m^1/2). The presence of carbon dioxide molecules leads to a 38% reduction in fracture toughness. Our simulation suggests that increasing C from zero to larger values continuously decreases the fracture toughness of ice, making the material more fragile. Based on the calculation result of the fracture toughness and the observation of the migration pattern of carbon dioxide molecules, we propose a model to summarize the mechanism that carbon dioxide exhibits and thereby enhances ice fragility as shown in the schematic figures plotted in figure 3(c). The carbon dioxide first moves closer to the crack tip in a flipping motion and then attempts to form hydrogen bonds by attacking water molecules under stress. This process repeats itself many times as shown, leaving broken hydrogen bonds behind and thereby effectively increasing the brittleness of the ice. Moreover, the continually increasing crack width prevents the ice from healing itself, making the bond breaking process irreversible. Similar chemically assisted fracture processes have been found in other crystalline materials [28]. For example, water interacts with mica and makes it brittle [29].

Zoom In Zoom Out Reset image size

Based on these results, we now ask the question how to quantitatively relate the fracture property of ice to the concentration of carbon dioxide. We use a bulk single crystal ice to obtain the basic mechanical properties by applying tensile stress on the (1 1 0) surface. The Young's modulus is measured as E = 4.2 GPa and the Poisson's ratio is ν = 0.3. According to linear elastic fracture mechanics, the fracture surface energy is $\gamma={K_{{\rm Ic}}^2(1-\nu^2)}/{2E}$. Therefore the energy required to break each hydrogen bond can be calculated from the fracture toughness as [30]

$$\begin{equation} E'=2\gamma A=\frac{AK_{{\rm Ic}}^2 (1-\nu^2)}{E}, \end{equation} \tag{ 3 }$$

where A = 15.6 Å² is the average surface area of a single hydrogen bond within the (1 1 0) plane.

We establish the connection between the carbon dioxide concentration and the material property of the ice via a simple statistical model based on Bell's equation [31]. The dissociation rate of hydrogen bonds at the crack tip without carbon dioxide can be estimated via [32, 33]

$$\begin{equation} \omega =\omega _0 \exp \left( {\frac{E_0 -Q_{{\rm HB}}}{K_{\rm B} T}} \right), \end{equation} \tag{ 4 }$$

where ω₀ is the natural frequency of bond vibration, E₀ is the external work before bond rupture without adding carbon dioxide and Q_HB is the hydrogen bonding energy. Adding the carbon dioxide around the water molecule increases the dissociation probability as

$$\begin{equation} \omega =\omega _0 (1+nC)\exp \left( {\frac{E'-Q_{{\rm HB}}}{K_{\rm B} T}} \right), \end{equation} \tag{ 5 }$$

where n is a constant factor, E' the biased external work before bond rupture with carbon dioxide. The similar expression has been used to analyse experimental results [34]. Combining equations (4) and (5) we obtain a E' ∼ C relationship as

$$\begin{equation} E'=E_0 -K_{\rm B} T\ln (1+nC). \end{equation} \tag{ 6 }$$

For each concentration C, we measure the fracture toughness in the simulation and obtain E' as given by equation (3). Therefore, the constant E₀ and C can be obtained by fitting the equation to the series of simulation results as shown in figure 4(a). Under the condition of T = 250 K, we obtain E₀ = 1.9 kcal mol⁻¹ and n = 548.5 according to the simulation results. Here the result given by E' in equation (6) is the most important result of this paper because we can use it to calculate the surface energy of ice under different conditions of C, and thereby larger scale simulations can be developed based on this result [35]. In the polar bulk ice, the average concentration of carbon dioxide is C ≈ 0.02% [17, 18], leading to E' = 1.85 kcal mol⁻¹, which means the surface energy decreases by ∼3% of its original value. This concentration value can reach C ≈ 0.1% for deep bulk ice [18], leading to E' = 1.68 kcal mol⁻¹, which means that the surface energy decreases by ∼12% compared with its original value. However, in the (naturally relevant) case the fracture toughness is affected by the local concentration of carbon dioxide around the crack tip area which is a high-stress zone. The migration of carbon dioxide as demonstrated in our study makes it clear that although the carbon dioxide has an overall low concentration of 0.02%, its local concentration around the crack tip can be significantly higher and may actually accumulate with crack propagation. Therefore, its actual effect can be larger than this estimate. For example, for C = 2%, we obtain E' = 0.66 kcal mol⁻¹, which means the surface energy decreases by ∼65% relative to its original value. We also repeat the set of simulations by only altering the temperature to T = 10 K and find good agreements between the simulation results for E' given by equation (3) and the theoretical results given by equation (6), as shown in figure 4(a). It is noted that E₀ is smaller than the hydrogen bonding energy (∼ 5.8 kcal mol⁻¹) than we measure for a single isolated bond because there is no confinement for the water molecule at the crack tip to break along a certain pathway. The hydrogen bonds can break and reform with other water molecules on the same fracture surface at the same time, decreasing the energy barrier for crack propagation. To make this argument more clear, we set up a simple simulation system with two water molecules and a single carbon dioxide molecule as illustrated by the inserted pictures in figure 4(b). The initial coordinates of the two water molecules and the carbon dioxide are obtained from the equilibrated ice crack tip structure with a C(in CO₂)–O(in H₂O) distance of 2.4 Å. We fix the coordinate of one of the water molecules and the carbon dioxide and use metadynamics to measure the free energy landscape of the system as a function of the displacement of the oxygen atom in the tensile and perpendicular directions. We find that the carbon dioxide close to the water molecule alters the hydrogen bond energy landscape by switching the conformation with minimum energy from the origin to the right region as shown in figure 4(b). For the case with no carbon dioxide the free energy of the conformation at the origin is measured to be 5.8 kcal mol⁻¹, while for the case with carbon dioxide the free energy valley is measured to be more shallow (4.2 kcal mol⁻¹ as shown), which makes the hydrogen bond more prone to rupture.

Zoom In Zoom Out Reset image size

The result given by equation (6) clearly shows that the effect of carbon dioxide on ice fracture is more than creating vacancies in the ice crystal. The equation gives an analytical relationship between the surface energy of ice and the combination of temperature and the concentration of carbon dioxide. For the non-carbon dioxide case, the estimated fracture toughness of pure ice (19.4 kPa m^1/2) agrees well with experimental measurements [27] (24.5 ± 5.3 kPa m^1/2). An increasing temperature and carbon dioxide concentration leads to a decreasing surface energy, which qualitatively agrees with the observation of the cause of deglaciation [13, 14]. Considering the difficulty of directly measuring the effect of carbon dioxide on ice fracture (in the existing literature), our study may provide a useful example to predict or analyse the ice fracture in different environments. We also hope that our work may inspire future, more detailed studies of this issue from an experimental perspective.

Using atomic simulations with the ReaxFF reactive force field, we investigated the effect of carbon dioxide on the fracture behaviour of ice. We find that the material strength and fracture toughness are decreased significantly under increasing concentrations of carbon dioxide molecules. This phenomenon is caused by the interaction between water and carbon dioxide molecules. Carbon dioxide molecules first adhere to the crack boundary by forming hydrogen bonds, and then migrate along the crack boundary towards the crack tip. The dissociation energy of hydrogen bonds at the crack tip is decreased under the attack by carbon dioxide. This migration-attacking process repeats itself and renders the ice crystal more brittle by mediating crack propagation. Our theoretical model quantitatively accounts for the effect of carbon dioxide on the surface energy and fracture toughness of ice. It could be instrumental for further studies of ice fracture in various chemical environments and may be scaled up by incorporating it into models of glacier dynamics.

This work was supported by AFOSR and DOD-PECASE.