The deformation and failure of strip footings on anisotropic cohesionless sloping grounds

Footing foundations are sometimes built on sloping grounds of natural sand which is highly anisotropic. The anisotropic mechanical behaviour of sand can significantly influence the bearing capacity of a foundation and the failure mechanism of its supporting slope. Neglecting sand anisotropy may lead to overestimated bearing capacity and under‐design of foundations. A numerical investigation on the response of a supporting slope under a strip footing is presented, placing a special focus on the effect of sand anisotropy. A critical state sand model accounting for fabric evolution is used. The nonlocal method has been used to regularize the mesh‐dependency of the numerical solutions. Predictions of the anisotropic model on the bearing capacity of strip footings on slopes are validated by centrifuge test data on Toyoura sand. Compared to the centrifuge test data, an isotropic model may overpredict the bearing capacity of the footing by up to 100% when the model parameters are determined based on test data on a horizontal bedding plane case. When the isotropic model parameters are determined based on test data where the bedding plane is vertical, the predictions of bearing capacity can be improved for some cases but the settlement at failure may be significantly overestimated. The soil body tends to move along the bedding plane upon the footing loading due to the non‐coaxial strain increment caused by fabric anisotropy. The slip surface appears to be deeper with lower bearing capacity when the preferred soil movement direction caused by bedding plane is towards the slope.

Bedding plane A F I G U R E 1 Schematic illustration of major principal stress directions of sand elements relative to the bedding plane axis (vertical) along the failure plane in an anisotropically deposited slope (after Uthayakumar and Vaid, 1998) on the shear strength and dilatancy of sand. [4][5][6][7][8] Under otherwise identical loading conditions (e.g., the same principal stress and void ratio), a sand sample shows the highest shear strength and most dilative response when the principal stress direction is perpendicular to the bedding plane. In the case of a strip footing sitting on a transversely isotropic slope, Figure 1 illustrates the deviation of the direction of major principal stress 1 with the bedding plane axis (denoted by the angle ) for sand elements along the slip surface in the slope. 9 The variation of with depth is apparent. It is typically negative (e.g., 1 tilts to the left of the vertical) at the top close to the footing, nearly vertical at a point under the footing (Element A) and almost horizontal at the toe of the slope (Element B). If the soil is assumed isotropic and testing conditions similar to that of Element A are used for determining its peak strength, the bearing capacity could be significantly overestimated. On the other hand, it may cause underpredictions of the bearing capacity if the peak friction angle for element B is used. 10 The importance of anisotropy on the performance of geotechnical foundations has indeed been recognised. [10][11][12][13] There have been experimental investigations of footings on anisotropic level ground. 4,[14][15][16][17] Both small-scale model tests and centrifuge physical tests have shown that the bearing capacity of a strip footing on sand is the highest when the deposition plane is horizontal, whereas the minimum bearing capacity, which can be 25% lower than the maximum one, is observed when the deposition plane is vertical. 12,15 The failure mechanism of footings on a sloping ground may differ substantially from that on a level ground (e.g., Meyerhof 1 ; Kimura et al. 12 ; Graham et al. 18 ;Narita,and Yamaguchi 19 ). Observations derived from the case of level grounds may not apply for footings built on slopes.
There have been a limited number of numerical studies offering an improved understanding of how fabric anisotropy affects the response of strip footings on sloping grounds. Huang and Tatsuoka 10 varied the friction angles with the direction of 1 to consider the fabric effect in their finite element study on the stability of sand slopes (see also Stockton et al. 20 ). Aghajani et al. 21 employed the artificial neural network to account for the strength anisotropy of soils in their analysis of the stability of sand slopes. Despite the improvemens they have brought on the topic, these existing numerical attempts cannot realistically capture the mechanical origins of anisotropic soil behaviour, such as the interplay between fabric and dilatancy, that underlies the deformation and failure of a footing foundation and dictates the changes of bearing capacity of a footing with bedding plane orientation and slope geometry. 8,17,22 Indeed, soil dilatancy has long been known to be closely related to soil strength and hence its engineering performance. Manzari and Nour 22 have shown that the stability of a slope can be significantly enhanced when the soil becomes more dilative, even though the peak friction angle does not change at all. Moreover, fabric anisotropy may cause non-coaxial strain increment with respect to the loading direction in sand. [23][24][25] Based on finite element modelling, Yuan et al. 25 have demonstrated that the non-coaxial strain increment has a detrimental effect on the bearing capacity of strip footings on level grounds. A first key step towards better understanding the role of fabric anisotropy is the adoption of a proper constitutive model in numerical simulation of the footing problem on sloping ground. Good candidates include early sand models accounting for the inherent fabric anisotropy (e.g., Pastor 26 ; Pastor et al. 27 ; Li and Dafalias 28 ; Dafalias et al. 29 ; Manzanal et al. 30 ) and recent models that account for fabric and fabric evolution. 23,31 The recent anisotropic models have been established within the anisotropic critical state theory, in which the effect of fabric and fabric evolution on sand behaviour is considered in conjunction with a well-defined critical state accounting for fabric anisotropy. 32 These models are theoretically more rigorous and practically more robust for predicting various soils behavior involving anisotropy. Indeed, they have also been successfully used in modelling the response of strip footings on level grounds with different bedding plane orientations. 17 This study presents a new finite element investigation on the effect of fabric anisotropy on predictions of strip footing on a sloping sand ground. The centrifuge tests on dry Toyoura sand slopes reported by Kimura et al. 12 will be used for model setup and model parameter calibration to facilitate better comparison. The critical state sand model accounting for fabric evolution proposed by the authors (Gao et al. 8 ; Gao et al. 23 ) is employed for the numerical study. At the material point level, the model has been demonstrated to provide satisfactory descriptions of the fabric effect on strength, dilatancy and the non-coaxial response of sand under proportional loading. Since this model considers strain-softening of dense sand, a nonlocal implementation of the model is employed in the finite element method (FEM) to regularize the meshdependency of the finite element solutions. [33][34][35] 2 CONSTITUTIVE MODEL AND NUMERICAL IMPLEMENTATION

