A universal hydro-mechanical coupled behavior model for clay-bearing strata—Molecular-level simulation approach

Clay minerals in soils and rocks exhibit large volume change upon interaction with water and this behavior becomes even more complex when the strata are being stressed by the engineering and environmental loads. Therefore, a realistic prediction of the hydro-mechanical behavior of the clay-bearing strata is always a challenge due to their coupled swelling-mechanical response in the cases of geotechnical and geoenvironmental engineering problems, nuclear waste storage in clay-bearing rock repositories, shale gas extraction, and other uses of clay in the manufacturing industry. All the existing behavior models have restricted applications in the engineering and other fields of practice mainly due to the partial consideration of the structure and fabric of clay-bearing strata in the model formulation. In this study, a hydro-mechanical behavior model has been formulated using the parameters acquired from the molecular-level simulations and modeling of the volume change and stress–strain behavior of the clay-bearing structure. The Molecular Mechanics and Molecular Dynamic simulations were performed on the natural structure of the clay-bearing strata formulated using Monte Carlo technique. The mathematical model, developed from the simulation results, can predict the overall hydro-mechanical behavior of clay-bearing strata for all possible combinations of clay minerals, non-clay minerals, salts causing cementation of the soil/rock structure, confining pressures, and the induced strain levels. The developed model has successfully been validated through laboratory and field testing on the clay-bearing strata in both the elastic and elasto-plastic regions of the stress–strain behavior and also from the data of two (02) swelling clays (MX-80 and FEBEX Bentonite) from the existing literature, supporting the universal nature of the developed behavior model.

Clay-bearing strata are characterized by the presence of the swelling clay minerals in the soil or rock structures.Clay minerals present in soils and rocks exhibit a high degree of volume change upon interaction with water and this behavior becomes even more complex when the strata are being stressed by the engineering and environmental loads.Clay-bearing strata are known to be problematic because of their volume change characteristics resulting from a continuous variation in the moisture regime of the unsaturated and partially saturated zones.These volume changes are accompanied by variations in the mechanical behavior of clay-bearing strata.The mechanical properties of natural and compacted clays are critical in the foundations' design and containment barriers for nuclear waste [1][2][3][4][5] .Foundations of the buildings are generally placed in partially saturated soil layers of the vadose zone.The increase in the moisture levels under the building foundations due to reduced evaporation, percolation of surface water(s), and the possibility for leakage from utilities will lead to significant changes in the mechanical properties of these soils.The hydro-mechanical behavior of compacted clays is also critical in the design of clay barriers for nuclear waste.The clay barrier is initially placed in an unsaturated state and is exposed to moisture fluctuations during its lifetime caused by heat-emitting waste and hydration from the parent/host rock.Similarly, slopes and excavation work in mudrocks/clayey rocks also require special attention to study the hydro-mechanical behavior 6 .Excavations in clayey rocks create unloading conditions (reduction in the confining stress) and expose the surrounding material to atmospheric interaction causing cyclic changes in the moisture levels.These changes in environmental and loading conditions lead to swelling of the exposed clayey rocks and the corresponding reduction in stiffness and strength.Further, the clays have applications in agriculture soils, landfill liners, and even in the field of medicine 7 .Therefore, understanding the hydro-mechanical behavior of clay-bearing strata is important for the application of clay-bearing soils/rocks in various disciplines and civil engineering.
The structure of the natural and compacted clay-bearing strata is very complex since they constitute different minerals with random arrangements and complex interactions of the particles ranging from millimeters to nanometers 8,9 .The overall material behavior is controlled by the structure and fabric of clay-bearing strata created by several factors such as clay-water-ion interactions, relative distribution and proportion of clay and non-clay particles, density, and moisture conditions 10 .Due to above-reported factors, the realistic modeling of clay-bearing soils/rocks has been very challenging for geotechnical engineers.Therefore, clay minerals are nano-materials and molecular-level simulations approach can be applied to determine the overall performance of clay-bearing soils/rocks.
Several macro-level reports 2,6,8,[11][12][13][14][15][16][17][18][19][20][21][22][23] , as well as molecular-level studies 7,[24][25][26][27][28][29][30][31][32][33][34] , have been conducted to model the volume change and mechanical performance of clay-bearing strata.All macro-level constitutive models are applicable to specific behavior(s) and cannot be generalized in the form of one universal model.The macrolevel models also do not explicitly incorporate realistic nano-level clay interactions with water and ions in their formulations.On the other hand, the existing molecular-level models for clay-bearing soils/rocks have mainly used a single clay crystallite in their simulations without incorporating the natural fabric and structure of clays.In addition, the existing molecular-level models also lack in formulating a mathematical model to assess the volume change and mechanical behavior.In the light of the above-described limitations of macro and molecular level models, it is imperative to perform realistic modeling of clay-bearing strata at the molecular level to predict the hydro-mechanical behavior with due consideration of the clay interactions with water and ions at the nano-level and linking it to the macroscopic level by considering all possible variations in the structure and fabric of clay-bearing soils/rocks.
In the current report, the authors present a comprehensive molecular-level simulations-based hydro-mechanical model for clay-bearing strata, which could be utilized for all possible combinations of clay and non-clay minerals with cementing agents.The model developed in this study is based on the basic intrinsic parameters utilized in various types and kinds of clay-bearing soils/rocks available worldwide.This research is comprised of three major levels: molecular-level simulations; model formulation; and model validation.Molecular-level simulations were carried out using material studio software followed by the model formulation using simulation results.Validation of the model was then performed through field and laboratory testing in the current study and utilizing the available database from the literature.All these phases of the study are discussed in the subsequent sections.

