Macroscopic Behavior and Microscopic Structure Evolution of Marine Clay in One-Dimensional Compression Revealed by Discrete Element Simulation

: Marine clay has been attracting in-depth research on its mechanical behavior and internal structure evolution, which are crucial to marine infrastructure safety. In the formation process of marine clay, including the sedimentation and consolidation stages, the compression behavior and internal structure evolution are highly dependent on the pore water salinity. Discrete element method (DEM) simulation is a powerful tool to study the microscopic mechanics behind the complicated macroscopic mechanical behavior of marine clay. In this study, a DEM simulation scheme is systematically proposed to numerically study the macroscopic beahvior and microscopic structure evolution of marine clay in one-dimensional compression that mimics the marine clay formation process. First, the proposed calculation scheme for double layer repulsive interaction and van der Waals interaction is introduced. Then, the developed DEM simulation scheme is validated by satisfactorily reproducing the experimentally observed one-dimensional compression curves and internal structure transition from an edge-to-edge/edge-to-face flocculated structure to a face-to-face dispersed structure. Finally, evolutions of coordinate number and fabric anisotropy are quantitatively evaluated in the microscopic view. The noticeable effects of ion concentration on the internal structure evlotion and mechanical behavior of marine clay have been examined and discussed.


Introduction
Marine clays are widely distributed along the coastal lines where large amounts of infrastructure are in service or under construction.A complete understanding of the physical and mechanical behavior of marine clay is crucial to infrastructure safety under varying conditions, such as pore water salinity change, depositional process alternation, wave loading, etc.The complicated behavior of marine clay is mainly attributed to its initial microscopic structure (fabric) and the evolution of its structure upon chemical, thermal, or mechanical disturbance.Fabric is found to be the dominant mechanism behind the distinct behavior between undisturbed and reconstituted clay samples [1,2].The microscopic structure has been found to be crucial to the development of a structure-informed constitutive model of marine clay [3].The pore electrolyte-clay interaction and its effects on the fabric and macroscopic behavior of marine clay have been under investigation for decades.Salinity content can increase shear modulus and alter the damping ratio of marine clay, thus imposing an influence on the site response under seismic loading [4].The development of extra quick clay in the Ariake Bay area is a result of low salinity (0.05~1.08 g/L), with a sensitivity from 50 to 1000 [5].The variations of geochemistry during deposition and after the deposition, play crucial roles in determining the fabric, bonding, and thus the mechanical behavior of marine clay [5].The data interpretation of an insitu cone penetration test needs to properly consider the pore-water salinity effects on soil type identification [6].Detailed investigation of the electrolyte-clay interaction leads to a modification of the well-established effective stress concept [7].However, most of previous research focuses on laboratory experiments, microscopic structure observations via scanning electron microscope, and field observations of the physical-mechanical behavior, while the internal structure evolution is simply qualitatively evaluated.The formation history of marine clay and its effects on clay structure are not fully examined due to the scarcity of deposition history records.
In recent decades, the numerical simulation of clay focusing on its fabric formation and evolution has gained great popularity.This is well supported by the rapid development of discrete element method (DEM) simulation techniques.In DEM, clay is viewed as an assembly of discrete objects and the fabric evolution can be captured by the displacements of and interactions between these objects.The discrete objects are line segments in two-dimensional analyses [8][9][10] and flat cuboids or flat ellipsoids in three-dimensional analyses [11][12][13][14][15][16].DEM simulation has been proven capable of capturing the major mechanical behavior of clay.However, the real clay platelets featured a complicated geometry, which has not been reproduced in previous DEM simulations.Apart from the geometry issue, the complicated interactions between clay platelets need to be captured for reliable DEM simulation.Different from coarse-grained soil, there exist three interactions between two clay platelets: a double layer repulsive interaction, van der Waals attractive interaction, and mechanical repulsive interaction.A systematic scheme to calculate the double layer repulsive interaction and mechanical repulsive interaction has been proposed in [11].The Gay-Berne potential function was used in [13] to describe the van der Waals interaction, but it is only applicable to clay platelets of flat ellipsoid shape.In [11], the van der Waals interaction between two finite flat cuboids was approximated by that between a finite flat cuboid and an infinitely extended flat cuboid.In [15], clay platelet is approximated by clump of balls arranged by row and column.So far, there has not been a complete DEM simulation scheme to capture the complicated clay platelet shape and the associated inter-platelet forces, which largely limits the study of marine clay in terms of its mechanical behavior and the associated microscopic mechanisms.
In this study, a complete scheme is proposed to calculate the interaction between clay platelets of irregular shapes.Then, the scheme is implemented into DEM to simulate the formation (sedimentation and consolidation) process of a marine clay by considering a one-dimensional compression process.The mechanical behavior and internal structure evolution during the compression are examined in detail.The implications of the simulation results are discussed.

