An anisotropic plasticity model incorporating fabric evolution for monotonic and cyclic behavior of sand

Fabric anisotropy has a significant influence on the mechanical behavior of sand. An anisotropic plasticity model incorporating fabric evolution is formulated in this study. Information on the overall stress–strain relationship and micromechanical fabric states from DEM numerical tests is used in the development of the constitutive model, overcoming the difficulties of fabric measurement in physical tests. The framework of the model and its formulations for fabric evolution, plasticity, and dilatancy enables it to capture the strength, shear modulus, and dilatancy of sand under both monotonic and cyclic loading. The model is validated against DEM numerical tests and physical laboratory tests on samples with different initial fabric, showing good agreement between the simulation and test results for the anisotropic stress–strain behavior of sand. The use of DEM test data also allows for the validation of the model on the micromechanical fabric level, showing that the model can reproduce the fabric evolution and its influence on key constitutive features reasonably well. The model is further applied to analyze the liquefaction behavior of sand, exhibiting the significant influence of fabric anisotropy on both liquefaction resistance and postliquefaction shear deformation.


Introduction
Anisotropy has long been acknowledged to be a salient feature of sand [8]. Sand deposited under gravity or subjected to anisotropic stress inevitably exhibits anisotropy in strength and deformation (e.g., [1,3,42,61,77,79,80]). As possible consequences of such anisotropic characteristics, the bearing capacity and earth pressure of sand can be significantly directional dependent [4,7,49]. Yoshimine et al. [88] showed that under undrained monotonic shearing in various directions, sand behavior can range from mostly dilative hardening to strongly contractive, approaching static liquefaction.
Due to its significance and unavoidable nature, a myriad of studies have been focused on various aspects of sand anisotropy. Research on strength anisotropy has been carried out to establish mathematical formulations for the anisotropic strength criterion of sand (e.g., [18,20,27,33,39,53]). The influence of anisotropy on dilatancy has been analyzed (e.g., [25,43,69,88]), with strong evidence indicating dilatancy to be directional dependent. Anisotropy in the small strain behavior of sand has also been investigated [21,24,32]. These important advances in different facets of anisotropic sand behavior undoubtedly contribute to the development of anisotropic constitutive models.
The origin of the anisotropic behavior of sand lies within the microstructural fabric of the material [13,46,56]. Two approaches for constitutive model development have been adopted to incorporate the anisotropic fabric and microstructure of soil. One approach introduces assumptions for the kinematic and mechanical characteristics of sand at the microlevel, such as the particle, particle group, or generically oriented plane levels [29,30,40]. The microlevel behavior is then integrated to obtain the macrolevel stress-strain relationship (e.g., [9,10,44,50]). Although this approach avoids most of the assumptions at the macrolevel, the assumptions made at the microlevel, which serves as an intermediary between actual granular sand and the continuum constitutive model, are nonetheless difficult to validate and calibrate. Another approach is to introduce a statistical representation of the microstructural fabric of sand to the continuum constitutive formulations (e.g., [15,35,47,64,65,79,85]), often making assumptions for the dependency of stress-strain relationship on the fabric tensor [57]. This approach requires quantification of the fabric tensor and its evolution and more importantly understanding of its role in various continuum constitutive components.
Although measuring fabric in physical tests is still technically challenging, virtual numerical tests through methods such as DEM (discrete element method [14]) have become an important means complementary to physical tests in the quantification and understanding of fabric in sand [17,22,31,34,45,62,71]. Aided by observations from DEM tests, Li and Dafalias [36] proposed the anisotropic critical state theory (ACST). ACST adds a fabric anisotropy condition to the classical critical state theory and introduces the influence of fabric anisotropy and its evolution on the state of sand. These advances have promoted rapid development in fabric tensor-based anisotropic plasticity models. Such models have been upgraded from only being able to consider inherent anisotropy with fixed fabric tensors [15,35,47] to incorporating the evolution of fabric [19,51,52,78,91].
Due to limitations in physical test measurement, fabric tensor-based anisotropic plasticity models have mostly only been validated indirectly against physical test results on the macrostress-strain relationship level, while quantitative validation of the micromechanical assumptions for fabric evolution and its influence on important constitutive components such as modulus and dilatancy has rarely been conducted [26]. DEM numerical test data can be adopted to validate the fabric anisotropy-related assumptions made for these models. On the other hand, existing anisotropic plasticity models have mostly been formulated and applied for monotonic loading, with little consideration for the influence of fabric anisotropy on the cyclic behavior of sand. However, evidence has suggested that fabric anisotropy can significantly affect the liquefaction resistance of sand and should be considered in the modeling of cyclic sand behavior [75,76,80,86,87].
This study aims to develop a plasticity model incorporating fabric evolution for monotonic and cyclic sand behavior and utilizes quantitative micromechanical information obtained through DEM numerical tests in the development and validation of the continuum-based constitutive model. Section 2 of this paper discusses the influence of fabric anisotropy and its evolution on the strength, shear modulus, and dilatancy of sand using DEM test data. Based on these understandings, the detailed multiaxial formulation of the proposed model is presented in Sect. 3. In Sect. 4, the model is validated for both macrolevel stress-strain relationships and microlevel fabric evolution against DEM numerical and laboratory physical test data and is then used to analyze the influence of fabric anisotropy on the liquefaction behavior of sand. The stress and strain in this study follow conventional soil mechanics sign conventions with compression being positive. Tensors are denoted by bold letters to distinguish them from plain scalars.