Model construction
The natural composition of clay-bearing strata is very complex owing to the presence of swelling clay minerals such as Na-montmorillonite (Na-Mt), non-swelling clay minerals such as kaolinite and interaction of these clay minerals with other non-clay minerals like quartz, gypsum, calcite and other salts.Cations and anions from these salts interact with the swelling clay minerals and cause cementation.Quartz acts as filler in the pore spaces and do not contribute to the volume change process.Furthermore, macro, micro and nano level pores exist in the clay-bearing strata with water and air phases.The overall volume change and mechanical behavior of the claybearing strata are thus controlled by the swelling clay minerals, density, degree of saturation and cementation between clay minerals.All such possible clay-water-ion compositions and interactions in the clay-bearing strata are simulated in this study through Materials Studio software 35 using molecular dynamics (MD), molecular mechanics (MM), and Monte Carlo (MC) techniques.
It is well established that the swelling mechanism of clays is governed by the interactions and processes among the nano-level clay particles, called as crystallites 31,[36][37][38][39] .The swelling response, in turn, defines the mechanical performance of these soils/rocks 7,[27][28][29]32 . In his study, the swelling-mechanical response of clay-bearing strata is modeled through molecular-level simulations of Na-montmorillonite with different cation exchange capacities (54, 90, and 144 meq/100 g).Na-montmorillonite has a tetrahedral-octahedral-tetrahedral (TOT) structure in which silica tetrahedral sheets are composed of one silicon atom surrounded by four oxygen atoms and the octahedral sheets contain aluminum or magnesium that are octahedrally connected by hydroxyls and oxygen atoms.The clay particles form by stacking of these sheets with dimensions of an order of nanometer and a length or width to thickness ratio of about 2000:1 9 .Based on the already published values in literature [40][41][42][43][44][45] and XRD findings by Ahmed and Abduljauwad 7 , a basic crystallite size of 2•6 × 5•4 × 2•0 nm was selected for the simulation of Na-montmorillonite (Fig. 1a).These individual crystallites were then grouped in a periodic unit cell, utilizing the MC technique to form bigger clay particles.The clay particles then combine to form aggregates that join with other non-clay particles to generate the overall fabric and structure of clay-bearing soils/rocks having micropores between clay particles, inter-aggregate and intra-aggregate pore spaces.Most of the unit molecules used in this study were obtained from the Nanoscale Simulation Laboratory at the University of Akron, USA 46 .These unit molecules consist of Na-montmorillonite with different cation exchange capacities (54, 90, and 144 meq/100 g), water, Calcite, Gypsum, and Pyrophyllite.Some other molecules were also drafted using Materials Studio software.All these molecules' charges were validated using the method of charge equilibration QEq 35 .

Simulation parameters
One of the basic steps in molecular-level simulations is the sorption of water molecules in the clay crystallites to simulate the swelling process.The sorption module of the software was used for the sorption of water on the selected clay crystallites using MC method (Fig. 1b).The MC method checks the fixing of the defined number