Methods
Marine clay is viewed as an assembly of clay platelets in this study.Each clay platelet is assumed to be extruded based on a planar shaped polygon with multiple edges.The extrusion distance is the thickness of the clay platelet.This assumption implies that the numerical clay in this study is applicable to clay minerals, such as kaolinite and illite, that do not present a significant spatially curving shape.The interactions between two clay platelets include double layer repulsive interaction, van der Waals attractive interaction, and mechanical repulsive interaction.The former two interactions are closely related to the shapes and relative alignment of two neighboring clay platelets of interest.Although the pore water is not explicitly simulated in the DEM model, the pore water effects are implicitly considered in the interaction law between clay platelets as outlined below.

Calculation Scheme of Double Layer Repulsive Interaction
Figure 1 presents an arbitrary alignment of two clay platelets.Without loss of generality, the planar shapes of clay platelet 1 and clay platelet 2 are assumed to have eight and six edges, respectively.As shown in Figure 1a, the normal directions of the two platelets are denoted by n1 and n2.A mid-plane is defined as the plane that bisects the acute angle made by the two platelets.A local coordinate system is defined with the normal direction of the mid-plane as z axis, and the direction of n2 × n1 as y axis, while the origin can be put anywhere in the mid-plane.Since the clay platelet surfaces are negatively charged, then in Figure 1a, the lower facet of clay platelet 1 and the upper facet of clay platelet 2 constitute a pair of inclinedly aligned, negatively charged facets.The electric interaction of the two clay platelets can be attributed to the interaction of the two facets.In Figure 1b, projecting the two facets on the mid-plane gives two intersecting polygons, as outlined by the thin dashed line segments.The intersection area is outlined by solid line segments.A side-view along the positive-y direction is presented in Figure 1c for easy understanding.The intersection area is the range of effective double layer repulsive interaction, while other parts of the two facets' projections contribute little to the interaction and are omitted.
Following the method in [11], the repulsive force can be calculated in a simplified way as follows.First, as presented in Figure 1b,c, the intersection area is divided into bands along the x-axis direction with a band width of w.The two facets (lower facet of clay platelet 1 and upper facet of clay platelet 2) are also divided in this way, as depicted Mid-plane in Figure 1c.The interaction between the two inclinedly aligned facets is now simplified as the summation of the interactions between a series of oppositely positioned band pairs with respect to the mid-plane, one from clay platelet 1 and the other from clay platelet 2 in a pair.Take band i in Figure 1b as example.The band has a length of li in the y-axis direction.The center distance of the two corresponding bands from the upper and lower platelet is hi in the z-axis direction, as shown in Figure 1c.According to the work in [11], the double layer repulsive force between the two inclinedly aligned bands, Fe,i, can be approximated by that between two parallelly aligned bands with the same distance hi, that is where n is ion concentration in the pore water, k = 1.38 × 10 −23 J/K is the Boltzmann constant, T is Kelvin temperature, and is the normalized mid-plane electric potential of two parallel plates at a distance of hi.In Equation (1), Fe,i's direction is normal to the mid-plane.The normalized mid-plane electric potential can be obtained by numerically solving the Poisson-Boltzmann equation and the solutions are presented in Figure 2. The normalized mid-plane electric potential is defined as where v is the ion valence, e * = 1.6 × 10 −19 C is the charge of an electron, and  is the mid- plane electric potential.In Figure 2, the normalized distance is defined as where  is the relative dielectric constant of pore water and is the dielectric constant in vacuum condition.
Figure 2 shows that the normalized mid-plane electric potential decreases non-linearly with the normalized distance.The curves also depend on the normalized surface electric potential of the clay platelet defined as ( )  being the surface electric potential.
By the summation operation, the double layer repulsive force, Fe,1, and moment, Me,1, applied on clay platelet 1 by clay platelet 2 can be expressed in the local coordinate system as e 1 e 1 = 0 0 e 1 e e 1 1 = 0 where Nx is the number of divided bands in Figure 1b, and xi and yi are the coordinates of band centers in the local coordinate system.The above derivations apply to any convex planar shape of clay platelet.It is noted that the core idea of the above simplified calculation scheme of a double layer repulsive interaction is to replace the inclinedly aligned facets with a series of staircase-like, parallelly aligned bands.The finer division of the bands leads to the higher accuracy.To address this issue, three clay platelets are used as illustrative examples.As shown in Figure 3, without loss of generality, the planar shapes of clay platelets A, B, and C are assumed to have four, seven, and five edges, respectively.The planar extent of the platelets is between 0.   When the inclination angle θ is changed, the net distance of the two clay platelets, dn, is kept constant at 1 nm. Figure 4 shows that the double layer repulsive force decays very quickly as the inclination angle increases.When θ is greater than 7°, the repulsive force is lower than 1‰ of that when θ = 0° and it can be omitted.This conclusion is intrinsically consistent with that in [16] for clay platelets with rectangular planar shape.Figure 5 presents the variations of double layer repulsive force with the relative inclination angle and the net distance between clay platelets A and B under various surface electric potentials.The division number Nx takes 8 here.Figure 5a shows that the repulsive force can be omitted when the inclination angle exceeds 10°, similar to conclusions from Figure 4. Figure 5b shows that the repulsive force decays very quickly as the net distance dn increases.When dn is greater than 10 nm, the repulsive force is lower than 1‰ of that when dn = 1 nm and can be omitted.