Observations on anisotropic behavior of sand in DEM
DEM is a powerful tool in assisting to bridge the gap between macroscopic anisotropic behavior of sand and its microscopic fabric origins. Results from 3D DEM triaxial tests conducted using the widely adopted open-source code Yade [60] are used to illustrate the influence of fabric and its evolution on several key constitutive components in this section, namely strength, modulus, and dilatancy. Four tests are conducted on samples with almost the same initial void ratio of 0.657 ± 0.005 but different bedding plane angles of 0°, 30°, 60°, and 90°, respectively. The bedding plane angle d here refers to the angle between the major principal stress in triaxial loading and the deposition direction, as illustrated in Fig. 1a. The vertical axis is the major principal stress axis. Over 200,000 particles (Table 1) are first randomly generated and deposited under gravity. Four samples consisting of approximately 40,000 particles each are cut out from the deposit after rotating it by 0°, 30°, 60°, and 90°angles, respectively, and are isotropically consolidated under mean effective stress p = 100 kPa. Constant p drained triaxial numerical test is then conducted on each sample through strain-controlled servomechanism loading using rigid walls. 3D clumped particles with aspect ratio of 1.5:1 formed by two rigidly connected identical overlapping spherical particles are used in the numerical tests. Linear elasticplastic contact is used with contact parameters in Table 1, and the loading walls are frictionless to limit boundary effects. The loading rate guarantees quasi-static conditions with the inertial number smaller than 1.0 9 10 -5 [41]. The homogeneity of the samples is evaluated by comparing the state of the sample at various subdomains following the procedures proposed by Fu and Dafalias [18]. Figure 1b plots the measured [5] deviatoric stress ratio g at both the peak and critical states for the four samples with different bedding plane angles. Here, the deviatoric stress ratio g = q/p, with q = r 1 -r 3 and p = (r 1 ? r 2 ? r 3 )/3, while r 1 , r 2 , and r 3 are the principal stress values. The peak g shows clear dependency on bedding plane angle, which reflects initial fabric orientation. The peak g for samples with lower d is significantly greater, agreeing with existing laboratory and DEM test results [23,46,64,83]. At the critical state, g becomes the same under triaxial compression regardless of the initial bedding plane angle, conforming to the critical state theory. The critical state is determined according to its definition, as the state when the sample approaches constant void ratio while continuing to deform in shear under constant stress [55,58]. The macroscopic stress-strain response of the tests is in more detail in Sect. 4 in comparison with constitutive model simulations.
The influence of fabric anisotropy on the shear modulus of sand can also be investigated using the numerical tests. An equivalent shear modulus, represented by the deviatoric stress increment divided by the deviatoric strain increment dq/de q , at different deviatoric strain values, is plotted in Fig. 2. The deviatoric strain e q = (2/9((e 1 -e 2 ) 2-? (e 2 -e 3 ) 2 ? (e 1 -e 3 ) 2 )) 1/2 , where e 1 , e 2 , and e 3 are the three principal strain values. At the initial state, the equivalent shear modulus decreases with increasing d, agreeing with the test observations by Yang et al. [82] on sand. As deviatoric strain increases, not only does the equivalent shear modulus decrease, the difference in equivalent shear modulus between the four samples also reduces. At e q = 0.1, the difference becomes almost indistinguishable. At the critical state, the equivalent shear modulus becomes zero, irrespective of the bedding plane angle.
The dilatancy D at different deviatoric strain values is plotted in Fig. 3. The dilatancy D = de v p /de q p , where the superscript p indicates plastic strain. The plastic strain is calculated by subtracting the elastic strain increment from the total strain in DEM, following Wan and Pinheiro [67]. The dilatancy at the initial state shows strong bedding plane angle dependency, with dilatancy increasing with increasing bedding plane angle, indicating stronger contraction. As loading progresses, the samples become dilative with negative dilatancy, and the samples with lower d become dilative at smaller deviatoric strain values. The and Nakata et al. [43] and qualitatively supports the dependency of dilatancy on fabric anisotropy proposed by Li and Dafalias [36]. The analysis on strength, shear modulus, and dilatancy indicates that these key constitutive components are significantly dependent on the bedding plane angle during the early stage of loading, but this dependency disappears as the samples approach the critical state. The origin of such behavior lies within the samples' fabric states and their evolution, the quantification of which is a major advantage of DEM numerical testing. Fabric tensors of sand can be defined based on different particle-scale features, including contact normal directions, particle orientations, and void orientations. In this study, we mainly focus on the contact normal fabric tensor, which has been found to be closely related to dilatancy anisotropy [69]. The deviatoric contact normal fabric tensor F is calculated based on Satake's formulation [57]: N is the number of contacts within the domain, v k is the unit norm vector in the normal direction of the kth contact, I is the 2nd-order identity tensor, and e is the void ratio. The term 1 ? e is introduced as a per-volume measure, achieving thermodynamic consistency with the continuum definition of fabric [37]. The anisotropic critical state theory suggests that a critical state fabric tensor exists as functions of the lode angle and mean effective stress, irrespective of the initial void ratio and fabric anisotropy [36]. This has been supported by evidence from DEM investigations [18,23,34,71,83]. Hence, F can be normalized by its norm ||F c || at the critical state, resulting in a unit norm tensor F n at the critical state: Figure 4 presents the fabric orientation h n , norm ||F n ||, and ||F n || -A n at different deviatoric strain values for samples with different d. A n is an invariant of the fabric tensor F n and the tensor-valued unit norm deviatoric loading direction tensor n in terms of A n = F n :n = tr(F n n); in triaxial loading, it is assumed n ¼ ffiffi  Fig. 4 Fabric at various states in constant p = 100 kPa drained triaxial 3D DEM numerical tests on for four samples with e in = 0.657 ± 0.005 and bedding plane angles of 0°, 30°, 60°, and 90°: a contact normal fabric orientation h n , i.e., the angle between the major principal contact normal fabric and the horizontal plane; b normalized contact normal fabric norm ||F n ||; c ||F n || -A n , A n = F n :n = tr(F n n) and the relative orientation of fabric and loading directions. Its value is between -||F n || and ||F n ||, and hence 0 B ||F n || -A n B 2||F n ||. h n is the angle between the major principal fabric orientation and the horizontal plane. The average ||F c || at the critical state is used to obtain the normalized fabric in Eq. (2). At the initial state, h n & d ? 90°, while ||F n || & 0.4. ||F n || -A n is distinctly different for samples with different bedding plane angles at the initial state (Fig. 4c). As deviatoric strain increases, the fabric orientation h n of each sample gradually evolves toward 90°, and the difference in h n between the four samples gradually reduces. Difference in norm ||F n || between the four samples is generated at low deviatoric strain levels, with ||F n || greater for samples with lower d. This difference in ||F n || reduces again at low deviatoric strain levels. The difference in ||F n || -A n between the four samples gradually reduces with increasing deviatoric strain. At the critical state, the fabric orientation becomes the same as the loading direction, while ||F n || reaches unity, and ||F n || -A n reaches zero, irrespective of the bedding plane angle. During its evolution, the norm of the fabric tensor ||F n ||, and subsequently A n , can reach a peak beyond its critical state value (e.g., samples with d = 0°and 30°in Fig. 4b), similar to the peak deviatoric stress ratio being greater than that at the critical state for dense sand. The evolution of fabric is the origin of the anisotropic behavior of sand observed in Figs. 1, 2, and 3. In particular, the evolution of ||F n || -A n for samples with different d shows strong resemblance to the strength and evolution of shear modulus and dilatancy. The initial difference in ||F n || -A n of the four samples is associated with the initial differences in stress-strain relationships. As fabric evolves toward the same critical state, the differences in stress-strain relationship diminish along with the difference in ||F n || -A n . Under the same conditions, a greater value of ||F n || -A n corresponds to smaller modulus and stronger contraction. This feature of ||F n || -A n makes it extremely useful in the formulation of constitutive models considering fabric anisotropy. ||F n || -A n adopted in this study is different from the term 1 -A n adopted to reflect anisotropy in many ACSTbased constitutive models. This choice is made with considerations of two very appealing characteristics of the term ||F n || -A n : (1) ||F n || -A n decouples the influence of fabric anisotropy intensity and fabric orientation on anisotropic behavior, as F n k k À A n ¼ F n k k 1 À n F : n ð Þ . This formulation allows for more control over anisotropic behavior when adopted in constitutive models; (2) the decoupling of the influence of fabric anisotropy intensity and fabric orientation on anisotropic behavior results in ||F n || -A n = 0 when the loading direction is the same as that of the fabric tensor (Fig. 4). This is very appealing for the calibration of constitutive models, where many parameters unrelated to anisotropy can be determined via traditional tests without the need to consider the coupled influence of fabric anisotropy, which will be discussed in more detail in Sect. 4. However, it should be noted that this choice neglects the influence of fabric anisotropy intensity on soil behavior when the loading direction is the same as that of the fabric tensor, which is a compromise.

