DEM study of fabric features governing undrained post-liquefaction shear deformation of sand

In an effort to study undrained post-liquefaction shear deformation of sand, the discrete element method (DEM) is adopted to conduct undrained cyclic biaxial compression simulations on granular assemblies consisting of 2D circular particles. The simulations are able to successfully reproduce the generation and eventual saturation of shear strain through the series of liquefaction states that the material experiences during cyclic loading after the initial liquefaction. DEM simulations with different deviatoric stress amplitudes and initial mean effective stresses on samples with different void ratios and loading histories are carried out to investigate the relationship between various mechanics- or fabric-related variables and post-liquefaction shear strain development. It is found that well-known metrics such as deviatoric stress amplitude, initial mean effective stress, void ratio, contact normal fabric anisotropy intensity, and coordination number, are not adequately correlated to the observed shear strain development and, therefore, could not possibly be used for its prediction. A new fabric entity, namely the Mean Neighboring Particle Distance (MNPD), is introduced to reflect the space arrangement of particles. It is found that the MNPD has an extremely strong and definitive relationship with the post-liquefaction shear strain development, showing MNPD’s potential role as a parameter governing post-liquefaction behavior of sand.


Introduction
Excessive deformation of liquefied sands is frequently observed in earthquakes and is one of the most damaging effects of soil liquefaction (e.g., [6,16,31,38]). The deformation of sand related to cyclic liquefaction can be reproduced in laboratories through undrained tests, and numerous studies have shown that significant shear deformation continues to be generated after ''initial liquefaction,'' i.e., in the ''post-liquefaction'' stage (e.g., [2,21,35,42]). In this work, ''initial liquefaction'' refers to the first occurrence of soil liquefaction (i.e., zero effective stress or excess pore pressure ratio of 100 %) during cyclic loading, following Seed and Lee's [28] definition. The conditions for triggering initial liquefaction have been the subject of extensive studies, and established criteria exist for both laboratory and field applications (e.g., [18,19,27,39]), which deal with ''pre-liquefaction'' (i.e., before initial liquefaction) behavior of sand. The work presented in this paper primarily focuses on sand behavior during continued post-liquefaction cyclic loading, which consists of alternating ''liquefaction states'' and ''non-liquefaction states.'' Undrained cyclic torsional laboratory tests (e.g., [40,42]) have produced a number of intriguing observations in the post-liquefaction stage: (1) The stress path of each load cycle follows almost exactly the same ''butterfly orbit,'' entering and exiting liquefaction two times per loading cycle (Fig. 1a). (2) Significant shear strain is generated within each liquefaction state after initial liquefaction (Fig. 1b). The amplitude of this strain is termed by Zhang and Wang [41] as the post-liquefaction shear strain at zero effective stress, denoted by c 0 . Figure 1c depicts the evolution of shear strain during the 11th load cycle of an undrained cyclic torsional test on Toyoura sand at a relative density D r of 70 % [40], where significant deformation occurs at the liquefaction state once during the first half load cycle in the positive loading direction and once during the second half cycle (denoted as cycle 11.5) in the negative direction. (3) Although the stress paths and the stressstrain relationship in the non-liquefaction states remain nearly identical among the loading cycles, c 0 generated within each liquefaction state increases as the cyclic loading continues (Fig. 1b, d). (4) The rate of the increase of c 0 gradually decreases until c 0 eventually saturates at a certain level c 0s (Fig. 1d).
Although sands' large post-liquefaction shear deformation and its progressive development have been observed in many undrained cyclic laboratory experiments (e.g., [5,17,35,42]), a widely accepted explanation for this phenomenon is still absent. Significant effort toward this end has been made in several constitutive studies through various assumptions associating c 0 with dilatancy and fabric history (e.g., [4,10,36,41]). However, as laboratory tests generally only provide macroscopic measurements of stress, strain, void ratio, pore pressure, etc., little progress has been made toward revealing the intrinsic microstructural and sand fabric evolution processes causing the accumulation and eventual saturation of c 0 at liquefaction. Understanding the aforementioned post-liquefaction behaviors of sand, particularly the microstructural (or fabric-related) mechanisms governing these behaviors, has significant implications for constitutive modeling of sand, and more importantly, it will push forward our overall understanding of granular material's mechanical behavior.
The current study attempts to solve a long-standing puzzle in soil mechanics concerning liquefaction: Sand in the liquefaction state has certain behaviors resembling those of fluids, with all particles ''semi-suspended'' without inter-particle effective stress, at least in the macroscopic scale. Intuitively, if liquefied sand is a simple fluid, all the history and states of solid particles' packing contacts (fabrics) would be ''erased'' upon entering liquefaction, and unbounded shear strain would develop without inducing effective stress. However, certain behavior of sand seems to imply that some memory of particle packing is preserved even as the material goes through the particle semi-suspension state. The fact that the amplitude of shear Fig. 1 Result of an undrained cyclic torsional test on Toyoura sand at D r = 70 %: a stress path; b stress-strain curve; c stress-strain curve during the 11th load cycle; d the evolution of the shear strain amplitude c 0 at liquefaction (one value for each half load cycle). (Data from Zhang, 1997 [40]. s is the torsional shear stress in torsional tests, p is the mean effective stress, c is the engineering torsional shear strain) strain c 0 progressively increases in the successive postliquefaction loading cycles and eventually saturates at a bounded value strongly suggests that certain characteristic(s) of sand particle arrangements at the moment of entering the liquefaction state must determine the amount of shear deformation required to take the material out of the liquefaction state. The main objective of this study is to identify a variable (or variables if necessary) to quantify the fabric (or microstructural) evolutions during post-liquefaction cyclic loading that explains the aforementioned intriguing post-liquefaction behaviors of sand. It is desired for this variable to have clear physical meanings, be directly measurable, and be able to predict the c 0 development in each loading cycle. The main challenge is that traditional fabric quantifications mostly, with the exception of the void ratio, rely on inter-particle relationships based on contacts, whereas the very signature of the liquefaction state is the ''absence'' of loading-bearing contacts, thus, conceptually suggesting a fabric entity related to such lack of particle contact. The current work will present a realization of this suggestion.
Direct visual observations and quantitative measurement of fabric characteristics of sand would be the most straightforward methodology to achieve these objectives. However, although recent developments in X-ray tomography technologies have made this class of methods possible (e.g., [1,34]), the spatial and temporal resolutions of the state-of-the-art technologies still fall short of tackling this specific problem. An alternative method is to use particle-based numerical simulation techniques, such as the discrete element method (DEM) [7]. Such methods can simulate granular materials and provide microscopic measurements that are yet very difficult or impossible to obtain in laboratory experiments, and have become an increasingly important tool in the study of the mechanical behavior of granular materials. Numerous applications of DEM in studying sand liquefaction have verified DEM's capability of reproducing the undrained cyclic behavior of sand (e.g., [8,15,20,23,29,37]). Ng and Dobry [23], Dabeet et al. [8], and Kuhn et al. [20] all successfully modeled the generation of large shear strain at liquefaction using relatively simple DEM models. Wei and Wang [37] confirmed the generation and saturation of shear strain at liquefaction using 2D DEM and introduced a new measurement of ''centroid distance'' which showed interesting relations to the post-liquefaction shear strain of sand; their work is conceptually closest to the present one. These DEM studies also validated the use of a constant volume constraint to simulate undrained loading, without having to incorporate an actual fluid phase due to the negligible influence of fluid flow during pseudo-static loading (e.g., [20,23]). However, if fluid flow is nontrivial to the problem under investigation, such as in partially drained loading conditions, coupled methods that can appropriately reflect particle-fluid interaction should be adopted [9,30,43].
The structure of the current paper is as follows. Section 2 presents the undrained cyclic biaxial 2D DEM simulation method and lays out the simulation program. Section 3 performs a comprehensive comparison between DEM simulation of sand's liquefaction-to-post-liquefaction behaviors and the counterpart observations in published laboratory work, thereby validating the adopted simulation method for the proposed objective. In Sect. 4, we prove that none of the existing conventional variables can adequately describe, explain, and predict the undrained postliquefaction deformation development through an exhaustive evaluation. Subsequently, in Sect. 5 a new fabric measurement is introduced to quantify the mechanicalgeometrical structure of granular assemblies at liquefaction states, which unravels micromechanical subtleties of undrained shear strain development. The evolution of the new fabric measurement and its complete interpretation within non-liquefaction states is studied and presented in Sect. 6, while its significance and its possible future extensions are discussed in Sect. 7.
2 Discrete element model 2.1 Model setup DEM simulations of 2D circular particles are conducted in this study to qualitatively investigate the undrained cyclic behavior of sand. 2D ideal particles are known to behave somewhat differently from real world 3D particles in a number of ways, having smaller void ratio values and coordination numbers, among others. However, 2D and 3D particles share many common microscopic mechanisms, particularly those related to particle interactions, which govern macroscopic behavior of granular materials in 2D and 3D alike. The fact that ''A-class predictions'' made using DEM [11] have been validated in later experimental observations [32] on actual sand samples provides strong confidence in 2D DEM's capabilities in qualitatively capturing behaviors of sand.
A DEM simulation package PPDEM that has been successfully applied in a number of sand fabric-related studies [11][12][13][14]33] is employed. Although PPDEM has the ability to model arbitrarily shaped 2D particles, circular particles are used to avoid the complication of strong anisotropic behavior and enable us focus solely on liquefaction. The samples consist of circular particles ranging from 0.30 to 1.00 mm in diameter, with a median particle size D 50 of 0.72 mm, and a uniformity coefficient D 60 /D 10 of 2.13. ''Master packs'' of particles are fabricated with different initial void ratios through a pluviation process described in detail by Fu and Dafalias [11]. Specimens for undrained cyclic biaxial simulations are then trimmed out of the master packs and isotropically consolidated under initial effective mean stress p in (Fig. 2a). After consolidation, undrained cyclic biaxial loading is achieved through controlling the velocities of the four walls enclosing the specimen. The two walls in the compression direction move toward the specimen at controlled velocities, and the velocities of the other two walls ensure constant enclosed area by the four walls. Constant velocities are used in the compression direction except for the ramping-up and ramping-down near loading direction reversal (i.e., the compression and extension directions are reversed) triggered at specified deviatoric stress values (i.e., q ¼ ðr y À r x Þ=2 ¼ q max , where q is the deviatoric stress and q max is the deviatoric stress amplitude during cyclic loading; r x and r y are the average normal stresses along the x and y directions of the sample, which are also the two principal stresses for 2D biaxial loading). The loading rate is sufficiently slow to guarantee pseudo-static states throughout the simulations.
A total of 17 undrained cyclic biaxial tests are simulated in this study. Variables including the void ratio, initial mean stress, deviatoric stress amplitude, pre-loading history, and specimen size are investigated (Table 1). Apart from simulation e21q20p100L, which consists of 12,000 particles for proving size-insensitivity, the other specimens all have approximately 5000 particles each. An inter-particle friction angle of 35°is used, and the friction angle for particle-loading wall contacts is 0.2°to minimize boundary constraint. The contact law used in the simulations was described by Fu and Dafalias [11], and the inter-particle contact normal and tangential stiffness parameters (K n and K s ) used in the simulations are 2100 and 700 GPa/m (the dimension is force/overlap area/unit thickness), respectively. For the pseudo-static simulations in this paper, parameters for the viscous damping components and the load rates are chosen so that the simulation results are insensitive to a moderate variation of their values. 20-120 loading cycles are performed for each simulation until the saturation of c 0 at liquefaction and each loading cycle consists of at least 40,000 time steps, ensuring sufficient temporal coverage and resolution. Drained biaxial compression and stress rotation simulations results for the same virtual material are available in Fu and Dafalias [13].

