A review of molecular dynamics simulations in the designing of effective shale inhibitors: application for drilling with water-based drilling fluids

When drilling for oil and gas, maintaining wellbore stability is of primary importance to reduce non-productive time and trouble cost. Shale swelling causes many problems related to stability when drilling with water-based drilling fluids (WBDF). For many years, it remains the number one cause of time wastage and well abandonment. Different shale samples have different chemical compositions that affect their behavior when in contact with water. This factor makes laboratory-based analysis and characterization of the swelling mechanisms and action of swelling inhibitors extremely challenging. Moreover, the need to replicate different conditions at which clay–water interactions might occur necessitates using a different technique. Molecular dynamics (MD) simulation can be used as a supplement technique to help interpret experimental studies, test and improve a theoretical model, and provide empirical data in high-pressure and high-temperature condition of the borehole. MD simulation applies Newton’s second law of motion to describe particles’ movement in a classical system. The technique can be performed on the time scale of nanoseconds, and in three dimensions, it is thus sufficient for the study of clay–water interaction at a molecular level. It provides a unique view of the clay mineral interlayer and surface activities. This work reviews the progress in MD simulations of clay swelling and its inhibition mechanisms for application in petroleum drilling operations.


Background
Shales are sedimentary rocks in marine basins, composed mainly of mud, silts, and clays. Although burial depth and tectonic stresses may alter their composition in different ways (Gholami et al. 2018), the term shale refers to all kinds of clay-rich sedimentary rocks in its simplicity. They are the most abundant kind of sedimentary rocks, accounting for well over 40% of the stratigraphic column (Al-Ani and Sarapaa 2008). Because of that, when drilling oil and gas wells, most of the time is spent penetrating through shales. Shale formations are responsible for over 90% of all wellbore instability cases in drilling operations (Chen et al. 2003;Gholami et al. 2018). The term clay refers to a natural material made up of powdery minerals and is usually plastic at appropriate hydration state and hardens when put to dry (Al-Ani and Sarapaa 2008;Ray 2013). Clay minerals are phyllosilicate minerals that impart plasticity to clay and account for about 50 to 60 wt.% of most shales. These minerals find use in a great range of applications and are likely the most utilized minerals.
In the petroleum industry, clay minerals are used together with other additives to prepare WBDF. However, drilling with WBDF comes with its share of challenges. It causes swelling of the drilled formation through the absorption of water by the clay minerals, necessitating inhibitor additives in the drilling fluid. Some clay minerals, such as the sodiummontmorillonite (Na-Mt), exhibit a high tendency to swell upon hydration (Salles et al. 2008). They can exceed up to ten times their original sizes (Fink 2012;Taylor and Smith 1986). While such characteristics might be of great value in applications such as a repository of radioactive wastes (Kaufhold et al. 2015), in the petroleum industry, the swelling of clay minerals brings wellbore stability problems with severe consequences (Anderson et al. 2010). The issues may include: below standard hole cleaning, stuck pipe, lost circulation, wellbore tightening, and well caving-in, naming just a few. Wellbore stability problems bring unwanted effects such as drilling delays, the occurrence of life-threatening accidents, loss of equipment, and even well abandonment (Moussa et al. 2017;Tan et al. 1996;van Oort 2003). The high cost of dealing with such situations makes research on wellbore stability an essential aspect of the petroleum industry. In addition to that, the increase in offshore activities in recent years has raised concerns over the ecological effects of the chemicals used in drilling operations (Fraser 2014;Pappworth and Caudle 2016;Zhang et al. 2016a). Worldwide, there are over 12,000 active offshore installations (Ars and Rios 2017). Across the globe, researchers are working on finding organic solutions to clay swelling and its related problems of wellbore instability.
Advancement in computer power and fast algorithms have given researchers an insight into the molecular world of particles in clay-water-inhibitor system interactions and the involved mechanisms behind clay swelling inhibition (Boek et al. 1995a; Moussa et al. 2017). This insight has enabled the identification of molecules that meet the requirements of environment conservation (Pappworth and Caudle 2016;Zhang et al. 2016a) and at the same time, capable of efficiently inhibit swelling, satisfying the requirements of the petroleum industry (Hodder et al. 2010;Suter et al. 2011;Zhong et al. 2013).

