The coal seam gas industry has raised public concerns about the potential risk of groundwater contamination, where gases leaked from coal seams are thought to pollute groundwater. However, the basic principles and controlling parameters for gas seepage from deep ground formations to the ground surface have not been fully understood. As a possible mechanism for gas transport in the subsurface environment, discrete bubble flow was previously investigated using laboratory experiments by Ma et al. (Water Resour. Res, 2015, 51 (6), 4359–4373). This study developed a multiphase computational fluid dynamic (CFD) model to simulate discrete bubbly flow in a two-dimensional granular porous media at the pore scale. Following the experimental setup from Ma et al. (Water Resour. Res, 2015, 51 (6), 4359–4373), a “point source” with preset bubble fluxes was specified in a simulating domain representing the flume size in the earlier experiments. There were around 7,000 granular particles within this domain to model the porous media. This numerical model was validated by conserving the gas mass in the simulating domain. The simulation results provide more physical insights into complex bubble transport behaviour in porous media through specific plume parameters. The breakthrough time of the bubble plume and the cross-sectional averaged velocity of ambient pore water flow were manifested to be proportional to the gas release rates in the logarithmic scales. Also, the bubble plume width was also observed to be proportional to the gas release rates. Moreover, the gas distribution on the top boundary could be observed. The outcomes were further tested against the scaling solutions proposed by Ma et al. (Water Resour. Res, 2015, 51 (6), 4359–4373) with disagreements. The limitations of this multiphase computational fluid dynamic model were finally discussed.