Calculation of stress, strain, and fabric quantities
The overall axial strains of the sample are calculated based on the displacement of the boundary walls as: where e x and e y are the normal strains along the x and y axes of the sample, which are also the two principal strains for 2D biaxial loading; L x0 and L y0 are the initial wall distance in x and y directions, respectively; DL x and DL y are the differences between current and initial wall distances in x and y directions, respectively. Following geomechanics conventions, compressive strain and stress are considered positive. c = e y -e x is referred to as the engineering biaxial shear strain, which is actually the maximum engineering shear strain on the 45°plane in 2D biaxial loading. The increment of c generated within a liquefaction state is defined in consistency with the experiment terminology in Fig. 1 as the post-liquefaction shear strain c 0 . Although the strain values in some simulations are too large to be considered within the ''small strain'' regime, the above ''engineering strain'' formulations are still used to be consistent with the prevailing result reporting conventions of laboratory testing. The average effective stresses of the samples are calculated based on a weighted summation of inter-particle forces [3]: where S is the area (2D) of the sample; k is an index that runs over all the inter-particle contact points, and N c is the total number of contacts; f x k and f y k are the x and y components of the kth contact force; l x k and l y k are the x and y components of the branch vector (connecting the centers of the two particles in contact) of the kth contact. The deviatoric stress is then defined as q = (r y -r x )/2, and mean stress p = (r y -r x )/2 for 2D biaxial loading. In the simulations conducted in this paper, the shear stress is observed to be less than 1 % of the mean stress.
In order to focus on the liquefaction state itself, we need a practical definition of the ''zero stress'' state for the DEM numerical simulations. Because particles continue to interact with each other through contacts and thereby generate small contact forces even in the liquefaction state, the mean effective stress calculated using Eq. (2) is never absolutely zero. We consider a liquefaction state to be reached when the mean effective stress of the sample is less than one percent of the initial mean effective stress. The value of c 0 is insensitive to a modest variation of this threshold value due to the much higher stiffness in nonliquefaction states compared to that in the liquefaction state.
The term ''fabric'' used in this study refers to features related to the spatial and geometrical configuration of grains in granular materials such as sand. Examples include void ratio, coordination number, and various fabric tensors that reflect the directional properties of the material, including particle orientation-, contact normal-, and voidbased fabric tensors, etc., [13,22,24]. However, other fabric features yet to be developed and fully appreciated could also be highly relevant to sand behavior under certain circumstances.
Void ratio e is an internal variable commonly used in soil mechanics that can be easily calculated in DEM. Another fabric quantity, the coordination number C, is defined as the average number of contacts per particle and can be easily obtained in DEM simulations as resolving inter-particle contacts is at the core of DEM. The contact normal fabric tensor F characterizes the spatial distribution features of inter-particle contact normal directions and is defined after Satake [26] to be: where n k is the unit vector representing the normal direction of the kth inter-particle contact. The intensity of contact normal fabric anisotropy was defined in Fu and Dafalias [12] as a c ¼ F I À F II with F I ; F II being the major and minor principal values of F, respectively. This definition yields always a positive value for a c , and it will be modified to a c ¼ F yy À F xx , with F xx and F yy being the normal components of F along x and y directions, respectively, in order to capture the sign change of fabric anisotropy intensity. For example, a vertical major principal fabric yields a positive intensity and a horizontal one produces a negative intensity measurement. This is necessary for the current study in order to differentiate processes taking place in the first half of each loading cycle and those in the second half. Notice that besides the sign issue the foregoing definitions are equivalent based on the fact that in bi-axial compression the off-diagonal components of F are zero, thus, F xx and F yy are in fact principal values of F (only in such case the definition of a c in terms of F xx and F yy makes sense). In the simulations conducted in this paper, the off-diagonal values of the fabric tensor are observed to be well less than 1 % of those in the x and y directions. Particle orientation-based fabric tensors are not applicable to the study of circular particles, and fabric tensors based on void shapes have been found to be strongly correlated with those based on inter-particle contact normal directions [13].
3 Validating DEM model's ability to capture liquefaction-related behaviors of sand Figure 3 shows typical stress and strain results from test e21q20p100, which has a 0.21 target void ratio, 20 kPa deviatoric stress amplitude, and 100 kPa initial mean stress. See the note beneath Table 1 for the naming convention of the simulations. The stress path in Fig. 3a shows that the DEM results closely resemble typical laboratory test results (e.g., Figure 1a) of sand in terms of progressive reduction of effective stress toward liquefaction and the ''butterfly orbit'' of stress path after initial liquefaction. More importantly to this study, the stress-strain relationship in Fig. 3b shows the generation of large but bounded shear strains each time the stress state goes through liquefaction, which occurs once during each monotonic half load cycle (from one load reversal instant to the next), echoing the observations of laboratory tests (Fig. 1b).
Similar to the results of laboratory undrained cyclic tests on sand, c 0 increases as the cyclic loading continues, and the increase slows down to eventually saturate at c 0s in the DEM tests, which is 0.223 for test e21q20p100 shown in Fig. 3c. Another test e21q20p100L on a larger sample of 12,000 particles with the same settings as test e21q20p100 is conducted for validation purposes, with the test results in Fig. 4 showing that the overall responses of the two samples are very similar. Note that the void ratio of the larger sample is slightly higher at 0.2101 than the smaller sample's 0.2077, which contributes to the larger sample's slightly stronger tendency to contract initially. The saturated c 0s for the larger sample is 0.229, only 2.7 % greater than the counterpart value from the smaller specimen, showing that the sample size of 5000 particles is sufficient and is used in the rest of the DEM tests.
The qualitative agreement between DEM and laboratory test results and that between numerical specimens of  An advantage of DEM simulation is its great convenience in quantifying fabric quantities. Figure 5 shows the evolution of the contact normal fabric anisotropy intensity a c in test e21q20p100. a c starts at close-to-zero after consolidation, indicating that the initial state is nearly isotropic. Prior to initial liquefaction, a c evolves cyclically with increasing amplitudes. After initial liquefaction, fabric anisotropy intensity a c is observed to change drastically at liquefaction states with a sign reversal within each liquefaction state. It tends to follow a consistent path after very few load cycles, with maximum and minimum values of 0.23 and -0.23 in test e21q20p100.
Another fabric quantity worth investigating is the coordination number C. Figure 6 shows the evolution of C for test e21q20p100. The coordination number is initially 2.9 and undergoes cyclic variations with an overall decreasing trend during pre-liquefaction loading. Once the coordination number drops below 2.2, the sample enters liquefaction, during which most of the shear strains occur (Fig. 6a), and does not regain effective stress until the coordination number exceeds 2.2 again (Fig. 6b). This threshold value of 2.2 not only remains valid among all post-liquefaction loading cycles, but also is shared by all the simulations performed. During liquefaction, the coordination number reaches a minimum of approximately 1.0 in this test. The results of fabric intensity and coordination number from the DEM simulations suggest that during liquefaction, sand samples can undergo a significant amount particle rearrangement, which could be related to the generation of post-liquefaction shear strain.