Overview of computer simulations
Computer simulation techniques have become extremely useful, if not essential, to understand the underlying principles behind clay swelling in clay-water systems and determine the mechanism behind its inhibition in clay-water-inhibitor systems. The simulations can reveal both the interlayer structure and the dynamics of clay minerals and the full swelling potential and its entropic and energetic components (Anderson et al. 2010). Thus, the computer simulations allow for a clear correlation between clay structure and swelling thermodynamics to be determined. Two families of computer simulations are employed in the studies, quantum mechanics and statistical mechanics, whereby the involved techniques are: computational quantum chemistry (CQC), molecular mechanics (MM), molecular dynamics (MD), and Monte Carlo (MC) simulation. Before looking at the MD simulation technique, it is vital first to understand the basics of clay swelling and have a general understanding. The remaining sections of this chapter are dedicated to the literature that covers just that.

Clay swelling
Swelling of clay occurs in the presence of water when interlayer cations acquire several hydration layers and force the clay layers apart. There is no single theory that describes the swelling of clay (Fig. 1). Over the years, researchers have come up with several ideas to explain the clay swelling mechanism, but only two have become widely accepted. Pioneering works in 1952 by Mooney et al. (1952) and in 1954by Norrish (1954 have shown that clay minerals swell via crystalline and osmotic patterns. (1) Crystalline swelling Crystalline swelling is a process whereby water layers are intercalated between individual 2:1 silicate layers (Laird 2006). The process by which the water layer accumulates in clay interlayer is discrete, which forms a one-, two-, and three-layer (sometimes up to four-and five-layer) hydrates. The forming of hydrate layers inside Mt interlayer space results in an increased interlayer separation (Boek et al. 1995b;Chavez-Paez and dePablo 2001;Norrish 1954). The range of typical interlayer separation in the crystalline swelling regime is from 9 to 20 Å (Nehdi 2014). The compaction pressures for dehydration for the last layer and the second layer of water molecules are about 70,000 and 25,000 psi (~ 480 to ~ 170 MPa) (Olphen 1962). Studies have shown, up to a depth of 10,000 ft (~ 3 km), the formation is dominated by clay-rich sediments with a fair composition of smectites but reduces rapidly past this value (Bjørlykke 1998;Brown et al. 2017;Hower et al. 1976), at such a depth, the significant stress is < 5000 psi (~ 30 MPa) (Wang 2010). However, it is crucial to have in mind that the non-swelling clay minerals, such as illites in deep intervals, will exhibit swelling if they were formed from the smectites transformation (Fink 2012). For drilling environments, in shale, the first several layers of water molecules always exist. Crystalline swelling is, therefore, insignificant to drilling and can occur in all types of clay minerals.
(2) Osmotic swelling Osmotic swelling is the second type of swelling where the concentration of cations between unit layers in a clay mineral is higher than that in the surrounding water. Water is drawn osmotically between the unit layers, and the interlayer spacing increases. This type of swelling can increase interlayer spacing from 20 up to 130 Å resulting in a considerable volume increase (Nehdi 2014). Osmotic swelling results in more massive overall volume increases than crystalline swelling, but only a few minerals, like Na-Mt, swell in this manner (Fink 2012). Smectite clay saturated with potassium cations (K + ) does not swell in this manner and forms only crystalline hydrates in aqueous suspension (Posner et al. 1963). K + saturated Mt's behavior explains why potassium chloride (KCl) has been referred to as the best inorganic salt shale inhibitor in petroleum drilling. It controls the swelling and hydration through the cation exchange reaction between the K + and the clay (Suter et al. 2011).

Montmorillonite clay minerals
Mt is a common clay mineral named after Montmorillon town in France, which was discovered in 1847 (Annabi-Bergaya 2006). It is the main constituent of bentonite that is usually used in the drilling fluid. The Mt is a member of the smectite group with two tetrahedral sheets sandwiching a central octahedral sheet (Wang 2010). It has a crystallographic structure based on the pyrophyllite model (Fig. 2). The isomorphic substitution of atoms in its lattice renders a permanent overall negative charge of the layers, which gives it the ability to absorb certain cations and retain them in an exchangeable state. This property makes it possible for the intercalated cations to be exchanged by other cations in a water solution.
For Mt's physical characteristics, one can visit the work of Ray (2013). This paper's structure is thus: the first chapter presented a background on shale, clay, and clay minerals and their associated problems when drilling with WBDF. Besides, it provided an overview of computer simulations and the roadmap for carrying out a typical MD simulation of clay hydration and swelling. The second chapter lays out the MD simulations requirements and reviews literature that describes the parameters involved in the technique and the progress made in the last two decades. In the third chapter, the different documented literature that suggests clay swelling mechanisms and its crystal structure-property are reviewed. Chapter four presents selected literature on the action of inhibitors in preventing clay swelling. Lastly, chapter five gives an outlook on the future of MD simulations in investigating clay-water-inhibitor interactions.

Molecular dynamics simulation
In any MD simulation system, interaction energies will depend on the potential function forms, the force field parameters, and the thermodynamic ensembles used. Figure 3 shows a typical flow plan for modeling and MD simulation of the clay-water-inhibitor interactions system. The ever-improving molecular simulation techniques and increasing accuracy of the potential energy parameters provide researchers a wide field of endeavor with limitless possibility. In the present section, we are going to introduce the MD simulation technique briefly.

Thermodynamics states
Any thermodynamics system has state variables that describe its macroscopic states, which are the number of particles (N), volume (V), temperature (T), pressure (P), total energy (E), and chemical potential (µ). In molecular simulations, some state variables are external parameters, while others are observables that need to be calculated. The following are the thermodynamic states commonly employed in MD simulations. (1)  can be used to test the adequacy of a time step by checking the conservation of total energy (Rappe et al. 1992).
(2) Canonical ensemble (NVT) (also called constant temperature molecular dynamics (CTMD)). In NVT, E and P are the observables to be calculated. The use of an NVT ensemble requires a thermostat (common ones are Nose-Hoover, the Berendsen, and the Andersen thermostats) to keep the temperature constant by adding or removing energy into the system.
(3) Isobaric-isothermal (NPT) ensemble differs from NVT dynamics in that it requires a barostat in addition to the thermostat. (4) Grand canonical (µVT) ensemble considers simulations with a varying number of particles. In µVT, the system is enclosed in a permeable container that allows adding or removing particles while keeping V, T, and µ constant.

Periodic boundary condition
MD simulations employ periodic boundary conditions (PBC) for a finite system to represent particles' interaction in an infinite real-life environment. PBC approximates an infinite system by replicating the simulation box in all directions, producing the original unit cell's images. A truncation is applied by introducing cut-off distance (r c ) to avoid particles interacting with their images. Usually, the r c is chosen that r c ≤ L 2 , where L is the diameter of the periodic box (Frenkel and Smit 2002). Additionally, all systems containing electrostatic interactions must be charge-neutral to avoid infinite charge generation when replicated.

Force field
A force field refers to the functional form and parameter sets describing the dependence of a system's energy on the coordinates of its particles (González 2011). Since the 1960s when we saw the first application of force fields, several force fields have been developed. The Universal force field (uff) is a generic force field by Rappe et al. (1992), the consistent valence force field (cvff) by Roberts et al. (1988), and the clay force field (Clayff) by Cygan et al. (2004). The Clayff was explicitly developed to describe clay-water interactions. These four are the most commonly used force fields in clay-water system interactions. It is important to remember that any MD simulation results are mostly dependent on the force field used. The availability of different force fields to choose from and the possibility to combine two or more guarantee the most reliable results.

Swelling mechanisms of Mt
When drilling with WBDF through shale formation, water penetrates the clay interlayer and causes it to swell by forming several water layers in a stepwise fashion. However, this way of swelling is not a concern for the petroleum drilling engineers; it only becomes a problem when excess water gets in the formation by the interlayer cations' osmotic action. The cations usually attract more water than the clay surface leading to swelling past the crystalline regime (> 20 Å basal spacing increase). In the following sub-sections, we review the MD simulations that explain the possible mechanisms of Mt swelling.

Effect of water content
It is known that clays swell through transitional phases between discrete layer hydration states. The water content is directly linked to the swelling of clay minerals as the equilibrium hydration states inside the interlayer space of clay minerals increase interlayer separation. With increasing uptake of water by the clay mineral, there will be the formation of water layers in the interlayer region (Berendsen et al. 1987;Georg et al. 2008;Katti et al. 2009;Shroll and Smith 1999;Yang et al. 2017;Zhang et al. 2014). Additionally, water content is related to forming two types of surface complexes, the inner-and outer-sphere complexes. As shown in Fig. 5, at low water content, the interlayer cations form inner-sphere surface complexes with the oxygen of the clay sheet and oxygen of the water molecule. However, at high water content, outer-sphere surface complexes are formed when a part of the innersphere surface complexes break from the clay surface and move in the interlayer's central plane (Hensen and Smit 2002). The loss of coordination to the surface increases the diffusion rates of the interlayer cations. This transition takes place because of the increase in interlayer separation and the number of intercalated water molecules. The inner-and outer-sphere surface complexes will co-exist in all the layer hydrates (Bourg and Sposito 2010;Chang et al. 1995Chang et al. , 1997Chang et al. , 1998. Therefore, as Sun et al. (2015) have shown it (Fig. 5), the increase in water content increases interlayer cations coordination, which adversely increases the Mt water uptake capability and swelling.

Effect of type of interlayer cation
Effect of charge of interlayer cation The valence of interlayer cation strongly influences clay mineral-water interactions (Na and Zhang 2006;Tao et al. 2010;Wu et al. 2015). The monovalent cations have low coulombic interactions with other interlayer particles compared to the higher valence cations (Zheng and Zaoui 2011). Zhang et al. (2016b) investigated the differences in the crystalline swelling regime of Na-and Ca-Mt. They showed that Ca-Mt had more substantial hydration swelling in the crystalline swelling than Na-Mt due to Ca 2+ forming large and stable hydration rings than Na + . However, Tao et al. (2010) showed that Na-Mt overall swelling was higher than that of Ca-Mt and K-Mt. Still, however, it is essential to note that the higher the valence, the higher the cation's hydration energy (Hensen et al. 2001;Li et al. 2019;Seppälä et al. 2016).
Accordingly, from these results, one can deduce that the crystalline swelling extent follows a well-known order for different Mt (Fig. 6). The hydration of the interlayer cation increases with the valence, so is the Mt swelling.

Effect of size of interlayer cation
Studies have shown that the size of the cation influences its residency. Hensen et al. (2001) studied adsorption isotherms of water in Li-, Na-and K-Mt. They found that the big sized K + mainly resided in the interlayer region, the small-sized Na + , and Li + organized in double layers close to the clay surface. However, the works of Moussa et al. 2017, Wu et al. 2015 contradicted this observation. Their results have shown that bigger sized cations such as K + and Cs + tended to sink into the hexagonal cavity and become immobilized. However, despite the contradiction, both simulations can be correct. The interlayer cation residency depends on the total simulation time and how long the system was allowed to stay in equilibrium before doing calculations. Additionally, from the production's trajectories, the mean residency time of the cations at each location in the interlayer space can be obtained and indicate the preferred residency. Whitley and Smith (2004) computed disjoining pressure for Na-, Sr-and Cs-Mt and revealed that at the smallest interlayer spacing, Cs-Mt pressure was greater than that of Na-and Sr-Mt because of the bigger sized Cs + . These variations in disjoining pressure that decomposed into energetic and entropic components led to different swelling free energies of the clays. In another study by Marry et al. (2002), Cs-Mt registered higher diffusion coefficients than Na-Mt because Cs + was less hydrated than Na + , which led to its high diffusion degree. The hydrated radii from their work agreed well with those reported by Zheng and Zaoui (2011) (Fig. 7). The reported MD simulations hydrated radii seem to propose that Mt saturated with smaller sized counterions will swell more than those saturated with larger ones.

Effect of migration and residency of interlayer cation
The mobility of interlayer particles, especially the counterions, impacts Mt's swelling pressure (Sun et al. 2016), an essential quantitative indication of expansive clays' swelling potential. The swelling pressure is the maximum normal force applied to soil in contact with an external solution to maintain a constant volume. The increase in interlayer particle mobility increases Mt interlayer's swelling pressure, leading to more swelling. The mobility of a particle in the interlayer region is dependent on the type of the particle and the confinement of the Mt surfaces (Marry et al. 2002;Ngouana et al. 2014;Zhang et al. 2014).
The residence of interlayer cation determines its degree of hydrations and the consequent swelling extent of Mt. The formation of outer-sphere complexes increases the hydration degree of interlayer cations. The outer-sphere surface complexes will start to form after two-layer hydrate when more water molecules are intercalated. In one-layer hydrate, interlayer cation will interact with the clay surface's oxygen to form inner-sphere surface complexes. Zhang et al. (2014) also found that water spends a long time in the cation's hydration shell than in bulk. Boek et al. (1995b) investigated Wyoming Mt with Li + , Na + , and K + counterions. Analyzing the interlayer structure, they found that both Li + and Na + detach from the clay surface and migrate to the interlayer region where they hydrate, agreeing with observation reported by Chávez-Páez et al. (2001). On the other hand, the K + migrates and binds to the clay surface, making it difficult to hydrate. The present literature suggests that when cations migrate to the interlayer's mid-plane, they tend to hydrate more and cause significant Mt swelling than when they stay close to the surface. However, as necessary as this might be to understand clay swelling, there are still contradictions on interlayer cations' actual residency from several reported findings (visit sub-section "Effect of size of interlayer cation"). Therefore, further MD simulations should be conducted to clarify where the cations are located and postulate why they choose those locations.
The hydration shell's size and stability around interlayer cation affect its ability to diffuse in the interlayer. Zhang et al. (2014) found that Ca 2+ self-diffusion coefficients were lower than that of Na + due to a larger and much stable hydration shell. These results indicate that hydration has an impact on particles' diffusion in the interlayer. Furthermore, larger hydration complexes in the interlayer have fewer diffusion pathways compared to the smaller ones. Figure 8 shows the calculated diffusion coefficients from various MD simulation works.

Effect of layer charge and charge substitution sites
Layer charge and charge substitution sites can cause different swelling behavior of clay minerals (Chávez-Páez et al. 2001;Kenji et al. 2000;Liu et al. 2008;Tambach et al. 2004). Liu et al. (2008) prepared two sets of Mt models by altering charge substitution sites in the tetrahedral sheet. Their simulation results revealed that the distribution of charge sites had a minimal effect on clay swelling, consistent with Ngouana et al. (2014) results.
In another work by Tambach et al. (2004), two sets of 8 unit cell Mt models were prepared. The first set had substitutions in both octahedral and tetrahedral charge sites with a total charge of 6e. The second set had substitutions only in octahedral with a total charge of 8e. The results showed that interlayer cations had more affinity to Mt's surface with both tetra-and octahedral substitution than the one with only octahedral substitution. Additionally, the model with an 8e layer charge had more water content as eight counterions attracted more water than six counterions in the other model. Similar swelling behavior with increasing layer charge due to the significant number of interlayer cations and greater ion hydration enthalpies was reported by Yi et al. (2018). Even so, opposite claims were made by Greathouse and Sposito (1998) and Seppälä et al. (2016). They indicated that Mt with lower layer charges experiences more significant swelling than one with a higher charge. The argument was made that coulombic attractions between the negatively charged tetrahedral-octahedral-tetrahedral (TOT) layers and positively charged interlayer cations adjoin the facing layers. So, the higher the charge, the more firmly the layers are held together. However, it should also be noted if the charge is excluded from the octahedral region, the force cannot be as strong as if it was due to the tetrahedral sites (Teich- McGoldrick et al. 2015). In another study by Sun et al. (2016), the influence of layer charge and charge location on the swelling Fig. 7 Ionic radius influence on cation hydration; values adapted from Zheng and Zaoui (2011) 1 3 pressure was found to be inversely related to the d-spacing (Fig. 9).
A 3D swelling pressure map in Fig. 10 shows the dependency of swelling pressure on the hydration state, charge magnitude, and charge location of Na-Mt.
Since atom substitution in Mt layers affects the number of balancing counter ions in the interlayer, it plays a significant role in Mt swelling (Figs. 9 and 10). As Yi et al. (2018) have reported, the interlayer exchangeable cations have more vital hydration ability than the Mt surface, indicating fewer cations will lead to less Mt water uptake and swelling. Hence, the need for more MD simulations to infer correct Mt swelling behavior and remove the prevailing contradictions on the effect of charge location and amount.

Effect of Mt structure
Suter and associates (2015) employed Ab initio MD simulation to investigate the interlayer and micropore structure of hydrated Mt clays. Their study showed that clay edge sites possess complex amphoteric behavior. Furthermore, the acid/base behavior was dependent on the considered face and the isomorphic substitution sites. They also found that in the tetrahedral substitutions of the clay edge (1 1 0), the adjacent exposed apical oxygens abstracts a proton from nearby water molecule behaving as a Bronsted base, whereas, with substitution in the octahedral layer of the clay edge, instead of drawing off a proton, the exposed apical oxygen increased  Zheng and Zaoui, 2011Chang et al., 1998Zheng and Zaoui, 2011Bourg and Sposito, 2010Zheng and Zaoui, 2011 Diffusion coefficient (10 -10 m 2 s -1 ) Layer hydrate Zheng and Zaoui, 2011Chang et al., 1997Zheng and Zaoui, 2011Chang et al., 1997Bourg and Sposito, 2010 Zheng and Zaoui, 2011Bourg and Sposito, 2010Zhang et al., 2014 Diffusion coefficient (10 -10 m 2 s -1 ) Layer hydrate Zheng and Zaoui, 2011Chang et al., 1997Chang et al., 1997Bourg and Sposito, 2010Zheng and Zaoui, 2011Chang et al., 1998 b) the number of hydrogen-bonded water molecules from one to two. Suter's observation has shown how complex interactions involving clay minerals can become. Therefore, more simulations that consider chemical reactions should be done to further our understanding of clay swelling. In 2017 Guomin Yang and associates studied cation hydration in Ca-and Na-Mt nanopores. They observed an increase in Mt hydrate stable states from one-to five-layer hydrates with increasing water content. In hydration states lower than the five-layer, the formation of a transient innersphere surface complex by Na + showed a strong affinity to clay ditrigonal cavities and the clay surface; however, it was not the case with Ca-Mt. Another study by Li et al. (2019) found that water mobility in the water-saturated mesopores was lower than water mobility in thin films adsorbed on the clay external basal surfaces. They thus revealed for given water content, the mobility of particles in the interlayer was lower than that from the exterior surface. Seeing that clay minerals are formed by the stacking together of clay sheets, many stacking arrangements will give rise to nanoand micropores. Hence, it is essential to understand cations' behavior in these different environments and see how they influence clay swelling.

