Compared to the continuum-based methods, the discrete element method (DEM) may reflect the effects of granular characteristics (such as the mean particle size d50, coefficient of uniformity Cu and coefficient of curvature Cc) on the behavior of soil mass. Hence, in this work, the rolling resistance linear (RRL) contact model in combination with the circular particles in the DEM, as a simplified approach of reflecting actual shapes of soil particles, is utilized to assess the effects of granular characteristics on the ultimate bearing capacity pu of the ground. Furthermore, a large amount of DEM simulations of ground with various granular characteristics are conducted, so that a dataset of relating the pu to the granular characteristics for machine learning (ML) can be achieved. Based on the dataset of DEM simulations, three ML algorithms, namely multiple linear regression (MLR), artificial neural networks (ANN) and extreme gradient boosting (XGBoost), are applied to training the prediction models of pu. The research results reveal that the footing rotations and the asymmetric failure pattern of ground in the DEM simulations are largely attributed to the granular characteristics of soil mass. For a ground composed of soil particles with larger d50 and larger Cu, the pu of the ground generally tends to be higher and the thickness of developed shear band tends to be greater. For a ground composed of soil particles with larger d50 but smaller Cu, the asymmetric failure pattern of ground is more legible and the predicted pu tends to be smaller. It was also found that among the investigated ML algorithms, the MLR algorithm may provide an explicit model for the pu, while ANN and XGBoost algorithms offer higher prediction accuracy.