Non-explosive volcanic eruptions can result in the formation of lava domes when highly viscous magmas extrude on the surface. To distinguish between different mechanisms during the evolution of lava domes, isothermal and non-isothermal numerical models were used to explore the processes behind the evolution of the internal structure and morphological shape. To model the lava dome growth, we assume that the lava behaves as a viscous fluid, and use the Navier-Stokes, continuity, and the heat equations, along with relevant initial and boundary conditions. Additionally, an advection equation is employed to describe the movement of the interface between air and lava. The lava viscosity depends on the melt viscosity and the volume fraction of crystals. The latter is determined from an evolutionary equation that describes simplified crystal growth kinetics of degassing-induced crystallization. Numerical modelling employs finite volume and volume of fluid methods. The models demonstrate how lava dome dynamics are influenced by crystallization, the radiative-convective heat transfer, and the rheological stiffening driven by cooling. If the lava viscosity mainly depends on the volume fraction of crystals (as in the isothermal case study), the lava dome grows upward initially and later advances horizontally. The non-isothermal model shows that a highly viscous layer (carapace) forms within the uppermost part of the lava dome due to cooling, which restricts rapid dome broadening by creating steep slopes on the sides during long periods of dome building. Both models provide insights into the underlying processes and mechanisms of the evolution of lava domes, including their growth, carapace formation, and morphological modifications.