Effect of hydrogen bonding
In 2006, Tambach and his fellows set up the simulation system to mimic clay mineral in contact with the water reservoir, imposing the water chemical potential and temperature. They found a free-energy barrier for the system to go from a one-to a two-and a three-layer hydrate that the system has first to overcome, which may cause it to stay in a meta-stable state. This finding led them to believe that the breaking and formation of hydrogen bonds within a water layer and between water molecules dominate the shrinking and swelling of clay minerals.

Effect of temperature and pressure
During petroleum drilling operations, the wellbore (Fig. 11) temperature and pressure increase at a gradient of 30 o C/km and 150 bar/km down the wellbore, respectively (Anderson et al. 2010;Bjørlykke 1998). By studying the interlayer spacing of Na-, Ca-, and Na/Ca-Mt, Sun et al. 2015 showed an increase in Mt swelling with temperature rise. The temperature effect is minimal at low water content and more pronounced at high water content. In another work, Zhang et al. (2019) observed insignificant basal spacing changes at the same water content, suggesting the basal spacing is independent of temperature.
The temperature increase does not influence the hydration energy of the interlayer cations. However, it reduces the size of hydration rings of interlayer cations by lowering the cation-water coordination (Fig. 12). The works of Moussa et al. (2017) and Sun et al. (2015) have shown that at high temperatures, water molecules gain enough energy to break free from hydration rings of cations.
Different from temperature, the swelling of Mt is weakly affected by pressure increase. However, Mt swelling response to pressure change is influenced by the type of interlayer cation and water content. Sun et al. (2015) showed that Na-Mt is stable to pressure change at one-and twolayer hydrates, but Ca-Mt begins to be stable from two-layer hydrate (Fig. 13).

