A previously developed semi-implicit method to solve a density dependent diffusion-reaction biofilm growth model on uniform Cartesian grids is extended to accommodate non-orthogonal grids in order to allow simulation on more complicated domains. The model shows two non-linear diffusion effects: it degenerates where the dependent solution vanishes, and a super-diffusion singularity where it approaches its upper bound. The governing equation is transformed to a general non-orthogonal ξ−η curvilinear coordinate system and then discretized spatially using a cell centered finite volume method. The nonlinear biomass fluxes at the faces of the control volume cell are split into orthogonal and non-orthogonal components. The orthogonal component is handled in a conventional manner, while the non-orthogonal component is treated as a part of the source term. Extensive tests showed that this treatment of the non-orthogonal flux component on the control volume face works well if the maximum deviation from orthogonality in the region of the grid where the biomass is growing is within 15-20 degrees. This range of validity is smaller than the one obtained with the same method for the simpler porous medium equation which is the standard test problem for degenerate diffusion equation but does not have all of the features of the biofilm model.