Stochastic Cellular Automata Modeling of CO2 Hydrate Growth and Morphology

Carbon dioxide (CO2) hydrates are important in a diverse range of applications and technologies in the environmental and energy fields. The development of such technologies relies on fundamental understanding, which necessitates not only experimental but also computational studies of the growth behavior of CO2 hydrates and the factors affecting their crystal morphology. As experimental observations show that the morphology of CO2 hydrate particles differs depending on growth conditions, a detailed understanding of the relation between the hydrate structure and growth conditions would be helpful. To this end, this work adopts a modeling approach based on hybrid probabilistic cellular automata to investigate variations in CO2 hydrate crystal morphology during hydrate growth from stagnant liquid water presaturated with CO2. The model, which uses free energy density profiles as inputs, correlates the variations in growth morphology to the system subcooling ΔT, i.e., the temperature deficiency from the triple CO2–hydrate–water equilibrium temperature under a given pressure, and properties of the growing hydrate-water interface, such as surface tension and curvature. The model predicts that when ΔT is large, parabolic needle-like or dendrite crystals emerge from planar fronts that deform and lose stability. In agreement with chemical diffusion-limited growth, the position of such planar fronts versus time follows a power law. In contrast, the tips of the emerging parabolic crystals steadily grow in proportion to time. The modeling framework is computationally fast and produces complex growth morphology phenomena under diffusion-controlled growth from simple, easy-to-implement rules, opening the way for employing it in multiscale modeling of gas hydrates.