The original constitutive model
The constitutive model has been developed based on the anisotropic critical state theory 32 in consideration of the fabric and fabric evolution during loading. 8,23 The model features with formulations incorporating the effect of evolving fabric on dilatancy, plastic hardening, and the non-coaxial response of sand under proportional loading. Details of the model formulations have been presented in previous studies. 8,23 Only key equations are presented in the following to facilitate the discussion of this study. The yield function of the model follows that used by Li and Dafalias 28 is the stress tensor, = ∕3 is the mean effective stress, is the Kronecker delta (= 1 for = , and = 0 otherwise), is the hardening parameter and ( ) is an interpolation function which describes the variation of critical state stress ration with the Lode angle of . 28 The plastic potential function in the space presents the following expression 8 where ℎ is a model parameter, is the anisotropic variable 8,23,32 and is calculated based on the current stress state and . is defined as a joint invariant of the loading direction tensor and fabric tensor The plastic potential function in Equation (2) is used to obtain the direction of plastic deviatoric strain increment as follows where is the loading index; ⟨⟩ are the Macaulay brackets which make ⟨ ⟩ = for > 0 and ⟨ ⟩ = 0 for ≤ 0. The term involving in Equation (3) enables the model to capture the non-coaxial response of sand due to fabric anisotropy. 23,36,37 The total plastic strain increment is where is the plastic volumetric strain increment and is the dilatancy equation expressed as follows where 1 , , and are model parameters, is the dilatancy state parameter 32 and is the state parameter. 38 The evolution of for is expressed as where ℎ 1 , ℎ 2 , and are model parameters and is the elastic shear modulus. 8 It is assumed that the fabric tensor evolves with the plastic deformation. It becomes co-directional with the loading direction and reaches a magnitude of one at the critical state 8,23,32 where is a model parameter. Note that the degree of anisotropy (= √ ) at the critical state is assumed to be one for simplicity. The fabric tensor used here is a phenomenological term that is not directly related to the fabric or particle characteristics of sand, such as contact normal distribution, or particle orientation. 39 More discussion on the fabric tensor for the anisotropic critical state theory can be found in Li and Dafalias. 39 The model can be readily recovered to an isotropic model as a special case by neglecting the fabric anisotropy in the entire model formulations. Specifically, it becomes an isotropic model when and the model parameters associated with fabric anisotropy ( ℎ , , ℎ 2 , and ) are set 0. The parameters for the anisotropic model are determined using the test data on sand with different bedding plane orientation with respect to the major principal stress direction. 8,23 The parameters for the isotropic model have been determined based on the test results in which the bedding plane is horizontal (perpendicular to the major principal stress direction) and vertical (parallel to the major principal stress direction), respectively. In the following sections, these two sets of isotropic model parameters will be called the horizontal and vertical parameters, respectively. Model parameters for both the anisotropic and isotropic models are summarized in Table 1