Fig. 11
A schematic illustration of the changing wellbore conditions with depths. Water invades the formation (a), low temperature and low pressure near the surface (b), and high temperature and high pressure deep in the wellbore (c) Most literature simulated clay behavior with temperature and/or pressure by fixing one while changing the other. However, as borehole drilling is concerned, they both increase with depth. Therefore, more MD simulations should be performed as it will be interesting to know how that will influence clay swelling and inhibition.

Mechanical properties of montmorillonite
Clay minerals, due to their structure, possess enormous elastic anisotropy. Some MD simulations have been conducted and provided interesting findings. For example, Sayers and Boer (2016) showed that clay's anisotropy was due to strong  (Moussa et al. 2017) covalent bonds within layers and much weaker electrostatic bonds in between. Because of the minerals' anisotropic nature, it is nearly impossible for experiments to be conducted on a single crystal but rather clay films. However, Carrier et al. (2017) showed a significant discrepancy in how thin clay films and clay particles behaved upon hydration. This sub-section introduces some of the MD simulation studies on the elastic behavior of clay crystals and their hydro-mechanical response with increasing water content.

Effect of water content
Different researchers (Carrier et al. 2014;Katti et al. 2007;Schmidt et al. 2005;Zhang et al. 2017;Zheng and Zaoui 2018) have shown that increasing water content decreased the elastic stability of Mt crystal (Fig. 14). Mt crystal configurations changes with the process of swelling and shrinking (Zheng and Zaoui 2018). Besides, Zheng and Zaoui (2018) observed hysteresis phenomena in bulk and shear moduli during both processes. Zhang et al. (2017) studied Na-Mt changes in mechanical properties with a variable component mixture of CO 2 and water. Hydrostatically testing the system at one-and twolayer hydrates, both volume and basal spacing decreased linearly with the pressure. The gradient of the one-layer hydrate was small compared to that of the two-layer hydrate indicating the crystal is much stable at low hydration.
The mechanical behavior of shale has concerned engineers for decades, especially when clay-water interactions are anticipated. When drilling shale formations with WBDF, the borehole structural integrity is significantly impaired. Therefore, these studies are essential not only for protecting the expensive equipment but also for people working in the fields. Carrier et al. (2014) studied elastic properties of swelling Mt and categorized their stiffness tensor results in in-plane coefficients (C 11 , C 22 , C 12 , and C 66 ) and out-of-plane coefficients (C 33 , C 13 , C 23 , C 44 , and C 55 ). They observed that Mt behaved differently depending on the crystal orientation; for instance, water content had a strong influence on the out-of-plane elastic properties than the in-plane properties. Fig. 13 The influence of pressure on the swelling of (a) Na-Mt and (b) Ca-Mt. The deviation in d-spacing is given relative to that at 0.1 MPa. The vertical lines indicate roughly the water contents of the 1-, 2-, and 3-layer hydrates (Sun et al. 2015) Ebrahimi et al. (2012) obtained similar results that led to the conclusion that clay mineral to be anisotropic in the transverse direction and only becomes isotropic after forming a one-layer hydrate. Additionally, Carrier et al. (2014) observed that in the dry state and at low water content, the Ca-Mt was found to be stronger than the Na-Mt along the z-direction.

