Mesomorphology of clathrate hydrates from molecular ordering

Clathrate hydrates are crystals formed from a guest–host solution—commonly a gas–water solution, thus the name clathrate hydrates of gas, or simply gas hydrates—where the guest molecules are imprisoned by polyhedral cages (cavities) formed of hydrogen-bonded water molecules.^1–6 The guest molecules do not bond with the water molecules and, therefore, can rotate and translate freely inside the cavities, leading to the entropic stabilization of the cage.⁷ The ratio of guest-to-cavity radii defines the possible geometry of the cages, and examples comprise pentagonal dodecahedra (5¹², that is, 12 pentagonal facets), oblate ellipsoidal tetrakaidecahedron (5¹²6²), spheroidal hexakaidecahedron (5¹²6⁴) tetrakaidecahedra, 5¹²6⁸ and 4³5⁶6³ cavities, and other unusual cages that might form due to geometrical frustration.⁸ The mixing of different cages can tile space, forming clathrate structure I (sI), structure II (sII), structure III (sIII),⁹ structure H (sH),¹⁰ other topological duals of Frank–Kasper phases,¹¹ and others.^12–15 Unoccupied cages also exist, stabilized by the neighbor-occupied cages, giving rise to the concept of occupancy.^16,17 Furthermore, more than one small guest molecule (e.g., H₂ and N₂ under high-pressure conditions) can stabilize one single cage.^18,19 Therefore, the ratio of guest-to-host molecules, known as hydration number, is not unique and depends on (i) the crystalline structure, (ii) the guest molecules, being a mixture of different molecules possible, and (iii) dynamic aspects, that is, how fast the guest molecules diffuse toward the crystalline growth front before the neighbor cages stabilize unoccupied cavities.

The complex nature of clathrate hydrates awakened the scientific curiosity of understanding the microstructures via molecular dynamics (MD) and Monte Carlo (MC) simulations.^20–22 As presented in Fig. 1(b), the microstructure covers the scales of Å to 10 nm and represents (i) the guest and host molecules, the cages, and the crystalline structure, (ii) the influence of ions (e.g., salts, with impacts in long-range interactions) and other molecules (e.g., alcohols and surfactants, implying the use of coarse-grained descriptions) as additives on displacing equilibrium or changing dynamics of molecular ordering, and (iii) strain accumulation via the formation of defects and grain boundaries. However, more than a scientific curiosity, clathrate hydrates also gained an important spotlight in engineering applications driven by the energy sector, especially the petroleum industry. The existence of water under high-pressure and low-temperature conditions in the presence of small hydrocarbons leads to the undesired formation of natural gas hydrates in offshore oil and gas production flowlines.²³ The development of thermodynamic models of phase equilibrium of clathrate hydrates based on van der Waals and Platteeuw ∼60 years ago²⁴ allowed offshore production of hydrocarbon in deep and cold waters by avoiding the formation of clathrate hydrates when using inhibitor additives (e.g., glycols and alcohols).

The more recent desire to reduce gas emissions augmented the attention to the potential use of clathrate hydrates in other applications, including (i) the storage of carbon dioxide in the form of clathrate hydrates in natural reservoirs of water under the ocean under pressure and temperature (P&T) conditions that stabilize clathrate crystals,^25–27 (ii) allowing clathrate hydrates to form without obstructing oil and gas flowlines, called gas hydrate management strategy, thus reducing the use of chemicals and their economical and environmental impacts,²³ (iii) the production of natural gas from marine hydrates as an alternative energy source,^28–30 (iv) the replacement of methane by carbon dioxide in marine hydrates,^31–33 and (v) the desalination³⁴ and recycling³⁵ of water. Of special attention in this study are applications (i) and (ii), which require the description of the dynamics (the timescales) of many length scales where clathrate hydrates manifest,^36–40 herein called the multiscales of clathrate hydrates.

Figure 1(a) shows the current understanding of the coupling between the multiscales of clathrate hydrates. Of extreme importance is distinguishing microstructure from mesomorphology and macromorphology. The microstructure deals with molecular ordering, linking potential energy at the molecule level to the formation of a crystalline structure [Fig. 1(b)]. The coupling of the molecular scale with smaller ones (quantum scale) is rather well described by feeding ab initio calculations to fit potential energy models. Methods to couple molecular to larger scales are, however, less defined. The competition between surface energy minimization and guest diffusion in a mesoscale environment (100 nm to 100 µm) leads to the formation of faceted crystallites, dendrites, and porous media [Fig. 1(c)]. The emergence of mesomorphology requires information about the free energy landscape coming from the molecular scale, as well as the competitive growth of clathrate crystals leading to heat and mass transfer-limited crystallization in the mesoscale.^37,41,42

Porous mesomorphologies^43–45 are undoubtedly the most important to macro-applications because of the impacts of mass transfer-limited crystallization when dealing with (i) low solubility guests in water (e.g., small hydrocarbons and carbon dioxide), (ii) conditions far from equilibrium, and (iii) sheared systems. After the induction time for the nucleation of clathrate hydrates,⁴⁶ the formation of porous patterns occurs in the timescale of a few seconds.³⁷ The remaining liquid water entrapped in the pores gets completely depleted of excess guest molecules in the same timescale, and further crystalline growth depends on the coupling between meso- and macroscales. Namely, (a) guest molecules convect from the macro-environment—the macromorphology, Fig. 1(d)—to the hydrate particles in the timescale of minutes to hours,⁴⁷ then (b) diffuse along the pores, and (c) crystallize to the pore walls.^37,40 As a consequence, the pores seal in time, blocking the pathway of the guest molecules from the macro-environment, and not all water in the system can be converted to clathrate hydrates. The sealing of the pores has a direct consequence (i) on obstructing (plugging) injection tips of carbon dioxide, limiting carbon storage operations, and (ii) on the surface state of the hydrate particles as wet (with a water layer) or dry (without a water layer). The latter impacts liquid-bridge agglomeration and deposition of clathrate hydrates,³⁸ changing the macromorphology of the multiphase flow pattern and leading to the obstruction of oil and gas flowlines.⁴⁰ 

Therefore, the mesoscale is key to bridging fundamental understanding to macro-applications. They, nevertheless, received much less attention in the literature compared to the molecular and macroscales. In this study, we discuss the emergence and evolution of mesomorphology and its coupling with molecular ordering, with a focus on porous patterns. The implications for the macroscale are only briefly discussed, and the reader can find more in previous studies.^{36–40,47,48} This is a multidisciplinary study covering the fundamentals of molecular dynamics, phase field modeling, and transport phenomena. We emphasize the necessity of simplified models that capture the main manifestations of the molecular scale, including diffusion, energetics, and orientation, to serve as a closure for mesoscale models. Several concepts imported from correlated fields (ice, polymers, alloys, and nanoparticles) are put in context to clathrate hydrates considering the anisotropic potential of water molecules, crystallization from a solution, and guest–host effects. We also discuss the open challenges to better understand each scale and their coupling.