Calculation Scheme of Van Der Waals Attractive Interaction
The van der Waals interaction between clay platelet 1 and clay platelet 2 is the summation of all interactions of molecule pairs respectively belonging to the two platelets.The van der Waals attractive energy, u, between two molecules are first given by in [17] as where B is the London constant, r is a vector connecting the centroids of two molecules, and thus r is the centroid distance of the two molecules.
The van der Waals attractive force between two molecules can be obtained by differentiation of Equation ( 6) with respect to r as In Figure 6a, the outlined solid volume of clay platelet 1 corresponds to the intersection area in the mid-plane (see Figure 1) and is denoted as clay platelet 1′.The same applies to clay platelet 2 and clay platelet 2′.Since the van der Waals attractive force decays with the distance very fast (to a power of −7), the intersecting clay platelets 1′ and 2′ would make the most contribution to the interaction while the other non-intersecting parts' contribution can be omitted.Similar to the band division in Figure 1b, the clay platelet 1′ is divided into bands and the van der Waals interaction between clay platelet 2′ and every divided band is calculated.In Figure 6b, taking band i of clay platelet 1′ as example, its van der Waals interaction with clay platelet 2′ can be approximated by the interaction with an in-plane infinitely extended platelet with the same thickness as platelet 2′, which has an analytical solution.The extended volume has little influence on the interaction since it is located far away from band i. Integrating Equation (7) over the volume of an in-plane infinitely large platelet with a thickness of t gives its van de Waals attractive force with a point as ( ) ( ) where h is the distance of the point to the upper surface of the platelet, and ρ1 is the molecule density in the platelet.Integrating Equation ( 8) over the volume of band i gives the van der Waals attractive force applied on band i as where A = 2 12 B    is the Hamaker constant, q is the acute angle made by the two clay platelets, and d1-d4 are the distances of the four vertices of band i to the upper surface of clay platelet 2′, as illustrated in Figure 6b.
Then, the van der Waals attractive force, Fv,1, and moment, Mv,1, applied on clay platelet 1 by platelet 2 can be determined by Similar to the results in Figure 4, although not given here for conciseness, it is found that Nx = 8 is fine enough to divide the clay platelet in Figure 6 for the accurate calculation of van der Waals attractive force.
When two clay platelets are aligned with the acute angle made by them being large, say greater than 70°, the calculation scheme presented above can be replaced by a simpler scheme as discussed below.As presented in Figure 7, the z-axis of the local coordinate system is aligned with the normal direction of clay platelet 2. Denote the closest point on clay platelet 1 to clay platelet 2 by E and EA is the extrusion edge of clay platelet 1.Among the two lateral facets sharing edge EA, i.e., facet EABF and facet EADH, facet EABF makes a smaller angle with plane oxy than that facet EADH makes.A thin cuboid, EABF-IJKL, is constructed by taking the facet EABF as one of its lateral facets and extending in the direction away from clay platelet 2. The extending distance takes the average edge length of the planar polygon of clay platelet 1.In fact, this average edge length is an ideal choice since the thus constructed cuboid would cover a large portion of the volume of clay platelet 1.A slight deviation from this extending distance has a very limited influence since the van der Waals interaction is mostly contributed by clay platelet 1′s molecules located near the clay platelet 2.An infinitely large clay platelet 2′ is constructed the same way as in Figure 6b.With the above considerations, the van der Waals interaction between clay platelet 1 and platelet 2 can be approximated by the interaction between the extended thin cuboid and the infinitely large clay platelet 2′, which has an analytical solution.
Integrating Equation ( 8) over the volume of the thin cuboid gives the van der Waals attractive force applied on clay platelet 1 as