Effect of Mt structure
As we have seen earlier, in any petroleum drilling operation with WBDF, we can expect a two-layer hydrate forming at the minimum. Additionally, clay minerals along the wall of a drill hole can have any orientation to the hole. Therefore, simulation studies like these are vital in finding molecules that will stabilize the wall and prevent it from collapsing. Carrier et al. (2014) employed an elastic bath method to investigate the elastic properties of swelling Na-and Ca-Mt at finite temperatures. The out-of-plane stiffness tensor coefficients were found to be very sensitive to temperature. When the temperature was increased from 0 to 300 K, the C 33 computed results differed by a factor of three, higher at 0 K and lower at 300 K. In contrast, the in-plane coefficients were not very sensitive to temperature and water content. They observed a small linear change that could be explained by the geometric effect due to swelling. Schmidt et al. (2005) and later-on Katti et al. (2007) found out that applying external stress compressed Na-Mt with more significant deformation in the interlayer than on the clay layers for stress ranging from 0 to 2.9 GPa (Fig. 15). Most of the stress-deformation response for the dry and onelayer hydrate was linear and became nonlinear with increasing water content. They observed that increasing stress resulted in insignificant one-layer deformation. Still, when two and three water layers were introduced in the interlayer, abrupt deformation was discovered at a stress range of about  (Schmidt et al. 2005) 0.5-1 GPa. The two-layer was compressed into monolayer thickness, and the three-layer was subsequently compressed to a two-layer thickness.