■ INTRODUCTION
Gas hydrates are solid structures formed at high pressure and low-temperature conditions from a solution of gas into water. 1,2 Such structures, which are classified as structure I (sI), structure II (sII), and structure H (sH), consist of large and small cages of hydrogen-bonded water molecules that encapsulate the gas molecules. Hydrate-forming gas species include, for instance, xenon (Xe), carbon dioxide (CO 2 ), and hydrocarbon molecules, such as ethane (C 2 H 6 ), propane (C 3 H 8 ), cyclopentane (C 5 H 10 ), and methane (CH 4 ). The type of structure formed depends primarily on the size of the gas molecule. Small gases such as CH 4 and CO 2 form sI hydrate structures, whereas large molecules such as C 3 H 8 form sII hydrate structures. Gas hydrates with even larger molecules like neohexane (C 6 H 14 ) together with a small one like CH 4 form sH hydrates. 3 Naturally occurring gas hydrates exist in sediments along the continental margins, deep lakes, and permafrost regions and are typically sI. 4 Due to their unique properties, including a large gas storage capacity, a large heat of formation and dissociation, and gas selectivity, gas hydrates have several possible industrial applications, 7,8 such as gas capture and storage, 9 refrigeration systems, 10 and water desalination. 11 Because of these applications, basic and applied research on CO 2 hydrate has been a topic of intense investigation. For instance, Zhang et al. 12 reported that an optimal amount of CO 2 liquefaction promotes CO 2 hydrate formation on porous media. Similarly, CO 2 capture micromechanisms and microstructures of stored forms together with hydrate growth rates were investigated by Zhang et al. 13 In their work, based on the large amount of dissolved CO 2 gas in water compared to other investigated gases, a mechanism for hydrate growth was proposed to promote CO 2 capture and storage. Recent advances in CO 2 sequestration and storage via gas hydrates have been recently summarized in a review article by Zheng et al. 14 In particular, various approaches of CO 2 sequestration via gas hydrates and their technical feasibility and potential storage capacity were discussed in that review. These approaches included storage in seawater, sediments under the sea floor, 15 permafrost regions, and CH 4 hydrate reservoirs via CO 2 −CH 4 exchange. However, to further develop hydrate-based technologies, it is imperative to investigate the growth behavior of gas hydrates and the main factors affecting their morphology. 5 Understanding hydrate crystal growth dynamics and morphology is necessary for instance in selecting their efficient formation method. Morphological studies could also deliver insights into the development of either innovative hydrate management alternatives (kinetic inhibitors and antiagglomerates) to mitigate the well-known oil and gas pipeline blockage due to internal hydrate-plug formation, 16 or hydrate promoters. 17 Concerning pipeline blockage, gas hydrate formation occurs as an interfacial phenomenon in which wettability and hydrate film growth and morphology within the pipeline dictate how easily the system may agglomerate due to capillary water bridging and sintering. 18,19 Thus, this work aims to study gas hydrate crystal growth and morphology computationally.
Morphological studies of gas hydrates involve investigating their characteristic appearance, such as the size, shape, and mode of growth. 5 In this respect, several experimental studies based on visual observations have revealed some features of the crystal growth dynamics and morphology of CO 2 hydrates 20−23 and hydrocarbon gas hydrates. 6,24−28 For example, crystal growth dynamics and morphology of hydrates were found to depend on the nature of the gas (i.e., hydrate-former) molecules and the thermodynamic conditions of crystal formation. The subcooling temperature ΔT has been used as the driving force for hydrate crystal growth in such experimental works. ΔT is the difference between the experimental temperature T ex of crystal formation and the triple gas-hydrate-water equilibrium temperature T eq under a given pressure (ΔT = T eq − T ex ). In almost all these experimental studies, the relation between ΔT and the crystal growth morphology of several hydrate structures was revealed. 5 A generic observation was that, as long as the hydrate-former molecules are initially in a gas state (see Figure 1a), 23,24 hydrate crystals preferentially form at the gas−water interface and grow in the form of a thin film covering the interface (see Figure 1b). Furthermore, 23,24 it is observed that hydrate crystals could grow in liquid water from the initial hydrate thin film, if the hydrate-former molecules had dissolved in liquid water to saturation, or a specific concentration lower than saturation, before the hydrate formation (see Figure 1c). Although such a two-stage process typically occurred throughout the entire range of subcooling covered in the experiments, the morphology of the hydrate crystal growing toward water depended very much on the degree of subcooling. 23,24 In particular, the experimental studies demonstrated that the morphology of the hydrate crystals that grew in the water from the hydrate thin film changed from skeletal polyhedral or columnar to dendritic depending on ΔT. 23,24 The mechanism behind this morphology change toward stagnant liquid water presaturated with a hydrateforming substance has been qualitatively rationalized based on the idea that such a variation depends on the growth rate of the hydrate crystal, which in turn is controlled by the diffusion of the hydrate-former molecules dissolved in the bulk liquid-water to the hydrate thin film surface. [23][24][25]29,30 This diffusion was proposed to be driven by the difference in hydrate-former (e.g., CO 2 and CH 4 ) concentration in liquid water adjacent to the hydrate thin film surface (i.e., the concentration of hydrateformer in liquid water in equilibrium with hydrate) and the hydrate-former concentration in the bulk liquid-water phase, given by the concentration in liquid water in equilibrium with gaseous hydrate-former (i.e., the solubility of hydrate-former in water). Such a difference in concentration increases with ΔT.
The CO 2 hydrate formation kinetics and morphology resulting from the interaction between liquid CO 2 and water has also been the subject of experimental investigations due to its potential impact on reducing CO 2 emissions. The motivation behind those studies is the possibility of injecting these emissions as compressed liquid CO 2 into the deep ocean sediments to store as gas hydrates. 20,21,31−34 A typical observation is the occurrence of different stages of hydrate Figure 1. Schematic of gas hydrate growth in a batch gas−liquid water system. 5,6 (a) Hydrate-former gas and water phases in equilibrium and separated by a flat interface. The concentration of hydrate-former in water in equilibrium with gas defines the solubility of the hydrateformer in water. In this panel, the system temperature, T ex , is slightly above the gas-hydrate-water equilibrium temperature, T eq . (b) Growth of an initial hydrate film covering the gas−water interface. (c) After the hydrate film covers the interface, the hydrate grows toward the water phase. Because it is assumed that the hydrate film blocks mass transfer across the interface, the hydrate growth morphology is only dependent on the solubility of hydrate-former in water. T ex is below T eq in (b) and (c). In this work, we are modeling the scenario of panel (c). formation, namely, hydrate nucleation, lateral and thickening hydrate film formation at the interface between liquid CO 2 and water, and hydrate growth. The CO 2 hydrate crystallization and growth from an aqueous CO 2 solution at pressures to 6.5 MPa have also been explored experimentally. 35 Most empirical morphological studies have been performed in batch systems, where the water is stagnant and presaturated with the hydrate-former substance (see Figure 1). However, a few studies have explored the growth morphology of gas hydrate crystals under flow conditions. 36−38 The crystal growth morphology of gas hydrates formed with multihydrate former gas mixtures has been reported. 26,39−42 The morphology of methane hydrate formation in porous media is also documented. 43 The morphological study of methane hydrate formation and dissociation in the presence of amino acids as kinetic promoters is also available in the literature. 44 Significant morphological differences using amino acid compared to the commonly used kinetic promoter, sodium dodecyl sulfate (SDS) surfactant, were found at the same concentrations under similar experimental conditions. The effect of salts, such as NaCl, on hydrate formation morphology has been studied as well. 45,46 The crystal growth morphology of structure-H hydrate has been experimentally explored. 28 More recent studies have reported the crystal growth morphology of gas hydrates with ozone, 47 CH 4 + tetrahydropyran (THP) hydrates, 48 and CO 2 hydrate in aqueous fructose solution. 22 Although the experimental studies of hydrate crystal growth morphology have achieved considerable progress during the last few decades, 5 the physical mechanisms leading to the morphological transition of gas hydrate crystals growing from a hydrate thin film into liquid water as a function of thermodynamic conditions remain unclear. The lack of quantitative information from current visual experimental studies could be the reason for such a lack of understanding. This is where computational modeling could provide significant insight. However, such computational modeling studies of gas hydrate growth morphology are scarce in the literature.
Molecular dynamics (MD) simulations are the main computational tools for understanding the fundamental molecular mechanisms behind gas hydrate growth. 49−53 However, MD simulations are limited to short simulation times and small domains and, therefore, cannot be used to perform large-scale computational studies of gas hydrate morphology. Yet, MD can be coupled with other mesoscopic theoretical frameworks, like the phase-field (PF) method, to model large-scale properties of gas hydrates. 54−59 This approach has been adopted by Tegze et al. 56 to predict the nucleation and growth rates of solid CO 2 hydrate in CO 2 aqueous solutions under conditions typical to oceanic natural gas hydrate reservoirs. Although PF theory is one of the most powerful methods that are used to describe solidification phenomena in undercooled liquids, 60 simulations based on it are generally very computationally demanding as well.
Motivated by those computational challenges, Buanes et al. 61,62 developed a simplified framework for simulating largescale hydrate kinetics based on a Monte Carlo (MC) cellular automaton (CA) approach to evaluate the most probable growth path. They also applied it to hydrate growth from a supersaturated aqueous solution of CO 2 under conditions typical to oceanic natural gas hydrate reservoirs. They concluded that the MC-CA approach has the benefit of being much more computationally efficient and is still giving results consistent with the PF theory. Although Buanes et al. 61,62 did not conduct a detailed sensitivity analysis of CO 2 hydrate growth morphology, the novelty of their work was combining the two numerical approaches to study gas hydrate kinetics.
The MC and CA modeling approaches are well-established mesoscopic procedures to simulate natural phenomena. In particular, CA was introduced by John von Neumann and Stan Ulam as a technique to model biological self-reproduction. 63 Since then, this approach has attracted significant attention in computer science, biology, and other fields, and it is, in fact, a well-established mesoscopic modeling approach in materials science. 64−66 A CA is a dynamic system composed of a lattice of cells. Each cell in the lattice can assume a state from a given set of states. The system evolves in time according to an updating rule which is a function of the cell state itself and the state of neighboring cells. All the cells are updated synchronously, and the updating rule can be either deterministic or probabilistic. A CA with probabilistic cell state change rules is called a probabilistic CA. A cellular automaton can also be combined with other computational methods to get what is known as hybrid CA modeling. The appealing features of the CA approach are its computational efficiency and its ability to capture complex real-world behavior from simple rules. In this respect, CA provides insights into the complex behavior of gas hydrate growth morphology and could even be fast enough to be implemented into flow modeling tools relevant to gas hydrate formation and flow assurance. 67 In this work, we build on and extend the MC-CA-type modeling approach by Buanes et al. 61,62 In particular, we implement a hybrid probabilistic CA modeling framework to explore qualitative and generic features of sI CO 2 hydrate crystal growth morphology toward liquid water presaturated with CO 2 . Thus, we focus on a batch system where the morphological hydrate growth proceeds toward a stagnant water phase from a pre-existing hydrate film. We assess computationally the proposed mechanism by which the hydrate crystal growth morphology is caused by the diffusion of CO 2 molecules to the hydrate film due to a CO 2 concentration gradient in the aqueous solution. To this end, we correlate the growth morphological change to the magnitude of ΔT and properties of the growing hydratewater interface, such as surface tension and curvature. Our findings could have implications for CO 2 storage technologies, since CO 2 hydrates have been suggested to be a stable storage form for this well-known greenhouse gas. 68−70 However, the fate of the stored CO 2 would strongly depend on the CO 2 hydrate formation pathway. Therefore, it is imperative to investigate the morphology of CO 2 hydrate crystals grown in well-defined systems and to know what factors control these morphological phenomena before moving toward more realistic hydrate conditions in the geological storage of CO 2 , where a continuous supply of CO 2 is expected.
The rest of the paper is organized as follows. In the next section, we present the hybrid probabilistic CA modeling approach and provide details about its numerical implementation. The two last sections are devoted to simulation results, discussion, and conclusions. that this hydrate film limits mass transport of CO 2 across the CO 2 -aqueous solution interface, 71 the hydrate formation and growth will happen only due to the CO 2 dissolved in water (see Figure 1). Furthermore, CO 2 hydrate formation will occur at the hydrate-water interface, and hydrate formation in the bulk of the water solution will be neglected. The CO 2 hydrate crystallization front morphology is explored as a function of the subcooling ΔT and hydrate-water interface properties.
As Figure 2 shows, the simulation domain of interest pertains to a sI CO 2 hydrate crystal region led by a phase boundary (i.e., hydrate-water interface) that propagates in time toward the water with dissolved CO 2 . This interface is the region of volume where the transition between the solid hydrate state and liquid state occurs and can be considered as an integration (solidification) layer that consists of a number of water cages under formation. The modeling approach's physical processes are (1) the thermodynamics of phase transformation between the two phases at the moving hydratewater interface, (2) the interfacial attachment (or integration) kinetics of CO 2 molecules at that interface, and (3) the diffusion of CO 2 molecules inside the stagnant liquid water region. 72 These processes are the centerpieces of this work, as we aim to simulate hydrate crystal growth and assess how the pertinent mechanisms give rise to morphological changes. The following sections describe how these processes are incorporated into the hybrid probabilistic CA framework implemented in this work. Following Ohmura et al., 23,24 the model does not consider the interface temperature increment due to heat release from the hydrate formation. This assumption is justified by the observation that the thermal transport of liquid-water and hydrate crystal is several orders of magnitude faster than the diffusion of CO 2 in liquid water. 23,24,73 Furthermore, CO 2 diffusion along the hydrate crystal is several orders of magnitude slower than in the liquid water phase and, therefore, will also be neglected. 23,24,71 With these aspects in mind, let us describe in the following section the computational simulation approach of this work.
Hybrid Probabilistic Cellular Automaton Framework. We introduce in this section the hybrid probabilistic CA approach implemented in this work. We will define the CA lattice of cells associated with the simulation domain presented in Figure 2, the possible states of the cells, and the probabilistic updating rule of those states. The latter will be based on the MC method. 61, 62 The interaction between cells due to mass transport by diffusion is performed using a MC implementation of Fick's law and the incorporation of CO 2 to the integration layer defined by the interface is modeled by Fick's law alone. 61,66 These implementations of diffusion and incorporation are combined with the CA probabilistic updating rule, giving rise to this work's hybrid probabilistic CA approach.
Cellular Automaton Lattice, State, and Composition. As Figure 3 shows, the CA lattice associated with our simulation domain is constructed by dividing it into cubic cells or voxels of equal and fixed volume, V = l 3 , where l will be of the same order as the hydrate-water interface thickness. 61, 62 The voxels of the hydrate crystal are portrayed with diagonal hash patterns, and the ones belonging to the hydrate-water interface have their borders highlighted by red dashed lines. The rest of the voxels are presented with border-lines in black. The bottom right of the figure shows a 2D representation of the voxel's neighborhood implemented in this work. It contains the so-called von Neumann neighborhood, in which a voxel is only connected to its first nearest neighbor (nn) voxels, and the Moore neighborhood, where the first and second nearest neighbors are included. 65 As Figure 3 also shows, there are no domain boundaries along the y direction because we assume periodic boundary conditions. However, closed boundary conditions are considered along the x direction. The latter defines two boundaries; the first involves voxels on the hydrate-water interface. The voxels at the end of the domain along the x axes provide the second boundary. The requirement for a voxel to be in the hydrate-water interface is having at least a hydrate voxel nearest neighbor. Note, however, that this hydrate-water interface boundary is moving toward the water region because, its voxels are allowed to perform phase transformations between a liquid-like state (water with dissolved CO 2 ) and a solid (or CO 2 hydrate) state.
The voxels are not only characterized by their neighborhood but also their state, denoted here by the Greek letter ϕ. For liquid water voxels with dissolved CO 2 (including the interface ones), ϕ takes the value of 0, and for a hydrate (i.e., solid) voxel, it takes the value of 1. The voxels are also characterized by their composition given by a CO 2 molar fraction of χ c = n c / Figure 2. Schematic representation of the simulated domain, showing a sI CO 2 hydrate crystal region led by a hydrate-water interface that propagates toward the stagnant water with dissolved CO 2 . This hydrate-water interface has a specific thickness. The modeling approach's physical processes are (1) the thermodynamics of phase transformation between liquid state (water with dissolved CO 2 ) and solid state (CO 2 hydrate), (2) the interface attachment (or integration) kinetics of CO 2 molecules, and (3) the mass transport by diffusion of CO 2 molecules inside the stagnant liquid water phase.
(n c + n w ) and a water molar fraction of χ w = 1 − χ c , with n c and n w been the moles of CO 2 and water, respectively, inside the voxel.
With the CA lattice, neighborhood, as well as voxel state and composition defined, we are ready to introduce the evolution rules of the hybrid probabilistic CA simulation framework, which are based on the three physical processes mentioned above, namely, phase transformation, attachment kinetics, and mass transport by diffusion. From a computational viewpoint, starting at a time, t, equal 0, these rules are applied at each iteration (after time Δt), to update the states and compositions of the voxels. The pertinent iterative process is explained in detail later, together with the definition of Δt. Let us first explain each of the rules.
Hydrate Phase Transformation. Let us describe the protocol for the transition from liquid water with dissolved CO 2 state and a CO 2 solid hydrate state, which pertains to the interfacial voxels (i.e., the voxels having at least one hydrate voxel as their neighbor, depicted in red dashed lines in Figure  3). The remaining liquid water voxels do not perform a phase transformation because hydrate formation in the bulk liquid water is not considered. In other words, nucleation is neglected in the current approach. At each time step, Δt, all interfacial voxels are allowed to attempt a phase transformation between liquid and solid hydrate. The basis of such a phase transformation is the Glauber transition probability for an interfacial liquid voxel to become a solid one, given as 74 (1) and the one for the remelting, where W h + W l = 1 (normalization), W h /W l = e −βΔF (detailed balance), and β = (k B T) −1 , with k B being the Boltzmann constant and T the temperature of the system, which is assumed to be constant during the simulation. 74 Thus, if a uniform random number (uniform deviate) r h is less than W h solidification takes place; however, if another uniform deviate r l is less than W l the solidified voxel returns to its original liquid state.
Following Buanes et al., 61,62 we assume that ΔF = Φ(ϕ, λ)Δf(χ c , T), where Δf(χ c , T) is the change in free energy when an interfacial voxel with CO 2 molar fraction χ c and temperature T changes its phase from liquid to solid, and Φ(ϕ, λ) is a term accounting for the solid−liquid interface effects, like for example surface tension and interfacial curvature. Thus, to obtain Δf(χ c , T) = f h − f l , we need expressions for the bulk free energy densities of the liquid water with dissolved CO 2 (f l ) and solid CO 2 hydrate (f h ) phases as a function of χ c and T. This work will use the free energy density expressions implemented in the literature [54][55][56][57]76 to obtain the equilibrium isobaric water-CO 2 phase diagram describing the zones of CO 2 hydrate stability in terms of temperature and composition. The free energy density for the solid hydrate state, f h , was constructed assuming a sI CO 2 hydrate structure with large cages fully occupied by CO 2 molecules and small ones empty. For further details regarding these free energy densities and the equilibrium phase diagram, refer to the Appendix.
The term Φ(ϕ, λ) defined as 61,62  (4) where the summation index k is over the first and second nearest neighbor voxels of the interfacial voxel in question (Moore neighborhood), ϕ k is their state (1 for solid hydrate and 0 for liquid), and w k is a weighting term equal to 2 for first nearest neighbors and 1 for second nearest neighbors. The value of 6 corresponds to a summation on a flat interface (i.e., κ = 0) and is taken equal to one-half the sum over the entire weighted neighborhood with all voxels in a solid hydrate state. We refer to the Appendix for an explanation of how this value is obtained.
Overall, ΔF (as it appears in the Glauber transition probabilities) considers the tendency to favor one phase or another based on bulk free energy minimization and the tendency to reduce the surface energy by reducing surface  Figure 2). The CA lattice consists of cubic voxels of volume, V = l 3 , containing CO 2 and water. There are periodic boundaries along the y axes and closed boundary conditions along the x direction. The interface is a moving boundary because their voxels, here represented in red dashed lines, can perform phase transformations between a liquid-like state (water with dissolved CO 2 ) and a solid (or CO 2 hydrate) state. The attachment kinetics and the diffusion mass transport between voxels are also highlighted. The bottom right of the figure shows a 2D representation of the voxel's neighborhood. The side length of the whole simulation domain in the x direction will be denoted by L x and the one for the y direction as L y . These two lengths are a multiple of the voxel side length l, whose size is chosen to be of the order of the interface thickness. Although not shown in the figure, the side length in the z direction is L z = l. The voxels are characterized by their composition, given by a CO 2 molar fraction of χ c and a water molar fraction of χ w = 1 − χ c , as defined in the main text. curvature (or surface area). In Figure 4, we present a schematic example in which, although the two liquid-to-solid transitions have an equal bulk free energy change, their Glauber transition probabilities are different because one leads to a flat interface and the other to a less favorable curved one. When λ = 0 (i.e., surface tension is ignored), Φ(ϕ, λ = 0) = 1, and the change in bulk free energy is the phase transformation determining factor. Thus, λ is a free parameter that is changed to explore the dynamics of this model.
Mass Transport. The mass transport mechanism considered in this work is diffusion. However, because the rate of CO 2 diffusion through the hydrate-water interface decreases rapidly relative to the diffusion in the liquid water region, and it is several orders of magnitude smaller in the solid hydrate structure, 78 mass transport via diffusion is considered only between voxels belonging to the water region. In other words, diffusion phenomena along the interfacial and solid hydrate regions and between such two regions are neglected.
Modeling diffusion is achieved using a probabilistic implementation of Fick's first law. 61, 62 We assume that voxelto-voxel diffusion occurs only between the first nearest neighbors voxels (von Neumann neighborhood). The voxels, however, can differ in composition. Thus, at each time step, Δt, one of the first nearest neighbors, i nn , is chosen randomly for each voxel, i, where i runs from 1 to the total number of water voxels in the system. Then, the net flux of CO 2 molecules from voxel i to voxel i nn is given as In this equation, c c i = n c i /V and c c i nn = n c i nn /V are the CO 2 concentrations in the voxels expressed in terms of the corresponding numbers of moles n c i and n c i nn and the voxel volume V = l 3 . The term A = l 2 is the area through the flux of CO 2 molecules pass during a time Δt, D c is the diffusion coefficient of CO 2 in water, l is the side length of the voxels, and Δn c i = n c i (t + Δt) − n c i (t). ϵ represents a random number from a Gaussian distribution with mean equal zero and standard deviation σ ≪ 1. If instead, c c i < c i nn , the flux is performed from voxel i nn to voxel i according to where Δn c i nn = n c where Δf l , which is obtained from the free energy density of the liquid water with dissolved CO 2 , f l (see Appendix), is the change of free energy of the system due to the diffusion process and β = (k B T) −1 . Figure 5a presents a schematic example of two voxels connected by a common area A through which the flux of CO 2 molecules and water occurs. The compositions of the voxels involved in the diffusion process are updated according to eqs 5 or 6, ensuring that χ c + χ w = 1 in each voxel.
Interfacial Attachment Kinetics. We assume a concentration-dependent driving force for the CO 2 uptake (i.e., adsorption) at the solid hydrate-liquid water interface (i.e., Figure 4. A liquid-to-solid transition as performed by the interfacial voxel in red dashed lines under two different neighborhood conditions. Because in both cases the voxel has the same initial CO 2 molar fraction, χ c , the value of Δf(χ c , T) is the same in both cases (as an example, we assume that Δf(χ c , T) > 0). However, Φ(ϕ, λ) is different in both cases because case (1) leads to a local interface with positive curvature, whereas case (2) leads to a flat interface. Φ(ϕ, λ) is calculated with eqs 3 and 4 using the Moore neighborhood represented in blue dotted lines. It is straightforward to verify that the Glauber transition probability for solidification (W h ) of case (1) is smaller than the one for case (2). If Δf(χ c , T) < 0, in order to ensure a similar outcome, Φ(ϕ, λ) = 1 + λκ has to be replaced by Φ(ϕ, λ) = 1 − λκ. Figure 5. Schematic of the diffusion process for CO 2 dissolved in voxels filled of water. Both voxels have the same constant volume V = l 3 , and the flux is through the area A = l 2 connecting both voxels. The side length of the voxels is given by l, where, as mentioned before, l is the interface thickness. The direction of the flux depends on the concentration difference between voxels and is accepted, provided the system's total energy is minimized. (b) Schematic of the interfacial attachment process (i.e., uptake of CO 2 at the interface). This is similar to (a) but with j being an interfacial voxel (represented by the red dashed line) with a neighbor j nn in the liquid phase. In this case, the net CO 2 flux is only in one direction, namely, from voxel j nn in the water region to the interfacial voxel, j.
integration layer). 23,24,79−88 Because the hydrate film is present at the start of the simulations, this driving force is assumed to be the difference between CO 2 concentration in liquid water near the interface and the CO 2 concentration in bulk liquid water. Thus, we consider a net molar flux of CO 2 molecules toward an interfacial voxel j as given by where n c j is the number of CO 2 moles in voxel j. As before, A = l 2 is the area through which the flux of CO 2 molecules pass during a small time interval Δt, D c is the diffusion coefficient of CO 2 in water, and l is the side length of the voxels. c c j nn = n c j nn /V denotes the CO 2 concentration of voxel j nn in the water region, with that voxel being a first nearest neighbor of interfacial voxel j (von Neumann neighborhood). Note that j runs from 1 to the total number of interfacial voxels. c c lh = n c lh /V is the concentration of CO 2 in liquid water near the interface that is in equilibrium with CO 2 hydrates as expressed in terms of the voxel volume and number of moles, In this work, the molar fraction χ c lh will be assumed to be given by the liquidus line of the equilibrium phase diagram presented in the Appendix. The implementation of eq 8 is justified from the steady state mass balance resulting from the equality between the net integration (or so-called enclathration) rate of CO 2 molecules into the partially formed cages at the interface, and the net molar flux toward that interface. 51, 52,87 For all ΔT values explored in this work, c c j nn ≥ c c lh , and the net flux of CO 2 molecule occurs from voxel j nn to the interfacial voxel j (see Figure 5b). In our simulations, the maximum value that c c j nn can have is that of the CO 2 -in-water solubility. Note nn . Because the volume is conserved, the CO 2 flux is compensated by a water flux in the opposite direction. The compositions of the voxels involved in the attachment process are updated according to eq 8, ensuring that χ c + χ w = 1 in each voxel.