Nonlocal formulation of the constitutive model
The model presented above employs a non-associated flow rule and considers the strain-softening response of sand. When this model is used in finite element modelling of a boundary value problem, the solution will be mesh-dependent, especially when localized deformation bands emerge in the post-peak loading stage. 8,17,34,35 Mesh dependency apparently affects the objectivity of predictions of the bearing capacity of foundations and shear band thickness and orientation that manifest the underlying failure mechanisms. The issue of mesh dependency can be addressed by various regularization approaches, including the nonlocal method to be adopted here. 33,40,41 In previous studies, the partially nonlocal approach has been widely used, wherein some of the state variables governing the strain-softening are assumed to be a function of the nonlocal plastic strain. 34,35 For the present model, the strain-softening is governed by several state variables, including , , and . Following Mallikarachchi and Soga, 41 the increment of void ratio is assumed to be nonlocal as below where positive means volume contraction and is the nonlocal volumetric strain increment expressed as where is the total number of integration points used for nonlocal averaging, , , and represent the weight function, volume, and local volumetric strain increment of the i-th integration point in a typical FEM discretization. The Galavi and Schweiger 34 weighting function expressed as below is used where is the internal length dependent on the soil particle size, is the distance between the current integration point and the -th integration point used for calculating the averaged value in Equation (12). It is found that the weighting function above gives better regularization results than the Gaussian normal distribution function. 35,41 Note that Equation (13) is the same as that in Summersgill et al. 35

Model implementation and calibration
The nonlocal model has been implemented in ABAQUS using the explicit stress integration method 42,43 via the user subroutine UMAT, where the large strain formulation proposed by Hughes and Winget 44 has been employed. In implementing the nonlocal model, is calculated using Equation (12) at the beginning of each increment and then used to update the void ratio at the end of the increment. 41 The remaining part of the UMAT is the same as that for the local model. Since the effect of the weight function at > 2 is very small, only the integration points within the radius of 2 are considered for Equation (12). 34 In the following simulations, the model prediction for the response of strip footings on a cohesionless soil slope will be validated first using the centrifuge test data on dry Toyoura sand in which the sand anisotropy has been studied. 12 The limitation of the isotropic model in predicting the bearing capacity will be discussed. The anisotropic model will be further used to examine how fabric anisotropy affects the strip footing behavior.

