# Lagrangian model with heat-carrying particles

### Bianca Tenti

### Stefano Alessandrini

The dispersion of pollutants in the atmosphere, whether from industrial emissions, wildfires, or other sources, poses significant challenges to air quality management and environmental protection. Understanding the behavior of plumes is crucial for predicting their dispersion patterns and potential impacts on human health and the environment. In this work, we present a new plume rise scheme based on heat transport. The idea at the basis of the new algorithm is the same as the actual scheme embedded in the Lagrangian Particle Model SPRAY-WEB. The temperature difference between the ambient and the plume and the vertical velocity of the plume are expressed on a fixed Eulerian grid. The particles are assigned with an equivalent momentum, temperature, mass, and density, transported as scalar quantities with the particles following the Lagrangian description of the motion. This allows us to reproduce the entrainment phenomenon as a mixing of two fluids (environmental air and plume) with different temperatures: the resulting temperature is given by Richmann's law.

The results obtained with the old plume rise algorithm and the new one are compared with Briggs (1975) analytical curve in the case of an idealized atmosphere with a neutral stratification and a constant horizontal wind and with experimental data. From the comparison in an ideal atmosphere, it emerged that with the new algorithm, the plume reaches a higher height than with the old one and the asymptotic trend obtained with both models follows the Briggs curve. As for the comparison with the measurements, the results obtained with the new algorithm are in better agreement with the experimental data than the old one.

- Article
(711 KB) - Full-text XML
- BibTeX
- EndNote

The computation of plume rise is one of the basic aspects for a correct estimation of the transport and dispersion of airborne pollutants and for the evaluation of ground level concentration. A forced buoyant plume rises under the action of its initial momentum and buoyancy. It experiences a shear force at its perimeter, where momentum is transferred from the plume to the surrounding air, and ambient air is entrained into the plume. This phenomenon, entrainment, is responsible for the plume diameter increase, of the decrease of its mean velocity and of the average temperature difference between air and plume as well. In the first stage the plume also spreads under the action of the buoyancy-generated turbulence but progressively the effect of ambient turbulence becomes predominant. In a calm atmosphere, plumes rise almost vertically, whereas in windy situations they bend over. In this case, the velocity of any plume parcel is the vector composition of horizontal wind velocity and vertical plume velocity in the first stage and then approaches the horizontal wind velocity. In the Eulerian dispersion models, the calculation of plume rise is based on the fluid dynamic equations, namely on the mass, momentum and energy conservation equations. A complete, exhaustive theory is not yet available. These equations are closed using the entrainment assumption proposed by Morton et al. (1956), which prescribes that the entrainment velocity, i.e., the rate at which ambient air is entrained into the plume, is proportional to the mean local rise velocity

In Lagrangian Particle Models (LPMs), plume rise can be dynamically computed, i.e., each particle, at each time step, can respond to local conditions: wind speed and direction, ambient stability and turbulence (both the self-generated and ambient ones). This allows a high degree of resolution to be obtained. Moreover, it allows simulating the interaction of a plume with a capping inversion layer and the mixture of different plumes in a “natural” way. However, introducing the plume rise computation in a LPM is not straightforward. In fact, the buoyant forces acting on the plume portions depend on the differences between their temperature and that of the background air. This 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 number of particles. This is practically impossible due to the huge amount of particles that it should be necessary to track in the computational domain. To overcome this problem, three different “hybrid” techniques (i.e., that consider together Lagrangian and Eulerian properties) have been proposed (Webster and Thomson, 2002; Anfossi and Physick, 2005): particles emission occurs at the final plume height computed by analytical models, such as the Briggs (1975) ones, an integral plume rise model is used at each time step for each particle and the derived velocities are added to the Lagrangian stochastic particle velocities (see, for instance, Anfossi et al., 1993) and a set of differential equations describing the time and space evolution of bulk plume quantities are solved at each time step (Webster and Thomson, 2002; Anfossi et al., 2010). An interesting method was proposed by van Dop (1992) in which a Langevin equation for the particle temperature is solved and the buoyancy of the particles is included in the Langevin equation for the particle velocity evolution. Though having advantages from the turbulence closure point of view, this approach has problems in the entrainment treatment. Another interesting approach comes from Bisignano and Devenish (2015), who have proposed a model with an equation for temperature fluctuations as well as for vertical velocity fluctuations.

