basis of similarity theory, which was previously used in centrifuge tests [30; 31; 32]. However, the grading of realistic sand usually spans many orders of magnitude, which means the voids created by the large particles need to be filled by massive volumes of fine particles [33; 34; 35]. Then, in the DEM models, a relatively narrow particle distribution has to be adopted, especially for three-dimensional simulations. Another limitation in the existing DEM simulations is that the realistic particle shapes were over-simplified as analytical shapes, spheres [36], clumps [37] and polyhedra [38]. These simplifications seem to be acceptable for quartz sand, but apparently not for calcareous sand [39; 40]. The more complex idealisation of particle shapes requires more computation efforts due to the complicated contact algorithm. Accurate representation of the particle size and shape is critical to the modelling results, and hence, needs to be calibrated prior to the simulations [41]. Unfortunately, a majority of the existing DEM simulations ignored this procedure and were based on assumptions of input parameters without justification.  
To promote the efficiency of the DEM and enlarge the modelling scales with higher realism of the particle shapes, parallel computing on the central processing unit (CPU) or graphics processing unit (GPU) has been explored by a limited number of researchers. References [42; 43; 44] developed single-CPU parallel algorithms using a loop-based parallel library OpenMP with a shared memory. [PERSON] and [PERSON] [45] presented a multiple-CPU parallel framework with a message passing interface on a distributed memory platform. In comparison with most commercially available CPUs, GPUs feature more dedicated cores in the processor, lower thread-scheduling cost and higher memory bandwidth; as such, they are ideal for arithmetically intensive and highly parallel computations. Currently, most GPU parallelisation is implemented using the 'compute unified device architecture' (CUDA) [46], which provides a friendly platform for coding in languages, such as C++ and Fortran [47; 48; 49; 50]. [PERSON] et al. [51] performed an industrial-scale analysis of mill charges with more than eight million polyhedron particles using a GPU-parallelised DEM code. Due to the differences in architectures between the CPU and GPU, DEM codes developed for the conventional CPU calculations need to be specially tuned to be efficiently parallelised on GPUs.  
The main bottleneck in parallelising the DEM is the potential data race between the computing cores writing to the common memory address concurrently for the creation of a contact pair list of particles. [PERSON] et al. [17] utilised atomic operations to update the particle list in the material point method by taking advantage of modern graphics hardware features, whose complication excludes most developers from their implementation. [PERSON] et al. [51] established a cell-based hash list for the neighbouring particles, in which extra overhead is taken on sorting the particles by their distances to the origin. [PERSON] et al. [52] placed each individual particle into a specific cell with very fine mesh, which often means a heavy memory requirement for a large-scale model along with its simplicity. The technique to update the particle list is vital to the acceleration effect of the parallel framework. In this paper, a new hybrid technique to establish the contact pair list is presented, which combines the advantages of the methods developed by References [51; 52]. A GPU-parallel framework of the DEM is presented with all the computations taken on the GPU cores. Three benchmark cases, the triaxial test of a sand specimen, the cone penetration test and sand flow due to dam break, are used to evaluate the performance of the GPU parallel strategy. Acceleration of the GPU parallel simulations over the CPU sequential counterparts is quantified in terms of speedup.