A Discrete Element Method(DEM) model is developed to study the particle breakage effect on the one-dimensional compression behavior of silica sands. The ‘maximum tensile stress’ breakage criterion considering multiple contacts is adopted to simulate the crushing of circular particles in the DEM. The model is compared with published experimental results. Comparison between the compression curves obtained from the numerical and experimental results shows that the proposed method is very effective in studying the compression behavior of silica sands considering particle breakage. The evolution of compression curves at different stress levels is extensively studied using contact force distribution, variation of contact number and particle size distribution curve with loading. It is found that particle breakage has great impact on compression behavior of sand, particularly after the yield stress is reached and particle breakage starts.The crushing probability of particles is found to be macroscopically affected by stress level and particle size distribution curve, and microscopically related to the evolutions of contact force and coordination number. Once the soil becomes well-graded and the average coordination number is greater than 4 in two-dimension, the crushing probability of parent particles can reduce by up to5/6. It is found that the average contact force does not always increase with loading, but increases to a peak value then decreases once the soil becomes more well-graded. It is found through the loading rate sensitivity analysis that the compression behavior of sand samples in the DEM is also affected by the loading rate. Higher yield stresses are obtained at higher loading rates. A Discrete Element Method (DEM) model is developed to study the particle breakage effect on the one-dimensional compression behavior of silica sands. The ’maximum tensile stress’ breakage criteria called multiple contacts is adapted to simulate the crushing of circular particles in the DEM Comparison between the compression curves obtained from the numerical and experimental results shows that the proposed method is very effective in studying the compression behavior of silica sands considering particle breakage. The evolution of compression curves at different stress it is found that particle breakage has great impact on compression behavior of sand, particularly after the yield stress is reached and particle size distribution curve with loading. probability of particles is found to be m acroscopically affected by stress level and particle size distribution curve, and microscopically related to the evolutions of contact force and coordination number. Once the soil becomes well-graded and the average coordination number is greater than 4 in two-dimension, the crushing probability of parent It is found that the average contact force does not always increase with loading, but increases to a peak value then reduce once the soil becomes more well-graded. It is found through the loading rate sensitivity analysis that the compression behavior of sand samples in the DEM is also affected by the loading rate. Higher yield stresses are obtained at loading rates.