In Alessandrini et al. (2013) a method for the buoyant plume rise computation is proposed and validated. It is based on the same idea described in papers related with chemical reactions in Lagrangian particle models (Alessandrini and Ferrero, 2009, 2011) 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. Consequently, 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. In this new approach, the Eulerian grid is used to compute the macroscopic properties of the plume locally, starting from the properties of the single particles within each cell. No equation of the bulk plume is solved on a fixed grid. Without considering the simplification of the model, this method has the advantage of being applicable also 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. A similar approach was later proposed by Kaplan (2015).

This plume rise scheme has been tested on some previous works (Alessandrini et al., 2013; Ferrero et al., 2019; Tenti and Ferrero, 2024b) and good results have been obtained. However, in some cases, this scheme produces increasing plume temperature which is clearly an unphysical situation. To overcome this problem, we present here a new plume rise algorithm for LPMs which is based on the transport of heat. When two fluids with different temperatures come into contact, the equation for reaching the equilibrium temperature through heat transfer is described by Richmann's law. The algorithm is based on the scheme proposed by Alessandrini et al. (2013), according to which each particle transports scalar quantities used to calculate the temperature difference between the plume and the ambient and the vertical velocity of the plume on a fixed Eulerian grid. The new model is compared with the Briggs analytical formula in an idealized neutral atmosphere, and it is tested against measurements taken during a water tank experiment in neutral and stable flow conditions (Huq and Stewart, 1996).

In this section, the new plume rise algorithm embedded in the LPM SPRAY-WEB is presented.

As in the previous algorithm (Alessandrini et al., 2013), the plume rise is calculated using an Eulerian grid. In Alessandrini et al. (2013), at the emission the particles are assigned with an equivalent mass of momentum (*m*_{w}) and of temperature difference (*m*_{T}) as follows

where *w*_{u} is the plume exit velocity, *A*_{s} is the source area, d$t={t}_{\mathrm{1}}-{t}_{\mathrm{0}}$ is the time step, ${N}_{{\mathrm{P}}_{\mathrm{step}}}$ is the number of particles emitted in the time step, *T*_{plume} is the temperature of the emitted plume and *T*_{a} is the ambient temperature at the source location. ${A}_{\mathrm{s}}\left|{w}_{u}\right|\mathrm{d}t$ is the volume occupied by the particles in the time step d*t*.

Afterwards, for the time *t*_{0}, the temperature difference between each cell of the Eulerian grid and the environment (Δ*T*_{c}), as well as its vertical velocity (*w*_{c}), are calculated by summing the contributions of particles contained within the cells themselves

where *M* is the number of particles within the cell and *V*_{c} is the cell volume.

After the lagrangian advancement of the particles at time *t*_{1}, the cell quantities are updated using the following equation for temperature difference

where *z*_{c} is the cell height, Γ(*z*_{c}) is the lapse rate of the ambient air at cell height *z*_{c} and 0.0098 is the adiabatic lapse rate. This temperature difference is used to retrieve the vertical velocity of the cell as follows

where *T*_{a} (*z*_{c}) is the ambient air temperature at height *z*_{c}, *g* is the gravitational acceleration, *S* is the cell horizontal cell cross section, *C*_{D} is the drag coefficient and *ρ*_{a} and *ρ*_{p} are the ambient air and plume density, respectively. To move from density to temperature (more easily handled) the ideal gas law was used. In previous works (Tenti and Ferrero, 2024a, b) the dependence of the model on the value chosen for the drag coefficient has been investigated: different formulations dependent on the Reynolds number have been tested, and most tend to asymptotically approach a value of 0.47. Here, we have chosen to use a constant value (*C*_{D}=0.5) because calculating this value at each iteration is computationally demanding for the model, and the result tends to approximate a value that is essentially constant anyway.

After the Lagrangian advancement, the temperature difference, and velocity masses of the particles are updated based on the temperature difference and velocity of the cell in which they are located at times *t*_{0} and *t*_{1}