h h h t h t h t h t A F F F F F n n n h h h h h t h t h t h t
where hA-hL are the distances from the thin cuboid's vertices to the upper facet of clay platelet 2, and n ， , and 3 EH n ， are the z components of the unit vectors along directions EA , EF , and EI .Similarly, the moment can be integrated.Figure 8 presents the variations of van der Waals attractive force with the relative inclination angle and the net distance between clay platelets A and B under various Hamaker constants.The division number Nx takes 8 here.In Figure 8a, when the inclination angle is less than 73°, Equation ( 9) is used, and otherwise Equation ( 12) is used.Here, 73° is just for illustration and other delimiting value around 70° can be used as well without changing the conclusion below.In Figure 8a, the smooth transition across 73° indicates that Equations ( 9) and ( 12) are both properly developed.Figure 8a shows that the van der Waals attractive force first quickly decreases as the inclination angle increases and then presents a slow decaying rate until the transitional angle about 45° is reached.After that, an increase of the attractive force is observed with the inclination angle and the increasing rate soars after the inclination angle exceeds Again, this conclusion is intrinsically consistent with that in [16] for clay platelets with rectangular planar shape.Figure 8a also shows that the van der Waals attractive force can be omitted when the inclination angle is between 10° and 80°, with the attractive force less than 1% of that in parallelly aligned platelets case.Figure 8b shows that the attractive force decays very quickly as the net distance dn increases.When dn is greater than 10 nm, the attractive force is lower than 1‰ of that when dn = 1 nm and can be omitted.

Total Interaction between Two Clay Platelets
In DEM simulations, a contact model is applied to describe the interaction between two discrete objects at the contact.The objects in three-dimensional DEM simulations can be balls, rigid blocks, ball clumps, walls, etc.In traditional DEM, two objects can interact only when they are in physical contact, i.e., dn = 0 nm.However, for the DEM simulation of clay, the total force/moment received by a clay platelet from surrounding platelets is a summation of traditional mechanical contact interaction (present when two objects being in physical touch) and the double layer repulsion and van der Waals attractive interactions (considered in certain net distance and relative inclination angle ranges as discussed above).
The double layer repulsion and van der Waals attractive interaction calculation scheme were implemented into DEM software PFC3D [18].The built-in linear rolling resistance contact model in PFC3D [18] is used as the physical touch interaction, where the normal, tangential, and rolling interactions can be expressed as where Fn and Fs are normal and tangential forces, respectively, Mr is the rolling resistance moment, Kn, Ks, Kr are the corresponding contact stiffnesses, un is the overlap, nc is the unit contact normal vector, Δ u s is the incremental tangential displacement, r Δq is the in- cremental rolling angle, μ is the frictional coefficient, μr is the rolling resistance coefficient, and R is an equivalent contact radius.
In Equations ( 13)-( 15), the contact stiffnesses are where E is the contact modulus and ξk is the normal-tangential stiffness ratio.Note that in Equations ( 9) and ( 12), when the net distance between two clay platelets is small, say less than 1 nm, the van der Waals attractive force would be very large.In fact, in net distance range smaller than 1 nm, the repulsive van der Waals interaction (not discussed above due to irrelevance) will start to dominate and the absorbed water molecule layers on clay platelet surfaces will start to overlap to generate repulsive force.The two issues are abstracted by setting a mechanical cut-off gap dm and when the net distance dn is less than dm, the following two aspects are considered: (1) van der Waals attractive and double layer repulsive interactions are kept constant, i.e., no change with decreasing dn; (2) the starting point of traditional mechanical contact interaction is shifted from dn = 0 to dn = dm [11].
Figure 9 visualizes the interactions between clay platelets A and C (see Figure 3), where the mechanical cut-off distance dm = 1.0 nm and the normal contact stiffness Kn = 60 N/m.In Figure 9, the net force is a summation of the double layer repulsive force, van der Waals attractive force and mechanical contact repulsive force, with positive value representing repulsion.It is interesting to note that the net force is repulsive when dn is less than 1.25 nm and a local maximum repulsion is observed at dn = 2.0 nm.The net force becomes attractive when dn is less than 1.25 nm and reaches a local maximum attraction at dn = dm (=1.0 nm here).When dn is less than dm, the mechanical contact repulsion starts to dominate.

