Biofilms and their formation dynamics are ubiquitous and complex in porous media. The mechanism of biofilm formation on solute transport behavior remains limited, which inhibits potential biofilm applications such as bioremediation. In this study, we present a new numerical solver, BioFOAM, based on the micro-continuum theory, to simulate the coupled pore-scale processes of biofilm formation, fluid flow and solute transport in heterogeneous porous media. The BioFOAM explicitly solves the Darcy-Brinkman-Stokes equation, the convection-diffusion equation, and Monod kinetics in an iterative way. Benchmark tests are conducted to validate and quantify regimes of biofilm formation. We find that the competition among diffusion, advection, and the growth kinetics controls biofilm formation patterns. This competition partially explains the emergence of anomalous transport features in the growth-clogging regime when the growth kinetics dominate over diffusion and advection. When the growth kinetics, diffusion, and advection are comparable, the growth and decay processes of biofilm reach equilibrium. When advection dominates other processes, biofilm formation could not occur. Finally, we apply our model to simulate biofilm formation in real quartz sand media. We observe strong velocity intermittency in the growth-clogging regime in quartz sand media. The velocity probability density function p(u(x)) for low velocities follows a power law (p(u(x))proportional to u(x)(alpha), with |alpha| increasing from |alpha| <0.05 to |alpha| > 1), which corresponds to the intermittency that enhances solute spreading in the breakthrough curves with typical anomalous features. These results indicate that the BioFOAM model is able to quantify biofilm formation patterns and simulate the growing interest in the effects of biofilm on solute transport behavior at the pore scale.