In the new scheme, we want to simulate the transport of heat by assigning to the particles a momentum mass (*m*_{w}) as described in Eq. (1), a temperature (*T*_{p}), a mass (*m*_{p}) and a denisty (*ρ*_{p}). The temperature and the density, computed using the ideal gasses law, are those of the plume

where *P*_{atm} is the atmospheric pressure and $\stackrel{\mathrm{\u203e}}{R}$ is the gas constant for dry air. The mass is given as

where *w*_{u}*A*_{s}d*t* is the volume occupied by the particles in the first temporal step (d*t*). Here we are considering the plume as a perfect gas with the same features as the air (same gas specific constant $\stackrel{\mathrm{\u203e}}{R}$) and specific heat *c*_{p}) and we are considering a constant pressure. These assumptions are reasonable approximations since the effects introduced by considering these differences are mostly negligible.

The temperature of the cell (*T*_{c}) derives from Richmann's equation for the mixing of two fluids with different temperatures

where *c*_{1} and *c*_{2} and *T*_{1} and *T*_{2} are the specific heats and the temperatures of the two fluids respectively. As previously mentioned we are considering ${c}_{\mathrm{1}}={c}_{\mathrm{2}}={c}_{p}$ so that Eq. (11) becomes

where *m*_{a} is the mass of the environmental air contained in the cell, *m*_{tot} is the total mass contained in the cell considering the air and the plume, ${m}_{\mathrm{plume}}={\mathrm{\Sigma}}_{i=\mathrm{1}}^{M}{m}_{\mathrm{pi}}$ is the total mass of the plume contained in the cell, obtained by summing the masses of all the particles within the cell and *ρ*_{c} is the cell density. The final density of the cell can be expressed as the total mass (air + plume particles) divided by the total volume (air + plume particles)

where *V*_{a} and *V*_{p} are the volumes occupied by the air and the plume particles rispectively.

Once the cell temperature has been computed, we estimate the temperature difference between the cell and the ambient substracting the air temperature

and this temperature difference is used in Eq. (6) to calculate the vertical velocity.

(Alessandrini et al., 2013)The vertical velocity of the cell is assigned to the particles inside that cell and it is used to calculate the vertical displacement (d*z*) of the particles as follows

According to this new algorithm, the particles move adiabatically, so that their temperature changes only due to the vertical lifting

The adiabatic assumption derives from the fact that the characteristic time of the particle displacement is less than the time necessary for a heat exchange between the particles and the ambient to occur.

Table 1 highlights and summarizes the differences between the two schemes.

Following the example of Alessandrini et al. (2013), the results obtained with the new plume rise scheme are compared with the Briggs (1975) formula for the plume height for the neutral case

where *F*_{M} and *F*_{B} are the momentum and the buoyancy fluxes respectively, *t* is the time and *u* is the ambient wind speed. In order to carry out this comparison, two simulations were performed in an idealized atmosphere, characterized by a constant horizontal wind of 3 and 5 m s^{−1}. The stack height is 120 m and the diameter is 4.5 m. The plume exit velocity is 21.4 m s^{−1}, the plume temperature is 439 K and the ambient air temperature at the stack height is 298 K, with an adiabatic constant vertical temperature profile. In this case the ambient turbulence is horizontally homogeneous and vertically distributed according to Hanna (1982) parametrization.

To compare the model with experimental data we considered measurements taken during the water tank experiment by Huq and Stewart (1996), in which they replicated a neutral stratified turbulent flow condition and conducted measurements and analysis on the dispersion of a buoyant plume originating from a single point source. Since the experiments are conducted in a water tank, in order to set up a simulation with a Lagrangian model, the scaling factors are calculated based on Froude number similitude and on the assumption $\frac{\mathrm{\Delta}\mathit{\theta}}{\mathit{\theta}}=-\frac{\mathrm{\Delta}\mathit{\rho}}{\mathit{\rho}}$, where *ρ* is the water density and *θ* is the potential temperature. The turbulent variables are estimated using the standard deviations of velocities and the dissipation rate of the turbulent kinetic energy. For all our simulations we set the drag coefficient to a constant value (*C*_{D}=0.5). In Fig. 1 the results obtained in an idealized atmosphere with a constant wind speed of 3 m s^{−1} (left) and 5 m s^{−1} (right) are presented: The red circles' lines are the results of the simulations obtained with the old plume rise scheme (Alessandrini et al., 2013), the green circles lines are the simulations results with the new one, and the continuous blue line is the height obtained with the Briggs equation (Eq. 17).