Clay Platelet Templates in DEM Simulation
Without loss of generality, it is simply assumed that the planar shape of the simulated clay platelet is polygon with eight edges.As illustrated in Figure 10, an initial polygon is constructed by connecting eight points located on a circle and thus it has an aspect ratio close to 1.0.The initial polygon is then extended in the x or y direction to an aspect ratio of 1.1, 1.2, 1.3, and 1.4.The nine polygons are scaled to have the same area of 0.64 μm 2 and then extruded to clay platelet thickness of 0.04 μm.The thus formed nine polyhedrons together constitute nine clay platelet templates for DEM model construction.The clay platelet dimensions are typical for kaolin platelets.

DEM Simulation of Sedimentation and Consolidation Processses
The commercial software PFC3D [18] was used for the numerical simulation.In PFC3D, clay platelets are represented by rigid blocks generated from the above constructed nine templates.Hence, 500 clay platelets are generated with random positions and directions in a slender box composed of six boundary walls, as shown in Figure 11a.The corresponding void ratio e is 21.0, avoiding any overlap between platelets.These clay platelets are randomly drawn from the nine templates.The solid density is set to 1650 kg/m 3 (considering the buoyancy).Note that 500 platelets are enough to capture the behavior focused on in this study, and the related computation cost is acceptable, which will be discussed later.The model parameters are listed in Table 1.A Hamaker constant of 1.0 × 10 −19 J is typical for a kaolinite-NaCl-water system.The surface electric potential may change in a wide range for clay minerals and a typical value of 200 mV is applied here.Numerous experimental data indicate that pore water chemistry can significantly change the physical and mechanical behavior of clay.In our simulations, the ion concentration is changed to investigate such effects.Salt concentration in marine clay pores poses large variation, depending on the specific location [19].For example, near river mouth where a river meets the ocean, the sea water can be diluted to a large extent.In Table 1, the ion concentration of 0.64 mol/L is a typical up-bound value for most ocean area while the reduced concentrations can be representative of diluted cases.The physical contact model parameters are selected based on common DEM simulation practice.The clay platelets in the box are allowed to settle down under gravity to simulate the sedimentation process under static water condition.The settled clay platelet assembly is then subjected to one-dimensional compression to simulate the consolidation process due to the weight of subsequent sediments lying above.One-dimensional compression is numerically achieved by moving the top boundary wall downward.The moving rate is very low to allow full arrangement adjustment of clay platelets.

Observations in Sedimentation Process
Under the action of gravity, buoyancy, double layer repulsive interaction, and van der Waals interaction, the platelet assembly arrives at a loosely organized structure, as shown in Figure 11b, resembling the so-called flocculated structure of marine clay commonly encountered at the sea bottom.The salinity condition at the time of sedimentation influences the spatial arrangement of the clay platelets: the stable void ratios after sedimentation are 4.05, 3.25, and 2.62 for ion concentrations of 0.64, 0.16, and 0.02 mol/L, respectively.