The molecular scale describes the structuring of guest and host (water) molecules to form the cages of the clathrate crystalline structure and, therefore, predict the molecular attachment rate—i.e., the growth rate, or the crystalline growth velocity—for a given concentration, pressure, and temperature (CP&T). Simplified (continuum) molecular attachment laws^49–54 are required for the boundary condition of crystallization–diffusion models that describe the mesoscale.³⁷ Such simplified laws, nevertheless, require (i) capturing the main trends of molecular attachment against the local CP&T in the bulk of the molecular dynamics domain and (ii) fitting against MD simulations, because the continuum approach cannot deal with configurational aspects of the structuring of molecules into crystals. Special attention must be given to the role of mass transfer-limited crystallization and guest–host cooperativity effects¹⁷ when dealing with guest molecules of low solubility in water that is of importance in macro-applications (e.g., small hydrocarbons and carbon dioxide).

We first discuss how to capture the dynamics of molecular attachment in a simplified, continuous model. Consider that the nucleation barrier was already overcome, and the existence of a solid–liquid layer of thickness d, where the liquid partially structures around the solid.⁵⁸ We first consider crystallization from a melt, e.g., (i) Lennard-Jones particles forming an fcc structure^50,51 or (ii) liquid water crystallizing into ice.^53,54,59 The partially structured liquid is not yet in the crystalline state, nor in the liquid phase. For crystallization to occur, molecules need to travel a distance similar to the thickness of the solid–liquid layer d ≈ λ ≈ ℓ_molecule that is in the order of magnitude of the mean free path λ, which is in the length scale of the molecule size ℓ_molecule when dealing with crystallization from a melt. Therefore, the rate of molecules that impinge this layer is proportional to D/λ², where D is the diffusivity.

Molecules have anisotropic potential energy and, therefore, require a given orientation to attach to the crystalline structure. In the case of ice, as shown in Fig. 2(e), the water molecules need to orient correctly to form the hydrogen bonds of the hexagonal rings, thus complying with the Bernal–Fowler ice rules. This means that only a fraction (f < 1) of the molecules that impinge the solid–liquid layer present the correct orientation to crystallize. For those that do, molecular attachment occurs at a rate proportional to the Boltzmann factor expressed in the semi-Gibbs ensemble,⁶⁰ exp[−Δμ/(k_BT)], where k_BT is the Boltzmann constant multiplied by the temperature at the solid–liquid layer and Δμ = μ_l − μ_s is the (average) variation of Gibbs free energy per unit mol (chemical potential) during crystallization. The latter is often simplified to the chemical potential variation from the free state (solution) to the solid (crystal) phase, neglecting chemical potential variations close to the solid–liquid interface.

Figure 2(a) shows the behavior of the WF law against subcooling, i.e., the distance from solid–liquid equilibrium simplified to a temperature difference. For low subcoolings, the crystalline growth velocity increases with subcooling due to a greater free energy minimization during phase change. Higher subcoolings, nevertheless, imply reduced temperature, which considerably decreases the diffusivity of molecules and consequently the rate of attempts of molecular attachment. The competition between driving force-limited and diffusion-limited regimes leads to a curved trend of the crystalline growth velocity against subcooling. Notice that diffusion can completely prevent crystallization, marked by a zero growth velocity, leading to the formation of a glass.

The fitting of the WF law to predict the growth velocities of different crystalline structures, at different crystalline orientations, and for different components, is, nevertheless, scarce in the literature. Fitting against MD simulations was developed for Lennard-Jones particles.^50,51 Early research applied the WF law to ice⁵⁹ and only recently was compared to experiments⁵³ and MD simulations^54,56 and expanded to clathrate hydrates.^55,56,62,63 Figure 2(b) compiles MD simulations and fitted WF law from the literature. The main fitting parameter is f, representing the impingement efficiency, but it also absorbs any deficiencies of the simplification hypotheses of the continuum approach. Growth velocities reported for argon modeled as Lennard-Jones particles^50,51 forming fcc achieve ≈10–80 m/s (!), that is, an almost instant crystallization in any macro-application. Ice^53,54 achieves peak velocities of ≈0.1–0.6 m/s, which evidences that the anisotropic potential energy of water molecules [the Bernal–Fowler ice rules, Fig. 2(e)] plays an important role in reducing the value of f. When considering clathrate hydrates formed by high solubility guests in water, the peak values of velocities are even smaller, of ≈0.04 m/s for ethylene oxide (EO) sI hydrates⁵⁵ and ≈10⁻³ m/s for tetrahydrofuran (THF) sII hydrates.⁵⁶ Growth velocities reported for low solubility guests, namely CO₂ and CH₄ sI hydrates,⁵⁷ and at high pressure (400 bars), stay in the same range as of high solubility guests. Lower growth rates are, however, expected at pressures in the range of macro-applications (40–100 bars) due to lower solubility. Furthermore, the MD setup and the crystalline orientation should be kept the same for a direct comparison between high (EO and THF) and low (CO₂ and CH₄) solubility guests, pointing out the necessity of further MD simulations to evaluate the role of guest solubility in water.

We now discuss the mechanisms leading to the variation in growth velocities reported in the literature [Fig. 2(b)] based on statistical arguments. If one considers that the water molecule forms hydrogen bonds with tetrahedron angles,⁶⁴ then each water molecule has ∼$14$ probability of orienting correctly to attach to the crystalline structure. When considering that half the water molecules of the ice hexagonal ring are stabilized by the crystal slab, then three more water molecules need to align correctly to complete the ring. This means that the growth velocity of ice Ih is, on average, only $143≈1.5%$ that of an isotropic molecule (e.g., Lennard-Jones) when considering only orientational configuration aspects. This explains the two orders of magnitude of the difference in peak velocity between ice Ih (basal plane) and fcc (111) of Fig. 2(b). If one further considers the connection of two adjacent hexagonal rings as roughly described by the aligning of one further water molecule, this value drops to $144≈0.4%$ and explains the three orders of magnitude of difference between peak velocities of ice Ih (basal plane) and fcc (100). Clathrate hydrates of high solubility guests as reported in Fig. 2(b) have a growth velocity of one to two orders of magnitude smaller than ice, explained by the larger number of water molecules that require alignment during the formation of clathrate hydrate cages.

