An accurate estimation of a buoyant plume rise is one of the basic aspects for a correct estimation of the ground level concentration emitted by an industrial stack. The introduction of the plume rise computation in a Lagrangian particle model is not straightforward. In fact the buoyant forces acting on the plume portions depend on the difference between their temperature and that of the background air that can be computed only considering the entrainment phenomenon. To do that, these models should take into account all the fluid simultaneously, i.e. by filling the whole domain with a large amount of particles. Several authors have suggested to overcome this problem by introducing a vertical velocity whose values come from analytical solution of Eulerian equations (Briggs,1975) or solving a set of differential equations describing the time and space evolution of bulk plume quantities (Hurley, 2005; Webster and Thomson, 2002; Anfossi et al, 2009). This approach is accurate and very efficient from the computational time point of view but some problems may arise in those situations where the analytical equation for the vertical velocity is difficult to define, i.e. where several plumes, released by different stacks located close to each other, mix themselves together. An alternative method for the buoyant plume rise computation is therefore proposed. It is based on the same idea described in a recent paper related with chemical reactions in Lagrangian particle models (Alessandrini and Ferrero, 2009) for simulating the ozone background concentration. A fictitious scalar transported by the particles, the temperature difference between the plume portions and the environment air temperature, is introduced. As a consequence no more particles than those inside the plume have to be released to simulate the entrainment of the background air temperature. In this way the entrainment is properly simulated and the plume rise is calculated from the local property of the flow.