Simulation methodology
A six (06) step simulation methodology comprising: (i) water molecules sorption on the individual clay crystallites, (ii) salts sorption, such as calcite and gypsum on clay crystallites, (iii) creation of soil fabric in the loose state by random sorption of four crystallites in the unit periodic cell, (iv) compaction of the soil to the required density, (v) stress relief of the compacted unit cell, and (vi) sorption of water to produce different swelling stages, was adopted using selected variations of cation exchange capacity (CEC), exchangeable cations, density, moisture content, and salts.These steps of the simulation program were repeated for three CEC changes in Na-montmorillonite (144 meq/100 g, and CECs of 90, 54) and subsequent variations in density, moisture condition, exchangeable cation, and salt.In the simulation for Na-montmorillonite at three different CECs (Na as 100% exchangeable cation), the clay particles were compacted at 0.01, 0.1, and 1.0 GPa at moisture contents of 10, 20, 30, and 40%.Cementation agents, such as Calcite (CaCO 3 ) and Gypsum (CaSO 4 .2H 2 O), were added to Na-montmorillonite at 10 to 20%.Simulation permutations with these salts were performed at 10% for Calcite and 10 and 30% for Gypsum.Exchangeable cations (i.e., Na + , K + , Mg +2 , and Ca +2 ) were also varied at 20, 40, and 60% with different combinations to investigate their effect on swelling and subsequent mechanical behavior.
The sorption of water molecules on the selected clay crystallites was achieved using the sorption module.After sorption of a maximum number of possible water molecules in every 25,000 steps, the system's energy was minimized using the software's molecular dynamics technique in the Forcite module.After the sorption of water molecules on single clay crystallites, the individual ions, e.g., Ca +2 , SO 4 -2 , etc., were sorbed on the waterbearing clay crystallites to simulate the process of cementation due to salts.A typical post-sorption particle of Na-montmorillonite is shown in Fig. 2a, and an enlarged view showing various inter-layer and intra-layer cations is provided in Fig. 2b.
During the deposition process, the clay particles first pack in a loose state and then get compacted under the consolidation process.Therefore, it was essential to simulate the clays' fabric at both stages.Using the Sorption module, simulation of the loose state is performed through random sorption of four crystallites with periodic boundary conditions in a 12•5 × 12•5 × 12•5 nm cubic unit cell (Fig. 3a).These loose particles were then compacted to the maximum achievable density using the NPT ensemble in the Forcite molecular dynamics module (Fig. 3b).Confining pressures of 0.01, 0.10, and 1.0 GPa were employed to model the different types and levels of laboratory compaction pressures and geological pressures.Each simulation was run for a varying period (ranging from 10 to 30 ns).The simulation time period was extended to a maximum limit where it was made sure that no more increase in density, under that confining pressure, is achievable and it has become constant.This is the stage where the particles have achieved the maximum possible packing by occupying all the possible empty spaces.This packing is governed by the size, shape and gradation of the clay particles.
Clay-bearing strata are subjected to the release of geological pressures (stress relief), making them overconsolidated in nature and prone to swelling.In the previous step, the Forcite MD module was run at a confining pressure of 0•001 GPa to simulate the stress relief of the compacted unit cell.The selected pressure of 0•001 GPa could be considered a representative of the over-consolidation pressure in most of the clay-bearing soils/rocks.The clay particles mixture after stress relief was subjected to sorption of water molecules to create different swelling stages.The process was repeated for different percentages of moisture content.
In this study, each absorption/adsorption run was carried out using Monte Carlo technique.In this simulation, the water molecules start absorbing/adsorbing to the highest energy points in the clay structure and continue till all the energy (water deficiency) points are covered by water molecules.Hence, in this manner, the simulation is continued till the adsorption/absorption process of the water molecules is completed.

Simulation outcomes
The mechanical behavior model formulated in this study considered cohesive energy density (CED), density, and moisture content as three (03) state parameters.The Forcite module of Material Studio software was used to determine CED and moduli values for all likely variations in the soil/rock structure.
CED is the energy required to separate molecules from each other and is given as: (1) where E coh is the cohesive energy of all molecules in Joules, E total is the total energy of the system, and E inter and E intra are the intermolecular and intramolecular energies, respectively in Joules.The cohesive energy per unit volume defines the CED and is generally represented in units of J/cm 3 .It is the measure of mutual attraction between the molecules 7 .
Both CED concept and volume change processes of clay-bearing strata are quite analogous to each other and, therefore, CED parameter is found to correlate well with all possible variations in the structure of clay-bearing soils/rocks resulting from changes in density, moisture content, CEC, and exchangeable cations 7 .
The CED value depends on the net result of all the interactions taking place in a clay-water-ion system including ionic hydration, adsorption, double layer repulsion, electrostatic attractions and van der Waals interactions 7 .The CED parameter varies with the changes in density, moisture content, cementation, exchangeable cations, and CEC and is directly/indirectly associated with all of these variables.In this study, the clay crystallites with high CEC have been observed to yield high CED values as compared to low CEC crystallites for the same density and moisture.This is because of the more number of hydrogen bonds in high CEC clays due to the higher charge deficiency centers which results in more electrostatic attractions in clay particles.The CED also shows an increasing trend with the increased bivalent cations and cementation.This is due to the extra bonding created by the cations and anions of cementing agents such as Calcite, Gypsum, etc., and/or exchangeable cations other than sodium.The clay crystallites with high densities have been found to achieve higher CED values due to close proximity of the crystallites or lesser moisture contents.The increase in CED owes to the additional attractions due to close range attractive forces on the crystallites.However, as the clay particles come closer due to the high www.nature.com/scientificreports/electrostatic attraction, the van der Waals forces become repulsive resulting in additional expansion.For the same CEC, the higher water content results in a decrease in the CED value due to balancing of charge deficiency centers by water molecules thereby lowering the electrostatic attractions.In conclusion, the CED parameter depends on density, moisture content, cementation, exchangeable cations, and CEC and controls the swelling and mechanical response of clay-bearing strata.The moduli values were determined using the Forcite module in the material studio at each swelling stage.The constant strain method was used for the modulus determination at strain levels varying from 0.001% to 10% in www.nature.com/scientificreports/four steps against each strain level.The software uses the Voigt-Reuss-Hill (VRH) approximation to determine the moduli values.The adopted simulation methodology in this study considered all possible variations in the strain levels (0.001% to 10%) that are generally observed in geotechnical engineering and, therefore, it is applicable to assess moduli values in the elastic and elasto-plastic regions of stress-strain behavior of clay-bearing strata.A three-dimensional view showing the clay-bearing strata structure with inter-aggregate and intra-aggregate pore spaces after CED and modulus simulations is shown in Fig. 4.
There could be directional arrangement of the clay particles during the compaction process.3D views of the repeated cells shown in Figs.3b and 4 reveal various nature of fabrics consisting of various directional arrangement of clay particles.The mechanical properties determined for these 3D fabrics are representative of the three-dimensional stress case and hence covers the effects of the anisotropy if there is any.