When dealing with clathrate hydrates of low solubility guests (crystallization from a solution of low concentration), the guest molecules need to diffuse a distance ℓ_diffusion ≫ d that is much longer than the thickness of the solid–liquid interface d. As depicted in Figs. 2(c) and 2(d), this distance is (i) proportional to the average distance between guest molecules and, therefore, (ii) inversely proportional to the guest concentration, $ℓdiffusion∝4π3NAC1/3$⁠, where N_A is Avogadro’s number with the guest concentration C in [mol/m³]. In other words, the lower the guest concentration, the larger the diffusion distance ℓ_diffusion, implying larger timescales for molecular attachment and smaller crystalline growth velocity, $v∝ℓdiffusion−1∝C1/3$⁠. The solubility of guest molecules in water limits the maximum concentration. Therefore, if one would consider a one-by-one guest molecule diffusion, then carbon dioxide, which is around one order of magnitude more soluble in water than methane,⁶⁵ implies growth velocities of clathrate hydrate crystals of $vsI,CO2vsI,CH4∝sCO2sCH41/3≈101/3≈2$⁠, where s is the solubility of guest molecules in water. The simulations reported in the literature [Fig. 2(b)], nevertheless, evidence that CO₂ sI hydrates grow one order of magnitude faster than CH₄ sI hydrates, both in the (100) crystallographic direction,⁵⁷ and therefore, the one-by-one guest molecule diffusion is not a good hypothesis for clathrate hydrates of low solubility guests.

Clathrate hydrates of low solubility guests present a guest–host cooperativity effect¹⁷ that stabilizes the cages, implying the diffusion of several guest molecules. As depicted in Fig. 2(f), the guest molecules stabilize the rings of water that form the faces of the cages,⁶⁶ giving the time window for the re-orientation of the water molecules at the solid–liquid interface to form hydrogen bonds without interacting with water molecules in the solution. Recognizing the role of this concomitant diffusion of several guest molecules, $vsI,CO2vsI,CH4∝2n$ (valid for saturation), where n is the number of guest molecules required to fulfill the guest–host cooperativity effect. The value of n was reported as 8 to 9 guest molecules per cage when dealing with the nucleation of sII hydrates,¹⁷ but this should be at least half for crystalline growth, because nearly half of the faces of the cages are already stabilized by the existing crystal slab, implying $vsI,CO2vsI,CH4≈10$⁠. This order of magnitude agrees with MD simulations from the literature⁵⁷ shown in Fig. 2(b) and also evidences how sensitive the dynamics of molecular attachment are regarding the solubility of guest molecules in water. Since the solubility of gases in liquids is sensitive to pressure, the growth velocity of clathrate hydrates is also expected to be sensitive to pressure. Furthermore, these estimations consider the saturation of guest molecules in water (that is, the concentration of guest molecules is the same as the solubility), whereas the concentration can be smaller in dynamic systems because of mass transfer-limited crystallization occurring at the meso-environment, to be discussed in Sec. IV A.

The derivation of the WF law of Eq. (1) considers only surface diffusion at the solid–liquid interface as it was developed for crystallization from a melt. Further development is required to consider (i) the diffusion of the solute from the bulk to the solid–liquid interface that is inherent in crystallization from a solution, and of extreme importance for low solubility guests, and (ii) guest–host cooperativity effects that are inherent in multi-component crystals as clathrate hydrates. The development of a continuum approach that captures such manifestations of molecular ordering is of extreme importance to serve as a closure for simulations of larger scales, to be discussed in Sec. IV A.

We now turn our attention to the setup of detailed MD simulations (aka computer experiments⁶⁷) to serve as a dataset to fit the WF law. First, one needs to correctly choose the potential energy model for the determination of dynamic properties. The choice of potential energy models for the guest molecules is rather straightforward, commonly coming from a coarse-grained (united-atom) approach that splits the contributions of (i) attraction due to London–van der Waals dispersion forces and repulsion due to the Pauli exclusion principle of different atoms, commonly represented by a Lennard-Jones potential, (ii) the harmonic bonds and their angles and torsion, and (iii) the distribution of electric charges over the geometry of the molecule. Several such approaches exist, and a common one for guest molecules in clathrate hydrates is the optimized potential for liquid simulations (OPLS).⁶⁸ 

When it comes to the description of the potential energy of water, the most common choices are⁶⁹ (i) the family of Bernard–Fowler/SPC/TIPnP models,^70,71 especially the TIP4P/Ice⁷² when it comes to clathrate hydrates, solved with Ewald summation⁶⁷ or other methods⁷³ to describe long-range electrostatic forces, and (ii) the coarse-grained model of water (mW)^64,74 based on the Stillinger–Weber potential⁷⁵ that reproduces many-body effects via energy minimization with angular orientation, capturing tetrahedron covalent bonds of group 14 of the periodic table (C, Si), and the tetrahedron hydrogen bonds of water molecules in ice and clathrate structures. The mW model is computationally faster as it (a) reduces the number of particles in the molecular dynamics domain by n when compared to a TIPnP model, (b) does not require rigid body algorithms⁷⁶ inherent of TIPnP models, and (c) avoids Ewald summation, at the cost of omitting long-range electrostatics.

Both water models (TIP4P/Ice and mW) reproduce well equilibrium properties of liquid water, ice, and clathrate structures. Dynamic properties are, however, predicted with orders of magnitude of difference.⁷⁷ For example, the mW model overestimates the crystalline growth velocity of ice by two orders of magnitude⁵⁴ when compared to experiments,⁵³ whereas the TIP4P/Ice captures the same order of magnitude. Coarse-grained models are known to overpredict dynamic properties (e.g., growth velocities) because of the reduced degrees of freedom. Two possible reasons for such discrepancies regard (i) the importance of long-range interactions in the estimation of dynamic properties and (ii) the angular energy minimization inherent in the mW model, which turns the molecular attachment rate of clathrate hydrates into a purely orientational problem following Bernal–Fowler ice rules, missing the guest–host cooperativity effect described from Fig. 2(f).

Further challenges that increase the computational cost of simulating clathrate systems come from (i) the large crystal lattices, requiring large simulation domains, and (ii) that major guest molecules of interest in macro-applications present low solubility and diffusivity in water, requiring simulation of long timescales (>μs).^20,22 The emergence of machine-learned models for water^79–81 proved to be a good strategy to deal with long-range interactions at a lower computational cost, capturing the orientational effects of the crystallization of ice.⁸² Understanding if machine-learned potentials capture guest–host cooperativity effects is still required when it comes to clathrate hydrates.

Both MD and MC simulations of clathrate hydrates can predict the stability of the crystalline structure and estimate the occupancy of the cages by the guest molecules.^16,20,74,83 Estimation of free energies by the Frenkel–Ladd method^84–88 provides a basis for testing and developing thermodynamic models (equation of state) that describe the chemical potential of clathrate hydrates in a continuous approach.^89–100 The latter serves to feed the driving force of meso- and macroscale crystallization models in a rigorous manner.^101,102 MD also provides insights into nucleation^{20,21,103,104} and crystalline growth.^{19,55,105–113}