Effect of the normal force
Hydrocarbons reserves can be located anywhere between 3 and 8 km in the ground (Titiloye and Skipper 2001). Before drilling, the rock strength and the in situ rock stresses are in equilibrium. However, the balance is disturbed by the drilling process resulting in a potential borehole instability. Therefore, more attention needs to be paid to study the clay minerals' structural stability change in these extreme conditions.

Inhibitor action on montmorillonite swelling
The use of WBDF necessitates the addition of shale inhibitors in the drilling fluid. MD simulation methods have proven to be a great tool in improving our understanding of different molecules' inhibition mechanisms and predicting inherent features for designing suitable ones. Bains et al. (2001) studied clay-water interaction with ethylene glycol (EG) and polyethylene glycol 300 (PEG300, relative molecular mass 300) each in a separate simulation. They found that swelling inhibition is a function of concentration, as the EG performed well only at a very high concentration (> 60%). In contrast, the PEG300 can perform relatively well at a low concentration of about 6%. Also, in their study, they have found that for improved and better performance, inhibitor molecules should possess well-defined hydrophobic regions with polar hydroxyl groups embedded in between. Additionally, the organic molecules should be a long linear chain with significant hydrophobic and hydrophilic areas. The hydrophobic areas will act as a seal to prevent water ingress during drilling operations, whereas the hydrophilic part will bind the Na + to the clay surface, preventing their hydration. Suter et al. (2011) studied swelling inhibition of the cationic and polyether molecules often used in WBDF. They substituted interlayer cation with monocationic swelling inhibitors, tetramethylammonium (TMA), and tetraethylammonium (TEA) cations. With these molecules, the hydration energies reached bulk water levels at lower hydration states than with Na + , indicating swelling inhibition. Additionally, when analyzing their interlayer density profiles in the dry state, both cations were in a monolayer formation and resided near the clay surface's negatively charged sites. As the water content increased, the molecules stayed above charge sites, and the water molecules were located between cations. Suter et al. (2011) also investigated the di-cationic amine hexamethylenediamine (HMDA) inhibition effect in Na-Mt interlayer, assuming a complete cation exchange. Basal spacing remained 13.0 Å from the dry state up to water content of 75 mg g −1 clay and increased sharply to 15.2 Å at a water content of 150 mg g −1 clay, which is a point where the hydration energies are similar to that of bulk water and above which there will be no more swelling. They also investigated polyether diamines PEG and PPG, and the hydration energy was found to be similar to that of the HMDA. However, the much lower water content of about 75 mg g −1 was required for the hydration energy to reach bulk water suggesting a better swelling inhibition than the HMDA.

