A newly developed two-phase mixture model is applied, in conjunction with a control-volume-based finite difference method, to numerically investigate boiling with thermal convection in a porous layer heated from below. The numerical procedure employs a fixed grid and avoids tracking explicitly the moving interface between the liquid and two-phase regions. Numerical results are obtained to shed light on the intricate interactions between boiling and natural convection as well as to explain experimental observations. Four distinct flow patterns that were observed in previous experiments are predicted. A quantitative comparison of the predicted and measured vapor volume fraction in the porous bed shows good agreement. The numerical results also agree with published linear stability results. In addition, the present study documents the effects of important parameters such as Rayleigh number, bottom heat flux, and aspect ratio.