Macroscopic Behavior in Consolidation Process
The obtained compression curves in e-logp plot are presented in Figure 12a, where p is the effective compression pressure.The one-dimensional compression responses resemble those of clays in experiments, with a linear trend when the compressive pressure exceeds 100 kPa.The ion concentration has a significant influence on the void ratio of clay under the same compression pressure.A higher ion concentration leads to a looser packing of clay platelets.For the studied range of ion concentration, a higher n leads to lower double layer repulsive force.Consequently, the van der Waals attractive force plays a more significant role to bond the platelets together to form a flocculated structure to resist the overall one-dimensional compression.This mechanism has been proposed and validated experimentally in the rich of literature, e.g., [20,21].The observation here is only valid for non-swelling kaolinic clays.For swelling clay such as bentonite, the microscopic mechanisms are different and are not to be discussed here.In [22], a void index Iv was introduced to normalize the intrinsic compression curve (ICL) of normally consolidated clay Iv is defined as The simulated one-dimensional compression curves in this study are replotted in the Iv-p plot with Equation ( 18) superimposed as the ICL curve in Figure 12b.It is found that the simulated data generally lie around the empirical ICL after p exceeds 100 kPa, validating the DEM simulation scheme.The deviation of data in the less than 100 kPa range could be attributed to the difference of sample preparation method between the lab experiments in [22] and the numerical simulation here.Such deviation is also observed in [1,21,23], among others, for marine kaolinite-rick clay.In lab experiments, the dry soil powder is mixed with water and stirred to a uniform state, while in numerical simulation, the sample is prepared by the sedimentation of clay platelets under gravity, thus forming a different fabric from that in experiments.This difference in initial fabric is gradually removed as the compression pressure increases.In all, the numerically simulated responses of marine clay are consistent with existing experimental results, thus validating the proposed DEM simulation scheme in this study.
Figure 13 presents the simulated spatial configurations of clay platelets at different compression stages for the case with n = 0.16 mol/L.It is well visualized that the clay platelets tend be aligned along the horizontal plane as the assembly becomes densely arranged due to one-dimensional compression.That is, a transition from an edge-toedge/edge-to-face flocculated structure to a face-to-face dispersed structure is observed.A quantitative evaluation of the microscopic structure evolution is presented below.

Microscopic Structure Evolution in Consolidation Process
First, the coordinate number of clay platelet assembly, CN, is analyzed, which is defined as where Np is the number of clay platelets and Nc is the number of interaction pairs between neighboring platelets.
Figure 14a presents the evolutions of e versus coordinate number CN during the simulated one-dimensional compression.The coordinate number increases in a highly nonlinear way as the void ratio e decreases, accompanied by approaching to an increasingly dense state.With the same e value, the assembly with higher ion concentration has a greater coordinate number, which can be attributed to the less double layer repulsion and the augmented flocculated structure.This is consistent with the conclusion drawn from experiments of kaolinite clay [24].Second, a fabric tensor Fij is defined to statistically describe the configuration of clay platelets in an assembly: where nk,i and nk,j are components of the unit normal vector of a clay platelet, and the indices i and j take values of 1, 2, and 3.
The fabric anisotropy can be assessed by fabric index With the same e value, the assembly with higher ion concentration exhibits lower anisotropy.This can be explained by the observation that less double layer repulsion due to higher ion concentration allows the van der Waals attraction to adhere clay platelets in a more stable structure to resist compressive deformation and platelet orientation adjustment.The well-known anisotropic behavior of marine clay in terms of hydraulic permeability, electrical conductivity, thermal deformation, modulus, and strength is closely associated with its internal anisotropic fabric.The relevance of fabric anisotropic evolution with void ratio and ion concentration revealed by the DEM simulations here can be utilized in future research to obtain a greater understanding of the physical, thermal, and mechanical behavior of marine clay.
The DEM simulation in this study is only applicable for the non-swelling type of clay, such as kaolinite clays.It will be extended to swelling type clays, such as bentonite, in future work.It should also be noted that there are various other platelet shapes in real kaolinic clay.The presented simulation results here should be conceived as an approximate and qualitative representation of the counterpart real clay.Compared with the cuboid in [11,12], disk in [13], and cuboid made of spheres in [15,16], the platelet shape in this manuscript is more flexible to replicate the geometry of clay platelets.Although the aspects of clay behavior focused on in [11][12][13]15,16] are well replicated, the platelet simulation scheme in this manuscript surely has a greater potential to reveal the effects of a rich range of factors, such as statistics of platelet shape, orientation, and alignment, as well as mineral combinations.
Due to the low computation efficiency, only 500 platelets were used in our numerical simulations.It is interesting to note that the one-dimensional compression behavior of marine clay and the effects of ion concentration were successfully replicated in this manuscript.Three further simulations with platelet numbers of 300, 1000, and 1500 were simulated under the condition of ion concentration of n = 0.16 mol/L, as presented in Figure 15.When the platelet number is low (say 300), the curve is noisy.When the platelet number increases from 500 to 1500, only a slight change can be observed.For the same computation expense reason, the clay platelet number was limited to 400 in [11,16], 500 in [15], 1000 in [13], and 1200 in [12], which were determined by trial and error.The aspects of clay behavior focused on in [11][12][13]15,16] are also well replicated.In conclusion, it seems to imply that 500-1000 platelets are enough for investigating the general trend of clay's macroscopic mechanical behavior under one-dimensional compression.However, if detailed quantifications of finer microscopic aspects are to be investigated, such as particle orientation and alignment, inter-platelet bond breakage in structured clay, or inter-platelet pore space evolution, a number greater than 500 should be required.