However, when it comes to predicting growth velocities, special attention must be given to heat and mass transfer limitations. That is, when clathrate hydrates form, (i) the released heat takes time to dissipate from the solid–liquid interface to the bulk and (ii) the depletion of guest molecules triggers a mass diffusion process from the bulk to the solid–liquid interface. In steady-state growth, temperature and concentrations at the solid–liquid interface are not the same as the bulk ones of the MD domain. The inherently small simulation domain of MD (<10 nm³) associated with periodic boundary conditions, closed systems (no insertion or removal of molecules), and thermostats that directly act on the overheated molecules close to the solid–liquid interface hinder such a steady-state growth mode and, therefore, are non-ideal for the estimation of molecular attachment rates. The use of heat baths to extract heat far from the solid–liquid interface shows that the solid–liquid interface locally overheats by ≈8 K for argon crystallizing into fcc.⁵⁰ This affects the energetic term of molecular attachment after efficient impingement, that is, the Boltzmann factor in the semi-Gibbs ensemble $exp−ΔμkBT$ of Eq. (1), a key factor for the molecular ordering. It still remains an open question whether such local overheating of the solid–liquid interface is of importance for clathrate hydrates or not, as all MD simulations so far consider thermostats that act in the entire MD domain.

Mass transfer limitation is most probably of major importance compared to heat transfer limitations, especially for guests of low solubility and diffusivity in water. Capturing mass transfer limitation at the molecular scale requires techniques that achieve steady-state growth with a constant guest concentration in the solution.^22,50,51,114 Because the nucleation of clathrate seeds from a bulk solution requires microsecond simulations,^17,20,66 the most common approach consists of defining a crystal slab in coexistence with the solution.^{22,50,51,58,114} By extracting heat, the slab grows until the entire MD simulation domain crystallizes, or until the solution completely depletes in guest molecules.⁵⁷ However, the front of crystalline growth requires a certain timescale to achieve the steady state and to allow for sampling of the growth velocity. For clathrate hydrates, this was achieved by an MD setup that considers a heat pulse that melts one of the solid–liquid interfaces of the crystal slab at the same velocity as the second interface grows.^22,114 In this sense, the velocity of the growing surface can be sampled in the steady state, in a closed simulation box, and with constant concentration.

Such a method, nevertheless, still implies the use of thermostats (e.g., Nosé–Hoover^115,116 or Bussi–Donadio–Parrinello¹¹⁷) to ensure a constant temperature of the system far from the heat pulse. As depicted in Fig. 3(a), fully predicting heat and mass transfer-limited crystallization at the molecular scale requires^50,51 (i) creating heat baths (i.e., velocity rescaling) encompassing regions far from the solid–liquid interface, (ii) solving the domain close to the solid–liquid interface via the microcanonical (NVE) or the isoenthalpic–isobaric (NPH) ensemble, thus ensuring accurate prediction of trajectories of the molecules that impinge the solid–liquid interface, (iii) translating the particles in the simulation domain, so the solid–liquid interface always stays far from the heat bath regions, which implies (iii.a) introducing molecules in the free state (solution) at the correct guest–host proportion to ensure a constant concentration far from the solid–liquid interface, and (iii.b) removing particles from the crystal phase. Aspect (iii) represents the MD simulation domain translating with the solid–liquid interface in the mesoscale, as presented in Figs. 3(b) and 3(c). That is, it represents an open system where the sampling is done in space (the traveled distance by the interface), and not in time as in the classic MD simulation domains that are closed and contain periodic boundary conditions.

The insertion and removal of molecules is nontrivial in MD, but it should be noticed, nevertheless, that such a simulation setup was already achieved in the 1980s for Lennard-Jones particles, and the important advance in computer technology should allow the application of the same method for clathrate hydrates. One important challenge was the large energy drift experienced by heat baths, incompatible with the large timescales required for steady-state sampling. This was, nevertheless, recently overcome by the use of a higher-order discretization of the Verlet algorithm in the heat bath region, called the extended heat exchange (eHEX) method.⁷⁸ 

The modeling of the mesoscale considers control volumes in the 100 nm that are much larger than any conventional MD domain. The solution of all mesoscale control volumes extends to 10 µm and solves temperature and concentration fields (that is, heat and mass transfer limitations in the mesoscale) in complex geometries (e.g., a porous medium). The coupling of molecular and mesoscales, therefore, requires feasible but sound simplification hypotheses. For the molecular scale, (i) any solid–liquid interface, no matter the complexity of the mesoscale geometry, is locally flat (unless dealing with pores smaller than 10 nm, out of context for clathrate hydrates), and (ii) the average temperature and concentration of the mesoscale control volume are the bulk property set in the MD simulation domain. The temperature and concentration of the mesoscale control volume change in the timescale of milliseconds to seconds and can be considered fairly constant for the molecular timescales (∼μs), allowing a quasi-static assumption. From the perspective of the mesoscale, (iii) the rate of consumption of guest molecules and their associated heat release (the latent heat of crystallization) are required as boundary conditions of the control volumes for solving concentration and temperature fields, but (iv) the continuum modeling approach of the mesoscale does not carry information on molecular ordering—namely, free energy minimization, orientational configuration, guest–host cooperativity effects, fluctuations, and heat and mass transfer limitations at the molecular scale, which is of importance when it comes to estimating the molecular attachment rate.

This is where the WF law comes into play. The WF law carries information on the molecular scale and serves as the boundary condition (closure model) for the larger scales. The current meso- and macroscale models of clathrate hydrates,^37,118–125 nevertheless, consider molecular attachment rates as a simple proportionality law with the chemical potential, $dngdt$ = k_iA_iΔμ, where $dngdt$ is the molar rate of guest consumption and A_i is the active surface of crystallization. Energetics-wise, the molecular attachment rate requires coupling with continuum thermodynamic models^101,102 to predict the chemical potential Δμ. Common approaches, however, simplify the problem to a fugacity-based driving force using Henry’s ideality,^{37,118–120,123} to the supersaturation in terms of guest concentration,^{121,122,124,125} or to the subcooling.^126–128 The range of subcoolings for macro-applications is often restricted up to 20 °C, and the linear relationship stops holding close to 20–30 °C for the EO, CO₂, and CH₄ sI hydrates [Fig. 2(b)] and much less for THF sII hydrates (around 3 °C based on the data points available). The use of a linear crystallization law, therefore, can only predict low subcoolings, as depicted by the dashed line in Fig. 2(b), but the definition of what a low subcooling means depends on the crystalline structure and composition.