SIMULATION OF STRIP FOOTINGS ON SLOPING GROUNDS
A group of centrifuge tests on strip footings placed near the crest of a sand slope have been reported by Kimura et al. 12 These tests were designed to investigate how sand density and anisotropy affect the bearing capacity of strip footings. They also performed footing tests on level ground with horizontal and vertical bedding planes for comparison. These tests are ideal for validating the model prediction of the problem of strip footing on sloping grounds. The same prototype of the footing ( = 1.2 m) and slope height ( = 2.4 m) is employed for the finite element simulation in this study ( Figure 2). The slope angle varies between 25 o and 35 o to explore its influence on the simulation results. The distance between the footing and slope crest is , with varying between 0 and 2. The effective unit weight of dry Toyoura sand is set to be ′ = 15.68 kN∕m 3 . The elastic stiffness and plastic modulus of the model is dependent on the mean effective stress. The model gives zero elastic stiffness and plastic modulus when mean effective stress is 0, indicating that the soil can collapse. It is thus difficult to simulate free-surface problems using this model. In most studies, either a small cohesion or confining pressure is applied on the free surface to prevent soil failure at low mean effective stress. In this study, a confining pressure of 2kPa is chosen to make sure that the simulations can continue after the bearing capacity is reached ( Figure 2). The initial stress state is generated using the gravitational loading method. Specifically, the soil is first assumed to be an elastic material with a large elastic modulus and Poisson's ratio of ν e =0.286. The gravity and surface loading (2 kPa) are then applied to generate the initial stress state. For the soil beneath the slope (y < 4.6 m in Figure 2), this gives a lateral earth pressure coefficient K 0 = ν e /(1-ν e )=0.4, which is close to the value reported in Okochi and Tatsuoka. 45 Vertical displacement is prescribed on the footing to simulate a rigid and rough footing case, while its horizontal movement is restrained. Standard boundary conditions in a typical slope stability analysis are applied to the simulated domain, that is, both the horizontal and vertical movements are fixed at the bottom of the domain, whereas only horizontal displacement is restricted along the left and right vertical boundaries. Figure 2 shows the simulated domain and mesh size for the slope with = 30 • . The same height and width of the simulation domain are used for two other cases of the slope angle with = 25 • and = 35 • . The adopted mesh has been based on a balanced consideration of computational efficiency and accuracy, but can definitely be improved, for example, through techniques such as mesh refinement or adaptive mesh. While we have chosen to use a uniform mesh to render better mesh quality, alternative ways include using finer mesh around the footings and coarser mesh at the boundaries to render more realistic capturing of shear band thickness. 34 Notable issues associated with mesh refinement techniques, including computational efficiency and convergence issues, have to be taken care of. Meanwhile, the extended finite element method (XFEM) method can be used to avoid mesh refinement for such simulations. 46 An initial degree of anisotropy 0 = 0.4 is adopted. The initial void ratio is 0 = 0.66 (relative density of = 90%). The internal length must be chosen for the simulations using the nonlocal model. Practical calibration of the internal length for a boundary value problem is a challenging task. 34,41 When the weighting function expressed by Equation (13) is used, the simulated shear band thickness is close to . 34 It has been observed that the real thickness of shear band in sand is about 10 times of the mean particle size 50 . Since the 50 of Toyoura sand is about 0.2 mm, = 2 mm should be used if realistic prediction of the shear band thickness is required. But this would lead to a very small mesh size because the maximum mesh size has to be smaller than ∕2 to guarantee that there are sufficient integration points for calculating the nonlocal variables. 35 Such a small mesh size will cause numerical convergence issues and significantly increase the computational time. 41 Therefore, the value of is typically chosen based on the size of the boundary value problem, which is much bigger than 10 50 . 35,41 The influence of on the − relationship predicted by the anisotropic model is shown in Figure 3, where is the vertical settlement and is the vertical loading on the footing. Eight-noded quadratic elements ) and the corresponding settlement is higher at bigger , due likely to relatively wider shear band and more even stress distribution caused in the soil by the use of bigger . 35 is chosen at 0.5 m to avoid using rather small meshes in this study. Indeed, the nonlocal anisotropic model gives a satisfactory prediction of the bearing capacity with = 0.5 m, which can be seen in Figure 5. Figure 4 shows the effect of mesh size on the original and nonlocal model predictions. The solution of the original model is highly mesh-dependent and does not converge as the mesh size decreases ( Figure 4A). Though the − relationship predicted by the nonlocal model is not completely mesh-independent, but the bearing capacity and settlement at failure which are of importance for a practical design show a reasonably small variation with the mesh size. Besides, the − relationship shows convergency as the mesh size decreases ( Figure 4B). The mesh sensitivity of the nonlocal model solution could be further reduced if nonlocal and are used, which may inevitably increase the complexity of the model formulations and its implementation. Notice that the averaging [Equation (12)] has to be done for all the nonlocal variables at the beginning of each increment, and therefore, more computational time is needed if more nonlocal variables are used. Besides, alternative regularization methods may be employed to attain more mesh-independent results. For instance, the XFEM has been specifically developed to ease the difficulties in solving problems with strong discontinuity in the displacement field (e.g., strain localization, and crack propagation), which can be effective for the problem under consideration. 46 Since the − relationship shows rather small change when the mesh size is smaller than 0.25 m, the mesh shown in Figure 2 (about 0.25 × 0.25 m) is thus chosen for all the simulations in this study.