Figure 1 .
Figure 1.Schematics of the calculation scheme of double-layer repulsive interaction between two clay platelets: (a) three-dimensional alignment of two clay platelets, (b) projections of clay platelets' facets on the mid-plane, (c) projections of clay platelets on the oxz plane.

Figure 2 .
Figure 2. Relationships between normalized mid-plane electric potential and normalized distance.
figures.If not stated otherwise, the parameters are T = 293 K,

Figure 3 .
Figure 3. Planar geometries of three clay platelets in analyses.

Figure 4
Figure4presents the variations of double layer repulsive force with the inclination angle θ between clay platelets B and C under different band division number Nx.When the inclination angle θ is changed, the net distance of the two clay platelets, dn, is kept constant at 1 nm.Figure4shows that the double layer repulsive force decays very quickly as the inclination angle increases.When θ is greater than 7°, the repulsive force is lower than 1‰ of that when θ = 0° and it can be omitted.This conclusion is intrinsically consistent with that in[16] for clay platelets with rectangular planar shape.Figure4also indicates that when Nx exceeds 8, the repulsive force will only change quite slightly as Nx further increases.It can be concluded that Nx = 8 is fine enough to divide the bind for the accurate calculation of double layer repulsive force.

Figure 4 6 NFigure 4 .
Figure4presents the variations of double layer repulsive force with the inclination angle θ between clay platelets B and C under different band division number Nx.When the inclination angle θ is changed, the net distance of the two clay platelets, dn, is kept constant at 1 nm.Figure4shows that the double layer repulsive force decays very quickly as the inclination angle increases.When θ is greater than 7°, the repulsive force is lower than 1‰ of that when θ = 0° and it can be omitted.This conclusion is intrinsically consistent with that in[16] for clay platelets with rectangular planar shape.Figure4also indicates that when Nx exceeds 8, the repulsive force will only change quite slightly as Nx further increases.It can be concluded that Nx = 8 is fine enough to divide the bind for the accurate calculation of double layer repulsive force.

Figure 5 .
Figure 5. Variations of double-layer repulsive force with (a) relative inclination angle and (b) net distance between two clay platelets under different surface electric potentials.

Figure 6 .
Figure 6.Schematics of the calculation scheme of van der Waals attractive interaction between two clay platelets: (a) three-dimensional alignment of two clay platelets, (b) projections of clay platelets on the oxz plane.

Figure 7 .
Figure 7. Schematics of the calculation scheme of van der Waals attractive interaction between clay platelets with large relative inclination angle.

Figure 8 .
Figure 8. Variations of Van der Waals attractive force with (a) relative inclination angle and (b) net distance between two clay platelets under different Hamaker constants.

Figure 9 .
Figure 9. Normal contact response curves between clay platelets.

Figure 10 .
Figure 10.Clay platelet templates used in DEM simulations.

e
 are the void ratios at consolidation pressures of 100 kPa and 1000 kPa, respectively, andC C  is the compression index.Based on one-dimensional compression of 26 types of normally consolidated clays, Iv is fitted in[22

2 Figure 13 .
Figure 13.Numerical samples of clay in DEM simulation with different void ratios (n = 0.16 mol/L).

Figure 14 .
Figure 14.Evolutions of microscopic structures: (a) e versus coordinate number CN and (b) e versus fabric anisotropy index IF. p the fabric tensor.Figure14bpresents the evolutions of e versus fabric anisotropy index IF during the simulated one-dimensional compression.The fabric anisotropy increases in a slightly non-linear way as the void ratio e drops, because of the clay platelet orientation adjustment during compression.

Table 1 .
Contact model parameters in DEM simulation.