Nano-level hydro-mechanical model
As explained earlier, CED can be effectively used in predicting the volume change and mechanical behavior of clay-bearing strata due to its sensitivity to changes in CEC, exchangeable cations, cementation, density conditions, and water.Therefore, this study has used it as a state parameter to predict the swelling behavior and mechanical properties of clay-bearing strata at different stages of swelling.
The Forcite module of Material Studio software was used to determine the CED and the corresponding modulus values.The simulation permutations were repeated for the adopted cases with varying density, cementation, moisture conditions, and exchangeable cations at strain levels of 0.001% to 10%.The outcomes of the simulations are plotted as a mechanical behavior model in Fig. 5.
The mechanical behavior model has been developed using a molecular level simulations approach, but the resulting structure and fabric of clay-bearing strata modeled in the software consists of both nano and micro-level pores, which are typically present in any overconsolidated hard clay/clayey rock structure.This is evident from the three-dimensional (3D) structure created using periodic boundary conditions (Fig. 4).Periodic boundary conditions create multiple unit cells of the resulting structure by replicating the model in all three dimensions infinitely.The amount of large-sized pores is almost minimal in the case of heavily overconsolidated clay-bearing soils/rocks because of the high densities, while there could be relatively high percentages of macropores in low density strata.
Two (02) types of trends can be observed from the mechanical behavior model plot (Fig. 5), one at small strain levels (0.01% and 0.001%) and the other at large strain levels (0.03% to 10%).Small strains range is applicable to foundations subjected to dynamic loads and corresponding field geophysical tests (i.e., in-situ cross-hole and downhole tests etc.).The large strain range is related to foundations, retaining walls, and conventional soil testing approaches 47,48 .
Three zones are obvious at a large strain level (Fig. 5).Zone I is bordered by CED values of 584 J/cm 3 and 1450 J/cm 3 .In this zone, the modulus values demonstrate an increasing trend with the raise in CED.This is mainly attributed to the increase in CEC with no contribution from the cementation between the particles.There is a tendency of reduction in the modulus values with an increase in CED in Zone II, bounded by CED values of 1450 J/cm 3 and 2640 J/cm 3 .Theoretically, the increase in CED in Zone II is only possible due to either the reduction in moisture content or an increase in cementation due to presence of salts such as Calcite, Gypsum, etc.It has been observed in molecular simulations that the presence of cementation in clay-bearing strata yields higher CED values than the limiting value of 2640 J/cm 3 for Zone II.Therefore, it can be interpreted that the increase in CED in Zone II is primarily due to a reduction in moisture content (dry nature of soil).The reduction in moisture content will cause the crumbling of dry clays and, ultimately, loss of strength and stiffness.In Zone III (CED > 2640 J/cm 3 ), there is again a trend in the rise of the modulus values with the increase in CED.
The increase in CED in this zone is caused by an increase in cementation due to the binding salts, which have resulted in an increase in modulus values.It can also be observed from Fig. 5 that keeping constant CED, the modulus values decay with an increase in strain level from 0.03 to 10%.The above discussion explained the material behavior at large strain levels, which is absent at small strain levels (0.01% and lesser).The modulus values at the small strain levels tend to increase at different rates in all zones.The behavior at small strain levels is essentially linearly elastic in nature, while it is elasto-plastic at higher strain levels.The modulus value at small strain level increases with the increase in CED value.The increase in CED value is associated with the increase in density and cementation and/or decrease in the degree of saturation.A similar trend is also observed by Pintado X. et al. 49 , where the modulus values at low strain levels increase with the dry density and at lower degrees of saturation.The small strain levels apply to foundations subjected to dynamic loads and low-strain geophysical tests.
Using the mechanical behavior model plots, fundamental equations are developed comprising the relationship between CED and modulus values, E, at various strain levels where E is the secant modulus.The relationships for all the three zones of the plots for large and small strain levels are represented below: Equations ( 2) to (4) represent relationships for large strain levels.However, for small strain levels, Eqs. ( 5) to ( 7) can be used: The numbers in the subscripts in Eqs. ( 2) to (7) represent the respective zone of the plot whereas large and small strain levels are shown by L and S, respectively.CED represents the cohesive energy density in J/cm 3 , and (2)  www.nature.com/scientificreports/ɛ is the strain level.As the simulated model has been constructed considering the natural composition of claybearing strata using a combination of clay and non-clay minerals with pore spaces, the strain, ɛ in Eqs. ( 2) to (7), represents the total strain of the clay-bearing strata.A linear interpolation could be performed for strain levels between 0.03% and 0.01%.The CED can be calculated by Eqs. ( 8) to (12) obtained by Ahmed and Abduljauwad 7 in their swelling behavior model as follows: where IWC is the initial water content.The correction due to CEC (ΔCED CEC ) is then applied using the following expressions: For CEC > 90 meq/100 g: For CEC < 90 meq/100 g: Equation ( 8) should also be corrected for binding/cementing salts and exchangeable cations (ΔCED cat ) as follows: where C is Calcite; L is hydrated lime, G is Gypsum; KCl is potassium chloride fraction of Na-montmorillonite (smectite) content, and Ca is exchangeable calcium cation; Mg is exchangeable magnesium cation; P is Palygorskite; K is potassium exchangeable cation fraction of the total cations.
For CECs other than 54 meq/100 g, Eq. ( 8) is normalized as: Furthermore, a correction for the confinement effect, p c, is also suggested in this study for the modulus values to account for the field stress conditions.This correction is given as: The CEC of clay-bearing soils/rocks is controlled by the relative amount of swelling-clay minerals in the sample.Therefore, in addition to the correction for the confinement effect to get the final modified modulus value (E m ), CEDm is also modified by the ratio of swelling-clay minerals to get the final CED value of the specimen for the mechanical behavior model.