Even in the cases where a linear crystallization law holds (and simplifies the mathematical efforts of modeling the larger scales), there is a major challenge in determining the kinetic constant of proportionality k_i. Literature data spread over six orders of magnitude even when dealing with similar systems and conditions.^{40,118,119,129–134} Isolating molecular attachment rates from mass transfer resistances in meso- and macroscale experiments is hard, if not impossible. Therefore, fitting such models over macroscale measurements—often in pressure cells/autoclaves/batch reactors,^{101,118–121,135,136} rock-flow cells,^137–143 and flowloops^144–156—requires (i) adopting a mass transfer resistance model and (ii) precisely measuring interfacial surfaces and the heterogeneous concentration and temperature delivered to it. These induce model and apparatus dependencies, which explain the limited success in measuring the constant of proportionality k_i.

The search for good k_i values to predict the molecular attachment of clathrate hydrates was abandoned in the early 2000s, on the grounds that (i) molecular attachment occurs in series to mass transfer resistances in a flat growth front (a facet), as described by the classic two-film crystallization model,^120,121 and (ii) molecular attachment is much faster. When dealing with porous clathrate hydrates, which are often the case in macro-applications, crystalline growth dictates how fast the pores seal.³⁷ Sealing the pores means blocking the pathway of guest and/or host molecules from the macro-environment to the active surface of crystallization. Therefore, precision on k_i values is of extreme importance when it comes to porous mesomorphologies,⁴⁰ which are discussed more thoroughly in Sec. III A.

The fitting of the WF law [Eq. (1)] against MD simulations is an avenue to generalize the kinetics of molecular attachment, thus avoiding the ill-estimated kinetic constants of proportionality k_i. The WF law depends on (i) dynamic and equilibrium aspects of the guest molecule by the diffusivity D and chemical potential variation Δμ, and (ii) different crystalline structures by the value of f relating to orientational and guest–host cooperativity effects necessary for cage stabilization, and (iii) the energetics of different guest molecules plays a role in the chemical potential variation Δμ. The impact of the local concentration of guest molecules (or in the case of a system close to equilibrium, the role of solubility) still needs to be included in the WF law to better predict crystallization from a solution. In addition, the chemical potential variation Δμ needs coupling with continuum thermodynamic models specific to clathrate hydrates.^89–102 Finally, the presence of ions (e.g., salts) and other additives (e.g., alcohol and surfactants) can further modulate the kinetics of molecular attachment, potentially (a) breaking the water ring motifs of the guest–host cooperativity, thus displacing equilibrium conditions and affecting the driving force (i.e., causing an inhibition effect^157–163), (b) capping active sites of crystallization,¹⁶⁴ decreasing the molecular attachment rate, or (c) changing the potential of mean force (PMF) of the solid–liquid surface, influencing the diffusion of guest molecules toward the solid–liquid interface.^165–167 All these aspects still require further evaluation.

We now focus on the emergence of mesomorphology of clathrate hydrate crystals at scales covering 100 nm to 100 µm. The mesoscale is too large to be modeled by discrete, particle-based approaches (MD, MC, and Brownian dynamics) and, therefore, requires continuum approaches, such as phase field models. The mesoscale, nevertheless, “feels” the influence of molecular ordering. Key phenomena comprise (i) surface energy minimization that explains the formation of facets (e.g., Wulff shapes¹⁶⁸), (ii) the Gibbs–Thomson effect, accounting for chemical potential variation with the curvature of the crystalline surface, thus leading to heterogeneous driving force distributions and formation of dendrites, and (iii) competitive crystallization of different surface growth fronts, causing both heat and mass transfer limitations at the mesoscale, which can lead to the formation of dendrites or, in more extreme cases, porous patterns.

In this section, we emphasize the emergence of porous patterns as the main mesomorphology for macro-applications of clathrate hydrates, which, however, received less attention compared to other mesomorphologies occurring under laboratory conditions.

Figure 4 shows mesomorphologies of clathrate hydrates collected from the literature. Faceted crystallites [Figs. 4(a) and 4(b)] can form in the laboratory when dealing with slow crystallization rates (close to equilibrium), quiescent systems, and high solubility guests in water. Facets form by minimizing the surface energy, which depends on the crystallographic direction of the facet. The crystal habit (the crystal shape) that minimizes the total surface energy is called a Wulff shape and is well described by geometric construction models.¹⁶⁸ The hexagon shown in Fig. 4(a) evidences that the crystallographic direction (111) limits growth, and the plate geometry is due to the crystallization assisted by an interface (in this case, a water droplet). The formation of rod-like crystal habits, as the one of Fig. 4(b), cannot be explained by surface energy minimization only, as rods break symmetry from the underlying crystalline structure. Symmetry breaking of crystal habits is an ongoing topic of research in the field of noble metal nanoparticles and relates to strain accumulation via defects and the competition between kinetics and equilibration during crystallization.^173,174 Formal models that capture symmetry breaking are, however, not yet established. Increasing driving forces incur in heat and mass transfer limitations that lead to the formation of dendrites, as depicted in Fig. 4(c), which further grows influenced by the Gibbs–Thomson effect. Dendrite growth is well captured by phase field and diffusion-limited aggregation models.

Effects due to solute depletion (mass transfer-limited crystallization) are strong when dealing with low solubility guests, leading to the formation of porous media^43–45 [Figs. 4(d)–4(h)]. The onset of hydrate formation (i.e., the nucleation visible at the macroscale) occurs after an induction time^46,176 that is highly stochastic and in the timescale of minutes to hours.^{46,137,145,177} Nevertheless, once clathrate hydrates form, they form fast, especially in sheared systems. Figures 4(e)–4(g) present macroscale observations of methane sI hydrates forming in an oil-continuous flow.³⁷ The water droplets, recognized as spherical bodies with reflective interfaces [Fig. 4(e)], are completely saturated with methane and persist flowing without any phase change for around 2 h in this specific experiment (the induction time). Macroscopically important nucleation events are observed by the turbidity of the interfaces [i.e., a change in the refractive index, Fig. 4(f)], followed by the formation of clathrate hydrate particles in a few seconds [Fig. 4(g)]. The fact that all interfaces become non-reflective and non-spherical evidence that liquid water vanishes as a free phase. The liquid water is, however, not fully converted to clathrate hydrates in these few seconds. If the entire amount of methane dissolved in water (∼160 mol/m³ at 278 K and 80 bars⁶⁵) converts to clathrate hydrates with an average hydration number of ∼6, then ∼98% of the particle is liquid water. Notice that any further water conversion requires the transport of guest molecules from the macro-environment to the active surface of crystallization, and this occurs in timescales of minutes to hours.^37,40,47

This evidences that the liquid water gets completely entrapped inside the hydrate particles in the first few seconds after the onset of hydrate formation in oil-continuous systems, which is of extreme importance for the coupling with agglomeration at the macroscale.³⁸ There is still a debate on the exact morphology of the hydrate particle formed in flowing systems as (i) a complete porous structure^{36–38,40,43–45,142} as depicted in Fig. 4(g), also called a spongeous particle because of the hydrophilic nature of clathrate hydrates,¹⁷⁸ or (ii) a porous shell formed around the water droplet,^23,124,125 as occurring in static systems as the one evidenced by Fig. 4(d). The computed tomography of clathrate hydrate particles so far assessed only quiescent systems and showed thin shells (∼50 μm). These shells are, nevertheless, too thin to stabilize mm-sized water droplets under shear stresses of turbulent multiphase flow; thus, sheared systems most probably present spongeous particles, although direct observation is still missing in the literature.