The relationships between conventional variables and post-liquefaction shear strain development
The variables evaluated in this section include loading parameters and fabric-related state variables commonly used in soil mechanics. The influence of loading parameters, which are external variables, is studied through conducting simulations with different deviatoric stress amplitude q max (20,25,30, and 35 kPa) and initial mean effective stress p in (60, 80, 100, 120 kPa), as listed in Table 1 on the same virtual sample. Since these two loading quantities remain constant for each test, they cannot be held accountable for the gradual increase in c 0 among loading cycles, but could be related to the final saturated post-liquefaction shear strain c 0s . Figure 7 shows the stress paths and stress-strain curves from two tests with different q max and p in , namely test e21q35p100 and e21q25p80. Although different loading conditions cause remarkable differences in stress paths and stress-strain curves, the eventual saturated shear strain amplitudes are almost identical. Figure 8 plots the results from all the tests with different deviatoric stress amplitudes and initial effective stresses, showing that c 0s is generally unaffected by the change in loading conditions, all falling  Fig. 6 Coordination number C in test e21q20p100: a C against shear strain c, b C against deviatoric stress q within the range of 0.215-0.225. This is a rather small range because the subsequent analysis will show that c 0s of the 17 simulations cover the range between 0.045 and 0.391. Moreover, this small variation is likely the consequence of the random error inevitable in any simulation or laboratory experiment since no systematic trend can be identified in Fig. 8a or b.
The void ratio e is one of the most basic internal fabric quantities of sand. Similar to the external loading variables studied, it does not change during undrained loading and thus does not account for the increase in post-liquefaction shear strain across loading cycles. Five DEM tests (e18q25p100, e19q25p100, e20q25p100, e21q25p100, and e22q25p100) on samples with different initial void ratios (0.1849, 0.1911, 0.2049, 0.2077, and 0.2235, respectively) are conducted under the same loading conditions (Table 1). Figure 9a, b depicts the stress paths and stress-strain curves from two tests with different void ratios of 0.2077 and 0.1849, respectively, where the looser sample experiences much greater shear strain than the denser sample does, with c 0s of 0.217 compared with 0.114. When the c 0s values of these five tests are plotted in Fig. 10 as hollow circular markers, a strong correlation can be observed between c 0s and the void ratio. The results so far (those in Fig. 9a, b and the hollow circles in Fig. 10) seem to suggest that the initial void ratio is the sole factor affecting c 0s , a simple, desirable relationship that makes intuitive sense. However, the next set of simulations will prove this is not the case.
After completing these five undrained cyclic biaxial tests, the specimens are reconsolidated under drained condition from a liquefaction state of zero shear strain to a mean effective stress of p = 100 kPa and then undergo the same undrained cyclic biaxial loading, which are identified as tests e18q25p100Re, e19q25p100Re, e20q25p100Re, e21q25p100Re, and e22q25p100Re. Note that the names reflect the void ratios of the original specimens. The reconsolidated specimens have void ratios of 0.1832, 0.1843, 0.1989, 0.1960, and 0.2068, respectively. Generally,  originally looser specimens experience greater void ratio reductions. Figure 9c shows the results of stress and strain from test e21q25p100Re. Compared with the results from its corresponding original test e21q25p100 in Fig. 9a, the test after reconsolidation (e21q25p100Re) results in a much smaller c 0s of 0.109, which is to be expected due to the decrease in void ratio during reconsolidation. However, comparison of the shear strain results between Fig. 9b and c shows that although these two test samples are significantly different in void ratio (0.1960 and 0.1849), they exhibit similar c 0s (0.109 and 0.114). The c 0s of the five tests conducted after reconsolidation is included in Fig. 10 as solid circular markers. The combined results in Fig. 10 clearly indicate that tests samples with very different void ratios could result in similar c 0s , and vice versa. Therefore, a unique relationship between the void ratio and saturated post-liquefaction shear strain c 0s does not hold and other fabric characteristics that are affected by the loading history must also have significant influence on c 0s . Figure 5 shows that the intensity of contact normal fabric anisotropy a c experiences a significant change in the amplitude of a c as the material evolved through liquefaction cycles. To further investigate the relationship between a c and post-liquefaction shear strain, a c0 (the difference between the a c value as the material enters and exits Hollow markers represent results for virginally consolidated specimens and solid markers represent results for re-consolidated specimens following the cyclic loading on the former liquefaction, respectively) and c 0 in test e21q25p100 are shown in Fig. 11. Although both quantities increase with increasing number of loading cycles and eventually saturate, a c0 saturates much sooner than c 0 does, and c 0 continues to increase substantially (from 0.133 to 0.217) after a c0 has already saturated. Figure 12 plots the saturated c 0s and a c0s in the aforementioned 10 different tests, showing that the a c0s is clearly unable to predict the saturated shear strain c 0s , which is understandable as the contact normal fabric tensor mainly provides a measurement for the orientation aspect of granular materials.
The simulation results show a significant drop in coordination number in each specimen during liquefaction (Fig. 6). Figure 13 compares the evolutions of the minimum coordination number at liquefaction (C min ) and c 0 in test e21q25p100, in a similar fashion to the comparison between a c0 and c 0 . Although C min decreases with increasing number of loading cycles, C min reaches a plateau much earlier than c 0 , and c 0 continues to increase after C min stabilizes. Saturated c 0s shown in Fig. 14 does not have a unique correspondence with C mins .
In summary, the results above show that the conventional external loading variables and internal fabric metrics investigated cannot provide a convincing explanation for and are inadequate in describing the generation and saturation of c 0 , and that other intrinsic quantities associated with this phenomenon must be sought after.

