In this paper, we develop a numerical model for simulating the solid-liquid-gas three-phase flow in unconsolidated particle layers. Based on the discrete element method (DEM) and the multiphase fluid model in the framework of the lattice Boltzmann method (LBM), a multiphase fluid-solid two-way coupling algorithm is proposed. In this model, the fluid-fluid interface is tracked using a phase-field method, and the multiphase fluid-particle interaction is tackled by a combination of the momentum exchange method for the flow field and the immersed boundary method for the phase field. We applied the method to simulate the upward migration of the leaked gas bubbles through a brine-filled sediment column at the seafloor, and investigated the influences of the leak flow rate and the interfacial tension on the bubble rising behavior. The results indicate three different flow regimes: connected finger flow, transition flow, and dispersed bubbly flow. These flow regimes can be distinguished by the dimensionless Eötvös and Weber numbers. The proposed numerical method can accurately characterize various multiphase interaction mechanisms at the mesoscopic scale and has powerful advantages in simulating complex fluid-particle coupling problems.