Despite the debate on the hydrate particle morphology, the porous medium formation is rather universal for low solubility guests.⁴³ The shear stresses of the environment have a role in propagating the emergence of the porous pattern to different scales, which can extend to (i) the entire water droplet in the mm-scale when hydrates form in sheared systems, forming spongeous particles, but (ii) only to the scale of 10 µm in quiescent systems, thus forming a porous shell. Even much larger bodies of clathrate hydrates present a porous nature, as presented in Fig. 4(h) for a sample of marine hydrates,⁴⁵ with pores in the scale of 100 nm.

The mesoscale received much less attention than the molecular scale (of fundamental interest to scientists) and the macroscale (of interest to engineers), but their understanding is imperative for bridging knowledge to application. That is, the mesoscale is the missing length scale to allow a complete depiction of clathrate hydrates. Of extreme importance is the determination of (i) the length scale of the propagation of the emergence of the porous medium, (ii) the length scale of the pores and their geometric disposition that dictates tortuosity and interconnectivity, and (iii) their links to thermodynamic conditions (CP&T) and shear stresses.

As shown in Fig. 5, Eq. (4) is a double-well function against the phase field ϕ, with wells located at ϕ = 0 (solution phase) and ϕ = 1 (crystal phase), and an energetic barrier in between both phases proportional to W_i. Parameter $−12<bi(T)<12$ represents the free energy difference between both wells (i.e., Δμ = μ_l − μ_s), where (i) b_i = 0 stands for equilibrium (or growth under nearly equilibrium conditions if the initial state of the system is not at equilibrium) with μ_s = μ_l, (ii) $−12<bi<0$ represents crystallization farther from equilibrium that implies a driving force μ_s < μ_l, and (iii) $0<bi(T)<12$ represents the dissolution of the crystal.

The coupled solution of Eqs. (5) and (6) gives the emergence of the mesoscale pattern and depends on the mobility coefficients M₁ (mobility of the phase field ϕ in the free energy landscape) and M₂ (mobility of the guest molecules in the host medium, related to the diffusivity of the guest molecule in water). Phase field models were extensively applied to clathrate hydrates^{186,190,193–195} under conditions close to equilibrium (b_i = 0) by coupling with free energy models specifically developed for clathrate hydrates.⁸⁹ They describe (i) the growth of crystal fronts¹⁹⁰ and their dissociation,¹⁹⁶ (ii) the emergence of dendrites when considering the Gibbs–Thomson effect to model the interfacial energy term,^41,190 and (iii) the nucleation and modes of growth of clathrate hydrates inside porous sediments.^197–199 That is, phase field models applied so far to clathrate hydrates can describe laboratory measurements for static systems and low driving forces,^{169–171,200–203} as already discussed in Fig. 4.

Systems of interest (low solubility guests, far from equilibrium, and sheared systems), however, lead to the emergence of a porous medium. Phase field modes can predict the emergence of porous patterns when the system concentration is such that it crosses the limit of metastability (spinodal decomposition) leading to a catastrophic phase transition. Clathrate hydrates in macro-applications do not trespass the limit of metastability; otherwise, there would be no induction time.⁴⁶ The porous medium, nevertheless, forms due to a competitive growth mechanism, as presented in Fig. 6. Because the system is far from equilibrium, several nuclei that are close to the critical size coexist. The nucleation of the first critical seed introduces a 3D perturbation that triggers nucleation elsewhere in the system [Figs. 6(a) and 6(b)]. The seeds then grow and deplete the guest molecules from the solution. This depletion is important for low solubility guests, triggering a complex concentration field that limits the growth of the seeds [Fig. 6(c)]. The further growth and merging of the seeds to minimize surface energy lead to the formation of the porous pattern [Fig. 6(d)]. The pores hold the remaining liquid water depleted of guest molecules.

Phase field models above the limit of metastability to predict the emergence of porous patterns of different geometries (e.g., bumps, stripes, and holes¹⁹²) has a long tradition in polymer science^205–208 and in solidification, phase transformation, and critical phase separation of alloys.^{182,191,209–213} The mass transfer-limited nature of clathrate hydrates of low solubility guests points that a porous medium can form even below the spinodal decomposition. Capturing the emergence of porous clathrate hydrates with phase field models still requires molecular aspects that include (i) the functional of the free energy density vs the phase field ϕ(r, t), which is not necessarily represented by the simplified double-well function of Fig. 5, (ii) the nucleation barrier at its links to W_i, (iii) the consideration of systems far from equilibrium and the computation of b_i according to the difference in chemical potential between solution and crystal phases via specific phase equilibrium models,^89–100 and (iv) the mobility over the free energy landscape M₁.

Another important challenge comes in determining the initial condition of the phase field model, where perturbations in the concentration field are required for the nucleation step. These fluctuations have ground in the molecular scale but are transported through larger distances in the existence of turbulent eddies in the Kolmogorov, Taylor, and integral scales in sheared systems. This is the reason why the porous medium formed in water droplets in flowing systems (sponges) propagates to much longer distances than the ones formed in static systems (porous shells). Understanding how the influence of shear stresses can be included in phase field models to change the amplitude of the fluctuations that serve as the initial condition to the solution of Eqs. (5) and (6) is a challenging and fundamental task to allow for understanding the emergence of porous mesomorphologies. Phase field models still require extensive sensitivity analysis to realize simplified expressions for the main features of the emerged porous pattern, including pore length scale, tortuosity, and interconnectivity. The latter are important inputs for crystallization–diffusion models that evaluate the pore sealing in time.

Whereas the formation of porous clathrate hydrates occurs in a few seconds after achieving the induction time, the evolution (sealing) of the porous medium can take minutes to hours because of the transport of guest molecules from the macro-environment to the active surface of crystallization. The prevailing description of transport phenomena in mesoscale growth kinetics of clathrate hydrates is the two-film model.^{120,121,124,125,214} This representation is, however, valid only for a smooth growth front of clathrate hydrates formed in the laboratory. The complexity rises in a porous medium, where concomitant diffusion and crystallization of guest molecules trigger a heterogeneous concentration field along the pores. As crystallization continues, the pores seal in time, eventually blocking further guest molecule diffusion from the macro-environment. In this section, we discuss (i) heat and mass transfer limitations at the mesoscale via crystallization–diffusion models, (ii) a simplified approach to estimate the timescale of pore sealing and the evolution of surface porosity, and (iii) a brief description of how macro-applications are affected by the understanding of pore sealing.

