Biofilms, which are nothing else than a cluster of bacteria, can be modeled using different approaches. They might be assumed to behave like an in-viscid fluid or a viscoelastic solid. Irrespective of the constitutive behavior, the spatiotemporal evolution of the biofilm structure is a very complex process due to its multiphysics nature. In fact, it involves mechanical deformation coupled with nutrient transport and biological growth. The problem becomes more sophisticated when it comes to multi-species colonies. The presence of biological growth poses serious challenges to the mathematical modeling of the problem in terms of thermodynamic consistency. This work presents a fully continuum-based and variational approach derived from thermodynamical principles, i.e. Hamilton's principle for dissipative processes. It ensures the compliance of the mathematical model with the thermodynamic laws. The growth of multi-species and multi-colonies biofilms is captured using the proposed method. The numerical implementation is carried out using the finite element method (FEM). Various test cases are provided to examine the impact of different conditions on the growth pattern. In particular, the availability and the gradient of nutrients, the proximity of different species, the relative motility of the species, and the geometrical constraints (mechanical obstacles) are thoroughly investigated to identify the most influential factors contributing to biofilm development. By that, the versatility and applicability of the model and the numerical tool are well established.