## 2 DEM Modelling  
The DEM simulation in this paper is carried out using the well-recognized Particle Flow Code(PFC\({}^{3\text{D}}\)) [32]. Approximately 27,000 particles in a cubic space of 110 mm\({}^{3}\) constitute the specimen (Figure 1), which is confined by three pairs of frictionless walls. The particle size distribution of the specimen mimics Toyoura sand (Figure 2) with a mean particle diameter of 0.214 mm and the coefficient of uniformity of 1.425. The nonlinear Hertz-Mindlin contact law [33] is used and the shear modulus and Poisson's ratio is 5.8 GPa and 0.15, respectively. The friction between the walls and particles is set to be zero. The parameters of the simulation in this paper are listed in Table 1.  
Figure 1: Cubic specimen with particles in the DEM simulation (walls are omitted). Figure 2: Particle size distribution in DEM simulation. \begin{table}
\begin{tabular}{c c} \hline \hline
**Specimen Dimensions (mm)** & **4.8 \(\times\) 4.8 \(\times\) 4.8** \\ \hline Particle density (kg/m\({}^{3}\)) & 2650 \\ No. of particles & 27,000 \\ Mean particle diameter (mm) & 0.214 \\ Inter-Particle friction coefficient & 0.5 \\ Wall-particle friction coefficient & 0 \\ Wall stiffness (N/m) & \(1\times 10^{10}\) \\ Contact law & Hertz-Mindlin \\ Shear modulus of particles (GPa) & 5.8 \\ Poison’s ratio of particle & 0.15 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Parameters in the DEM simulation. The clump of aspect ratio 1.7 (Figure 3) is adopted to produce inherent anisotropy. Initially, the anisotropic specimen is prepared with reference to the method adopted by [PERSON] and [PERSON] [34]. The particles are generated with a coincident orientation and undergo two-stage isotropic consolidation with 10 kPa consolidation pressure. The orientations of the particles are fixed in the first stage and are then set to be free in the second stage. As illustrated in Figure 4, the angle \(\delta\) between the long axis orientation of clumps and horizontal plane represents the inherent anisotropy principal direction. Five specimens of different anisotropy principal direction (\(\delta\) = 90\({}^{\circ}\), 60\({}^{\circ}\), 45\({}^{\circ}\), 30\({}^{\circ}\), and 0\({}^{\circ}\)) but the same anisotropy degree is prepared. The contact normal distributions of these five specimens are plotted by 3D rose graphs in Figure 5. For specimens with angle \(\delta\) equaling 90\({}^{\circ}\), 60\({}^{\circ}\), 45\({}^{\circ}\), 30\({}^{\circ}\), and 0\({}^{\circ}\), the scalar anisotropy parameter \(a_{r}\) are 1.291, 1.336, 1.354, 1.361, and 1.311, respectively, representing the nearly identical anisotropy degree. From the rose graph, specimens with different initial anisotropic principal directions can be produced effectively by this specimen preparation method. As a reference object, an isotropic specimen is made, with random orientation in the generation stage and no restraint of orientation at any stage. In the isotropic compression stage, the inter-particle friction angle is 0.15 while in the \(K_{0}\) loading stage the angle is 0.5. After consolidation, the void ratios are 0.623, 0.622, 0.626, 0.626, 0.627, and 0.624 respectively, for the isotropic specimen and the other five ones of \(\delta\) = 90\({}^{\circ}\), 60\({}^{\circ}\), 45\({}^{\circ}\), 30\({}^{\circ}\), and 0\({}^{\circ}\). The effect of initial void ratio has been eliminated as much as possible since their values are roughly equal. Figure 4: Generating initial anisotropic specimens. Figure 3: Three kind of particle shapes in the simulation. Next, for evaluating \(K_{0}\)-value, one-dimensional compression is performed in the DEM simulation with the stress path, similar to that in laboratory triaxial tests. Due to the difficulty of measuring lateral stress in oedometer tests, triaxial tests are usually performed and an initial isotropic stress is applied to stand the sample before the \(K_{0}\) tests. Therefore, all the samples are initially isotropic consolidated, after which \(K_{0}\) loading is conducted.