Consider that the porous hydrate medium already formed, Fig. 7, and presents an average pore size r. In addition, consider that any excess guest molecules in the host medium crystallized during the few seconds of the emergence of the porous pattern. The larger environments of macro-applications, however, continue furnishing guest molecules that diffuse in water along the pores and crystallize to the pore walls. The guest concentration C(z = 0) = C₀ at the pore opening to the outer surface (pore tip) serves as a boundary condition to the pore scale. This boundary condition changes in time C₀(t) according to heat and mass transfer limitations occurring at the macroscale, which is the link between meso- and macroscales.³⁷ 

The competitive crystallization–diffusion of guest molecules along the pore leads to a guest concentration distribution, as presented in Fig. 7(a). That is, the driving force Δμ is not the same over the entire active surface of crystallization A_i (the pore walls). Every differential surface dA_i of the pore walls is, nevertheless, submitted to a given thermodynamic state (CP&T) that dictates the molecular attachment rate, which is the link between molecular and mesoscales.

The use of linear crystallization laws gives analytical solutions to Eq. (7) for a variety of boundary conditions.³⁷ These solutions are similar to heat transfer fins,²¹⁵ as expected from the analogy between heat and mass transfer.²¹⁶ Fins are commonly employed in engines and microchips to extend the surface and enhance the net heat transfer. Analogously, the pores of clathrate hydrates extend the active surface of crystallization and enhance the net mass transfer rate. Therefore, the pores of clathrate hydrates can be regarded as “mass transfer fins” that increase the net crystallization rate. Because the clathrate hydrate pores seal in time, the enhancement in mass transfer is not constant.

The analytical solution of Eq. (7) presents an exponential decay of the guest concentration along the pore, as sketched in Fig. 7(a). The total consumption of guest molecules per pore comes by integrating the (linear) crystallization law along the entire (cylindrical) pore, $ng(t)=∫0t∫0∞2πrkiHg/w(C−Hg/wfg/w/h)dzdt$⁠. Upscaling the consumption of guest molecules to the macro-environment requires (i) multiplication by the number of pores per particle, which depends on the average surface porosity of the particle that evolves in time (to be discussed in Sec. IV B), and (ii) the number of particles in the system that can change with agglomeration in sheared systems, discussed elsewhere.³⁸ The estimation of consumption of guest molecules is of extreme importance to macro-applications, relating to (a) the amount of CO₂ stored in the form of clathrate hydrates and (b) flow deceleration⁴⁷ and further slurry destabilization^39,217–219 in the context of gas hydrate management strategies of hydrocarbon production.

The crystallization–diffusion equation, Eq. (7), captures mass transfer-limited crystallization at the pore scale. Nevertheless, the current approach applied to clathrate hydrates still lacks descriptions of (i) a chemical potential-based driving force²²⁰ instead of a fugacity-based one, with further coupling to equations of state specifically developed for clathrate hydrates^89–100 for the estimation of the chemical potential, and (ii) the use of the WF law of Eq. (1) instead of a linear proportionality law, with all implications already discussed in Sec. II. Extension to transient terms in 3D porous geometries is possible by coupling the lattice Boltzmann method (LBM) to the volume of fluid (VoF) method to treat the growing surface, recently developed in the field of hydrated cement paste microstructures.²²¹ 

Crystallization–diffusion models also need expansion to guest mixtures to predict the concomitant growth of different clathrate structures evidenced in experiments.^222–227 Consider, for example, a mixture of methane and propane that commonly forms sII hydrates. Because propane diffuses slower than methane, the ratio of methane/propane concentration increases along the pore. Therefore, certain depths within the pore are achieved by methane only. If methane arrives with sufficient concentration to produce a driving force of crystallization, and in the absence of propane, it forms sI hydrates on top of the already-formed sII hydrate during the onset of formation.

P&T can also present a heterogeneous distribution along the pore, changing the local driving force for the evaluation of molecular attachment rate, and this still requires further evaluation in the literature. Latent heat is released upon the formation of clathrate hydrates, which then requires to be diffused away. Preliminary estimations⁴⁰ show that temperature variations along the pore should not exceed 0.05 K, and therefore, heat transfer limitations seem rather unimportant at the pore scale. The distribution of temperature of a population of clathrate hydrate particles in a sheared system was estimated to achieve up to 0.2 K of heterogeneity,²²⁸ and therefore, heat transfer limitations are of secondary influence in the m-scale. Drastic temperature variations can, however, occur in hydrocarbon production flowlines at the km-scale.⁴⁷ Therefore, it is fair to model the pores as quasi-static in temperature, that is, (i) all pores of all particles in a macroscale control volume of ∼1 m in size present a nearly constant, average temperature, but (ii) this average temperature needs updating with time as clathrate hydrates keep forming in the timescale of minutes to hours. Notice that such conclusions on the importance of heat transfer limitation could, nevertheless, change given the high uncertainty of the kinetic constant of proportionality k_i used in these estimations. We emphasize once again the importance of inputting a better description of molecular attachment models [i.e., the WF law, Eq. (1)] to feed meso- and macroscale clathrate hydrate models.

Finally, pressure can also vary along the pores because (i) complex porous patterns present the distribution of radii of curvature of the solid–liquid interface, leading to different driving forces for nucleation and growth via the Gibbs–Thomson effect,¹⁹⁷ (ii) the number of particles in the system that can change with agglomeration in sheared systems, discussed elsewhere.^{37,204,229–231} and (iii) the pressure of the external medium can vary due to (iii.a) the formation of liquid bridges upon collisions between clathrate hydrate particles, whose curvature implies a variation in the local pressure via the Young–Laplace effect, (iii.b) the fluctuating nature of the dynamic pressure of sheared systems (i.e., turbulent multiphase flow), and (iii.c) the pressure drop of the system as a whole, e.g., as particles flow along a flowline. An estimation of these effects requires coupling of the crystallization–diffusion models with permeability models (e.g., pore network models²³²) and multiphase flow models.^47,233

There are several challenges to better estimate the evolution of the surface porosity of hydrate particles. A 3D description of the porous medium can be achieved by using the lattice Boltzmann method^234–239 or pore network models.^240–244 In the field of clathrate hydrates, these methods are commonly applied to estimate the permeability of hydrate-bearing sediments, that is, hydrates formed inside a porous rock, but could be extended to estimate the permeability of the porous pattern of clathrate hydrates themselves, of importance to agglomeration in the macroscale.³⁸ Fluid–structure interactions also occur when clathrate hydrates form in a sheared environment (e.g., a turbulent multiphase flow), leading to particle squeezing and reopening of the pores.³⁶ The curve-fitted parameter λ of Eq. (9) represents such a phenomenon³⁷ but does not bring any information on the mechanical strength of the porous structure. In addition, the development of phase field models to capture the emergence of porous clathrate hydrates, as discussed in Sec. III B, will feed the initial porosity ɛ_in and the length scale of the pore radius r. Finally, the value of the local concentration C₀ at the pore tip requires coupling to larger scales via the mass transfer resistance system that depends on mass transfer coefficients (shear rate) and interfacial surfaces (i.e., the multiphase flow pattern^245–247) dictated by the sheared conditions. The latter was already accomplished for steady-state multiphase flow with dilute suspension hypotheses^40,47,248 but still requires evaluation using more sophisticated flow models.^{233,249–254}