■ NUMERICAL IMPLEMENTATION
The dimensionless parameter, Figure 6. Flowchart for the main steps involved in the hybrid probabilistic CA simulation. At each step, Δt, the attachment, diffusion, and phase transformation processes proceed sequentially.   (9) from eqs 5, 6, and 8 defines the length and time scales of the system. We will use D c = 10 −9 m 2 /s. 62,89 By choosing the iteration time as Δt = 1/9 × 10 −9 s and the side length of the voxels as l = 10 −9 m, the dimensionless diffusion coefficient is fixed at D c ′ = 1/9. 61, 62 The simulation starts at t = 0 s with a given domain configuration and composition of the associated voxels. Then within each iteration time Δt, the phase transformation, diffusion, and attachment processes are implemented to update the state and composition of the voxels. A flowchart of the steps involved in the simulation is shown in Figure 6.
A simulation domain of total volume V = Nl 3 is defined, where N is the total number of voxels in the domain. As presented in Figure 3, these voxels have an equal volume of V = l 3 which remains constant during our simulations. Initially, the whole domain is divided into three regions. The first one represents the CO 2 hydrate thin film and has a volume of V h = N h l 3 , where N h is the number of solid hydrate voxels. The second region corresponds to the water with dissolved CO 2 voxels having at least one nearest neighbor solid hydrate voxel (i.e., hydrate-water interface) and has a total volume of V in = N in l 3 , where N in is the total number of interfacial voxels. The third region is the one containing the rest of the water with dissolved CO 2 voxels and has a total volume of V w = N w l 3 , in such a way that Although N is constant, N h , N in , and N w evolve in time due to the hydrate phase transformation of the interfacial voxels. However, initially we assume that N in < N h < N w (note that in Figure 3, N in = N h < N w ).
After setting up the simulation domain and initializing the composition of the voxels, time is advanced by t + Δt. Then, the interfacial voxels are identified, and the interfacial attachment process is performed for each interfacial voxel. The mass transport between the water region voxels proceeds by selecting one first nearest neighbor for each voxel and realizing the diffusion of CO 2 and water molecules. The current simulation step ends with the solidification attempt of all interfacial voxels. After that, the time step increases again by Δt. As specified in Figure 6, the procedure is repeated until a final stipulated time is reached.