3.1
Predictions by the anisotropic model Figure 5 shows a comparison between the predicted and measured bearing capacity factor for the tests on slopes with different slope angle . The model prediction is in good agreement with the test data when ≥ 1. However, the predicted is higher than the measured one at = 0 and = 0.5. Possible reasons accounting for such discrepancy include: (a) the rough boundary condition applied to the footing imposes more constraints than in the real test condition, leading to higher predicted . The rough boundary condition allows no relative movement between footing bottom and the sand immediately underneath, causing overprediction of . Indeed, as mentioned by Kimura et al., 12 any relative movement between sand and footing would reduce the . It is expected that the overprediction caused by this boundary condition is more significant when the footing is closer to the edge of the slope crest, where soil may be more easily mobilized due to less lateral support. A proper simulation of the relative movement between sand and footing needs the determination of the friction coefficient between them which is not available in the literature. (b) Real sand samples in the tests of Kimura et al. 12 may have not been uniform. There could be some locations with larger void ratio which has a detrimental effect on the bearing capacity. 47 If the samples were uniform with the same void ratio, the test data could be closer to the model predictions. (c) The initial effective stress state generated using the gravity loading method may not be realistically close to that in the tests. The difference between the real and assumed initial stress states is expected to be bigger near the slope, as the slope was formed by excavation in the tests, 12 whereas the gravity loading method used in this study may be more appropriate when the slope is formed by vertical deposition. 48 Such a difference in the initial stress state may have more influence on the behavior footings which are close to the slope crest. (d) The constitutive model may need to be improved. The initial stress state is anisotropic in this problem. But it is shown in Gao and Zhao 49 that this model gives less satisfactory prediction for the stress-strain relationship of sand when the initial stress state is anisotropic. The plastic flow rule used in the general plasticity models can be employed to improve the model prediction of sand behavior with initially anisotropic stress state. 26,50 The − curves may partially reflect the load-displacement response of the footing. Predictions of the − relationship for all the three slopes are shown in Figure 6. No full − curves for these footings are available from the study Kimura et al. 12 for comparison. The settlement at failure has been used as an important parameter for geotechnical design. 51 The predicted is 0.075 − 0.09 (Figure 6), which is within the typical range reported in existing studies. 12,51,52 Moreover, at the same , is higher when the slope angle is smaller, because a steeper slope may offer less horizontal confinement to the movements of the soil beneath the footing (Figures 5 and 6).
In the following discussion of slope failure, the commonly used term slip surface in slope stability analysis is employed to denote a localized shear band breaking through from the crest to the slope surface or toe. Shear bands, in terms of cumulative shear strain, initiating from the right side of the footing and propagating through the slope surface or toe over the post-peak loading stage, have been observed in most of the simulations. Full slip surfaces with significant strain localization have not been observed in most cases because the simulations were terminated at a relatively moderate settlement since the force equilibrium cannot be achieved due to severe distortion of elements at the footing edges. Figures. 7(A-C) shows the contour of cumulative shear strain and slip surfaces after the state for three simulations in the case of = 30 • . Clear slip surface is observed at =0 only. The slip surface for the other case indicate those can be expected at a larger footing settlement. 17 A trapped wedge with small shear strain under the footing can be observed when > 0, which is consistent with previous studies. 1,53,54 Graham et al. 1 have shown that there could be two slip surfaces beneath the footing. But only one slip surface is observed for each simulation, which starts at the footing edge and ends at the slope toe. It has been found that Figure 7D shows the increment of second-order work for the simulation with = 30 •  Figure 7A), which is in agreement with the findings of previous studies. 46,[55][56][57] No obvious localization of negative second-order work increment has been observed in other cases because full slip surfaces have not developed at the maximum settlement.