Another important challenge remains in validating porous medium models. Visualization is easier in systems containing high solubility guests in water because they can form at ambient pressure, but those do not form porous patterns. Expanding visualization to systems of interest in macro-applications requires visualization techniques in pressurized systems. Common techniques include (i) particle video microscopy (PVM)^149,153,255 to visualize scales covering 10 to 100 µm (i.e., the smaller droplets and particles of the macromorphology) and (ii) high-speed imaging^{37,148,152,155,156} to visualize the scales covering mm to cm (the larger droplets and particles, and agglomerates and deposits). Unfortunately, these visualization techniques cannot handle the observation of the porous medium.

X-ray computed tomography was applied for hydrate formation over static droplets,^172,256 but the sapphire window thickness required for the high-pressure system does not allow for positioning the beam close enough to the sample to achieve resolution at the pore scale. Another technique is inserting clathrate hydrate samples in a scanning electron microscope (SEM),^43–45,257 giving an estimate of the surface porosity, but not able to fully reconstruct the 3D porous medium. It is also important to notice that the use of tomography and microscopy techniques misses the dynamic evolution of the porous medium, which is intrinsically related to the shear stresses of the flow. One could extract samples at different time lapses after the onset of formation in an experiment in a sheared system (e.g., a rock-flow cell^137–143 or a flowloop^144–156), cryogenize them to avoid any further diffusion and evolution of the porous medium, and then analyze them with tomography or microscopy. This is nontrivial as it requires dealing with samples at high pressure; otherwise, the results would be biased by a quick dissociation of the clathrate hydrates. A direct validation of phase field models to understand the emergence of porous patterns, and of crystallization-diffusion models to understand pore sealing, can only be achieved with the observation of the porous medium and are key to bridging the molecular scale (fundamental science) to the macroscale (engineering applications).

We briefly discuss the implications of the description of the mesoscale to the macroscale. As already discussed in the Introduction, one macro-application of clathrate hydrate is in storing carbon dioxide in water reservoirs under the ocean.^25–27 In this case, the macro-environment consists of a tip drilled into the reservoir. The reservoir is a porous bed that contains water (brine). The tip injects carbon dioxide, which solubilizes in water and diffuses through the porous sediments. Clathrate hydrates nucleate and grow until blocking the pores, causing a reduction in the permeability of the pore sediments.²⁶ The diffusion of guest molecules can further occur through water entrapped in the pores of the clathrate hydrates. Nevertheless, there is a timescale for the complete sealing of the pores, Eq. (8). Once the pores get sealed (both the pores of the sediments and the pores of the clathrate hydrates), further carbon dioxide injection from the same injection tip is not possible. Whereas sequestration and storage of ∼250 GtC (gigatonnes of carbon) are required²⁵⁸ in a period of 50 years to limit global warming to 3 K, the water reservoirs under the ocean discovered so far²⁶ can theoretically store 10³–10⁵ GtC, which is promising. Nevertheless, the mesoscale prevents the complete conversion of the water reservoirs into clathrate hydrates because of pore sealing. Notice that drilling several injection tips threatens the economic feasibility of such a hydrate-based application, but hopefully, a better description of the mesoscale in the future could bring together the necessary simulation tools to allow for creative and innovative engineering solutions.

Another important application is allowing clathrate hydrate formation without important agglomeration and deposition in hydrocarbon production flowlines as a means of reducing the use of chemical additives and the consequent economic and environmental impact of the production operation.²³ In this case, the macro-environment is a kilometer-long flowline linking the oil well to a platform or the coast. Clathrate hydrates can form because water is always produced alongside oil and natural gas, in a high-pressure and low-temperature environment. It has recently been shown that the agglomeration of hydrate particles in oil-dominant flows can be prevented by quickly sealing the hydrate pores, thus avoiding the formation of water bridges between colliding particles.³⁸ That is, by “plugging” the clathrate hydrate pores at the mesoscale, one can prevent plugging the flowlines at the kilometer scale. Further coupling with multiphase flow models to predict slurry suspension^{217–219,259–261} then makes possible hydrocarbon production with a smaller amount of additives.³⁹ 

This study gave a perspective on the emergence and evolution of mesomorphologies of clathrate hydrates and their links to molecular ordering. Unlike classic, isotropic (e.g., Lennard-Jones) particles, molecular ordering into clathrate hydrates deals with anisotropic potential energies (water) that induce orientational effects (the Bernal–Fowler ice rules) and guest–host cooperativity effects to stabilize the clathrate structure. We highlighted the use of the Wilson–Frenkel law as a means to generalize crystalline growth velocities instead of the current linear crystallization law that prevails in the literature. We discussed the necessary molecular dynamics setup to fit the structural information into this simplified, continuum approach.

We then discussed mesoscale phase field models to predict the emergence of mesomorphologies, with an emphasis on porous patterns that form in low solubility guests and conditions far from equilibrium that are common in macro-applications. Phase field models couple the free energy landscape from the molecular scale with mesoscale phenomena comprising interfacial energy minimization and competitive growth and merging of seeds. Whereas the estimation of driving forces is well-established based on the van der Waals and Platteeuw model, the shape of the free energy landscape and the mobility of the phase field over it are not. We also discussed the role of perturbations in sheared systems in propagating the formation of the porous pattern to longer length scales.

The porous medium further evolves depending on transport phenomena from the macro-environment, diffusion of guest molecules along the pores, and crystallization to the pore walls. The Wilson–Frenkel law is a potential closure relationship to enhance the capabilities of crystallization–diffusion models considering averaged molecular ordering effects. Understanding the timescale of pore sealing considering these multiscale manifestations is an avenue to allow for macro-applications of clathrate hydrates.

C.L.B. acknowledges the sponsorship of the Alexander von Humboldt Foundation through the Humboldt Research Fellowship for postdoctoral researchers and the Emerging Talents Initiative (ETI) of the Friedrich-Alexander-Universität Erlangen-Nürnberg.

The authors have no conflicts to disclose.

Carlos L. Bassani: Conceptualization (equal); Formal analysis (lead); Funding acquisition (lead); Writing – original draft (lead). Michael Engel: Writing – review & editing (equal). Amadeu K. Sum: Conceptualization (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article.