Basic equations and elasticity
Based on the anisotropic behavior of sand and its fabric origins observed in the previous section, an anisotropic plasticity model is proposed. The model builds on the framework proposed by Wang et al. [73], which has been proven to be highly effective in providing a unified description of sand of different conditions from pre-to postliquefaction under both monotonic and cyclic loading [11,68,70], but does not consider anisotropy. Here, the evolution of the fabric tensor and its influence on the strength, shear modulus, and dilatancy of sand is incorporated into the model. The basic equations in incremental form for the model in multiaxial stress space are: e is the strain tensor, the volumetric strain is denoted by e v ¼ tr(eÞ, and the deviatoric strain tensor is e ¼ e À e v =3I. Superscripts e and p represent elastic and plastic, respectively. r is the effective stress tensor, p ¼ tr(rÞ=3 is the mean effective stress, and s ¼ r À pI is the deviatoric stress. The deviatoric stress ratio tensor is here defined as K and G are the elastic bulk and shear moduli, respectively, L is the loading index, m is the deviatoric strain flow direction, and D is the dilatancy. In this model, the dilatancy is decomposed into a reversible and an irreversible part, following Zhang and Wang [87]. h i are the Macaulay brackets with hxi ¼ x for x [ 0 and K and G are defined following Richart et al. [54]: Acta Geotechnica (2021) 16:43- 65 47 G o and j are material constants, and p a = 101 kPa is the atmospheric pressure. Anisotropy in elastic modulus is not considered in this study.

Plasticity
Plastic loading occurs when the load index L is greater than zero: n is a unit norm deviatoric tensor representing the loading direction, which can be determined following the mapping rule in Fig. 5. H is the plastic modulus. Figure 5 shows the three surfaces and corresponding mapping rule of the proposed model in the p plane of the stress ratio space. The maximum stress ratio surface serves as the bounding surface of the model. The critical state surface determines the ultimate critical state stress ratio. The reversible dilatancy surface determines the stress ratio beyond which reversible dilatancy occurs. The shape of these three homothetic surfaces is determined by a function g(h) of the load angle h, for which the formulation of Zhang and Wang [90] is adopted, defined as: where M c = M and M o are the critical state deviatoric stress ratio at triaxial compression and simple shear, which follow the relationship between these two values in the Mohr-Coulomb criteria. The formulation of g(h) can be replaced based on the needs to reflect different features of the influence of intermediate stress coefficient b, which is not the primary focus of this study, and does not affect the overall concept of the model. The size of each surface at h = -30°is determined by the maximum, critical state, and reversible dilatancy stress ratios, respectively, at the triaxial compression state, which will be presented in the following formulations. For the mapping rule, n is the unit norm deviatoric tensor normal to the maximum stress ratio surface at r, which is the projection of r À a in on the surface (Fig. 5).
Here, a in is the previous load reversal point and is updated upon each load reversal to the current r, i.e., when L in Eq. (8) becomes negative.
The plastic modulus H is defined using bounding surface plasticity [16]: h is a model parameter for plastic modulus, M is the stress ratio at critical state in triaxial compression, and is also a model parameter. The state parameter W = ee c proposed by Been and Jefferies [6] is used for state dependency behavior. The critical state void ratio e c is a function of the mean effective stress, which is determined as e c ¼ e 0 À k c ðp c =p at Þ n following Li and Wang [38], where e 0 , k c , and n are model parameters. n p is a model parameter dictating the sensitivity of the plasticity modulus to changes in void ratio. M m is the equivalent maximum stress ratio at the triaxial compression state. When the current equivalent stress ratio g/g(h) is smaller than Mexp(-n p W), M m is the maximum stress ratio that has occurred during loading.

Influence of fabric anisotropy on plastic modulus
The influence of fabric anisotropy on both the strength and shear modulus anisotropy of sand is introduced in Eq. (10) through the term exp(-D 1 (||F n || -A n ) 1/2 ). As discussed in Fig. 4, 0 B ||F n || -A n B 2||F n ||, with A n = F n :n. In Eq. (10), greater ||F n || -A n leads to a smaller exp(-D 1 (||F n || -A n ) 1/2 ), which results in smaller peak deviatoric stress ratio and also smaller H, in accordance with the DEM results shown in Figs. 1, 2, and 4. In the constitutive model, F n is a normalized deviatoric fabric (2), the norm of which is 1 at the critical state. The fabric evolution equation for the model builds on the incremental formulation by Li and Dafalias [36]: This incremental fabric evolution formulation guarantees that at the critical state, the normalized fabric tensor F n becomes the same as n, as suggested by the results in Fig. 4 and many other studies [18,23,34,71,83]. c is a fabric evolution rate model parameter. The incorporation of the dilatancy D following Yang et al. [84] allows ||F n || to reach beyond 1, as observed for dense sand and in the DEM results in Fig. 4b. The use of D to capture the fabric evolution of sand with different densities has been shown to be effective when simulating DEM test results on samples with various densities [75]. For dense sand, peak ||F n || exceeds its critical state value, while for loose sand, ||F n || gradually increases toward its critical state value. Only one model parameter is introduced in Eq. (11) to try and keep the fabric evolution formulation as simple as possible for the current model.
In the proposed model, it is assumed that the deviatoric strain flow direction m = n, for simplicity. If the influence of fabric anisotropy on the noncoaxiality between the strain increment and stress is to be considered, a linear combination between n and the direction of F n in the form of m ¼ ðn À xF n =jjF n jjÞ=jjðn À xF n =jjF n jjÞjj may serve as possible solution, with x being a positive parameter less than 1. Calibration and validation of such noncoaxial behavior should ideally be conducted for radial loading in the p plane in different directions with respect to the fabric orientation and more importantly, for loading with rotation of the principal stress axes, which is not within the scope of the current study.

Dilatancy
The dilatancy's D is a combination of a reversible part D re and an irreversible part D ir following the findings of Shamoto et al. [59] and Zhang [89]: This decomposition has two major benefits: (1) It allows for separate control over contraction and dilation, which is very important when working with cyclic loading; (2) using the formulation of the model, it allows for the generation and eventual saturation of the postliquefaction shear strain during undrained cyclic loading. According to its definition, the volumetric strain e p v;re caused by reversible dilatancy is always nonpositive [87], generating during loading and releasing after load reversal. Hence, D re is determined separately for generation and release following Wang et al. [73]: r d is the intersection of r and the reversible dilatancy surface, i.e., is the deviatoric stress ratio beyond which reversible dilatancy is generated in triaxial compression. n d is a model parameter for statedependent dilatancy. The reversible dilatancy generation D re,gen follows: d re;1 is a reversible dilatancy parameter. After load reversal, reversible is released until e p v;re becomes zero again. The reversible dilatancy release D re,rel is: By definition, irreversible dilatancy induces nonnegative contractive volumetric strain e p v;ir only. Zhang and Wang [90] showed that e p v;ir accumulates asymptotically during cyclic loading, while its accumulation rate also decreases during every single monotonic shearing event. Based on these observations, Wang et al. [73] proposed an isotropic formulation for the irreversible dilatancy D ir,iso : a is another irreversible dilatancy model parameter. The accumulation rate during a monotonic shearing event. c mono is the accumulated shear strain since the last stress reversal. c d;r is a model parameter that can be assumed to be 0.05 by default. v is introduced to enhance contraction upon load reversal. Figure 3 indicates that the dilatancy of sand is significantly anisotropic, especially during contraction. To take dilatancy anisotropy into consideration, the irreversible dilatancy D ir in this model is proposed as:

Influence of fabric anisotropy on dilatancy
D 2 is the positive-valued irreversible dilatancy anisotropy model parameter. By adding the fabric anisotropy term (||F n || -A n )D 2 to the isotropic irreversible dilatancy D ir,iso in Eq. (13), the irreversible dilatancy D ir in Eq. (14) increases in contraction with increasing ||F n || -A n , as indicated by the results in Figs. 3 and 4. Figures 3 and 4 also show that the bedding plane angle has only little influence on the dilatancy during the dilative phase, which has also been discussed in Wang et al. [75]. These observations suggest that dilatancy anisotropy should be able to be reflected by the anisotropic irreversible dilatancy alone in the proposed model. Hence, the influence of fabric anisotropy on D re is neglected.
The development of total, irreversible, and reversible dilatancy during a typical undrained triaxial loading and reverse loading path is plotted in Fig. 6, to visualize the performance of the dilatancy formulation. The model parameters for Toyoura sand, discussed later in Sect. 4, are used here, under the undrained stress path shown in Fig. 6a. In the model, the total dilatancy D is the sum of the irreversible dilatancy D ir and the reversible dilatancy D re . At the initiation of loading, D re remains inactive at zero and D = D ir , being positive, i.e., contractive (Fig. 6b). When the stress ratio exceeds the reversible dilatancy surface in Fig. 5, reversible dilatancy D re is generated following Eq. (14). Also, as loading continues, D ir decreases gradually according to Eq. (16). When the sum of D ir and D re becomes negative, the model exhibits overall dilatancy, causing the undrained stress path to show increase in mean effective stress p in Fig. 6a. Upon reverse loading, D ir increases significantly due to the combined effects of Eq. (16) and fabric anisotropy in Eq. (17), while D re follows Eq. (15) and becomes strongly contractive. D re quickly falls back to zero during reverse loading, and D once again becomes the same as D ir until the stress state reaches beyond the reversible dilatancy surface again.

Liquefaction state
A distinct feature of the models by Zhang and Wang [90] and Wang et al. [73] is the formulation for the generation of large yet bounded shear strain at the liquefaction state of zero effective stress, which uniquely captures the cyclic liquefaction behavior of sand. This formulation was proposed based on the insights obtained from test observations of cyclic drained and undrained loading of sand by Shamoto et al. [59] and Zhang [89], and is adopted here. At the liquefaction state, as the increments of p and s are zero, no elastic strains are assumed to be generated. However, dr is not zero and the dilatancy equations are still assumed functional. As the model reaches liquefaction state, it is initially contractive. During the contractive phase, a ''virtual'' elastic volumetric strain is generated due to ðde e v Þ virtual ¼ de v À de p v . This is to offset the inconsistency between the plastic volumetric strain with the volumetric strain boundary condition, which has been shown to be caused by the change in internal structure of the material [70]. As shear strain continues to occur at the liquefaction state, r increases and the model eventually becomes dilative. During the dilative phase, the ''virtual'' elastic volumetric strain must be fully erased before the model exists liquefaction. Hence, according to Eq. (4), sufficient shear strain must occur at the liquefaction state to facilitate the dilatancy process, which accumulates asymptotically with increasing number of load cycles due to the asymptotic accumulation of irreversible dilatancy-induced volumetric strain e p v;ir . This ''virtual'' elastic volumetric strain is fully released when sand leaves liquefaction state, and does not come into the model in nonliquefaction states, where the sum of the elastic and plastic volumetric strain is consistent with the volumetric strain boundary condition.

Model parameters
The proposed model has a total of 17 parameters, as listed in Table 2. In total, 14 of these parameters are identical to those in the isotropic model by Wang et al. [73], the calibration of which can also follow the same procedure as stated in the previous work. Two of these 14 parameters, namely d re,2 and c d,r , can generally use default values of 30 and 0.05, respectively. Three new fabric anisotropy-related parameters, D 1 , D 2 , and c, are introduced in the proposed anisotropic model. Two of the three new parameters, D 1 and D 2 , can be calibrated based on results from triaxial or torsional shear tests on samples with different bedding plane angles, such as the tests by Oda [46] and Yoshimine et al. [88]. The fabric evolution parameter c and the normalized initial fabric tensor should ideally be obtained by measuring the fabric of sand during loading. Though measuring the microstructure of sand is now possible (e.g., [2,66]), the application of such technology is still not yet common. For now, the determination of the initial fabric tensor and the calibration of all three fabric anisotropy parameters for actual sand without microscopic measurement would still need to rely on an indirect trial-and-error process using tests on samples with different bedding plane angles. In this study, the initial fabric tensor F n,in is not considered as a model parameter along the same logic as that of the initial void ratio e not being considered a model parameter, although many existing ACST models use initial fabric as a model parameter. The initial fabric of soil is its inherent property and is independent of the constitutive model. Although direct measurement of fabric in laboratory tests is still difficult compared with void ratio, DEM already allows ''real measurements'' of initial fabric. The current development of DEM can provide a valuable means to evaluate ACST-based models' simulation of the true material fabric, rather than treating fabric as an idealized concept. The detailed model parameter calibration process is described as follows: 1. Monotonic drained tests on samples with bedding plane angle of 0°and different densities and mean effective stress values are first used to calibrate 10 model parameters G 0 , j, h, M, d ir , n p , n d , k c , e 0 , and n following the conventional method described in Wang et al. [73]. Here, the advantage of using ||F n || -A n to reflect anisotropic behavior is highlighted. When the bedding plane angle is 0°, ||F n || -A n is always 0 and does not affect the calibration of these 10 parameters. 2. Monotonic drained tests on dense sand samples with bedding plane angle of 90°are then used to determine the two fabric anisotropy influence parameters D 1 and D 2 . In the proposed model, the peak deviatoric stress ratio g f ;90 and maximum volumetric strain in contraction e v;max;90 for dense sand samples with bedding plane angle of 90°are influenced by D 1 and D 2 , respectively. Therefore, D 1 and D 2 can be determined via g f ;90 and e v;max;90 through an iterative approach: where n is an iteration index. The initial iterative values can generally be taken as . In most cases, adequate convergence is achieved after two to three iterations, i.e.,   The original isotropic model has been validated against a wide range of tests on sand of different void ratios and under different confining stresses [70,73], allowing for the evaluation of the proposed anisotropic model performance in this study to focus on the influence of fabric anisotropy. The model is first scrutinized against DEM numerical monotonic drained triaxial tests on samples with different initial void ratio and bedding plane angle, including the four presented in Sect. 2, and then against undrained triaxial tests on samples with different bedding plane angles and undrained cyclic simple shear tests, on both macro-and microlevels. Though several plasticity models incorporating fabric tensors have been developed, quantitative validation of their fabric-related formulations on the micromechanical level has rarely been reported, especially under cyclic loading. Undrained and drained monotonic tests [46,88] and undrained cyclic tests [12,28] on various sand samples are also simulated, validating the model against physical tests. Numerical analysis is also performed to investigate the influence of fabric anisotropy on the liquefaction behavior of sand.

Validation against DEM tests
The four DEM constant p = 100 kPa triaxial tests on samples with e in = 0.657 ± 0.005 and d = 0°, 30°, 60°, and 90°, discussed in Sect. 2, along with two other constant p triaxial tests on samples with e in = 0.693 and 0.690 and d = 0°and 90°, two undrained conventional triaxial compression tests on samples with e in = 0.655 and 0.661 and d = 0°and 90°, and four undrained cyclic simple shear test on samples with e in = 0.655 and d = 0°, and e in = 0.693 and d = 0°are used to validate the model. In the simulations for these tests with various stress paths on samples with various void ratio and initial fabric using the proposed constitutive model, a single set of model parameters is used, while the initial normalized fabric tensor F n of each sample obtained from the DEM tests is used directly as the input of the simulations. The 17 model parameter values used for these simulations are provided in Table 2. Figure 8 shows the overall deviatoric stress ratio-deviatoric strain and void ratio-deviatoric strain relationships from both the DEM tests and the constitutive model simulations for the constant p triaxial tests on samples with e in = 0.657 ± 0.005. The proposed model is able to reproduce the anisotropy in both strength and void ratio evolution observed in the DEM tests. The peak g is greater for samples with lower d. Samples with greater d initially experience stronger contraction. Both the stress and void ratio converge to each respective unique value at the critical state.
The equivalent shear modulus dq/de q and dilatancy D for these DEM tests and constitutive model simulations on the four samples with e in = 0.657 ± 0.005 are plotted against the deviatoric stress ratio g in Fig. 9. For the equivalent shear modulus, the constitutive model simulations agree with the DEM test results in terms of samples with higher d having lower dq/de q . The model simulations show faster degradation of equivalent shear modulus with respect to the g, which is due to the specific formulation of the plastic modulus in Eq. (10). For the dilatancy, the simulated results agree well with the DEM results. Both the simulated and DEM results show greater initial contraction for samples with higher d. Also, as dilation occurs, the    Fig. 9 Shear modulus dq/de q and dilatancy D results from DEM numerical tests and constitutive model simulations of constant p = 100 kPa triaxial loading on four samples with e in = 0.657 ± 0.005 and bedding plane angles of 0°, 30°, 60°, and 90°: DEM numerical test results for a shear modulus dq/de q versus deviatoric stress ratio g; b dilatancy D versus deviatoric stress ratio g; constitutive model simulation results for c dq/de q versus g; d D versus g constitutive model is able to reflect the peak of ||F n || reaching beyond 1 for dense sand. It should be acknowledged that the simulation of fabric evolution based on the relatively simple form of Eq. (11) is far from perfect, especially for the results in Fig. 10b, d. Recently, through detailed analysis of DEM test data, Wang et al. [74] showed that the contact normal fabric evolution is not only dependent on its own state but also is affected by other particle-scale features, including particle and void orientations. More realistic simulation of fabric evolution can be achieved with formulations incorporating a few extra model parameters. However, considering the current state of the art for fabric measurement and model calibration, the current model still adopts the simple formulation in Eq. (11), to limit the number of model parameters.
Corresponding to the tests with e in = 0.657 ± 0.005 and d = 0°, 30°, 60°, and 90°, two other constant p triaxial tests on looser samples with e in = 0.693 and 0.690 and d = 0°a nd 90°are also simulated using the model. Using the same set of parameters, the model is also able to simulate the behavior of the looser samples under triaxial loading, which exhibit lower peak g and stronger initial contraction (Fig. 11). The anisotropic behavior due to different initial d is again well reflected.
Two undrained conventional triaxial compression DEM tests on the dense samples, with e in = 0.655 and 0.661 and d = 0°and 90°, are also conducted and simulated, again with the same set of model parameters. Undrained loading is achieved in DEM through applying constant volume constraint. The simulation from the constitutive model can capture the difference in stress path and stress-strain relationships for the two samples with different bedding plane angles to a good extent (Fig. 12). Stronger initial contraction is observed for the 90°sample, evident from the more significant initial decrease in p. The simulations show slight underestimation of the deviatoric stress compared with the DEM results. These results show that the proposed model is not only valid for drained triaxial tests on samples with various void ratio and initial fabric, but also for undrained tests as well.
Few anisotropic models have been validated against cases of undrained cyclic loading. Here, an undrained cyclic simple shear DEM test on a sample with e in = 0.655, d = 0°, and ||F n,in || = 0.43 is conducted and used to validate the capability of the model in capturing the cyclic liquefaction behavior of anisotropic sand, on both macroscopic and microscopic levels. The cyclic simple shear DEM tests are conducted on samples isotropically consolidated at p = 100 kPa, and the bedding plane angle d refers to the angle between the deposition direction and the vertical direction, and shear is then applied to the horizontal plane. Figure 13 shows the stress path, stress-strain relationship, and fabric evolution during undrained cyclic loading in both DEM test and constitutive simulation. The proposed constitutive model is able to capture the decrease in effective stress during undrained cyclic loading, eventually leading to initial liquefaction and the ''butterfly orbit'' of the stress path after similar number of load cycles as that in DEM (Fig. 13a, e). The accumulation of large shear strain at liquefaction during each load cycle is well reproduced by the constitutive model, though generation of shear strain in the constitutive simulations is slower (Fig. 13b, f). Figure 13c, g also zooms in on the shear strain in the first six load cycles. The key features of fabric evolution are also surprisingly well simulated by the constitutive model (Fig. 13d, h), considering that the fabric evolution formulation has yet to be directly compared to that in actual physical or numerical tests.
In Fig. 13d, h, the contact normal fabric norm ||F n || is multiplied by the sign of F zz -F xx -F xz to reflect not only the value of ||F n || but also the direction of the fabric tensor, for visual clarity, especially after initial liquefaction. Prior to initial liquefaction, relatively mild evolution of fabric is observed in both DEM test and constitutive simulation, where the oscillation of fabric in simulation results within each cycle is less significant than that in the DEM results. After initial liquefaction, the fabric tensor is found to evolve abruptly at the liquefaction state within each loading cycle, from one direction to another, reflected by the oscillation of ||F n ||sign(F zz -F xx -F xz ) from positive to negative and vice versa. As with the development of shear strain, the evolution of fabric in cycles after initial liquefaction predicted by the constitutive model is slower than that in the DEM test and shows slower achievement of the saturated level. These quantitative differences between Further undrained cyclic simple shear DEM test is conducted on the sample with e in = 0.655, d = 0°, and ||F n,in || = 0.43 for two other cyclic stress ratios (CSR = s max /p 0 0 ), and on another sample with e in = 0.693, d = 0°, and ||F n,in || = 0.41, and are simulated using the proposed model. Figure 14 shows the stress path (sp) and stress-strain relationship (s -c) during undrained cyclic loading in both DEM test and constitutive simulation. The increase in number of load cycles needed to reach initial liquefaction with decreasing CSR is well reflected, while the amplitude of postliquefaction shear strain is unaffected by CSR, as was found in Wang et al. [72]. The number of load cycles required to reach initial liquefaction for the sample with e in = 0.655, d = 0°, and ||F n,in || = 0.43 under CSR of 0.1, 0.2, and 0.3 is 46, 8, and 4 in DEM, and 40, 8, 4 in constitutive modeling, respectively. Decrease in liquefaction resistance and increase in postliquefaction shear strain amplitude with increase in void ratio is also well modeled.

Validation against laboratory tests
Although DEM test results may show quantitative differences to real sand due to the difference in particle shape and contact, they can be qualitatively representative of sand and offer the benefit of direct measurement of fabric quantities. This benefit allows for the direct validation of fabric anisotropy-based constitutive models. With the aid of DEM tests, direct evaluation of the dependency of shear modulus and dilatancy on fabric anisotropy, and the evolution of fabric can be achieved, as in this study. Nonetheless, it is still very important to apply the model to real sand. In the case of physical tests, there are not enough data corresponding to all the DEM test stress paths for the same sand and no corresponding measurement of microscopic fabric quantities. Consequently, the macroscopic results of several classical test results on two types of sand are simulated, using model parameters in Table 2. Both undrained monotonic and cyclic torsional shear tests on the widely used Toyoura sand are first simulated using the same set of model parameters.
The simulations for undrained monotonic loading on Toyoura sand with e = 0.825 ± 0.004 and bedding plane angle d of 15°, 30°, 60°, and 75°are presented in Fig. 15. The tests were conducted by Yoshimine et al. [88] using a hollow cylinder torsional apparatus, with p in = 100 kPa and b = (r 2 -r 3 )/(r 1 -r 3 ) = 0.5. The model parameters for the simulations are listed in Table 2. As the fabric information is not available, the initial fabric tensor orientation is assumed to be d ? 90°, and the initial fabric norm is assumed as ||F nin || = 0.4, similar to that those in the DEM tests. The three fabric-related model parameters are estimated. Both the simulated stress paths and stress-strain curves for the different samples agree well with the test data. The sample with d = 15°shows mostly dilative behavior, while the d = 75°shows strong initial contraction, almost reaching static liquefaction. Figure 15e, f also plots the evolution of dilatancy D and contact normal fabric norm ||F n || obtained from the constitutive simulations. As the samples in these tests are relatively loose, the samples are initially strongly contractive, and dilatancy D is initially greater for samples with greater bedding plane angles. The evolution of contact normal fabric of the different samples also shows different features, ||F nin || gradually increases toward 1 for the sample with d = 15°, while it initially decreases and then increases toward 1 for the sample with d = 75°. It should be noted that by introducing the dilatancy in the fabric evolution Eq. (11), for a strongly contractive loose sample, the initial fabric increment direction may become opposite to the current loading direction. Whether this feature is realistic still requires justification from further test evidence.
Using the exact same model parameters and initial fabric state, the undrained monotonic triaxial tests reported by Yoshimine et al. [85] are also simulated. These tests are conducted under different initial confining pressures, under both triaxial compression (TC) and extension (TE). Overall good agreement between simulation and test results is achieved for both initial confining pressures. Stronger contraction under TE is observed in both simulation and test results, with the TE tests reaching static liquefaction, and generating significant shear strain at liquefaction state (Fig. 16c, d). The TC results show initial contraction followed by subsequent dilatancy and increase in p (Fig. 16a,  b). The constitutive model simulations overestimate the increase in deviatoric stress with increasing deviatoric strain (Fig. 16c, d), which may be due to the simple g(h) function in Eq. (9), which does not introduce any parameters other than M. If necessary, an alternative function form may be used for g(h).
The simulations for undrained cyclic hollow cylinder torsional shear tests on Toyoura sand with e = 0.846 under CSR = 0.16 and e = 0.768 under CSR = 0.182 are presented in Fig. 17, and the tests were conducted by Chiaro et al. [12] and Koseki et al. [28], respectively. The model parameters are the same as those used to simulate the monotonic tests by Yoshimine et al. [88], listed in Table 2. The initial fabric norm is again assumed as ||F nin || = 0.4, as the sample preparation technique was largely consistent for these tests. The good agreement between simulation and test results under undrained cyclic loading again highlights the model's ability to capture both monotonic and cyclic behavior of sand, reproducing the decrease in effective stress during undrained cyclic loading and also the postliquefaction shear strain generation. Both test and simulation results exhibit similar increase in number of load cycles to liquefaction with decreasing void ratio and decrease in postliquefaction shear strain with decreasing void ratio.
A drained cyclic torsional test on Toyoura sand with e = 0.730 reported by Zhang and Wang [90] is also simulated using the model with the same set of model parameters listed in Table 2. The initial fabric norm is assumed to be the same as those of the previous tests as ||F nin || = 0.4. As shown in Fig. 18, the stress-strain loop is well captured using the model and its parameters and the accumulation of contractive volumetric strain from constitutive simulation agree well with the laboratory test results.
It should be noted that different batches of even the standard Toyoura sand may show large inconsistencies in terms of material behavior and different preparation methods would result in significantly different initial state      [88] of the samples. Therefore, it should not be expected that this exact same set of parameters and initial state values will always be applicable to all existing tests using Toyoura sand.
The drained triaxial compression tests on a quartz sand at 65% relative density and p in = 100 kPa with bedding plane angle d of 0°, 30°, 60°, and 90°by Oda [46] are also simulated using the proposed model (Fig. 19), to evaluate  the model's performance for drained loading, using the parameters in Table 2. Again, the initial fabric tensor orientation is assumed as d ? 90°, while ||F nin || = 0.4. Similar to the simulations of the DEM triaxial tests, the model is able to capture the anisotropic deviatoric stress and volumetric strain development for the four different samples.

Influence of fabric anisotropy on liquefaction behavior
As exhibited in the validations against DEM data, one of the main strength of the model framework adopted in this study is its capability in reflecting the cyclic and especially liquefaction behavior of sand. Here, undrained cyclic simple shear simulations are conducted to investigate the influence of fabric anisotropy on liquefaction resistance and liquefaction-induced shear strain. Simulations on samples of e in = 0.655 with bedding plane angle d = 0°a nd initial fabric norm ||F nin || = 0.01 and with d = 45°and ||F nin || = 0.43 are carried out and compared with the results for the simulation in Sect. 4.1 on the sample with e in-= 0.655, d = 0°, and ||F n,in || = 0.43.
The shear stress s versus mean effective stress p and shear stress s versus shear strain c results for the simulations are plotted in Fig. 20. These simulation results, along with that in Fig. 13, clearly depict the significant influence of fabric anisotropy on the liquefaction behavior of sand. Comparison of Fig. 20a, b with Fig. 13d, e indicates that greater fabric anisotropy can lead to a significant decrease in liquefaction resistance. The number of load cycles needed to reach initial liquefaction reduces from 12 to 7 as ||F nin || increases from 0.01 to 0.43. These results corroborate the observations from laboratory tests [48,81,86,87] and DEM numerical tests [72,74,76]. The results from the simulations on samples with the same ||F nin || but different d show that the fabric anisotropy orientation significantly affects the shear strain accumulation (Fig. 20c, d and  Fig. 13d, e). When the initial fabric orientation is symmetric with respect to the cyclic shearing direction, the postliquefaction shear strain accumulation is also symmetric with respect to zero (Fig. 13e). When the initial fabric orientation is asymmetric with respect to the cyclic shearing direction, accumulation of postliquefaction shear strain in one direction occurs (Fig. 20d).
The influence of initial fabric on liquefaction resistance is a result of its influence on dilatancy. When the fabric tensor F n and loading direction n are not proportional, greater ||F n || can lead to greater irreversible dilatancy, according to Eq. 17. Figure 21 shows the dilatancy for the two simulations with the same d = 0 but different ||F nin ||. Increase in ||F nin || causes the positive dilatancy to increase, indicating stronger contraction tendency.

Conclusions
Based on existing understanding of the anisotropic behavior of sand from laboratory tests and micromechanical DEM numerical test data, an anisotropic plasticity model incorporating fabric evolution is proposed in this ig. 19 Stress and strain results from laboratory tests and constitutive model simulations of conventional drained triaxial compression loading on for four samples with bedding plane angles of 0°, 30°, 60°, and 90°: laboratory test results for a deviatoric stress r 1 -r 3 versus axial strain e 1 ; b volumetric strain e v versus axial strain e 1 ; constitutive model simulation results for c r 1 -r 3 versus e 1 ; d e v versus e 1 . Laboratory test results from Oda [46] study. The development and validation of the proposed model using macro-and microlevel data from DEM numerical tests on granular materials provides a paradigm for the formulation and evaluation of models incorporating micromechanical features of sand, aiding the endeavor to bridge the gap between micro and macro.
The proposed model builds on an isotropic model that can achieve a unified description of sand of different conditions from pre-to postliquefaction under both monotonic and cyclic loading. A normalized fabric tensor and its evolution rule are introduced in the model. By incorporating the fabric tensor and its invariant with the loading direction in the formulation for plastic modulus and dilatancy, the model is able to reflect the influence of fabric on the strength, shear modulus, and dilatancy of sand.
The model is directly validated against DEM test results, using the fabric information obtained from the numerical tests, exhibiting good performance for both the macrolevel stress-strain relationships and the microlevel fabric features. This validation process using DEM numerical test results provides a new means of validating micromechanical-based constitutive models on both micro-and macroscales. Validation of the model against a wide range of different laboratory tests is also carried out, showing good agreement with the laboratory test results, highlighting the model's adaptability to different stress paths.
The influence of fabric anisotropy on the undrained cyclic behavior of sand is investigated using the proposed model. Initial fabric tensor norm and orientation are found to significantly affect the liquefaction resistance and postliquefaction shear strain, respectively. Greater initial fabric norm leads to weaker liquefaction resistance, which is due to the influence of fabric anisotropy on dilatancy, agreeing with observations from laboratory and DEM tests.
In the proposed model, a relatively simple formulation for the evolution of fabric is currently used. Though the simple formulation is able to capture the overall fabric evolution observed in DEM tests, distinct difference still exists between the model simulation and DEM test fabric evolution results. Further studies should be carried out to investigate the evolution of and the relationship between fabric tensors based on various particle-level characteristics, such as the contact normal, particle orientation, and void orientation, under both monotonic and cyclic loading. More accurate formulations for fabric evolution and its influence on the mechanical behavior can be developed based on such analysis. Also, the influence of intermediate stress coefficient b on the evolution and the effects of fabric should be further investigated and considered in constitutive model formulations. Although the current model mainly focuses on the influence of fabric anisotropy on the shear strength, shear modulus, and dilatancy of sand, it should be acknowledged that fabric anisotropy can also affect other behaviors as well, which should also be explored in future studies. In applications where the small strain behavior of sand is important, anisotropic elasticity should be considered. For loading with a significant rotation of principal stress axes, fabric anisotropy induced noncoaxiality between strain increment and total stress may also be important.