The neighboring particle distance
The previous analysis has shown that the coordination number C is a strong indicator of whether the material is in a liquefaction state or not. Once C of the virtual material    drops below 2.2, the sand grains lose the ability to form an effective load-bearing skeleton. It is common sense knowledge that it is possible for 2D grains to stack on each other to form a 1D column with a C value of two. Although this column can bear load, it is unstable and susceptible to ''buckling'' [25]. In a stable 2D structure, each particle needs to have three or more contacts. This argument seems to contradict the 2.2 threshold coordination number. The reason for this discrepancy is that we do not need all the particles to participate in the load-bearing skeleton and it is well known that a force chain involving only fraction of particles can do the job. The reasoning here also implies that the threshold value, although should certainly be smaller than three, is likely to depend on the grain size distribution. However, the present work did not quantitatively investigate this aspect.
Once the material is in a liquefaction state, the material is overall in a ''semi-suspended particle'' regime, in which the coordination number is associated with random transient contact between particles instead of load-bearing structures. The current coordination number tells how many more contacts, on a per-particle basis, the material needs to build a stable load-bearing structure. More directly related to c 0 is the distance between each particle and its neighboring particles with which the contacts are necessary for constructing the skeleton. The microstructural essence of c 0 is nothing but the inter-particle displacement or macroscopic deformation necessary for bringing these neighbor distances to zero.
To test this hypothesis, the concept of ''Neighboring Particle Distance'' (NPD) is proposed to quantify the aforementioned particle distances in the ''semi-suspended particle'' regime. The NPD of an individual particle is the mean surface-to-surface distance between this particle and its n closest neighbor particles, with n being the number of contacts needed to support a stable load-bearing structure; this n can be therefore defined as the integer immediately greater to the smallest coordination number supporting a stable non-liquefaction state. n = 3 in the 2D case as discussed in the previous paragraphs. Figure 15 shows a conceptual illustration of surface-to-surface distances between a 2D particle and its 3 closest neighbor particles, where the NPD for the particle at the center is the mean of D1, D2, and D3. In this example case D1 = 0 as the two particles are in contact while D2 and D3 assume positive values. In DEM, the possible geometric intrusion of two contacting particles may result in a negative net surface-tosurface distance. However, a zero distance is used for the calculation of NPD in such situations. For a material consisting of an assembly of particles, the mean value of NPD over all particles is a fabric characteristic or internal variable of the assembly that can be called MNPD and defined analytically by where N is the number of particles in the assembly and k is an index that runs over all the particles within the scope of analysis. MNPD is proposed with the intention to reflect the amount of rearrangement needed for a granular assembly at liquefaction to reach a stable load-bearing state, which would in turn be associated with post-liquefaction deformation and provide a logical explanation for the phenomenon observed during undrained cyclic tests on sand. Like other indices that measure the extent of contacts, such as the coordination number, the MNPD measures the extent of contact loss. Figure 16 shows the evolutions of MNPD and mean effective stress (p) during two half loading cycles (the 8th and 15th cycle) in test e21q25p100. Upon entering liquefaction, MNPD first increases to reach a maximum value (MNPD max ) and then gradually decreases as the shear strain increases until the sample eventually regains effective stress. By comparing Fig. 16a, b, we see that the later loading cycles are associated with both larger shear strain amplitudes and larger MNPD max values than the earlier cycles are. This suggests that the proposed fabric measurement MNPD does reflect the amount of particle rearrangement during each occurrence of liquefaction, the full potential of which is critically evaluated in the rest of this section. Note that we only study the responses of MNPD in liquefaction states in the present section, and the next section focuses on the interpretation of MNPD in nonliquefaction states.
Similar to the analyses of the contact normal fabric intensity and minimum coordination number, the developments of MNPD max and c 0 in test e21q25p100 are Fig. 15 Conceptual illustration of the surface-to-surface distance between a 2D particle and its three closest neighboring particles. The fourth and fifth closest particles are also shown but they do not participate in the calculation of NPD. The distance between the particles have been exaggerated for visualization compared in Fig. 17. MNPD max and c 0 exhibit extremely similar patterns of increase and saturation with respect to loading cycles and have a strong correlation with each other, unlike the previously analyzed fabric quantities which reach saturation levels much earlier (Figs. 11, 13). Figure 18 plots MNPD max against c 0 for all the post-liquefaction cycles of all the 17 simulations conducted in this study, which has a total of 1172 data points. Despite the inevitable random variation, these data suggest a very strong correlation between MNPD max and c 0 in all postliquefaction cycles. This means that not only can MNPD max uniquely reflect the post-liquefaction shear strain c 0 in a particular test (Fig. 17), but it is also able to uniquely determine the c 0 for any test on the same material regardless of specimen initial states, loading conditions, and loading histories. Figure 19 shows the relationship between the final saturated c 0s and MNPD maxs values in the 17 different tests. With the correlation between MNPD max and c 0 of each load cycle in mind, it is then not surprising to find that an even stronger correlation exists between saturated MNPD maxs and c 0s for all the void ratios and loading histories included. One might argue that relatively strong correlations also exist between void ratio e and c 0s , between a c0s and c 0s , and between C and c 0s , as suggested by Figs. 10, 12, and 14, respectively. However, when compared with the much stronger and much more definitive relationship between MNPD maxs and c 0s shown in Fig. 19, those correlations are at best circumstantial and the scattering of the data indeed reflect void ratio, contact normal fabric anisotropy intensity, and coordination number's inability to predict post-liquefaction deformation development. A comparison between Fig. 19 and Fig. 10 indicates that while the overall void ratio is the same, internal particle arrangement and structure of sand could change and have significant influences on the behavior of sand, making MNPD a more appropriate fabric measurement than void ratio e in describing the undrained postliquefaction cyclic shear strain.   Correlation between c 0 and maximum mean neighboring particle distance at liquefaction (MNPD max ) in each half loading cycle after initial liquefaction in 17 different tests. Each of the 1172 data points represents a half load cycle after liquefaction. When fitted with a second order polynomial, the data yielded a coefficient of determination of 0.89 6 NPD in non-liquefaction states Figure 16 shows that as the sample leaves liquefaction and starts to bear an increasing amount of stress, MNPD actually tends to increase. This is somewhat peculiar since intuitively a non-liquefaction state should have smaller MNPD than a liquefaction state and higher mean effective stress should be associated with tighter particle packing and thus smaller MNPD. To explore this discrepancy, Fig. 20 depicts in detail the development of MNPD and coordination number C in various loading segments for simulation e21q25p100. We pick the data of 2000 time steps from the 17 to 18th loading cycles at which c 0 has already saturated. Three distinct segments can be identified: the loading segment (with the absolute value of q increasing), the unloading segment (|q| decreasing), and the liquefaction segment (q & 0) as denoted in Fig. 20. In addition to the peculiar increase of MNPD during loading, a close examination of Fig. 20b, d reveals a decreasing trend of MNPD during unloading, which is also counterintuitive for the same reasons.
We discover that the counterintuitive trends of MNPD in both the loading and unloading segments are caused by biased statistics: The particles that have high NPD values dominate the overall MNPD results, but they do not participate in load-bearing in any substantial way. To demonstrate this, we sort the particles based on individual particles' NPD values at each time step. We plot the evolution of MNPD for of the top 25 % particles in Fig. 20f and that for the lower 75 % in Fig. 20h. We calculate the coordination numbers for the top 25 % particles (based on NPD ranking, not coordination number ranking) and the lower 75 % separately and plot them in Fig. 20g, i, respectively. It becomes apparent that within the whole particle assembly, the top 25 % particles and the lower 75 % behave in very different ways. In the loading segment, MNPD for the top 25 % particles increases, while that for the lower 75 % actually shows a clear declining trend. Because the MNPD value of the top 25 % particles is approximately 50 times greater than that of the lower 75 %, the overall MNPD is dominated by the former. The coordination number of the top 25 % is below 1.0 in all three loading segments, indicating that these particles remain in the semi-suspended state even when the overall specimen is in a non-liquefaction state. On the other hand, the lower 75 % particles have a coordination number greater than 2.5 during loading, meaning that these particles form the load-bearing skeleton. The same observations apply to the unloading segments, where the load-bearing lower 75 % particles have an increasing MNPD, but it is obscured by the decreasing MNPD of the semi-suspended top 25 % particles in the overall statistics. This choice of the 25-75 % division is somewhat arbitrary and was made with a trial-and-error procedure.
These observations suggest that in the non-liquefaction states, the overall MNPD is statistically limited as an indicator of material state or behavior. The MNPD values are dominated by particles that are in the semi-suspension state, whereas the contribution of the particles actually participating in the load-bearing skeleton is severely obscured. However, we argue that this does not undermine the MNPD's role as a strong indicator and predictor of material behavior in the liquefaction state since all particles are in the semi-suspension state in liquefaction. It is possible that MNPD for a specific group of particles might provide a greater prediction power than the simple overall MNPD does, but this possibility is not pursued in the current study because the relationship obtained from the simple overall average (such as that in Fig. 19) is already satisfactorily definitive. Fig. 19 Relationship between saturated post-liquefaction shear strain c 0s and saturated maximum mean neighboring particle distance at liquefaction (MNPD maxs ) in 17 different tests. When fitted with a second-order polynomial, the data yielded a coefficient of determination of 0.98 Acta Geotechnica (2016) 11:1321-1337 1333 The great relevancy of MNPD to sand's liquefaction behavior is also demonstrated by how it provides insight into the different behaviors between virginally consolidated specimens and re-consolidated (after liquefaction) specimens, which has been puzzling in our previous analysis. The stress paths in Fig. 9c, a show that compared with the virgin sample, the reconsolidated, denser sample requires less load cycles to reach initial liquefaction and exhibits a stronger tendency to contract. This is at odds with the patterns observed for the tests on the original samples, where denser specimens always had higher cyclic liquefaction resistance. The decrease in liquefaction resistance of samples reconsolidated after previous post-liquefaction undrained cyclic loading has also been reported by Wahyudi et al. [35] in their cyclic simple shear test results. In our DEM simulations of e21q25p100 and e21q25p100Re, although the reconsolidated sample has a lower void ratio than the original sample, it has a higher initial MNPD value, which is 3.53 9 10 -3 mm compared to 3.25 9 10 -3 mm for the original sample, thereby Fig. 20 A part of the loading process in test e21q25p100 focusing on the development of MNPD and coordination number C: a shear strain c showing the load step frame from which the excerpt was taken; b mean effective stress p; c deviatoric stress q; d MNPD of the entire assembly; e coordination number C of the entire assembly; f MNPD of the top 25 % of particles based on NPD ranking; g C of the top 25 % particles based on NPD ranking; h MNPD of the lower 75 % particles based on NPD ranking; i C of the lower 75 % particles based on NPD ranking providing a possible explanation to the discrepancy. This shows that the new fabric measurement is not only effective in describing the post-liquefaction shearing of sand, but could also be useful in non-liquefaction states. However, as we have pointed out in analyzing the MNPD of different portions particles, realization of the full potential of MNPD in characterizing non-liquefaction states requires a comprehensive investigation beyond the scope of the current study.