Predictions by the isotropic model
The parameters for isotropic models are conventionally determined using the results of tests with horizontal bedding. 8,17 The bearing capacity predicted by the isotropic model (parameters being determined using horizontal bedding) for strip footings on slopes with different is shown in Figure 8. The same FEM mesh in the anisotropic model case is used. It is evident that the predicted by the isotropic model is always higher than the measured value. The maximum and minimum overprediction is about 100% ( = 25 o and 35 o with = 0) and 20% ( = 25 o with = 1.5), respectively. This is consistent with the findings in Chaloulos et al. 17 and Gao et al. 8 The isotropic model can only describe the stress-strain relationship of a typical soil element underneath the footing with a vertical major principal stress direction (e.g., = 0 • in Figure 1), due to the previously mentioned procedure for parameter determination. However, the actual value of varies inside the soil slope. Figure 9 shows the predicted variation of by the anisotropic model for soil elements along the slip surface under a slope with = 35 • ( = 0.5 and = 0.1 ). The value of changes from −18 • at the bottom of the footing to 70 • at the toe of the slope. For all the soil elements with ≠ 0, the isotropic model offers predictions of higher peak shear strength and more dilative response. The overprediction is higher at bigger . One example is given in Figure 10, which shows that the isotropic model can reasonably capture the stress-strain relationship of the sand sample with horizontal bedding plane ( = 0). But it gives the same prediction for sand samples with both horizontal and vertical bedding planes. It is evident that the isotropic model overpredicts the peak deviator stress and volume expansion for the sample with vertical bedding plane ( = 90 • ), and hence offers an overestimation of the bearing capacity of strip footing.
Since the isotropic model does not account for the effect of fabric anisotropy on sand behavior, different sets of model parameters can be calibrated when the data with different bedding plane orientations are used. In this study, tests data for the case of vertical bedding plane orientation have also been used to determine model parameters for the isotropic model ( Table 1). The model prediction for the soil element behaviour is shown in Figure 11. Since the soil shows the lowest shear strength and least dilative responses when the bedding plane is vertical under the same condition of stress state and void ratio, the isotropic model underpredicts the peak deviator stress and volume expansion for the samples once the bedding plane is not vertical (Figure 11). F I G U R E 7 Contours of total shear strain for FEM simulations with = 30 • (A-C) and the increment of the second-order work for the simulation with = 30 • and = 0 (D). The dashed lines indicative of potential slip surfaces are added for better illustration Figure 12 shows the − relationship predicted by the isotropic model with the horizontal and vertical parameters. While − curves predicted using the horizontal parameters show an obvious peak of for all cases, keeps increasing when the vertical parameters are used. The − curves shown in Figure 12B and D are typically observed for loose sand, 12,52 though the sand used in the simulations has a relative density ≈ 90%. Following the previous studies, 17,51 the value of at = 0.1 has been used to determine the bearing capacity predicted by the isotropic model with vertical parameters (Figure 13). When the bearing capacity is determined this way, the isotropic model prediction becomes very F I G U R E 8 Comparison between the bearing capacity predicted by the isotropic model (parameters from horizontal bedding) and centrifuge test results on dry Toyoura sand with different slope angles and distance parameters (data from Kimura et al., 1985): close to the test results in some cases ( Figure 13B). Obvious underprediction of the bearing capacity can be observed at = 25 • with > 1 and = 35 • with = 2.
The above comparison highlights both the limitation of an isotropic model in predicting the behavior of strip footing on an anisotropic slope and the consequence of data selections on model calibration. If the isotropic model is to be calibrated by a horizontal bedding plane case, the overestimated bearing capacity may render a design risker in strength. Calibrations of the isotropic model based on other bedding plane cases may help to improve the predictions on bearing capacity but cannot essentially avoid excessive overestimations on the ultimate settlement at failure. For the dense Toyoura sand studied here, the − curves predicted by the isotropic model with vertical parameters does not show an obvious peak