Up-scaling of the model
Based on the level of compaction or consolidation, natural or engineered clay-bearing strata contain varying proportion of micro to macropores.Percentages of macropores are almost minimal in the case of heavily overconsolidated/compacted clay-bearing soils/rocks owing to their high densities, while there could be relatively high percentages of macropores in low density strata.Therefore, the models formulated using molecular-level simulations may require an upscaling for the cases where macropores are predominantly present.In this study, therefore, a correction, ΔIDD in the form of upscaling has been applied to account for the macropores in the clay-bearing strata.The factor, ΔIDD represents the difference in densities of clay-bearing soils/rocks at the molecular and macroscopic levels due to the presence of macropores and is found to be ranging from about 1.0 to 3.0 depending upon the proportion of macropores in the strata.The upscaling mechanism is described below in the form of equations.
The initial dry density (IDD) of the clay-bearing strata resulting from molecular-level simulations is related to normalized CED value, CEDm, through the following correlations 7 : The factor, ΔIDD is defined as follows: where IDD CEDm is the initial dry density from molecular-level simulations, and IDD specimen is the initial dry density of strata.The modulus values obtained from the mechanical behavior model are reduced by the factor ΔIDD to account for macropores in the strata and to get final modulus values.

Model validation
The model validation is mainly comprised of a testing program consisting of basic characterization tests, mineralogical study, and deformation properties, followed by the determination of CED and modulus values from the test results.Model validation has been carried out in two steps: (i) E-moduli validation at various strain levels from field and laboratory tests performed in this study and from an existing database, and (ii) validation through complete stress-strain plots.

E-moduli validation from current study
The mechanical behavior model is validated through laboratory testing on natural swelling soil samples collected from clay-bearing strata and reconstituted laboratory specimens (control specimens) utilizing several proportions of soil constituents.Many soil constituents were selected to identify the utmost feasible composition of the natural clay-bearing strata.Thus, individual soil constituents were attained from various known sources to prepare control samples.Natural samples from clay-bearing strata were acquired from the Qatif and Hofuf areas of Saudi Arabia, known for their shallow subsurface heavily over consolidated clay/claystone deposits [50][51][52][53] .The large intact lumps and blocks of heavily over consolidated clay/claystone samples were obtained from the existing test pits and were wrapped with plastic sheets to maintain the natural moisture content of these samples.Laboratory compacted bentonite samples were prepared from the Bentonite collected from a known source in Saudi Arabia (Kanoo Est., KSA).Bentonite was taken out of bags, oven-dried for 24 h, and then put into air-tight containers for further sampling and testing.The control specimens were then prepared by static compaction in a compression machine.
Calcium carbonate was collected from Techno Pharmchem, India, and Gypsum samples were taken from Phosphate Plant, RIC, KSA.Sand samples were collected from the dunes in Jubail, KSA (Table 1) and mixed with Bentonite to formulate the control specimens with different CED.
Basic characterization tests consisting of the grain size distribution (ASTM D6913), moisture content determination (ASTM D2216), and Atterberg limits (ASTM D4318) were performed on field-collected swelling clays and control specimens.Cation Exchange Capacity (CEC) tests were also carried out on swelling soils using Rayment & Higginson (1992) Method 15A2 54 .The results of basic characterization tests are summarized in Table 2.
X-ray diffraction (XRD) tests were performed using an Ultima IV X-Ray Diffractometer with copper radiation generated at 40 kV and 40 mA at the Research Institute (RI) in KFUPM.In this study, both quantitative and qualitative evaluations were performed to determine the type and amount of clay and non-clay minerals present in the samples.Randomly orientated samples (powdered specimens) were prepared by drying the sample in the oven and pressing it against the metal surface.Relative percentages of the different clay and non-clay minerals in natural clay-bearing strata of Hofuf and Qatif area and Bentonite are presented in Table 3.
To study the stress-strain response for the possible changes in fabric and structure, moisture-density relationships were established for the control specimens using the modified Proctor test procedures as per ASTM D1557.The results of the modified Proctor test are provided in Fig. 6.The control samples with recognized amounts of clay and non-clay minerals were arranged by static compaction in a compression machine at the moisture content and density on the "dry side of optimum," "optimum," and "wet side of optimum" using the results of the modified Proctor test to create different forms of fabric and structure.The undisturbed natural samples from clay-bearing strata were tested at the natural moisture contents and densities.The testing matrix for the experimental laboratory program carried out in this study is summarized in Table 4.These control and natural undisturbed specimens were then subjected to Unconsolidated Undrained (UU) Triaxial Compression tests as per ASTM D2850 to determine the stress-deformation properties.The stress-strain plots of UU Triaxial www.nature.com/scientificreports/Compression tests for the controlled specimens and natural samples at a cell pressure of 300 kPa (the stress conditions that the material will face in the field under most project conditions) are shown in Fig. 7. Crosshole seismic tests were performed as per ASTM D4428 using Olson CS Equipment on clay-bearing strata of the Hofuf area to validate the mechanical behavior model at low strain levels of 0.01% and lesser.In this test, three boreholes (spaced 3.0 m from center to center) were drilled in a straight line.The seismic waves were generated in one of the boreholes (called a "shot hole"), and its horizontal travel time was recorded in the other two boreholes (referred to as "receiver holes").The holes were drilled and prepared for the tests before the testing by placing PVC pipes and grouting the annular space.The field procedure consisted of the generation of seismic waves in the shot hole and recording the time of arrival of waves in the receiver holes.The readings were recorded at 1.50 m, 2.0 m, 6.50 m, 7.0 m, and 7.50 m in swelling clay/claystone layers at two different locations.www.nature.com/scientificreports/Based on the velocity of propagation of the waves (shear and compression) in the medium, the modulus values of the sub-surface layers were determined.A summary of the results of cross-hole tests is tabulated in Table 5.