Action of organic inhibitors
WBDF, in the past few decades, has increasingly been used in drilling for oil and gas because of its low cost and being less environmentally pollutant. Therefore, more simulations should be done to improve the efficiency of swelling inhibitors. It is useful for reducing the drilling costs and the environment, especially when organic inhibitors are used.

Role of the cation
The addition of inorganic molecules in hydrated clay mineral interlayer lowers cation coordination to water reducing the mineral swelling due to cations hydration (Jiafang et al. 2014). The inorganic salts addition promotes cation exchange by introducing foreign cations in the interlayer region of hydrated Mt. For instance, Na + 's exchange by K + from KCl addition will introduce K + that has lower hydration energy than Na + (Fig. 16).
Similar results were obtained by Camara et al. (2019), where they studied the influence of different inorganic salts addition in Na-Mt interlayer. NH 4 Cl addition performed better in reducing basal spacing increase compared to the addition of AlCl 3 , MgCl 2 , and FeCl 3 (Fig. 17). Additionally, inorganic salt addition into Mt interlayer increases ionic interactions between interlayer particles themselves and the clay sheet's minerals. This increase in coulombic interactions improves the crystal's elastic strength (Jiafang et al. 2014).
Despite being less pollutant compared to other drilling fluids, high concentration of inhibitors used with WBDF (e.g., K + from KCl) has raised concern, particularly now that we see an increase in offshore operations (Fraser 2014;Olsgard and Gray 1995;Pappworth and Caudle 2016;Zhang et al. 2016a). Therefore, we need more MD simulations that focus on finding organic cations. Fig. 16 The d-spacings (solid lines, left vertical axis) and discrete hydration energies (dashed lines, right vertical axis) of Wyoming montmorillonite with Na + as the majority counterion (a) and with K + as the only cation (b) (Suter et al. 2011) Fig. 17 Basal spacing changes with temperature for different layered Na-Mt after adding inorganic salts; a-c represent one, two, and three hydration layers, respectively (Camara et al. 2019)