NUMERICAL INVESTIGATION OF FABRIC EFFECT ON THE BEHAVIOUR OF STRIP FOOTINGS NEAR SLOPES
Predictions of anisotropic and isotropic models on the strip footing on a slope with horizontal bedding plane have been presented in last section. In real cases, the bedding plane can be inclined in any direction. This causes non-alignment of F I G U R E 1 3 Comparison between the bearing capacity predicted by the isotropic model (parameters from vertical bedding) and centrifuge test results on dry Toyoura sand with different slope angles and distance parameters (data from Kimura et al., 1985): fabric anisotropy with the loading direction and hence non-coaxial strain increments, which collectively affect the strain localization in sand and failure mechanism of strip footings. 25,36 In this section, the anisotropic model will be used to analyses how bedding plane orientation and non-coaxial strain affects the behavior of strip footings. orientation is given in Gao et al. 8 At the same , smaller is observed when is bigger ( Figure 14A). When the bedding plane is vertical ( = 90 • ), the predicted is 24% and 19% lower than that for the corresponding case with horizontal bedding at = 0 and = 1.5, respectively. The difference in caused by bedding plane orientation is much smaller when the ground surface is levelled. For the same soil with similar density, this difference is about 10%. 8,12 The − curves for strip footings on the slope with = 45 • and = 90 • are shown in Figure 13B and C. It is evident that the initial slope of the − curve is smaller when is bigger. Similar to the simulations with horizontal bedding, slip surfaces have been observed in some cases. Two examples are given in Figure 15 ( = 0.5, = 45 • , and 90 • ). Comparing to the corresponding case with horizontal bedding (Figure 6A), the depth of the slip surface is shallower when is bigger. This is consistent with the observations on level grounds in Kimura et al. 12 The difference in the and − relationship is caused by the initial fabric and fabric evolution of the soil. Figure. 16 shows the fabric evolution of one element (Element A) which lies on the slip surfaces of all cases with different bedding plane orientations ( = 30; = 0 • , 45 • , and 90 • ). The initial value of the anisotropic variable is smaller when is bigger ( Figure 16A). is dependent on the degree of anisotropy and relative orientation between the soil fabric and the loading direction. Since the initial is the same for these simulations, the difference in the initial is caused by the relative orientation between the soil fabric and the loading direction. The constitutive model gives higher shear strength and more dilative response for bigger when the other conditions are the same. Both higher soil strength and more soil volume expansion have beneficial effect on the bearing capacity of footings. 22 From the perspective of the micromechanics on sand, higher means that the soil fabric is at more optimum state for sustaining the external loading. 23  highest at = 0 • throughout the tests ( Figure 16A), which makes the highest when the bedding plane is horizontal. keeps increasing at = 0 • and = 45 • , because there is no significant rotation of the principal fabric directions. Increase in is primarily due to the increase in for these two cases. When = 90 • , the fabric changes the principal directions to become co-directional with the stress, which leads to decrease in and increase in . 23 If the soil was loaded to larger strain level, would increase again. 23 In all the simulations, and have not reached the critical state value of one.

Effect of initial bedding plane orientation
The simulations above have used positive (major principal stress direction tilting right to the vertical). However, real bedding plane orientations can be negative . The response of strip footings with negative may be different when the bedding plane is not vertical. To demonstrate this, Figure 17 further shows the simulations for footings on sand with = 45 • and = −45 • (slope angle = 30 • ). The value of is higher when is negative. The initial slope of the − curve is similar for positive and negative cases. But the decrease of with is slower after the peak at negative . Comparing to the case with positive , the slip surface is shallower when is negative ( Figure 15A and Figure 18). These differences in the footing response can be explained with the assistance of Figure 19. It has been shown by Gao and Zhao 36 F I G U R E 1 9 The preferred sliding direction of soil with different bedding plane orientations that the soil tends to move along the bedding plane when it is loaded due to the non-coaxial strain increment caused by fabric anisotropy (Figure 19). More discussion on the effect of non-coaxial strain increment on footing behaviour will be provided in the following subsection. When is positive, the bedding plane promotes more soil movement towards the slope, leading to detrimental effect on the bearing capacity and makes the slip surfaces deeper (Figure 19). But the soil tends to be pushed to move away from the slope when is negative ( Figure 19). This makes smaller soil body move towards the slope (shallower slip surface) and increases the bearing capacity of footings. Indeed, the preferred soil movement can be inferred from the magnitude of horizontal reaction force on the strip footings ( Figure 20). Since the horizontal movement of the footings is constraint in the simulations, there will be horizontal reaction force on the footing due to soil movement towards the slope. The direction of the horizontal reaction force is shown in Figure 19 ( ℎ1 and ℎ2 ). Smaller horizontal reaction force indicates that smaller soil body tends to move towards the slope. Figure 20 shows that the maximum value of ℎ2 ( = −45 • ) is about 1/3 of ℎ1 ( = 45 • ), which means that larger soil body tends to move towards the slope at = 45 • .