■ SIMULATION RESULTS
The hybrid probabilistic CA modeling approach is developed to model the crystal growth morphology of gas hydrate during its formation in stagnant water with a dissolved gas hydrateforming substance. As a case study, we consider the variations in the morphology of sI CO 2 hydrate. The primary intent is to correlate those variations with the degree of subcooling, ΔT, and the local interface curvature, quantified by the dimensionless parameter λ and the variable term κ.
The simulations start from a preformed CO 2 hydrate thin film exposed to water with dissolved CO 2 . We will assume that for all ΔT values explored in this work, the initial CO 2 molar fraction in all voxels of the interface and water regions is around the constant value χ c o > χ c lh , where, χ c lh will be a decreasing function of ΔT (see the liquidus line of the isobaric phase diagram in the appendix section). In the following, all the results will be expressed in molar fractions, and in all figures the time will be presented in units of microseconds (μs). Following Tegze et al., 56 we will perform simulations for conditions typical for seabed reservoirs. Thus, the pressure will be P = 6.2 MPa, and the triple equilibrium temperature will be assumed to be T eq = 282.5 K (see Figure 15 in the appendix section).
With the length and time scales of the model defined implicitly by the parameter D c ′, the only parameters left for varying are the subcooling, ΔT, and the strength of the local curvature effect, λ. Increasing ΔT leads to an increment of the CO 2 uptake at the interface and the subsequent hydrate growth rate, while λ tends to reduce the surface area of the interface by forcing it to take a flat shape. In our simulations, each voxel of the water and interfacial regions are randomly initialized with a molar fraction of χ c o = 0.029 + δ, where δ is a random number from a Gaussian distribution with mean equal zero and standard deviation σ = 10 −4 . 61,62 This initial molar fraction represents the amount CO 2 in water in equilibrium with the vapor or gas phase (or solubility of CO 2 in water). As Figure  17 in the Appendix shows, this quantity impacts the hydrate crystal formation, growth and morphology. 90,91 In all cases, we consider an initial hydrate thin film with a thickness of 2 × 10 −3 μm in the x direction, extended all along the y direction.
Morphology. In Figure 7, we present the spatiotemporal evolution of a CO 2 hydrate planar front, for ΔT = 1.1 K and λ = 0.5. Although not shown in the figure, the planar front moves slowly beyond 44 μs but eventually stops when there is no driving force to continue growing. To better appreciate the mechanism behind the CO 2 hydrate planar front growing presented in Figure 7, we plot in Figure 8a the CO 2 molar fraction profile in the water region at y = 0.64 μm and along the x axes, for the four panels of Figure 7. As the figure shows, at t = 6 μs (red dotted line), the CO 2 molar fraction in the liquid water equals the initial average value of χ c o = 0.029, at the far right; the CO 2 molecules have not yet diffused from those regions toward the interface. Near the interface, a minimum rapidly evolves. Thus, the planar front motion is viewed as a moving interface determined by the diffusion of CO 2 molecules toward that interface from the bulk liquid phase. This diffusion is driven by the difference in CO 2 molar fractions between the bulk and the region near the interface. As the diffusion continues, the CO 2 molar fraction at the far right decreases, and the growth eventually stops once the concentration gradient is zero. In Figure 8b, the position of the front is tracked by following the minimum CO 2 molar fraction value near the interface. The position of interface follows a power law (length ∝ t 1/2 ), indicating a planar crystal diffusion-controlled growth in which the planar crystals do not grow steadily. Thus, its growing velocity decreases proportionately to t −1/2 . 72 In Figure 9, we present the spatiotemporal evolution of CO 2 hydrate from the hydrate film, but now for ΔT = 9.1 K, keeping the same value for λ = 0.5. As the figure shows, initially, a planar front evolves toward the liquid water region. However, needle-like crystals eventually form due to local random fluctuations at the interface and advance deeper into the liquid water phase with a faster local growth rate than the neighboring planar regions. This behavior indicates that the planar front is unstable for such a large degree of subcooling.
While some crystals grow with a parabolic needle geometry, others develop side branches forming what are known as parabolic dendrites. 72 As Figure 9e shows, the difference in growth rate between the planar regions and the growing crystals is because the CO 2 molar fraction gradient is more prominent in the water region in front of the crystal tips. A stronger gradient implies a more significant flux of CO 2 molecules toward the interface, resulting in a larger local hydrate phase transformation probability and hence a faster growth rate. This positive feedback mechanism is known as Mullins-Sekerka instability. 92 It often leads to pattern formation, one of the most common ones being the dendrite shape pattern in which a tree-like structure develops. This instability also causes columnar, needle, and sword-like patterns observed during the growth of gas hydrate crystals. 23−25 The simulations above suggest that large enough values of ΔT destabilize, through the Mullins-Sekerka mechanism, the sI CO 2 hydrate planar fronts that move toward liquid water with dissolved CO 2 . On the other hand, large values of parameter λ are expected to act as a stabilizing factor for the planar interface because its goal is to reduce the total interfacial area to its minimum (i.e., the interfacial area of a planar front). Thus, to better quantify the interplay between these two competing model parameters, let us introduce the following quantity, (10) where N in is, as mentioned in previous sections, the total number of interfacial voxels at a given time step, and N in o is a constant expressing the total number of interfacial voxels for the special case of a planar interface. ⟨N in ⟩ is the average value of N in over the number of time steps. Thus, for a growing planar interface, M = 1 (i.e., N in = N in o ), while for an irregular one, M > 1 (i.e., N in > N in o ). From Figure 10, where M (see eq 10) is plotted as a function of ΔT for several λ values, it is evident that the extent of ΔT for which R is equal or close to the ideal value of unity increases with λ. This behavior indicates that, indeed, the stability of the planar front increases with λ and decreases with ΔT. In Figure 11, we present the spatiotemporal evolution of sI CO 2 hydrate growth for some of the data points of Figure  10. For λ = 0.01, the planar front is unstable for all simulated values of ΔT. In particular, for ΔT = 2.1 K, parabolic needlelike crystals emerge from the initial hydrate thin film. For ΔT = 5.1 K, the needle-like crystals penetrate deeply into the water region and become parabolic dendrites by developing side branches. Finally, for ΔT = 7.1 K, the side branches of the parabolic dendrites grow larger, and the crystals penetrate further into the water region. The figure also shows that the number of crystals that grow into the water region increases with ΔT for all values of λ except for λ = 1, where only planar fronts are observed. The transition from irregular crystal growth to planar crystal growth as a function of λ, for a given ΔT, occurs gradually.
Hydrate Conversion. We now explore how the fraction of hydrate within the simulation domain evolves in time. This quantity is obtained by counting the number of solid hydrate voxels at each time step, including the initial hydrate thin film voxels, and dividing it by the total number of voxels. Thus, a hydrate fraction of 1 means the whole simulation domain is CO 2 hydrate. Figure 12 presents this conversion rate as a function of time for three representative values of λ and ΔT. For λ = 1 (see Figure 12a), the hydrate grows as a planar front for all the chosen ΔT values (see bottom panels in Figure 11), and hydrate fraction is in proportion to t 1/2 . In other words, we get a planar crystal diffusion-controlled growth. 72 For λ = 0.5, with ΔT = 2.1 and 5.1 K (see Figure 12b), the hydrate fraction also exhibits a planar crystal diffusion-controlled growth, as Figure 11g and Figure 11h show. Although in Figure 11h a protrusion emerges at the end of the simulation, the hydrate grows as a planar front during most of the simulation.
For ΔT = 7.1 K, the dynamic evolution of the hydrate fraction exhibits a different behavior. Early on diffusioncontrolled planar crystal growth is dominant. However, after a certain simulation time, the hydrate fraction suddenly increases steadily (in proportion to t). This behavior is typical of diffusion-controlled growth of needle crystals with parabolic tips, where the growth occurs mainly at the tip region. 72 Eventually, the hydrate fraction reaches a plateau because the driving force to growth diminishes. For λ = 0.01 and ΔT = 7.1 K (see Figure 12c), the hydrate fraction increases steadily with time until a plateau starts to emerge. As Figure 11c shows, needle-like and dendrite parabolic crystals form for such values of subcooling. The hydrate fractions for ΔT = 5.1 and 2.1 K also increase steadily with time after an initial short time, during which the hydrate fraction evolution resembles the one of a planar front. It is interesting to point out that the curve for ΔT = 5.1 K reaches a plateau, while the one for ΔT = 2.1 K does not. This behavior occurs because the one with ΔT = 5.1 approaches a zero driving force within the time scales considered, whereas the other still does not. Panels (a) and (b) in Figure 11 show how the crystals grow for the two subcooling values just mentioned. As expected, the parabolic crystals for ΔT = 2.1 K grow more slowly because of a smaller associated driving force.
Larger Simulation Domain. Finally, in Figure 13, we plot the spatiotemporal evolution of sI CO 2 hydrate toward liquid water for a larger simulation domain. As in previous cases, the stability of the planar front is determined by the competition between the Mullins-Sekerka destabilization mechanism for large ΔT and the surface tension stabilization factor quantified by the parameter λ. However, in this case, the number of available CO 2 molecules in the water domain is higher, and the hydrate growth driving force lasts longer. Thus, the model yields larger parabolic dendrites, with side branches that develop even smaller branches, separated by barely developed Figure 10. Quantity M, quantifying front morphology (see eq 10), as a function of the degree of subcooling ΔT, for several values of λ. Each simulation point is an average over four independent realizations, and the error bars denote ± one standard deviation around that average. The total simulation time for each data point is t = 24 μs. In this case, L x = L y = 0.128 μm. Other parameters are as indicated in the main text.
needle-like parabolic crystals. The density of crystals increases with ΔT and decreases with λ. ■ DISCUSSION AND CONCLUSIONS CO 2 hydrate technologies offer potential alternative solutions to energy and climate related challenges. 69 However, a clear understanding and prediction of morphological variations of CO 2 hydrate crystals are necessary to support the further development of such technologies. Hydrate crystal shape should strongly influence the kinetics of the hydrate formation process, and it should be taken into account in designing equipment, such as transportation devices and reactors.
In this work, we implemented a hybrid probabilistic CA modeling framework to explore morphology variations of CO 2 hydrates growing in stagnant liquid water presaturated with CO 2 . The framework offers an alternative to the PF modeling approach for gas hydrate growth in terms of implementation simplicity and computational efficiency. [54][55][56][57]76 These two aspects are essential when modeling gas hydrate formation in multiphase flow systems where modeling frameworks of gas hydrate growth kinetics must be coupled to computational flow modeling tools. 67 The modeling approach developed and implemented in this work constitutes a first step toward an accurate, efficient, and easy-to-implement computational framework for investigating gas hydrate crystal morphology. Figure 11. (a−l) Spatiotemporal evolution of the simulated sI CO 2 hydrate growing from a seeded hydrate thin film toward liquid water with dissolved CO 2 , for several values of ΔT and λ. All parameters are as in Figure 10, and all panels are at t = 24 μs. After that time, the hydrate barely grew because the driving force for hydrate growth is minimal in some cases or even zero in others. The color bars of the most right panels apply to all panels.
Our simulations focused on the sI CO 2 hydrate crystal growth and morphology variations toward aqueous CO 2 solution from a preformed sI CO 2 hydrate thin film and under conditions typical to oceanic natural gas hydrate reservoirs. Under these conditions, the natural gas hydrates can be converted into the significantly more stable CO 2 hydrate in the presence of liquid CO 2 or aqueous CO 2 solution while natural gas is released. 14,56 Phase transformation and attachment kinetics at the hydrate growing interface and mass transport by diffusion in liquid water are the main physical processes accounted for in our modeling approach. The framework is developed based on the assumption that gas hydrate morphology is controlled by the mass transfer of the CO 2 molecules dissolved in the bulk of the liquid water toward the growing hydrate crystal surfaces, and that the corresponding rate of mass transport increases with the level of subcooling, ΔT. 23−25 This mass transport results from the difference in CO 2 concentration between the bulk of the aqueous solution and that near the hydrate-water interface. 23−25 We explored the morphological growth changes as a function of ΔT and properties of the hydrate-water interface, such as surface tension and curvature. The modeling approach allowed us to explore the stability of the moving sI CO 2 hydrate planar fronts as a function of ΔT and parameter λ, which quantifies, in a qualitative manner, the impact of surface tension and curvature. A higher value of parameter λ acts as a stabilizing factor of the flat interface, while a higher ΔT tends to destabilize it. We found that the planar interfaces follow a power law, length ∝ t 1/2 . We also found that the planar fronts become destabilized for a particular combination of ΔT and λ, giving rise to parabolic needle-like and parabolic dendrite crystals. As expected for parabolic crystals, the crystal growth proceeds steadily with time (i.e., in proportion to t). 72 In many applications, a planar front could be more desired, and in this respect, thorough experimental and computational investigations of the impact of ΔT and surface tension on the stability of gas hydrate planar fronts could be relevant.
An appealing feature of our CA model is that it reproduces complex hydrate morphologies under diffusion-controlled growth. It is possible to extend the framework to study the problem of CO 2 hydrate melting. Although diffusion across and along the interface is neglected, the model can be readily incorporated into the framework. A further extension to a whole 3D simulation domain is possible, and the domain boundary conditions can be easily adapted to the problem under consideration. Other developments of our CA modeling framework could also incorporate the impact of the temperature rise at the interface due to phase transformation and its removal by thermal diffusion. 61,62 For a fixed temperature, the impact of increasing pressure on the morphology variation can also be studied within this framework. 24 Furthermore, the extension of the framework to study gas hydrate morphological variations in oil−water systems instead of the liquid CO 2 − water one studied in this work is an exciting avenue for future research as it is considering other hydrate formers. 93,94 Applying our simulation approach to the case of CH 4 hydrate crystal growth and morphology is a topic for future studies. In this case, bulk-free energy densities entering the transition probabilities of the hydrate phase transformation must be considered for solid CH 4 hydrate and water with dissolved CH 4 . Another aspect to consider is that, given the differences in solubility of CO 2 and CH 4 in water, the kinetic gas hydrate growth in these two cases will differ. 55,95 The solubility of CO 2 in water is higher than the one of CH 4 ; 55,95 therefore, the kinetic rates of CO 2 hydrate growth should be more prominent than the one of CH 4 hydrate under the same operating conditions. 70,95 A first step in comparing these two important hydrate formers has already been performed by Kvamme et al. by using the computationally costly PF approach combined with MD. 54 One can also apply our computationally efficient framework to Xe, which has emerged as attractive laboratory alternative for studying hydrate formation and dissociation. 58,59 Although this work's hybrid probabilistic CA approach reproduces features of diffusion controlled gas hydrate growth, it uses lumped parameters to describe the impact of complex phenomena on that growth. [54][55][56][57]61,62,76 It is however possible to derive the parameter λ governing the strength of the curvature effect on solidification and the noise terms ϵ added to the diffusion fluxes from first-principles. 61,62 Another point that needs attention is a more rigorous treatment of the phase transformation along the water-hydrate interface concerning the time scale of solidification and hydrate stability in terms of cages occupancy. 61,62 Addressing these issues will make the model more predictive. However, as mentioned above and already noted by Buanes et al.,61,62 the computational efficiency of the easy-to-implement CA modeling approach makes it a valuable supplemental tool to more rigorous approaches. Thus, for instance, PF theory could be used to calibrate the free parameters of the hybrid probabilistic CA model by doing simulations in small systems. Then the CA model can be used to perform simulations of larger systems in scales relevant to real experiments, which are not currently attainable using the PF framework.
Although CO 2 hydrate is at the center point in hydrate technologies of CO 2 capture and separation, the number of experimental studies on the crystal growth and morphology of CO 2 is relatively scarce. We still can, however, compare our simulation results to published literature. 5,20,21,24,[31][32][33][34][35]96,97 For example, our simulations agree with the observation of CO 2 hydrate growth and dendrite formation as reported by Tohidi et al. for pressures to 6.5 MPa. 35 Our results are also consistent with the experimental observation of the transition from skeletal polyhedral or columnar to dendritic on CO 2 hydrate crystals in liquid water in contact with gaseous CO 2 , as a function of the subcooling parameter. 5,23,24 However, in contrast to the water−gaseous CO 2 system, under our simulation conditions in which CO 2 is expected to be in the liquid state, the driving force for CO 2 hydrate formation may be higher, 32 meaning that a low degree of subcooling will be needed to destabilize the planar hydrate front and produce dendrite-like morphologies.

