When we construct continuous and/or discrete mathematical models in order to describe a real-life problem, these models should possess various qualitative properties, which typically arise from some basic principles of the modelled phenomenon. In this paper we investigate this question for the numerical solution of initial-boundary value problems for parabolic equations with nonzero convection and reaction terms with function coefficients in higher dimensions. The Dirichlet boundary condition will be imposed, and we will solve the problem by using linear finite elements and the theta-method. The principally important qualitative properties for this problem are the non-negativity preservation and different maximum principles. We give the conditions for the geometry of the mesh and for the choice of the discretization parameters, i.e., for theta and the time-step sizes, under which these discrete qualitative properties hold. Finally, we give numerical examples to investigate how sharp our conditions are.