Effect of non-coaxial strain increment
Fabric anisotropy can cause non-coaxial strain increment in sand. 24,25,58 The non-coaxial strain increment can occur when fabric and loading direction are not coaxial, and the degree of non-coaxiality decreases gradually with shear strain. 24,25 The present model can describe the non-coaxial sand behavior in a natural way by employing a fabric-dependent plastic potential, which has been discussed in previous studies. 23,36,37 Specifically, the model parameter ℎ has an influence on the non-coaxial strain increment. Since the degree of non-coaxiality is typically small in monotonic loading with small or no change in principal stress directions, the value of ℎ is small. The model response becomes coaxial when ℎ = 0. Figure 21 shows the effect of parameter ℎ on a sand element response in drained plain strain tests. It is evident that there is small increase in the peak deviator stress and small decrease in volume expansion for both the horizontal and vertical bedding planes when ℎ is set 0. But the behaviour of strip footings is significantly affected by the non-coaxial response ( Figure 22). When the non-coaxial behaviour is neglected, the predicted becomes about 20% higher ( Figure 22). This is consistent with the results obtained by Yuan et al. 25 showing that the non-coaxial strain increment has a detrimental effect on the bearing capacity. The effect of non-coaxiality is twofold. First, the soil shows higher deviator stress when there is no non-coaxial strain (Figure 21), which increases the bearing capacity. Secondly, non-coaxial strain causes earlier F I G U R E 2 2 Effect of non-coaxial sand response on the behaviour of strip footings: (A) the relationship between distance parameter and bearing capacity factor; (B) the load-displacement relationship F I G U R E 2 3 The contour of incremental shear strain at the peak state for the strip footing with λ=0.5 and = 30 • : (A) non-coaxial strain increment considered ( ∕ = 0.09); (B) no non-coaxial strain increment ( ∕ = 0.09) development of slip surfaces. Figure 23 shows the contour of incremental shear strain at the peak state for the strip footing with = 0.5. The incremental shear strain is expressed as where is the deviatoric strain increment for each time step. It is evident that the region with high concentration of incremental shear strain has extended to the slope toe when the non-coaxial strain increment is considered ( Figure 23A). In contrast, if non-coaxial strain increment is not considered, the localization of incremental shear strain may incept from beneath the footing ( Figure 23B). More localized incremental shear strain shown in Figure 23A leads to early development of the slip surface, which has a detrimental effect on the bearing capacity. This phenomenon has also been observed in Zhou et al. 59 and Yuan et al. 25

CONCLUSIONS
A numerical study of the fabric effect on deformation and failure of strip footings on a sloping ground has been presented. The numerical simulations are compared with a group of centrifuge tests, in which sand anisotropy has been studied. A critical state constitutive model for sand accounting for fabric evolution has been used for the study. The nonlocal theory has been used to regularize the mesh-dependency of the numerical solutions. The main conclusions are: 1. The anisotropic model offers satisfactory prediction for the bearing capacity of strip footings on sand slopes with different slope angles. Single slip surfaces have been observed in most simulations. 2. Satisfactory predictions of the bearing capacity and settlement at failure cannot be achieved simultaneously by an isotropic model. When the isotropic model parameters are calibrated by test data on a horizontal bedding plane case, the model offers significant overestimations of bearing capacity. If the model is calibrated by tests with other bedding plane orientations, the prediction of the bearing capacity can be improved. But it may yield significant overprediction of settlement at failure. 3. The bearing capacity and depth of the slip surfaces decrease when the bedding plane is inclined. The minimum bearing capacity which is observed at vertical bedding plane can be 19% to 24% lower than the maximum one with horizontal bedding. The soil tends to move along the bedding plane when it is loaded due to the non-coaxial strain increment caused by fabric anisotropy. The slip surface is deeper and bearing capacity is lower when the preferred soil movement direction caused by bedding plane orientation is towards the slope. 4. The non-coaxial strain increment caused by fabric anisotropy has a detrimental effect on the bearing capacity of strip footings near slopes, which is consistent with existing research findings. 25 Non-coaxial strain increment may decrease the peak deviator stress of soil elements and causes earlier development of slip surfaces.

A C K N O W L E D G M E N T S
The authors would like to thank Dr Hansini Mallikarachchi at Ramboll UK Limited and Dr Peter Grassl at University of Glasgow for their support in the implementation of the nonlocal constitutive model.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.