In this work we consider a mathematical and computational model for biofilm growth and nutrient utilization. In particular, we are interested in a model appropriate at a scale of interface. The model is a system of two coupled nonlinear diffusion--reaction partial differential equations (PDEs). One of these PDEs is subject to a constraint, which can be characterized as a parabolic variational inequality (PVI). Solutions to PVI have low regularity which limits the numerical scheme to low order. We analyze the numerical approximations of this model by implementing two numerical methods: (i) low order Galerkin finite element method, and (ii) mixed finite element method with the lowest order Raviart-Thomas elements. We show the well-posedness of the approximate problems, derive rigorous error estimates, and present numerical experiments in 1D, 2D, and 3D that confirm the predicted estimates. We also illustrate the behavior of biofilm and nutrient dynamics in simple and in complex porescale geometries.