E-moduli validation from existing literature
The molecular-level simulation results were validated from the existing literature using data from two (02) swelling clays (MX-80 and FEBEX Bentonite).The FEBEX and MX-80 are the bentonites used as Engineered Barrier systems in Finland, Sweden, and Spain.The mineralogical composition and material properties of both FEBEX (calcium bentonite) and MX-80 (sodium bentonite) samples have been presented in Tables 6 and 7.The information about MX-80 is taken from Kiviranta and Kumpulainen 55 and Kiviranta et al. 56 , whereas the information about FEBEX properties is collected from Enresa 57 .www.nature.com/scientificreports/Data on the mechanical properties of MX-80 and FEBEX bentonites were collected from Pintado X. et al. 49 .The resonant column tests on these bentonites were performed on different dry densities and degrees of saturation, and shear modulus G max or G 0 values were determined (Tables 8 and 9).

Stress-strain behavior model validation
The model developed in this study can be applied to predict stress-strain curves for clay-bearing soils/rocks.It can be used to predict the complete stress-strain behavior including elastic, elastoplastic, peak strength, and post peak behavior.Since the model is developed from the realistic simulations of the interaction between clay and non-clay minerals, pore water, and dissolved salts, it results in the stress-strain plots which constitute by default the elastic and elasto-plastic behaviors.
The following step-wise procedure is adopted to generate the stress-strain curves from the mechanical behavior model:  www.nature.com/scientificreports/ Step 1 The nano-level model parameter, cohesive energy density (CED), at initial conditions of moisture content and dry density is determined using Eq. ( 8).
Step 2 Eqs. ( 5) to ( 7) are used to find the modulus value for the initial sample conditions at low strain level of 0.001%.
Step 4 During the sample shearing stage, the CED and density of the sample change in response to various strain levels.This change in CED and density for each strain increment can be determined by a factor equal to 10 -4 × (E 2 -E 1 ).The E-value for the next step is then determined against new CED and density values after each strain increment.
Step 5 The deviator stress can then be obtained from E-values after each strain increment using Eq. ( 18), and results can be plotted as stress-strain curves.
where σ 1 and σ 3 are the major and minor principal stresses, ν is Poisson's ratio, and ɛ is the strain.

