Abstract: To better understand the catastrophic mechanism of mud inrush disasters of tunnels, a numerical computing platform based on the coupled lattice Boltzmann method-discrete element method (LBM-DEM) is developed and used to simulate the evolution process of mud inrush of tunnels. According to the simulated results, the effects of particle bonding strengths of disaster-causing media, groundwater pressures and sizes of mud inrush holes on the characteristics of mud inrush of tunnels are analyzed. The results show that coupled LBM-DEM simulation can well reproduce the evolution process of four successive stages of mud inrush of tunnels: "starting, accelerating, decelerating and stabilizing". The failure form of unbonded disaster-causing media after mud inrush is approximately straight, whereas the failure form of disaster-causing media with certain bond strength is generally arc or parabolic. The expanded range of failure zone and mud inrush mass both decrease with the increase of the inter-particle bond strength of disaster-causing media. The higher the water pressure is, the faster the mud inrush mass increases after the occurrence of mud inrush disasters, and the more the final mud inrush mass is, which is more remarkable when the inter-particle bond is much stronger. When there is no bond between particles of the disaster-causing media, the models with different sizes of mud inrush holes have basically the same failure zone and mud inrush mass after stabilization. However, when a certain strength of bond is formed between particles, the mud inrush mass increases faster and the final mud inrush mass is more with the increase of the sizes of mud inrush holes. The mud inrush disasters are the result of the combined action of disaster-causing geo-materials, groundwater pressure and tunnel excavation.