We study the dynamics of hydrothermal plumes with the 3D time-dependent, Eulerian, adaptive mesh refinement code GERRIS, which solves the equations of viscous, incompressible hydrodynamics. We have implemented a new module into Gerris that treats buoyancy-driven turbulence by means of a subgrid mode. Our model is validated in numerical experiment and applied to the dynamics of a rising plume. First we simulate hydrothermal plumes in a static environment and compare our results to the widely used integral models (MTT or Briggs' model). The entrainment coefficient that we deduce from simulations falls into the range of the experimentally determined values. We also investigate the ratio between the level of the neutral-buoyancy layer and the maximum plume height. This ratio is frequently used to estimate plume heat flux via the measured level of neutral buoyancy. Although the ratio is only moderately (less than 10%) higher than the one predicted by the integral model, heat flux estimations can be substantially different. Finally, we explore the importance of background currents. We find that the simulated trajectories agree with integral models in the rising stage but the subsequent oscillations around the neutral-buoyancy layer are damped much more quickly and the level of the neutral buoyancy is also higher, same as the calm environment cases. By simulating the oscillation of a plume with suppressed transported turbulence and find a stronger oscillation than the original simulation, we suggest that a significant fraction of the difference between our model and the integral model can be explained by the absence of the turbulent transport of the latter.