In general, with the new scheme, the plume rises more, especially when the wind speed is lower (Fig. 1, left). This could be linked to the fact that we are considering the rise of the particles to be perfectly adiabatic when instead minimal heat exchange occurs. Although the plume with the new algorithm is higher than with the old one, the asymptotic trend is very similar to the Briggs curve. When the wind speed is higher (Fig. 1, right) the height of the particles simulated with the new scheme overlaps perfectly with the Briggs curve (green circles), while with the old scheme (red circles), the plume is slightly lower.

Figure 2 shows the comparison with the experimental data (Huq and Stewart, 1996) with the old scheme (left) and the new scheme (right). The three different runs are characterized by different momentum (*F*_{M}) and buoyancy fluxes (*F*_{B}) as stated in Table 2.

The black lines in the plots represent the trend of the plume height according to the Briggs curve (Eq. 17). Close to the source, the mechanical momentum prevails and the plume should follow the $\mathrm{1}/\mathrm{3}$ power law trend, while moving further away the buoyancy flux begins to prevail and therefore the plume's trend should follow the $\mathrm{2}/\mathrm{3}$ power law trend. In the plots, both the results obtained with the simulations and the measurements are scaled with the buoyancy length scale (${l}_{\mathrm{B}}=\frac{{F}_{\mathrm{B}}}{{u}^{\mathrm{3}}}$), to be able to compare the results of runs with different characteristics.

As for the near-source field we found similar behaviors with both the schemes, while they differ far from the source, in the $\mathrm{2}/\mathrm{3}$ power-law slope regime. With the old scheme (Fig. 2, left), similar results were obtained for the runs with similar buoyancy flux (3 and 9) and these, although showing a correct asymptotic trend, differ from the experimental measurements. Run 6, having a different value of *F*_{B}, tends to align with the measurements. On the other hand, with the new scheme (Fig. 2, right), all runs tend to collapse onto the same values, resulting in better agreement with the measurements.

After comparisons with the Briggs curve and experimental data, the new plume rise scheme gives similar results to the old one but it seems to better reproduce the experimental data.

In this study, we introduced a novel plume rise scheme based on heat transport by particles. This method simulates entrainment through the mixing of two fluids – air and plume particles – at different temperatures, with the resultant temperature calculated using Richmann's law. We compared this new algorithm with the existing one in SPRAY-WEB, developed by Alessandrini et al. (2013). Both approaches employ an Eulerian grid to determine plume characteristics, particularly vertical velocity and temperature differential with the external environment. The primary distinction lies in the treatment of particles: while Alessandrini et al.'s (2013) model assigns particles with momentum and temperature difference, our new algorithm also includes mass and density to better simulate heat transport. We evaluated the performance of both algorithms against the analytical Briggs curve in an idealized atmosphere and against experimental data from a water tank experiment. Our findings indicate that the new scheme performs better at higher ambient wind speeds, while maintaining correct asymptotic behavior at lower speeds. Comparisons with experimental measurements reveal that the new algorithm enhances model performance, accurately replicating experimental data regardless of the buoyancy flux involved.

The software is available through a repository on GitLAB upon request to the following address e.ferrero@isac.cnr.it.

The data are not publicly accessible because we are not the owners.

AS, FE and TB were responsible for the development of the algorithm, TB ran the simulations and FE and TB drafted the manuscript.

The contact author has declared that none of the authors has any competing interests.

This article is part of the special issue “EMS Annual Meeting: European Conference for Applied Meteorology and Climatology 2023”. It is a result of the EMS Annual Meeting: European Conference for Applied Meteorology and Climatology 2023, Bratislava, Slovakia, 3–8 September 2023.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

This paper was edited by Gert-Jan Steeneveld and reviewed by two anonymous referees.

Alessandrini, S. and Ferrero, E.: A hybrid Lagrangian-Eulerian particle model for reacting pollutant dispersion in non-homogeneous non-isotropic turbulence, Physica A, 388, 1375–1387, https://doi.org/10.1016/j.physa.2008.12.015, 2009. a