Concluding remarks
This study investigates at grain level the phenomenon of the buildup of large but bounded shear deformation of sand at the liquefaction state during undrained cyclic loading. The most significant contribution of the current study is the establishment of a definitive connection between post-liquefaction shear deformation development, which is previously well-known yet poorly understood, and a new, theoretically measurable intrinsic fabric metric with a clear physical interpretation. This new fabric measurement, Mean Neighboring Particle Distance (MNPD), is formulated to capture the microstructural features of granular materials that govern deformation behavior in the liquefaction state. This work unveils a new aspect of mechanical behavior of granular soils that traditional fabric measurements have largely overlooked. While the void ratio reflects the overall packing density of particles and various fabric tensors reflect the orientation aspects of soil, the proposed MNPD provides an effective scalar fabric quantity that reflects the internal arrangement of the particle and void system, especially in the particle semi-suspension regime, and can have different values for the same void ratio. This new metric is shown to have significant influences on the mechanical behavior of soil and would be an important addition to existing fabric descriptions of granular soils. The findings in this study are vital to developing macroscopic and practical models, based on physical evidence, that can reasonably reproduce and predict the deformation of sand during liquefaction and have the potential to provide better guidance to engineering practice. The fact that we have proved that conventional fabric measurements simply are unable to achieve our intended purpose, namely providing a definitive correlation with shear strain development during cyclic liquefaction, suggests that the new quantity is not only novel, but it is also essential for characterizing a fabric feature related to sand stress-strain behavior in liquefaction. This is achieved by correlating shear strain, a non-state variable (i.e., a variable which cannot be defined solely by the current state of a sample) to MNPD which is a state variable (i.e., a variable that in principle can be measured solely by the current state of a sample). Many aspects of this new quantity and potentially derived quantities should be further studied to enhance our understanding of sand behavior.
As the first step into the research of a new class of fabric measurements, the DEM simulations are based on 2D circular particles to capture the general response of cohesionless granular materials under undrained cyclic loading. Many observations and conclusions are expected to still hold true for more complex shaped particles and also in 3D, as the core concept behind NPD, inter-particle distances, determine the amount of deformation that can take place within the liquefaction state regardless of particle shapes and dimensionalities. For 3D cases, the number of contacts n needed to support stable load-bearing structure would be different from that in 2D and can be determined in 3D following the same procedure as was done for 2D in this study. The most likely value of n should be 4 in 3D as tetrahedron in 3D is directly analogous to triangles in 2D. This speculation makes a plausible suggestion that the smallest coordination number that can support a stable nonliquefaction state in 3D is between 3 and 4. Although the choice of n = 3 seems to have yielded an MNPD metric with satisfactory prediction power in 2D, it is also possible that a smarter choice of n (e.g., particle-dependent) could further improve this fabric metric.
NPD and MNPD can be easily measured in any particlebased simulations and are theoretically measurable for real world sand specimens using advanced imaging techniques. However, Fig. 20 shows that for the particle sizes simulated, a resolution finer than 0.1 lm is necessary to meaningfully resolve the evolution of MNPD. Considering that the specimens used are 40 mm in each dimension, this means that a spatial resolution to specimen size ratio of 1:400,000 would be necessary, which is approximately two orders of magnitude more precise than the state-of-the-art techniques. It therefore turned out to be a wise and fortunate decision to have based the current study on DEM simulations. Nevertheless, the direct measurability, even though currently only theoretically possible, of MNPD is an important advantage of this new fabric measurement over other empirical variables used in many constitutive models, which can only be determined by curve fitting.
The analysis of the relationship between MNPD and a specimen's contraction (or dilatancy) tendency prior to liquefaction at the end of Sect. 6 indicates that the internal arrangement of particles reflected by MNPD, which is not captured by the void ratio, affects not only post-liquefaction sand behavior, but also certain pre-liquefaction behaviors. However, the roles of different statistical portions of MNPD are shown to be quite different. Thus, the relationship between NPD and the dilatancy of granular materials should be investigated in depth. Acquiring a better understanding of NPD not only in liquefaction state but also in non-liquefaction state would be of great significance in further investigation of the role of fabric in the behavior of sand and provide a possible path to incorporating such fabric measurements into a continuum constitutive framework for practical purposes.
NPD is defined in this study as the mean surface-tosurface distance between a particle and its n closest neighboring particles, which does not take into consideration the possible effect of particle sizes. For example, the fourth closest neighboring particle that is large in size may be more important than the third closest neighboring particle that is very small, as in the end the center particle is likely to form a stable load-bearing configuration with the large neighbor particle. In this case, a weighted definition for neighboring particle distance based on relative particle sizes may be useful pending further investigation. Normalization of the NPD and MNPD measurements is not considered in this study as all the simulations used the same virtual material, although a normalization of MNPD into a dimensionless form would facilitate constitutive modeling. However, normalizing MNPD against the mean particle size or a representative particle size is a straightforward modification and could be appropriately adopted in future quantitative studies, especially when multiple material types are involved.