Hydro-mechanical behavior
Based on the mineralogical composition, CEC, initial moisture content, and density values, CED values were determined for the samples used in this study and also for MX-80 and FEBEX bentonites using the relationships presented in Eqs. ( 8) to (12) and the results are provided in Table 10.The modulus (E i ) values for the natural undisturbed and control specimens of clay-bearing strata at any given strain level were estimated using stress-strain plots of UU triaxial compression tests (Fig. 7) using the hyperbolic function 58 59 was used to estimate the variation of G and the corresponding E with the strain level.
In Eq. ( 19), γ 0.7 corresponds to the strain at which the modulus G is 0.722G o .The following relationship, given by Pintado et al. 49 , was used for γ 0.7 and shear modulus at a small strain.
The modulus values were determined for all the tested specimens of natural and compacted clays, and the results are summarized in Table 11.Table 11 also includes the modulus values obtained from crosshole seismic tests, represented by sample IDs H-2 to H-6, to validate the mechanical behavior model at low strain levels.The estimated modulus values for two (02) swelling clays (MX-80 and FEBEX bentonite) against varying CED values are also presented in Table 12.The calculated modulus values were then plotted on the mechanical behavior model plots developed in this study (Fig. 8).
Based on these results, the CED values range from 1130 to 9543 J/cm 3 , with moisture contents varying from 25.0 to 50.6% for both natural and control samples.The test results presented in Table 10 dictate that the final CED values of the tested samples vary from 633 to 4103 J/cm 3 , and the modulus values range from 0.021 to 1.0 GPa at strain levels varying from 0.10% to 10.0%.The modulus values from low-strain crosshole tests were found to be 5.77 to 6.09 GPa against the final CED values of 3995 to 4137 J/cm 3 .
As discussed earlier, three zones could be identified from the mechanical behavior model plot, (i) Zone I between a CED value of 584 J/cm 3 and 1450 J/cm 3 in which the modulus values increase with an increase in CED, (ii) Zone II from a CED value of 1450 to 2640 J/cm 3 where there is a tendency of reduction in the modulus values with the increase in CED, and (iii) Zone III having CED of more than 2640 J/cm 3 and showing an increasing trend in the modulus values with the increase in CED due to cementation.8 that all the tested bentonite samples fall in Zone I.The final CED values of the tested bentonite samples range from 791 to 931 J/cm 3 .The modulus value at a 1.0% strain level and a CED of 791 J/cm 3 is estimated as 0.106 GPa, whereas it is 0.236 GPa at a CED of 931 J/cm 3 .This confirms an increasing trend in the modulus values in Zone I with an increase in CED.The FEBEX bentonite samples are lying in Zone II of the mechanical model.The decrease in modulus values is observed in FEBEX bentonite samples with an increase in CED, which can be very well correlated with the lower moisture content of 4.70%, 3.70%, and 3.70% in samples F-3, F-4, and F-7, respectively.The undisturbed natural samples of swelling clays of Hofuf and Qatif areas fall in Zone III.The mineralogical composition of Hofuf and Qatif clays shows the presence of Gypsum and Calcite (Table 6) in the samples.The CED values of Hofuf and Qatif clays are calculated as 4103 J/cm 3 and 3167 J/cm 3 , respectively.The higher CED values in Zone III are caused by the increase in cementation due to binding salts, which increases the stiffness values with CED.The modulus values of Hofuf and Qatif clays are estimated to be 0.891 GPa and 0.387 GPa, respectively, at a 1.0% strain level, validating an increasing trend of the modulus values with CED in Zone III.
A comparison is also provided in Tables 11 and 12 between the laboratory and field-determined modulus values from crosshole tests and the ones predicted by the mechanical behavior model through Eqs. ( 2) to (7).It can be seen from the results that the modulus values of the tested samples at strain levels of 0.001 to 10% are in very close agreement with the model predictions.However, a major discrepancy was seen for the bentonite samples with Calcite as a cementing agent, which can be ascribed to the relatively low solubility of Calcite in distilled water.www.nature.com/scientificreports/Typically, foundations of civil engineering structures are placed in partially saturated soil layers.The moisture levels of the subsurface strata often increase after construction due to the reduced evaporation, percolation of surface water(s), and the potential of leakage from utilities.An increase in moisture level in the clay-bearing strata below foundations results in the swelling strains and associated reduction in the CED value depending upon the exchangeable cations, CECs, and binding/cementing agents present in the strata.In response to the swelling strains, the modulus value also decreases as a function of the difference between initial and final densities with respect to changes in CED value.This reduction in modulus value causes additional compressive strains under the constant superimposed loads from the structure.The net result of these coupled swelling and compressive strains defines the overall hydro-mechanical response of clay-bearing strata under the foundations.The overall hydro-mechanical response is formulated through the coupling of the mechanical behavior model developed in this study with the existing nano-level swell behavior model by Ahmed and Abduljauwad 7 .The swell potential can be determined using the swell behavior model 7 from the following Equation: The initial dry densities (IDD) are calculated using Eq. ( 13) to Eq. ( 15).The final dry densities (FDD) are related to CEDm through the following correlations 7 : where FDD is in units of g/cm 3 , CED is in J/cm 3 , and P, seating/confining pressure, is in kPa.
Incremental swelling can be estimated on each incremental intake of water content using the slope = (IDD-FDD)/(FWC-IWC).Final water content can also be assessed from the swelling behavior model as: where FWC, final water content, is in the units of percent and CED vw is in J/cm 3 and can be calculated from Eq. ( 27) given by Ahmed and Abduljauwad 7 .
An illustration of the coupled hydro-mechanical behavior is made using a general case of a shallow foundation placed in a partially-saturated, clay-bearing strata (Bentonite in this case) with initial conditions of IWC = 25% and IDD = 1.327 g/cm 3 .The saturation level in the clay-bearing strata is considered to increase from an initial water content of 25% to a final water content value of 35% after construction.The swell potential calculated through the swelling behavior model in ten equal increments of moisture content is presented in Table 13.A typical characteristic strain of 1.0% is considered for the modulus determination through Eq. (2) using the mechanical behavior model, and values are provided in Table 13 corresponding to each increment in water content.The results are also presented in Fig. 9, showing a decrease in CED values upon swelling and the corresponding reduction in the modulus values.The reduction in the modulus values from 0.206 GPa to 0.073 GPa against 5.13% swelling will cause additional compressive strains under the constant load from the structure.The foundation will finally be subjected to the net result of these coupled swelling and compressive strains.
The validation of coupled behavior through experimental data is also presented in Fig. 9.The UU Triaxial Compression and swelling-potential tests (ASTM D4546) were performed on the samples prepared at an initial moisture content of 25% and then allowed to swell by adding water in increments of 1% to reach a final moisture content of 35%.The predictions of swelling percent and modulus values are in well agreement with the values obtained from test results.The swelling behavior model proposed by Ahmed and Abduljauwad 7 has been verified using several laboratory control samples and the undisturbed samples acquired from the field in various studies carried out by Ahmed and Abduljauwad 7,60,61 .This model has also been successfully used by Abduljauwad and Ahmed 62 for the prediction of swelling of field clay deposits.

