Boiling heat transfer is widely used in industry and daily life and can be enhanced by porous media structure. Previously HPM have been shown to significantly improve heat transfer performance, but how the porosity gradient affects boiling heat transfer remains unknown. In this paper, the subcooled flow boiling in vertical square channels partially filled with gradient porous media was numerically investigated by using a newly developed solver named GPMPhaseChangeFoam embedded in the open-source platform OpenFOAM. The results indicate that in positive GPM inserted channels, when z ≤ 0.16 m Nuz increases with increasing Hp*, and when z ≥ 0.16 m Nuz increases with increasing Hp* from 0.2 to 0.8 but deceases with Hp* rises from 0.8 to 1.0, and the effect of Hp* on Nuz is more pronounced in the entrance region over the exit region. In negative GPM inserted channels, Nuz increases with increasing Hp* throughout the channel length in positive GPM inserted channels, and the effect of Hp* on Nuz is less pronounced in the entrance region than the exit region. Flow resistance increases with Hp*, but changes little with Δε, and Nufd increases with Hp* in both positive and negative GPM inserted channels. However, Nufd decreases with an increase in ∆ε from 0 to 0.3 in positive GPM inserted channels, from -0.3 to 0 in negative GPM inserted channels, and ∆ε = -0.3 at Hp* = 1.0 is 13.4 % higher than that of ∆ε = 0.3. Furthermore, the negative GPM inserted channel behaves higher flow boiling heat transfer performance than the HPM and positive GPM.