■ APPENDIX Free Energy Densities
In this work, we adopt the free energy densities for liquid water with dissolved CO 2 and CO 2 hydrate implemented in refs 54−57, and 76. The free energy densities are parameterized by the CO 2 molar fraction, χ c , the temperature, T, and the pressure, P. Thus, the free energy density of the liquid water with dissolved CO 2 is given as where v m is the average molar volume (assumed independent of either phase or the concentrations), and f l,w and f l,c are the molar free energies of water and CO 2 in the liquid solution, given as where Figure 13. (a−f) As in Figure 12, but for larger scale, L x = L y = 0.5 μm. All panels are at t = 44 μs. The color bars of the most right panels apply to all other panels. is the molar free energy of pure water in kJ/mol, and R is the universal gas constant, v w is the molar volume of water, v c is the molar value of CO 2 , γ l,w (χ c ) is the activity coefficient of water in the solution, γ l,c (χ c ) activity coefficient of CO 2 in the solution, and P o is a reference pressure (normally 1 bar where the molar free energy of water in hydrate in kJ/mol is (19) and that of CO 2 is, is the occupancy, with ν L = 3/23. In this case, a sI CO 2 hydrate structure is assumed, with large cages fully occupied by CO 2 molecules and small ones empty. In the previous expression, is the molar free energy of the empty hydrate in kJ/mol, and (23) is the molar free energy of the CO 2 inclusion. 54

Common-Tangent Construction and Isobaric Phase Diagram
Given the functional forms of the free energies densities, f l and f h , the equilibrium compositions, χ c lh and χ c hl are obtained from (24) and  (25) eqs 24 and 25 express equality of the chemical potentials at equilibrium. Figure 14 shows the free energy densities from eqs 11 and 18 for two different temperatures as a function of the CO 2 molar fraction. Given these free energy densities, the CO 2 equilibrium composition of the coexisting hydrate (χ c hl ) and water (χ c lh ) phases can be obtained by the common tangent points as represented by the dashed lines and stars in those figures. These common tangents are obtained by solving the previous two equations for χ c lh and χ c hl . Figure 14. (a, b) Bulk free energy density as a function of the CO 2 molar fraction for CO 2 sI hydrate (solid blue lines, f h ) and water with dissolved CO 2 solution (dashed black lines, f l ) as specified by eqs 11 and 18 free energy densities are for the case in which the CO 2 molecules only occupy the large cages, with the small ones empty. The dotted line in both figures is the common tangent which expresses that the chemical potentials are equal in the two phases. The stars are the corresponding equilibrium CO 2 compositions in the coexisting hydrate (χ c hl ) and water (χ c lh ) phases. χ c hl is the CO 2 molar fraction in the hydrate in equilibrium with water with dissolved CO 2 , while χ c lh is the CO 2 molar fraction in water in equilibrium with CO 2 hydrate.
In Figure 15, we present the equilibrium isobaric phase diagram for the water-CO 2 system. The two black dashed lines are obtained from the common-tangent method applied to the free energy densities presented in Figure 14. They represent the CO 2 molar fraction in water in equilibrium with sI CO 2 hydrate, χ c lh (liquidus line), and the one in sI CO 2 hydrate in equilibrium with water, χ c hl (solidus line). These two lines, together with the triple equilibrium temperature T eq (solid blue line), define the equilibrium stability region of sI CO 2 hydrate coexisting with liquid water with dissolved CO 2 , denoted as L w + H, in terms of temperature and CO 2 molar fraction, χ c , for the given pressure. We assume that T eq = 282.5 K for such pressure (see solid blue line). Inside the L w + H region, sI CO 2 hydrate, denoted by H, coexists with liquid water with dissolved CO 2 , denoted by L w . Hydrate exposed to water within areas of stability, with respect to pressure and temperature, and with a higher molar fraction than χ c lh , will grow. Hydrate under the same condition but exposed to water with a lower molar fraction than χ c lh will melt. Figure 15 also shows a schematic representation of the CO 2 -in-liquid-water solubility line in a blue dash-dotted style. This line gives the CO 2 molar fraction in liquid water in equilibrium with vapor CO 2 . A real calculation of the triple temperature and CO 2 -in liquid water solubility line requires a free energy expression for the vapor CO 2 and is out of the scope of this work. Details on such calculations, can be found in refs 54−57, and 76. Above the triple point, the hydrate is unstable, and liquid water with dissolved CO 2 coexist with vapor CO 2 (L w + V region). The only stable phase in the region denoted by L w is liquid water with dissolved CO 2 .

Curvature of a Flat Surface
The value 6 in the curvature expression (see eq 4 in main text) is an approximation value for the weighted summation of that expression on a flat interface. It corresponds to κ = 0 and is taken to equal one-half the sum over the entire weighted neighborhood with all voxels in a solid hydrate state 66 (see Figure 16).

Impact of the Initial Concentration of CO 2 in the Water Region
In this section of the Appendix, we briefly explore the impact of the initial amount of CO 2 in the water domain on the hydrate growth. Each voxel of the water simulation domain is started with a CO 2 molar fraction given by χ c o = χ + δ, where δ is a random number from a Gaussian distribution with a mean equal to zero and standard deviation σ = 10 −4 . Here, we will consider the cases, χ = 0.020, 0.025, 0.029, and 0.033.
As Figure 17 shows, the amount of hydrate covering the water region increases as a function of χ at each simulation Figure 15. Equilibrium phase diagram for the liquid water-CO 2 system at a pressure of P = 6.2 MPa. In the region denoted by L w , liquid water with dissolved CO 2 is the stable phase. Liquid water with dissolved CO 2 coexist with vapor CO 2 in the region, L w + V, whereas in the region denoted by L w + H, water with dissolved CO 2 can coexist with a sI CO 2 hydrate with large cages fully occupied and small ones empty. The two dashed lines are obtained from the common-tangent method using the free energy densities as presented in Figure 14 and outlined in this Appendix. The liquidus line (χ c lh ) is the CO 2 molar fraction in water in equilibrium with CO 2 hydrate (left black dashed line). The solidus line (χ c hl ) is the CO 2 molar fraction in hydrate in equilibrium with the water with dissolved CO 2 (right black dashed line). The assumed triple vapor-hydrate-water equilibrium temperature T eq = 282.5 K is presented as a solid blue line, while the CO 2 -in-water solubility line is shown as a blue dash-dotted line. An exact calculation of the latter two lines requires a free energy expression for the vapor CO 2 . Although not shown in this figure, after the solidus line, there is a small region where only solid hydrate, H, is stable. After that region, a coexisting region between hydrate and vapor CO 2 exists (H + V). For details on this phase diagram, we refer to refs 54−57 and 76. Figure 16. Template used in order to calculate the curvature of a flat surface. This curvature is taken to equal one-half the sum over the entire weighted neighborhood with all voxels in a solid hydrate state. Figure 17. Hydrate fraction as a function of time for several initial molar fractions of CO 2 in the water domain. The hydrate fraction is calculated as the number of voxels in the solid phase (hydrate) divided by the total number of voxels in the simulation domain. Thus, a hydrate fraction equal to 1 means the whole simulation domain is CO 2 hydrate. These simulations are for ΔT = 5.1 K, and, as mentioned in the main text, λ, is a dimensionless parameter quantifying the impact of the curvature of the hydrate-water interphase. Each voxel of the water simulation domain is started with a CO 2 molar fraction given by χ c o = χ + δ, where δ is a random number from a Gaussian distribution with a mean equal zero and standard deviation σ = 10 −4 and χ = 0.020, 0.025, 0.029, and 0.033. The corresponding molar fraction of water is time. For χ = 0.020 (solid black line), no hydrate growth is observed for the simulation time considered. For χ = 0.025 (red dotted line) and 0.029 (blue dotted dashed line), a diffusion-controlled growth dominates. Finally, for χ = 0.033, early on, diffusion-controlled planar crystal growth is dominant. However, the hydrate fraction suddenly increases steadily after a specific simulation time. As mentioned in the main text, this later stage corresponds to the diffusioncontrolled growth of needle crystals with parabolic tips, where the growth occurs mainly at the tip region.