Stress-strain behavior prediction
The model has been successfully applied to predict stress-strain curves for swelling soils.Figure 10 shows the stress-strain curves predicted from the model compared with the experimental curves generated from UU Triaxial Compression Tests on both natural and controlled samples.The stress-strain curves predicted from the mechanical behavior model captured successfully the main features of the stress-strain response of clay-bearing strata, including the initial elastic response, plasticity, peak strength, and post peak behavior.The experimental curves were obtained for three (03) replicated specimens of each set of conditions to observe the variation in the stress-strain behavior among the tested specimens.The tests on reconstituted replicated specimens have resulted in similar stress-strain plots due to the controlled sample parameters.The natural clay/claystone samples, however, have shown little variation in the peak stress and the corresponding strain values.This variation in peak stress is associated with the natural variation of several minerals quantified in the field samples of the Qatif clay/ claystone.The variation in mineralogy is presented by upper and lower bound model predictions against the maximum and minimum gypsum content of 12% and 9% in the field samples, respectively.

Conclusions
In molecular-level simulations, a matrix of clay-bearing strata constituting both clay and non-clay particles, pore water, and dissolved salts was created to model interactions among these particles.The results of molecular-level simulations were compiled to develop a hydro-mechanical model that can be utilized to predict the  swelling-mechanical behavior of the clay-bearing strata.The model can be applied to all possible combinations of clay and non-clay minerals with cementing agents.
The study further concludes that the coupled hydro-mechanical response under given loading conditions can be well assessed through the coupling of the mechanical behavior model developed in this study with the existing nano-level swell behavior model by Ahmed and Abduljauwad 7 .
The proposed model can be effectively used to predict the overall material behavior of clay-bearing strata, covering all possible variations in the structure and fabric of clay-bearing soils/rocks having nano, micro and macropores.The model is characterized by a set of few parameters only, including cation exchange capacity, density, moisture content, and cohesive energy density, along with field stress conditions.These parameters can be determined from the basic characterization tests on the representative samples with great accuracy.
The predictions of stiffness moduli and stress-strain plots from the proposed model have been found to be in well agreement with the results obtained from the macro-level tests on the field and laboratory samples.Since the model is developed from the realistic simulations of the interaction between clay and non-clay minerals, pore water, and dissolved salts, it results in the stress-strain plots which constitute the elastic and elasto-plastic behaviors.The developed model has also been successfully validated from the existing literature using the data from two (02) swelling clays (MX-80 and FEBEX Bentonite).This suggests that the model can be used in all fields, such as civil engineering, agriculture, the petroleum industry, and waste management.

Figure 2 .
Figure 2. (a) Water sorbed particles of Na-montmorillonite, (b) A close up view of Na-Montmorillonite particle after water and cations sorption.

Figure 3 .
Figure 3. (a) Cubic unit periodic cell representation of four water sorbed crystallites, (b) Dispersed fabric at 0.001 GPa pressure and 10% moisture content.

Figure 4 .
Figure 4. Three-dimensional view of multiple unit cells after CED and modulus simulations of Na-montmorillonite with non-clay particles (water content = 10%, strain level = 1.0%).

Figure 5 .
Figure 5. Mechanical behavior model plot for swelling soils.

Figure 9 .
Figure 9. Coupled hydro-mechanical behavior and its validation.

Figure 10 .
Figure 10.Model predictions of stress-strain plots; (a) Bentonite samples on dry side of optimum, (b) Bentonite samples mixed with 20% sand at optimum moisture content, (c) Bentonite samples on wet side of optimum, (d) Clay/claystone samples of Qatif at natural moisture content.

Table 1 .
Gradation analysis of sand used in this study.

Table 3 .
Relative amount of minerals in clay-bearing soil samples.

Table 4 .
Testing matrix for laboratory program carried out in this study.

Table 6 .
Relative amount of minerals in FEBEX and MX-80.

Table 11 .
Modulus determination from test results and comparison with mechanical model predictions.E i , modulus value from field and laboratory tests; E m , final modified modulus value after corrections and upscaling.It can be observed from Fig.

Table 12 .
Modulus determination for FEBEX and MX-80 and comparison with mechanical model predictions.