In this study, three-dimensional simulations of cone penetration tests (CPT) in cemented sandstone have been investigated using the Discrete Element Method (DEM). We generated the realistic grain size distribution of sand- stone in the numerical simulations using the Ustyurt-Buzachi Sedimentary Basin samples. Numbers of numerical CPT tests at di erent bond strength values were performed with the penetrometer vertically moved down at a constant rate. The results of the real-world particle size distribution show that the cone penetration resistance and side friction increase with increas- ing depths; and with increase in bond strength the cone resistance and side friction decrease, while the friction ratio increases. The result of numerical CPT tests in cemented sandstone was found to be in good agreement with the Soil Behaviour Type (SBT) classi cation system from CPT data. Keywords: Discrete Element Method; cone penetration; cemented sandstone 1. Introduction The main applications of the cone penetration test (CPT) in geotechni- cal engineering are for classi cation and identi cation of soil types, and for stratigraphic pro ling. In a eld cone penetration test, a cone with its sleeve Corresponding authors Email address: aigerim.rakhimzhanova@nu.edu.kz (Aigerim Rakhimzhanova) is moved down into the ground, where the cone resistance is the force re- quired pushing the tip of the cone and the friction sleeve is the force required to push the sleeve through the ground. Then the friction ratio is calculated by the ratio between sleeve friction and cone penetration resistance. The classi cation of soil type can be identi ed from such basic CPT data. Here, we provide a brief summary of the literature relating to the cone penetration test in granular soils that has been numerically investigated us- ing the Discrete Element Method (DEM) proposed by Cundall and Strack (1979). Huang and Ma (1994) and Jiang et al. (2006) numerically (2D) stud- ied the mechanism of deep penetration in granular soils; Jiang et al. (2014) numerically (2D) examined the mechanism of inclined CPT tests in granular ground; Jiang et al. (2015) numerically (2D) studied the CPT mechanism un- der di erent gravity conditions; Arroyo et al. (2011), Butlanska et al. (2010, 2013) and Ciantia et al. (2016) numerically (3D) modelled the CPT test in a calibration chamber taking into account the particle crushing behaviour. The paper presents results obtained from 3D DEM simulations of cone penetration tests on medium dense samples of cemented sandstone. This research extends recent development of a simple 3D bond contact model for cemented sandstone material in Rakhimzhanova et al. (2019, in press). This study was successfully applied to the triaxial shear tests of cemented sandstone, where the macro responses of the numerical samples were a good match for the experimental results of cemented sandstone. This study was successfully applied to the triaxial shear tests of cemented sandstone, where the macro responses of the numerical samples were found to be in good agreement with the experimental results of cemented sandstone. This paper is organized as follows. The numerical contact model will be described in Section 2. In Section 3, the numerical simulation of cone penetration test details and procedure are presented. The numerical results are described in Section 4 and conclusions are reported in Section 5. 2. Numerical methods 2.1. The Discrete Element Method In the theoretical model integrated in the DEM model, all motions of particles are guided by Newton's equations of motion and the components 2 (i= 1;2;3) inx,y, andz- directions of the translational and rotational accelerations of each particle are given by Eq.(1) and Eq.(2), respectively: midvi dt=Fci+migi (1) Idwi dt=FtiR (2) where,mi;viand!iare the mass, linear and angular velocities of particle i, respectively; gis the gravity, Fciis the solid interaction force between particles,FtiRare the moments due to the tangential components of the contact forces, and Iis the moment of inertia. 2.2. Contact models 2.2.1. No adhesion In all simulations, particles are considered to be elastic and frictional, where the particle-particle and wall-particle contact forces without adhesion are calculated using the theories of Hertz to determine the normal force, and Mindlin (1949) to calculate the tangential force. According to the Hertz theory, the normal force Fndue to relative approach of the spheres is Fn=4 3Ep R 3 (3) where, 1=Rthe relative curvature of the surface of contacting spheres of radiiRi(i= 1;2) and is expressed as 1 R=1 R1+1 R2(4) andEi(i= 1;2) are the relative Young's moduli and it is given below 1 E=(12 1) E1+(12 2) E2(5) whereiare the Poisson ratios of two contacting spheres. From Eq.(3) the normal contact sti ness is de ned as kn= 2Ep R (6) 3 The tangential force depends on the increasing or decreasing of the normal force and it is calculated incrementally. For the increasing normal force Fn0, the tangential force at the time step ithis Fi t=Fi1 t+ki t (7) and for the decreasing normal force Fn<0, it is de ned as Fi t=Fi1 tki t ki1 t +ki t (8) whereis the rigid displacement of the sphere centers. For FnFn, the tangential force is given as Fn=Fn (9) whereis the interparticle friction coecient. The tangential contact sti - ness is calculated by kt= 8Gp R (10) and it is de ned by the Mindlin (1949) `no-slip' solution, where Gis the complex shear modulus de ned by 1 G=(21) G1+(22) G2(11) 2.2.2. With adhesion A simple 3D bond contact model for cemented sandstone material de- veloped by Rakhimzhanova et al. (2019, in press) was used to calculate the interparticle normal contact force with adhesion, which simply modi ed the existing JKR model for auto-adhesive silt sized particles provided by John- son et al. (1971), and no-slip model proposed by Mindlin (1949) to calculate the tangential contact force. Johnson (1977) provided the relationship be- tween normal interaction force at contact, Fn, and and the relative approach between two particles, (see Fig. 1). Initially, when two particles come into contact ( = 0) because of van der Waals attractive forces, the normal contact force, Fn, instantly drops down to point A and at this point the normal contact force equals Fn=8Fnc=9. The normal contact force moves from Point A from point B in case of loading (compression) and during the unloading (decompression) it returns to point A due to the elastic response and the relative approach is equal to 0 (there 4 Figure 1: JKR theory is still a nite area of contact) and all the work done during the compression will be recovered. The tensile force increases during the extra work from point A to point C and it decreases until Point D when the contact between two particles breaks, where the normal contact force Fn=5Fnc=9 and the negative relative approach = f, where f=3F2 nc 16RE21=3 (12) Nevertheless, the cement bonds of oil led sandstone are not as stretchable and break in a brittle manner. Therefore, in the simple 3D contact model for cemented sandstone material developed by Rakhimzhanova et al. (2019, in press) the bond breaks at point C, where the normal contact force Fn=Fnc and the relative approach = f=32=3. The maximum tensile force required to break the contact (at point C) is Fnc= 1:5R(13) where  is the work of adhesion and  = 1+ 2= 2 are the surface energies of the two solids. In the simple 3D contact bond model, any subsequent new contact force is calculated by the Hertzian theory. 5 3. Numerical simulations In this study, the cone penetration tests have four stages: particle gener- ation, pluvial deposition, compression, cone penetration (see Fig. 3). Three- dimensional DEM simulations of cone penetration tests have been performed using 10000 frictional elastic spheres bonded together over a range of bond strengths (interface energy of adhesion) of  = 0 ;5;10;20 J/m2 (where  = 0 J/m2 is an uncemented sample) and compressed at 1 MPa of overburden pressure, where in all cases particle rotations were prevented. 3.1. Particle generation and deposition In this study, the same eight di erent sand grain sizes and material prop- erties used by Rakhimzhanova et al. (2019, in press) were utilized for the CPT test simulations, where the grain sizes were initially selected from the Qicpic dynamics image analyzer measurements of Shabdirova et al. (2016) for the Ustyurt-Buzachi Sedimentary Basin samples (Fig. 2). Particles were randomly generated (Fig. 3a) within a thin rectangular parallelepiped of di- mension 8.4 mm in width, 21 mm in height and 1.4 mm in thickness: 0.355 mm (306), 0.3 mm (610), 0.275 mm (1119), 0.25 mm (1973), 0.22 mm (1558), 0.2 mm (1505), 0.18 mm (1750), and 0.15 mm (1179). On the basis of the Uni ed Soil Classi cation System (USCS) these grain sizes pass through the sieve sizes # 45-100, which indicates that the soil type is S { sand with the subdivision of ne sand. The thin parallelepiped workspace was divided into 720 small DEM boxes (box dimension of 0.0007 mm): 12 boxes in width, 30 in height and 2 in thickness. The following material properties were used for all sand particles: exper- imental particle density = 2605 kg/m3, Poisson ratio = 0:3 and Young's modulusE= 70 GPa. Further, a vertical gravity eld was applied to the system and all particles were pluvially deposited to form the initial homoge- neous bed (Fig. 3b). One nite wall was created at the base of workspace that one may avoid/prevent the periodic boundary condition at this side. The following material properties of iron were used for all walls: density = 7900 kg/m3, Poisson ratio = 0:29 and Young's modulus E= 210 GPa. 3.2. Compression and cone penetration Three top nite walls were created at the top of the deposited specimen and pushed down at a wall velocity of 0 :1 m/s in order to vertically compress the sample (Fig. 3c). When a vertical stress level (overburden pressure) 6 Figure 2: Grain Size Distribution reached 500 kPa, then the bond strength values  = 0 ;5;10;20 J/m2 were set. When the overburden pressure reached 1 MPa, the top three wall veloc- ities were reset to 0. After the compression of DEM cemented sandstone, a standard cone penetrometer was rst created above the bed and then pushed down at a high constant rate of 0 :1 m/s in order to avoid computational cost. The penetrometer had a \radius" Rof 0:7 mm and an apex angle of 60(Fig. 3d). It was de ned as 2 rigid frictional walls inclined to X-axis at 60to simulate the cone/tip penetrometer and two rigid vertical frictional walls to simulate the penetrometer sleeve (Fig. 4. During the simulations the cone was always in contact with about 33 35 particles, which give the resistance value. The surface of initial ground is treated as the origin and the distance between the tip point (lowest point of tip) and the origin is de ned as the penetration depth y. The cone penetrometer was removed from the penetrated specimen once the cone tip had reached 2/3 of the depth. 4. Numerical results Figs. 5 a, 5 b and 5 c show the evolution of cone resistance, friction sleeve and friction ratio data at di erent penetration depths for  = 0 ;5;10;20 7 (a) Particle generation (b) Pluvial deposition (c) Compression (d) Cone perforation Figure 3: 3D DEM cone penetration simulation stages 8 Figure 4: Cone penetrometer (modi ed from Jiang et al. (2006) J/m2, respectively. The cone resistance is de ned as the vertical pressure measured on the cone qc=Pfy 2Rt(14) wherefyis the contact force vertically acting on inclined wall (cone); Ris the cone radius; tis the thickness of the bed ( t= 2R), and the frictional sleeve,fc, is de ned as fs=Pfy 2Ht(15) wherefyis the vertically acting contact force on the frictional sleeve; His the frictional sleeve height in the ground. It can be seen that the cone penetration resistance and the friction sleeve increase with increasing depths. Similar behavior was obtained in 2D DEM CPT simulations reported by Jiang et al. (2006), and in the real eld CPT data reported by Robertson (1990). The friction ratio, Rf, is calculates as Rf=fc qc100% (16) 9 Figure 5: The numerical CPT data 10 whereqcis from Eq.(14) and fsfrom Eq.(15). The results show that the cone resistance and side friction decrease, and the friction ratio increases with increase in bond strength. Robertson et al. (1986) and Robertson (1990) proposed CPT-based charts as a guide to predict the Soil Behavior Type (SBT) based on the basic eld CPT data where only the total cone penetration resistance, qt, and the sleeve friction,fc, are available (Fig. 6). Campanella and Robertson (1982) pro- posed that the cone resistance, qc, can be corrected to the total cone resis- tance,qt, as qt=qc+u(1a) (17) whereuis the pore pressure measured between the tip of cone and the sleeve friction and ais the net area ratio calculated as a=d2=D2, where D is the diameter of the cone and dis the cone diameter of the friction sleeve. This correction should be made to account for unequal end area e ects because the water pressure acts on an area immediately behind the cone tip. Since the numerical simulations of cone penetration tests conducted on dry particles the net area, a, was equal to 1 and the assumption of qc=qthas been made. The data points of 3D DEM simulations of CPT tests of cemented sand- stone have been added to this classical soil classi cation chart. The average values of cone resistance and friction sleeve was taken starting from the depth 0:0015 m, once they became more stable. It is shown that the CPT data of samples with higher bond strengths of  = 10 and 20 J/m2 locates in SBT zone 12, which is identi ed as overconsolidated or cemented sand to clayey sand; and  = 0 and 5 J/m2 locate in SBT zone 11, which behaves as very sti , overconsolidated ne grained sand. The SBT zone 11 and 12 for all samples can be explained as at the end of compression stage, the material that contains only ne sand particles (Fig. 2) had a porosity of about 0.38, which are similar to those porosities obtained for 3D DEM simulations of triaxial shear tests of cemented sandstone on medium dense samples compressed at 1MPa and  = 0 ;5;10;20 J/m (Rakhimzhanova et al. (2019, in press)). For  = 5 J/m2 it can be explained that, at the lowest value of  = 5 J/m2, the bond strength is soft and more bonds were broken, compared with higher values of  = 10 and 20 J/m2, with a tendency to behave like an uncemented sand. The SBT zone 11 for uncemented sand ( = 0 J/m2) can be explained in that the sand becomes sti er as initial density and/or con ning pressure increases (frictional response). Similar behavior has been experimentally obtained for compressive behavior of ne 11 Figure 6: SBT classi cation system from CPT data, as proposed by Robertson et al. (1986) sand by Kabir et al. (2010) and numerically studied by Rakhimzhanova et al. (2019, in press). 5. Conclusion In this paper, the 3D DEM simulations of cone penetration test of ce- mented sandstone have been investigated at various values of bond strength. Material properties of reservoir equivalent sandstone the Ustyurt-Buzachi Sedimentary Basin samples were reproduced for the numerical experiments. We have extended the previous research on 3D bond contact model of ce- mented sandstone material to the numerical simulations of CPT tests and the identi ed bond strength values to be used for veri ed with experimental results by Rakhimzhanova et al. (2019) were used for the simulations of CPT test. The results show that with increase in bond strength the cone resis- tance and side friction decrease, while the friction ratio increase; the cone resistance and side friction increase with increase in penetration depths. It has been demonstrated that the numerical CPT tests result of cemented sandstone are in good agreement with the Soil Behavior Type (SBT) classi - cation system from CPT data. The ndings suggest that this method could also be helpful for practicing engineers to control perforation technique in the sandstone eld. 12 Acknowledgements This research was supported by the Nazarbayev University Faculty De- velopment Competitive Research Grant No 110119FD4502 and by the EU Project GEO-RAMP (H2020-MSCA-RISE-2014 GA 645665). The authors gratefully acknowledge Prof. Stefano Utili (Newcastle University) and Dr. Helen Cheng (University College London) for their help during the rst au- thor's visits to the UK.