Role of the anion
The negatively charged particles are not likely to belong in the interlayer space of clay minerals. The work of Tournassat et al. (2016) has shown that the Cl − is going to be excluded from the interlayer and will be partially banned from entering nanopores of clay minerals. When the water content was increased gradually to five water layers, Cl − ions had access to the interlayer only from three to five-layer hydrates. Another simulation that reported similar results was by Hedström and Karnland (2012). These results seem to suggest that anions have limited access to the interlayer space. Therefore, at low water contents, the anions in WBDF simply serve as carriers of the cations, and once the dissociation occurs, the anions remain in the bulk phase. More MD simulations can be conducted on anion's role in swelling inhibition when adsorbed on the surface and edges of clays as well as when they are (even if it is for a short time) in the clay nanopores.

Conclusions
This article reviewed MD simulation works on clay-water interactions to investigate clay swelling mechanisms and the inhibitors' actions. The review presented the key insights needed to quantify the observations and understand the mechanisms involved in clay swelling and its inhibition for the petroleum drilling application. In this final section, some of the things not addressed by documented literature are pointed out.
1. In a petroleum drilling operation, a mixture of bentonite and different chemicals are usually used to prepare WBDF. Since mixing different chemicals can have a synergistic effect, it will be interesting to simulate inhibitor and other vital ingredients of the WBDF interactions with clays. Some time ago, due to less computer power, such computationally expensive systems could not be studied. However, there have been considerable advancements in computer hardware recently. Therefore, systems that consider the impact of the presence of other chemicals can be investigated. 2. Different MD simulations in the literature use different models to represent water in clay-water interaction. Some of these works have reported discrepancies when other water-based models are used. However, no clear explanations are given on the decision criteria for choosing any particular model. Following the work of Ferrage et al. (2011), it will be much proper for future simulations to adopt models that best fit the experimental values. Apart from simplifying computations, there is no other meaningful reason given for making a specific choice. 3. There are many contrasting findings on the influence of layer charge and charge location on the swelling of Mt. These differences necessitate a more in-depth analysis of the atomic-scale interactions to quantify the observations and explore the involved mechanisms more accurately. 4. Most MD simulation literature does not analyze Mt swelling as a function of inhibitor concentration. Therefore, it will be captivating to know the optimum concentration for the best inhibitor performance. Besides, knowing this is extremely helpful as it can reduce the costs of drilling and save the environment if it is found that low concentration can achieve a similar effect as high concentration. 5. Finally, to make the MD simulations closely resemble clay-water interactions in real life. Instead of using the non-reactive force fields (e.g., Clayff), it will be compelling to investigate the clay-water-inhibitor interaction using a reactive force field, such as that by Chenoweth et al. (2009).