Alessandrini, S. and Ferrero, E.: A Lagrangian particle model with chemical reactions: application in real atmosphere, Int. J. Environ. Pollut., 47, 97–107, 2011. a

Alessandrini, S., Ferrero, E., and Anfossi, D.: A new Lagrangian method for modelling the buoyant plume rise, Atmos. Environ., 77, 239–249, 2013. a, b, c, d, e, f, g, h, i, j

Anfossi, D. and Physick, W.: Lagrangian particle models, in: Air Quality Modeling Theories, Methodologies, Computational Techniques, and Available Databases and Software, Fundamentals, vol. II, edited by: Zannetti, P., Chap. 11, The EnviroComp Institute and the Air & Waste Management Association, The ISBN 0-923204-86-5, 2005. a

Anfossi, D., Ferrero, E., Brusasca, G., Marzorati, A., and Tinarelli, G.: A simple way of computing buoyant plume rise in a lagrangian stochastic model for airborne dispersion, Atmos. Environ., 27, 1443–1451, 1993. a

Anfossi, D., Tinarelli, G., Trini Castelli, S., Nibart, M., Olry, C., and Commanay, J.: A new Lagrangian particle model for the simulation of dense gas dispersion, Atmos. Environ., 44, 753–762, 2010. a

Bisignano, A. and Devenish, B.: A model for temperature fluctuations in a buoyant plume, Bound.-Lay. Meteorol., 157, 157–172, 2015. a

Briggs, G.: Plume rise predictions, in: Lectures on Air Pollution and Environmental Impact Analyses, edited by: Haugen, D., Workshop Proceedings, 29 September–3 October, Boston, Massachusetts, American Meteorological Society, Boston, Massachusetts, USA, 36, 59–111, ISBN 978-1-93-57-04-23-2, 1975. a, b, c

Ferrero, E., Alessandrini, S., Anderson, B., Tomasi, E., Jimenez, P., and Meech, S.: Lagrangian simulation of smoke plume from fire and validation using ground-based lidar and aircraft measurements, Atmos. Environ., 213, 659–674, 2019. a

Hanna, S. R.: Applications in Air Pollution Modeling, Springer Netherlands, Dordrecht, 275–310, ISBN 978-94-010-9112-1, 1982. a

Huq, P. and Stewart, E.: A laboratory study of buoyant plumes in laminar and turbulent crossflows, Atmos. Environ., 30, 1125–1135, https://doi.org/10.1016/1352-2310(95)00335-5, 1996. a, b, c, d

Kaplan, H.: A Lagrangian Particle Transport and Diffusion Model for Non-passive Scalars, Bound.-Lay. Meteorol., 156, 37–52, 2015. a

Morton, B. R., Sir Geoffrey Taylor, F. R. S., and Turner, J. S.: Turbulent gravitational convection from maintained and instantaneous sources, Proceeding Royal Society of London A, 234, 1–23, https://doi.org/10.1098/rspa.1956.0011, 1956. a

Tenti, B. and Ferrero, E.: Evaluation of Turbulence Depending Drag Coefficient in Plume Rise Model for Fire Smoke Dispersion, Atmos. Environ., 323, 120411, https://doi.org/10.1016/j.atmosenv.2024.120411, 2024a. a

Tenti, B. and Ferrero, E.: Smoke plume from fire Lagrangian simulation: dependence on drag coefficient and resolution, Air Quality, Atmosphere & Health, 17, 1–10, 2024b. a, b

van Dop, H.: Buoyant plume rise in a Lagrangian framework, Atmos. Environ., 26, 335–1346, 1992. a

Webster, H. and Thomson, D.: Validation of a Lagrangian model plume rise scheme using the Kincaid data set, Atmos. Environ., 36, 5031–5042, 2002. a, b

A new plume rise scheme based on heat transport by particles was presented: the entrainment is simulated by the mixing of two fluids (air and plume particles) with different temperatures and the resulting temperature is given by Richmann's law. The new algorithm is compared with the one that is currently included in SPRAY-WEB by Alessandrini et al. (2013). The new scheme seems to behave better when the ambient wind speeds are higher, but the asymptotic behavior is correct even with lower speeds.

A new plume rise scheme based on heat transport by particles was presented: the...