An Analytic Model for Dispersion of Rocket Exhaust Clouds: Specifications and Analysis in Different Atmospheric Stability Conditions

ABSTRACT:

This paper brings the first results concerning a new analytic model to evaluate atmospheric dispersion in rocket launch scenarios, namely Generalized Integral Laplace Transform Technique for Rocket effluent dispersion model. The model is constituted by three different modules, being two pre-processors (micrometeorological parameters and deposition parameters) and the dispersion program. The dispersion calculations are made through the Generalized Integral Laplace Transform Technique solution of the two-dimensional, time-dependant, advectiondiffusion-deposition equation. The results show a set of simulations using data from the Chuva Project from the Centro de Lançamento de Alcântara for both stable and unstable planetary boundary layers, in order to evaluate model performance and illustration.

KEYWORDS:
Rocket effluent; Pollution dispersion; Analytical model

INTRODUCTION

The advance of science and technology over the last decades has, undoubtedly, led to countless benefits to mankind. However, the residuals generated in many processes relative to these improvements are worldwide distributed, affecting the environment as a whole. In order to prevent and remedy environmental damages due to human activity, it is necessary that science also provides tools to assess and manage the possible negative effects of its residuals on the Earth's system. One example of these residuals is the air pollution, which has a great anthropogenic component. We may be aware of some human sources of air pollution such as automotive and industrial releases, biomass burn and so on. There are also other sources that may not be not so well known or studied, but are significant. One example is the spacecraft, which emits a large amount of pollution to the atmosphere in a few seconds during its launch. Some of the contaminants released (in a huge quantity) during a rocket launch are considered extremely harmful to human health even at low concentrations, such as the HCl. On the other hand, there are pollutants that become dangerous due to the enormous concentrations that might remain up to 30 minutes in the lower troposphere, as CO and Al₂O₃.

The rocket launch process includes the immediate, previous and ahead instants to the launch. Within this time interval, combustion occurs, in which a large hot and buoyant cloud is formed and ascends until it reaches thermal equilibrium with the atmosphere. This referred cloud is called 'ground-cloud', because it is generated near the ground level (Nyman 2009), and in normal launch conditions (excepting explosions and engine trials) it concerns a short-term release, in the order of 10 s (Bjorklund et al. 1982).

The modelling of rocket exhaust generally concerns this cloud of pollutants. It encompasses a wide range of complex atmospheric phenomena, such as turbulence and buoyancy. In the USA, the Rocket Exhaust Effluent Diffusion Model (REEDM) is traditionally used in this task. REEDM is a Gaussian dispersion model constituted by several algorithms designed to assess dosage, peak concentration and deposition. It is also used for calculating cloud and plume rise and turbulence (Bjorklund et al. 1982). In Brazil, a model called Modelo Simulador de Efluentes de Foguetes (MSDEF) or Model for Simulation of Rocket Effluents (in a free translation) was recently developed, which is based on the mathematical procedure Advection Diffusion Multilayer Model (ADMM) to solve the advectiondispersion equation in a semi-analytic fashion, also incorporating some of the North American assumptions (Moreira et al. 2011). Bianconi and Tamponi (1993) developed an analytical model for the unsteady three-dimensional advection, diffusion equation for short-term releases, and tested its sensitiveness in different stability conditions, mixing layer heights and wind speeds for different release durations.

This work uses the Generalized Integral Laplace Transform Technique (GILTT) instead of the ADMM method in order to simulate rocket exhaust dispersion. The GILTT method consists on some few standard steps (Moreira et al. 2009): expand the concentration in series based on the eigenfunctions of an auxiliary problem, using cosine functions; replace this expansion on the advection-diffusion equation and take moments, leading to a matrix ordinary differential equation; solve the equation by the Laplace transform, analytically. In 2012, Buske et al. formulated a three-dimensional solution for the stationary advection-diffusion equation using GILTT, for general vertical profiles of wind and eddy diffusivity, which proved valid against Kinkaid and Copenhagen experimental datasets. In 2013, Gonçalves et al. developed a solution for the two-dimensional stationary advection-diffusion equation though GILTT, using Bessel functions instead of cosine functions, as base to expand the concentrations in series, and also the source condition was different. The solution was tested against Copenhagen, Praire-Grass and Hanford experiments and against the results obtained with the cosine function base. The Bessel function approach showed very good results against experimental data and proved to numerically converge faster than the cosine approach.

Here, we display the solution, via GILTT, of the transient two-dimensional advection-diffusion equation for a finite (shortterm) release. It has already been tested against the experimental Hanford-83 and Copenhagen datasets, with success, for a first test and demonstration of the gravitational settling importance to an accurate representation of physical processes, such as dry deposition (Bainy et al. 2015). Due to its application to rocket lauches, it was named Generalized Integral Transform Technique for Rocket effluent dispersion model (GILTTR). This paper presents the mathematical model for dispersion, deposition and micrometeorological parameters, and the application of the model in different atmospheric stability conditions.

METHODOLOGY

THE MATHEMATICAL MODEL

The model considers the dispersion of a pollutant according to Eq. 1.

(1) $\frac{\partial \mathit{c}}{\partial \mathit{t}}+\mathit{u}\frac{\partial \mathit{c}}{\partial \mathit{x}}-{\mathit{v}}_{\mathit{g}}\frac{\partial \mathit{c}}{\partial \mathit{z}}=\frac{\partial }{\partial \mathit{x}}\left({\mathit{K}}_{\mathit{x}}\frac{\partial \mathit{c}}{\partial \mathit{x}}\right)+\frac{\partial }{\partial \mathit{z}}\left({\mathit{K}}_{\mathit{z}}\frac{\partial \mathit{c}}{\partial \mathit{z}}\right)-\mathit{\lambda }\mathit{c}-{\wedge }_{\mathit{s}}\mathit{c}$

where:

c = c(x,z,t) is the crosswind integrated contaminant concentration; u is the mean wind speed; v_g is the gravitational settling of the specie; κ_x and κ_z are, respectively, the longitudinal and vertical eddy coefficients (function of z alone); λ is the physical-chemical decay coefficient; λs is the scavenging coefficient. The tri-dimensional concentration (c*) profile is assumed to be Gaussian-shaped, as in:

(2) $\mathit{c}*\left(\mathit{x},\mathit{y},\mathit{x},\mathit{t}\right)=\mathit{c}\left(\mathit{x},\mathit{z},\mathit{t}\right)·\frac{1}{\sqrt{2{\mathit{\pi }\mathit{\sigma }}_{\mathit{y}}}}{\mathit{e}}^{-\left[\frac{{\left(\mathit{y}-{\mathit{y}}_0\right)}^2}{2{\mathit{\sigma }}_{\mathit{y}}^2}\right]}$

The initial and boundary conditions are, respectively, given by:

(3) $\mathit{c}\left(\mathit{x},\mathit{z},0\right)=0$

(4) $\begin{array}{c}{\mathit{K}}_{\mathit{z}}\frac{\partial \mathit{c}}{\partial \mathit{z}}={\mathit{V}}_{\mathit{d}}\mathit{c}\mathrm{at}z=z_0\hfill \\ {\mathit{K}}_{\mathit{z}}\frac{\partial \mathit{c}}{\partial \mathit{z}}=0\mathrm{at}z=h,L_{*}\hfill \end{array}$

where:

V_d is the deposition velocity at the surface; z0 is the roughness length; h is the planetary boundary layers (PBL) height; L_* is a distance far away from the source.

The source is punctual and established in a height H_s, and its condition is given in Eq. 5.

(5) $\mathit{c}\left(0,\mathit{z},\mathit{t}\right)=\frac{\mathit{Q}}{\mathit{u}}\left[\mathit{\eta }\left(\mathit{t}\right)-\mathit{\eta }\left(\mathit{t}-{\mathit{t}}_{\mathit{r}}\right)\mathit{\delta }\left(\mathit{z}-{\mathit{H}}_{\mathit{s}}\right)\right]$

where:

Q is the emission rate; t_r means the duration of release; H_s is the source height; η is the Heaviside step function; δ is the Dirac delta function.

The advection-diffusion equation, given by Eq. 1, is solved by GILTT (Moreira et al. 2009). The first step is to apply the Laplace transform in Eq. 1, resulting in:

(6) $\begin{array}{c}\mathit{r}\overline{\mathit{c}}-\mathit{c}\left(\mathit{x},\mathit{z},0\right)+\mathit{u}\frac{\partial \overline{\mathit{c}}}{\partial \mathit{x}}-{\mathit{v}}_{\mathit{g}}\frac{\partial \overline{\mathit{c}}}{\partial \mathit{z}}=\hfill \\ \frac{\partial }{\partial \mathit{x}}\left(K_{\mathit{x}}\frac{\partial \overline{\mathit{c}}}{\partial \mathit{x}}\right)+\frac{\partial }{\partial \mathit{z}}\left(K_{\mathit{z}}\frac{\partial \overline{\mathit{c}}}{\partial \mathit{z}}\right)-\mathit{\lambda }\overline{\mathit{c}}-{\wedge }_{\mathit{s}}\overline{\mathit{c}}\hfill \end{array}$

where:

x, z, r) = L{c (x, z, t); t → r}). The next step is to define the Sturm-Liouville auxiliary problem, and then expand the concentration in a sum. This auxiliary problem is chosen as function of the boundary conditions of the original problem. In this case it will be: is the Laplace transform of concentration in the variable t ( (

(7) ${{\Psi }^{″}}_{\mathit{n}}\left(\mathit{z}\right)+{\mathit{\lambda }}_{\mathit{n}}^2{\Psi }_{\mathit{n}}\left(\mathit{z}\right)=0\mathrm{for}{z}_0<z<h,$

With the same boundary conditions of the original problem, it is:

(8) $\begin{array}{c}{{\Psi }^{\prime }}_{\mathit{n}}\left(\mathit{z}\right)=0\mathrm{for}z=h\hfill \\ {\mathit{K}}_{\mathit{z}}{{\Psi }^{\prime }}_{\mathit{n}}\left(\mathit{z}\right)-{\mathit{V}}_{\mathit{d}}{\Psi }_{\mathit{n}}\left(\mathit{z}\right)=0\mathrm{for}z=z_0\hfill \end{array}$

The solution of the auxiliary problem serves as base to the expansion of the concentration in series and is given by:

(9) ${\Psi }_{\mathit{n}}\left(\mathit{z}\right)=\mathit{cos}{\mathit{\lambda }}_{\mathit{n}}\left(\mathit{z}-\mathit{h}\right)$

where:

λ_n refers to the eingenvalues of the auxiliary problem and satisfies the transcendental equation:

(10) $\frac{{\mathit{V}}_{\mathit{d}}}{{\mathit{K}}_{\mathit{z}}\left({\mathit{z}}_0\right)}={\mathit{\lambda }}_{\mathit{n}}\mathit{tan}\left[{\mathit{\lambda }}_{\mathit{n}}\left(\mathit{h}-{\mathit{z}}_0\right)\right]$

that is solved by the Newton-Raphson method. This way, it is possible to expand the concentration as:

(11) $\overline{\mathit{c}}\left(\mathit{x},\mathit{z},\mathit{r}\right)=\sum _{\mathit{n}=0}^{\infty }{\overline{\mathit{c}}}_{\mathit{n}}\left(\mathit{x},\mathit{r}\right){\Psi }_{\mathit{n}}\left(\mathit{z}\right)$

The problem now consists in determining the values of _n(x, r). Substituting this equivalence in Eq. 6, taking moments (applying the operator ∫₀^h () Ψ_m dz, where Ψ_m is a function orthogonal to Ψ_n), truncating the sums to N terms, applying the derivation rules and rewriting the terms, we obtain:

(12) $\begin{array}{c}\sum _{\mathit{n}=0}^{\mathit{N}}\frac{\partial {\mathit{c}}_{\mathit{n}}\left(\mathit{x},\mathit{r}\right)}{\partial \mathit{x}}\underset{0}{\overset{\mathit{h}}{\int }}\mathit{u}\left(\mathit{z}\right){\Psi }_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}-\hfill \\ {\mathit{v}}_{\mathit{g}}\sum _{\mathit{n}=0}^{\mathit{N}}{\mathit{c}}_{\mathit{n}}\left(\mathit{x},\mathit{r}\right)\underset{0}{\overset{\mathit{h}}{\int }}{{\Psi }^{\prime }}_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}=\hfill \\ \sum _{\mathit{n}=0}^{\mathit{N}}\frac{{\partial }^2{\mathit{c}}_{\mathit{n}}\left(\mathit{x},\mathit{r}\right)}{\partial {\mathit{x}}^2}\underset{0}{\overset{\mathit{h}}{\int }}{\mathit{K}}_{\mathit{x}}\left(\mathit{z}\right){\Psi }_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}+\hfill \\ \sum _{\mathit{n}=0}^{\mathit{N}}{\mathit{c}}_{\mathit{n}}\left(\mathit{x},\mathit{r}\right)\underset{0}{\overset{\mathit{h}}{\int }}{{\mathit{K}}^{\prime }}_{\mathit{z}}\left(\mathit{z}\right){{\Psi }^{\prime }}_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}+\hfill \\ \sum _{\mathit{n}=0}^{\mathit{N}}{\mathit{c}}_{\mathit{n}}\left(\mathit{x},\mathit{r}\right)\underset{0}{\overset{\mathit{h}}{\int }}{\mathit{K}}_{\mathit{z}}\left(\mathit{z}\right){{\Psi }^{″}}_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}+\hfill \\ \left(\mathit{\lambda }+{\wedge }_{\mathit{s}}+\mathit{r}\right)\sum _{\mathit{n}=0}^{\mathit{N}}{\mathit{c}}_{\mathit{n}}\left(\mathit{x},\mathit{r}\right)\underset{0}{\overset{\mathit{h}}{\int }}{\Psi }_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}\hfill \end{array}$

This equation may be written in the matrix fashion, such as:

(13) ${\mathit{Y}}^{″}\left(\mathit{x},\mathit{r}\right)+\mathit{F}{\mathit{Y}}^{\prime }\left(\mathit{x},\mathit{r}\right)+\mathit{GY}\left(\mathit{x},\mathit{r}\right)=0$

where the unknowns are vectors whose elements are Y(x,r) = c_n(x,r), Y'(x,r) = c_n'(x,r), Y'(x,r) = c_n'(x,r), and F = B^-1D and G = B^-1E are matrixes constituted by:

(14) $\begin{array}{c}B=\left\{b_{n,m}\right\}={\int }_0^{\mathit{h}}{\mathit{K}}_{\mathit{x}}\left(\mathit{z}\right){\Psi }_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}\hfill \\ D=\left\{d_{n,m}\right\}=-{\int }_0^{\mathit{h}}\mathit{u}\left(\mathit{z}\right){\Psi }_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}\hfill \\ E=\left\{e_{n,m}\right\}={\mathit{v}}_{\mathit{g}}{\int }_0^{\mathit{h}}{{\Psi }^{\prime }}_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}-\hfill \\ {\mathit{\lambda }}_{\mathit{n}}^2{\int }_0^{\mathit{h}}{\mathit{K}}_{\mathit{z}}\left(\mathit{z}\right){\Psi }_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}-\hfill \\ \left(\mathit{\lambda }+{\wedge }_{\mathit{s}}+\mathit{r}\right){\int }_0^{\mathit{h}}{\Psi }_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}+{\int }_0^{\mathit{h}}{{\mathit{K}}^{\prime }}_{\mathit{z}}\left(\mathit{z}\right){{\Psi }^{\prime }}_{\mathit{n}}\left(\mathit{z}\right){\Psi }_{\mathit{m}}\left(\mathit{z}\right)\mathit{dz}\hfill \end{array}$

reminding that Ψ_n"(z) = λ_n² Ψ_n (z).

Following the work of Moreira et al. (2009), the next step is to reduce the order of Eq. 13, leading to:

(15) ${\mathit{Z}}^{\prime }\left(\mathit{x},\mathit{r}\right)+\mathit{H}·\mathit{Z}\left(\mathit{x},\mathit{r}\right)=0$

where:

Z(x, r) = col Yx, r), Y'(x, r)], and the matrix H has a block form H = .

Applying the same procedures to the source condition (Eq. 5), we obtain:

(16) $\mathit{Y}\left(0,\mathit{r}\right)=\mathit{Q}\left[\frac{1}{\mathit{r}}-\frac{{\mathit{e}}^{-\mathit{r}.{\mathit{t}}_{\mathit{r}}}}{\mathit{r}}\right]\frac{\Psi \left({\mathit{H}}_{\mathit{s}}\right)}{\mathit{u}\left({\mathit{H}}_{\mathit{s}}\right)}{\mathit{A}}^{-1}=\mathit{Z}\left(0,\mathit{r}\right)$

where:

A is a diagonal matrix given by A = a_n,m = ∫₀^hΨ_nΨ_m dz.

Applying the Laplace transform in Eq. 15, we obtain:

(17) $\mathit{s}\overline{\mathit{Z}}\left(\mathit{s},\mathit{r}\right)-\mathit{Z}\left(0,\mathit{r}\right)+\mathit{H}·\overline{\mathit{Z}}\left(\mathit{s},\mathit{r}\right)=0$

where:

s, r) is the Laplace transform of the vector x, r). Assuming that the matrix H has distinct eigenvalues, we may write H = X.D.X^-1, where D is the diagonal matrix of the eigenvalues and X is the respective matrix of eigenfunctions. Replacing this relation in Eq. 17 and rearranging the terms, one leads to: ( (

(18) $\left(\mathit{sI}+\mathit{X}·\mathit{D}·{\mathit{X}}^{-1}\right)·\overline{\mathit{Z}}\left(\mathit{s},\mathit{r}\right)=\mathit{Z}\left(0,\mathit{r}\right)$

where:

I is the identity matrix (since X.X^-1 = 1).

In order to solve Eq. 18 for Z(x, r), we must apply the inverse Laplace transform, as in:

(19) $\mathit{Z}\left(\mathit{x},\mathit{r}\right)=\mathit{X}·{ℒ}^{-1}\left\{{\left(\mathit{sI}+\mathit{D}\right)}^{-1}\right\}·{\mathit{X}}^{-1}·\mathit{Z}\left(0,\mathit{r}\right).$

The term to which the inverse transform was applied is:

(20) ${ℒ}^{-1}\left\{{\left(\mathit{sI}+\mathit{D}\right)}^{-1}\right\}=\mathit{P}\left(\mathit{x},\mathit{y}\right)=\left(\begin{array}{cccc}{\mathit{e}}^{-{\mathit{d}}_1\mathit{x}}& 0& \cdots & 0\\ 0& {\mathit{e}}^{-{\mathit{d}}_2\mathit{x}}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\mathit{e}}^{-{\mathit{d}}_{\mathit{N}}\mathit{x}}\end{array}\right)$

The solution is, then, given by:

(21) $\mathit{Z}\left(\mathit{x},\mathit{r}\right)=\mathit{X}·\mathit{P}\left(\mathit{x},\mathit{r}\right)·{\mathit{X}}^{-1}·\mathit{Z}\left(0,\mathit{r}\right)=\mathit{M}\left(\mathit{x},\mathit{r}\right)$

where:

M(x,r) = X.P(x,r) and ξ = X^-1Z(0,r).

Alternatively, we may rewrite Eq. 21 as:

(22) $\left(\begin{array}{c}\mathit{Y}\left(\mathit{x},\mathit{y}\right)\\ {\mathit{Y}}^{\prime }\left(\mathit{x},\mathit{y}\right)\end{array}\right)=\left(\begin{array}{cc}{\mathit{M}}_{11}\left(\mathit{x},\mathit{r}\right)& {\mathit{M}}_{12}\left(\mathit{x},\mathit{r}\right)\\ {\mathit{M}}_{21}\left(\mathit{x},\mathit{r}\right)& {\mathit{M}}_{22}\left(\mathit{x},\mathit{r}\right)\end{array}\right)\left(\begin{array}{c}{\mathit{\xi }}_1\\ {\mathit{\xi }}_2\end{array}\right)$

And to obtain the unknown vector ξ, we must solve the linear system:

(23) $\left(\begin{array}{cc}{\mathit{M}}_{11}\left(0,\mathit{r}\right)& {\mathit{M}}_{12}\left(0,\mathit{r}\right)\\ {\mathit{M}}_{21}\left({\mathit{L}}_{*},\mathit{r}\right)& {\mathit{M}}_{22}\left({\mathit{L}}_{*},\mathit{r}\right)\end{array}\right)\left(\begin{array}{c}{\mathit{\xi }}_1\\ {\mathit{\xi }}_2\end{array}\right)=\left(\begin{array}{c}\mathit{Y}\left(0,\mathit{r}\right)\\ {\mathit{Y}}^{\prime }\left({\mathit{L}}_{*},\mathit{r}\right)\end{array}\right)$

So far, the only numeric approximation concerns the truncation of the summation in Eq. 11. In this work, N was established as 90. The solution given in Eq. 22 leads to Y(x, r), or _n (x,r). In order to obtain c (x, z, t), it is necessary to multiply the solution by the inverse GILTT (Eq. 11) and invert the Laplace transform in the time variable. This, here, was done through the Gaussian quadrature scheme, as:

(24) $\mathit{c}\left(\mathit{x},\mathit{z},\mathit{t}\right)=\sum _{\mathit{k}=1}^{\mathit{M}}\frac{{\mathit{P}}_{\mathit{k}}}{\mathit{t}}{\mathit{A}}_{\mathit{k}}\sum _{\mathit{n}=0}^{\mathit{N}}{\overline{\mathit{c}}}_{\mathit{n}}\left(\mathit{x},\frac{{\mathit{P}}_{\mathit{k}}}{\mathit{t}}\right){\Psi }_{\mathit{n}}\left(\mathit{z}\right)$

where: M is the order of the quadrature; A_k and P_k are respectively the weights and roots of the tabulated quadrature scheme (Stroud and Secrest 1966).

More details on the solution and the inversion of the Laplace transform procedure may be found in the work of Moreira et al. (2009).

PARAMETERIZATIONS

As mentioned before, the model is constituted by three different modules: micrometeorology, deposition and dispersion. The first two are pre-processors, which receive information from radiosondes and generate the necessary input for the dispersion model. A simplified flow-chart is shown in Fig. 1, highlighting the main calculations (from several ones, as will be seen) of each module from the model.

Figure 1
Flow-chart of the schematic model with the main features.

Micrometeorological Parameters

This module calculates stability and scale parameters, requiring for such wind and temperature profile (from two different surface boundary layer levels) data and PBL height. Ideally, this basic data may be obtained on a meteorological sounding or even an output from a weather numeric model, and all the other necessary atmospheric parameters are obtainable from them. The quantities searched here are: Monin-Obukhov length, Richardson number, friction velocity, convection velocity and temperature scale. The Richardson number is calculated as:

(25) $\mathit{Ri}=\frac{\mathit{g}}{\mathit{T}}\frac{\frac{\Delta \mathit{T}}{\Delta \mathit{Z}}}{{\left(\frac{\Delta \mathit{u}}{\Delta \mathit{Z}}\right)}^2}$

where:

g is the gravitational acceleration (9.81 m/s²); T is the potential temperature and the gradients are calculated as finite differences in a linear approximation.

The stability parameter ζ is given by:

(26) $\mathit{\zeta }\equiv \frac{\mathit{z}}{\mathit{L}}=\left\{\begin{array}{c}\mathit{Ri},\mathit{if}\mathit{Ri}<0\\ \frac{\mathit{Ri}}{1-5\mathit{Ri}},\mathit{if}0\le \mathit{Ri}\le 0.2\end{array}\right\}$

where L is the Monin-Obukhov length, calculated as:

(27) $\mathit{L}=\left\{\begin{array}{c}\frac{\overline{\mathit{z}}}{\mathit{Ri}},\mathit{if}\mathit{Ri}<0\\ \frac{\overline{\mathit{z}}\left(1-5\mathit{Ri}\right)}{\mathit{Ri}},\mathit{if}0\le \mathit{Ri}\le 0.2\end{array}\right\}$

where:

). is the geometric mean of the heights from which wind and temperature data were collected ( =

The non-dimensional gradients of heat and momentum (Ф_h and Ф _m, respectively) are calculated through similarity relations such as:

(28) $\begin{array}{c}{\Phi }_{\mathit{m}}={\Phi }_{\mathit{h}}=\left(1+4.7\frac{\mathit{z}}{\mathit{L}}\right)\hfill \\ \mathrm{for}\mathrm{stable}\mathrm{conditions}\left(z/L\ge 0\right),\mathrm{or}\hfill \\ \Phi _{\mathit{m}}{}^2={\Phi }_{\mathit{h}}=\left(1-15\frac{\mathit{z}}{\mathit{L}}\right)\hfill \end{array}$

for unstable conditions (z/L < 0).

Finally, the friction, convective velocity and temperature scales are computed respectively by:

(29) ${\mathit{u}}_{*}=\frac{\mathit{k}·\mathit{z}}{{\Phi }_{\mathit{m}}}\frac{\Delta \mathit{u}}{\Delta {\mathit{z}}^{\prime }}$

(30) ${\mathit{w}}_{*}={\left(\frac{-{\mathit{gu}}_{*}{\mathit{\theta }}_{*}\mathit{h}}{\mathit{T}}\right)}^{1/3}$

(31) ${\mathit{\theta }}_{*}=\frac{\mathit{k}·\mathit{z}}{{\Phi }_{\mathit{h}}}\frac{\Delta \mathit{T}}{\Delta {\mathit{z}}^{\prime }}$

where:

k is the Von Karmán constant (0.4).

Deposition Parameters

In order to provide the necessary data of gravitational settling for Eq. 1 and to satisfy the boundary conditions of the original and auxiliary problems, a deposition module was built according to a resistance analogy (Seinfeld and Pandis 2006). In this analogy, the deposition velocity is inversely proportional to the sum of all resistances offered by the Earth-atmosphere system; it is: V_d = (r_a + r_b + r_c)^-1 + V_g. Generally, the resistances are functions of PBL stability and physical state (particulate or gaseous) of the pollutant. This module is based on the works of Zhang et al. (2001) and Seinfeld and Pandis (2006), and, despite the several calculations this module encompasses, the ultimately necessary ones are the deposition velocity (V_d) and gravitational settling (V_g).

The aerodynamic resistance (r_a) is formulated by:

(32) ${\mathit{r}}_{\mathit{a}}=\frac{\mathit{ln}\left({\mathit{Z}}_{\mathit{ref}}/{\mathit{Z}}_0\right)-{\Psi }_{\mathit{h}}}{{\mathit{ku}}_{*}}$

where:

Z₀ is the roughness (seasonal, here assumed 0.06 and 0.19 for the dry and rainy seasons, respectively; Roballo and Fisch 2008); Z_ref is the reference height where the parameter is being evaluated; Ψh is a stability function given by:

(33) $\begin{array}{c}{\Psi }_{\mathit{h}}=\left\{\begin{array}{c}\frac{-5{\mathit{Z}}_{\mathit{ref}}}{\mathit{L}}\mathit{if}\mathit{\zeta }\ge 0\\ \mathit{exp}[0.598+0.309\mathit{ln}\left(-{\mathit{Z}}_{\mathit{ref}}/\mathit{L}\right)-\end{array}\right.\hfill \\ 0.09{\left(\mathit{ln}\left(-{\mathit{Z}}_{\mathit{ref}}/\mathit{L}\right)\right)}^2]\mathit{if}\mathit{\zeta }<0\hfill \end{array}$

The quasi-laminar resistance (r_b) for gaseous contaminants is written as:

(34) ${\mathit{r}}_{\mathit{b}}=\frac{5{\mathit{Sc}}^{2/3}}{{\mathit{u}}_{*}}$

where:

Sc = ν/D is the Schidt number; ν = 1.5 x 10^-5 m²/s is the viscosity of air; D is the species molecular diffusion (1.38 x 10^-5 m²/s for CO₂, 1.81 x 10^-5 m²/s for CO, and 1.5 x 10^-5 m²/s for the other gases, according to Massman (1998) and ADMS-CERC (2012). For particulates, it is expressed by:

(35) ${\mathit{r}}_{\mathit{b}}=\frac{1}{{\epsilon }_0{\mathit{u}}_{*}\left({\mathit{E}}_{\mathit{B}}+{\mathit{E}}_{\mathit{IM}}+{\mathit{E}}_{\mathit{IN}}\right){\mathit{R}}_1}$

where: ε0 is a constant equal to 3; E_B, E_IM and E_IN are the collection efficiencies for the mechanisms of Brownian diffusion, impaction and interception, respectively; R₁ is a correction factor, here considered as 1. The efficiencies are computed as:

(36) ${\mathit{E}}_{\mathit{B}}={\mathit{Sc}}^{-\mathit{\gamma }}$

where:

γ is a constant between 1/2 and 2/3, here adopted as 0.56 as function of Alcântara soil use and vegetation (Zhang et al. 2001); Sc = ν/D, where D is now the Brownian diffusivity, given by:

(37) $\mathit{D}=\frac{{\mathit{k}}_{\mathit{b}}{\mathit{TC}}_{\mathit{c}}}{3\mathit{\pi }\mathit{\mu }{\mathit{D}}_{\mathit{p}}}$

where:

k_b = 1.38 x 10^-23 J/K (Boltzmann constant); μ = 1.8 x 10^-5 kg/ms (air dynamic viscosity); T is the temperature (K); D_p is the particle diameter; C_c is the Cunningham correction factor, function of the molecules mean free path in air (λ_a = 6.7 x 10^-8 m) and particle diameter:

(38) ${\mathit{C}}_{\mathit{c}}=1+\frac{2{\mathit{\lambda }}_{\mathit{a}}}{{\mathit{D}}_{\mathit{p}}}\left(1.257+0.4{\mathit{e}}^{-\frac{0.55{\mathit{D}}_{\mathit{p}}}{{\mathit{\lambda }}_{\mathit{a}}}}\right)$

The impaction efficiency coefficient is computed by:

(39) ${\mathit{E}}_{\mathit{IM}}=\frac{{\mathit{St}}^2}{1+{\mathit{St}}^2}$

where:

St = V_gu_* /gA is the Stokes number for vegetated surfaces as function of the settling velocity (V_g), to be presented, gravitational acceleration, friction velocity and a characteristic radius, A, tabulated according to the land use and seasonal category (Zhang et al. 2001). For Centro de Lançamento de Alcântara (CLA), this coeficient was estimated as 2 mm. The gravitational settling (settling velocity) for small particles is defined by the Stokes' law:

(40) ${\mathit{V}}_{\mathit{g}}=\frac{\mathit{\rho }{\mathit{gD}}_{\mathit{p}}^2{\mathit{C}}_{\mathit{c}}}{18\mathit{\mu }}$

The interception collection efficiency is computed by:

(41) ${\mathit{E}}_{\mathit{IN}}=\frac{1}{2}{\left(\frac{{\mathit{D}}_{\mathit{p}}}{\mathit{A}}\right)}^2$

The surface resistance (r_c), a complex issue, was treated in a simple way in this work. For gases poorly reactive to the surface, it was assumed r_c = 1,000 s/m; for reactive gases, r_c = 30 s/m; and for parti culate matter, r_c = r_a × r_b × V_g. A more detailed discussion may be found in Seinfeld and Pandis (2006).

Turbulent Parameterizations

The simulations were carried out using the Degrazia formulation for the turbulent diffusion. To stable conditions the chosen vertical eddy diffusivity is given by (Degrazia et al 2000).

(42) ${\mathit{K}}_{\mathit{z}}=\frac{0.3\left(1-\frac{\mathit{z}}{\mathit{h}}\right){\mathit{u}}_{*}\mathit{z}}{1+3.7\left(\frac{\mathit{Z}}{\wedge }\right)}$

where:

λ = L(1 - z/h)^5/4 ; h is the PBL height; L is the Monin-Obukhov length.

To an unstable PBL, it is given by (Degrazia et al. 1997):

(43) $\begin{array}{c}{\mathit{K}}_{\mathit{z}}=0.22\left({\mathit{w}}_{*}\mathit{h}\right){\left(\frac{\mathit{z}}{\mathit{h}}\right)}^{1/3}{\left(1-\frac{\mathit{z}}{\mathit{h}}\right)}^{1/3}\hfill \\ \left[1-{\mathit{e}}^{-4\mathit{z}/\mathit{h}}-0.0003{\mathit{e}}^{8\mathit{z}/\mathit{h}}\right]\hfill \end{array}$

The horizontal eddy diffusivity (K_x) was neglected in this work, being considered of small contribution in comparison to mean wind horizontal advection.

The mean wind vertical profile is given by the power law that is formulated as (Panofsky and Dutton 1984):

(44) $\frac{\overline{\mathit{u}}}{{\overline{\mathit{u}}}_1}={\left(\frac{\mathit{z}}{{\mathit{z}}_1}\right)}^{\mathit{\alpha }}$

where:

₁ are the mean wind speed at the levels z and z₁; α is a coefficient related to atmospheric turbulence. For CLA, α was computed as a value between 0.19 and 0.27 for dry and rainy seasons, respectively (Roballo and Fisch 2008). In this paper, α was admitted as 0.2. and

The lateral dispersion coefficient was computed as:

(45) ${\mathit{\sigma }}_{\mathit{y}}=\frac{{\mathit{\sigma }}_{\mathit{v}}{\mathit{xS}}_{\mathit{y}}}{\mathit{u}}$

where:

x is the source distance; u is the mean wind speed; S_y(x) = (1 + 0.0308x^0.4548)^-1 and the horizontal wind fluctuation is:

${\mathit{\sigma }}_{\mathit{v}}=\left\{\begin{array}{c}1.92{\mathit{u}}_{*}\mathit{if}\mathit{\zeta }\ge 0\\ {\mathit{u}}_{*}{\left(12-\frac{0.5\mathit{h}}{\mathit{L}}\right)}^{1/3}\mathit{if}\mathit{\zeta }<0\end{array}\right.$

RESULTS AND DISCUSSION

For test cases it was chosen an emission of the particulate Al₂O₃ (that corresponds to 28.2% of total mass of pollution and whose density is 3,950 kg/m³), with assumed diameter of 2.5 μm. The source was established at a height of 150 m with and emission rate of Q = 5.2 x 10⁵ g/s¹ and a release duration of t_r = 15 s. There were four cases selected for the test simulations, as presented in Table 1, in which two are convective and two are stable cases. All the necessary data (radiosonde vertical profiling, precisely wind speed and temperature) were taken under subscription from the Chuva Project website (http://chuvaproject.cptec.inpe.br). Temperature profile, changed into potential temperature, was used in order to obtain PBL height (visual estimative) and, with the wind speed profile, was used to obtain other micrometeorological parameters necessary to run the model. It is worth noticing that local time equals to UTC time minus 3 h, meaning that the soundings were obtained during middle afternoon (Cases 1 and 4) or early night (Cases 2 and 3).

Table 1
Selected Cases: date, time, Obhukov’s length, Richardson number, friction velocity, convective velocity and velocities of gravitational settling and ground deposition, respectively.

Case	Date	Time (UTC)	PBL height (m)	L	R i	u *	w *	Vg (cm/s)	Vd (cm/s)
1	01/03/2010	18:00	410	-36.02	-0.38	0.39	1.30	0.08	0.13
2	02/03/2010	23:45	400	266.05	0.07	0.47	NA	0.08	0.16
3	04/03/2010	00:02	550	172.60	0.07	0.32	NA	0.08	0.17
4	24/03/2010	18:05	560	-108.62	-0.19	0.40	1.26	0.08	0.17

• NA: not applicable.

Figure 2 shows the vertical eddy diffusivity profiles for the selected cases. The vertical axes scale is dimensionless in order to eliminate the effects due to different PBL heights; however, the horizontal axes scale was kept unchanged, exposing the magnitudes of turbulent dispersion. As expected, the unstable cases presented turbulent diffusivities approximately 10 to 20 times greater than stable situations. Also, it may be important to point out that vertical turbulent diffusivity in unstable PBL is stronger on the upper-centre of the boundary layer, differently from the stable PBL, in which this occurs on the lower portion of the layer.

Figure 2
Vertical eddy diffusivities profile for (a) Stable cases – Case 2 (dotted line) and Case 3 (continuous line) — and for (b) Unstable cases – Case 1 (dotted line) and Case 4 (continuous line).

Figure 3 shows the wind speed profiles simulated by the model and the observed case. It is noticeable that the used wind profile formulation has a better fit on the unstable PBL and it causes a significant underestimation on the stable ones, mainly on the residual layer. It is not clear if this dichotomy is related to atmospheric stability or the magnitude of wind speed itself, since stable cases have also greater wind speeds. Simple statistics (correlation - COR, and normalized mean square error - NMSE) between simulated and observed cases leads to: COR = 0.954 and NMSE = 0.004 for Case 1 (C1); 0.907 and 0.080 for Case 2 (C2); 0.992 and 0.137 for Case 3 (C3); 0.713 and 0.033 for Case 4 (C4). The cases will be discussed with respect to the PBL stability.

Figure 3
Vertical profiles of observed (green lines) and calculated (blue lines) wind profiles for (a) Case 1, (b) Case 2, (c) Case 3, and (d) Case 4.

UNSTABLE CASES

Cases 1 and 4 are related to the (slightly) unstable PBL. The PBL winds on C1 are classified as light to moderate breeze (from about 2 to 8 m/s), while on C4 they are light to gentle breeze (about 2 to 4 m/s), although the vertical eddy diffusivity on C4 reaches values of almost 50% greater than the ones in C1. The concentration isolines for Cases 1 and 4 are displayed on Figs. 4 and 5, respectively, for 5, 10, 15 and 20 min after the release and at the height of 1.5 m above ground level. As one might notice, for the assumed source conditions (height and time-release), the same for all cases, C4 presents greater concentrations than C1: the simulations made resulted in maximum concentration of 12 mg/m³ for C1 and about 20 mg/m³ for C4, both after 300 s after engine burn, which can be better observed in Fig. 6. We may conclude that the greater values of eddy diffusivity allied to smaller wind velocities should have led to this situation. Those factors are responsible for a small advection (horizontal spread) and a high vertical transport (turbulent diffusion), increasing the effectiveness of the gradient transport: moving pollution from the buoyant cloud downwards ground level. C1, though, has a greater horizontal (downwind) reach than C4, due to major wind speeds.

Figure 4
Isolines of concentration (mg/m³) at z = 1.5 m for Case 1 simulation.

Figure 5
Isolines of concentration (mg/m³) at z = 1.5 m for Case 4 simulation.

Figure 6
Concentration versus distance from the source at y = 0 (axis of the greatest concentrations) at different times after release for (a) Case 1 and (b) Case 4.

STABLE CASES

Cases 2 and 3 are constraints of stable PBL. For these cases, wind speeds range is from 2 to 8 m/s for C2 and from 2 to 6 m/s for C3. Concerning the vertical eddy diffusivity, there is not a significant relative difference, except for the greater magnitude reached in C2, as well as for the lower relative height in which this peak occurred. The maximum concentration simulated in both cases is about 0.6 mg/m³, after 600 s for C2 and after 900 s for C3, as in Figs. 7 (C2) and 8 (C3), and in more details in Fig. 9. We may consider this time difference between the peaks as a result of the more intense turbulent diffusivity in C2. The wind speed seems also to play an important role: since this quantity is greater in Cases 2 and 3, along with a much smaller turbulent diffusivity, horizontal transport of contaminant gains crucial importance against the almost non-effective vertical transport, allowing the wind to carry pollution further before it reaches ground level. The two main consequences are the small concentration in the lower PBL and the further potential distance from the source it may reach (when compared with the unstable cases).

Figure 7
Isolines of concentration (mg/m³) at z = 1.5 m for Case 2 simulation.

Figure 8
Isolines of concentration (mg/m³) at z = 1.5 m for Case 3 simulation.

Figure 9
Concentration versus distance from the source at y = 0 (axis of the greatest concentrations) at different times after release for (a) Case 2 and (b) Case 3.

DISCUSSION

After examining the four simulated cases, there are some considerations to be made, concerning both the model itself and the dispersion on the different atmospheric conditions. The model presents some important limitations. First, the source is considered punctual (instead of a more accurate volume source like a sphere or a cilinder), which may lead to an overprediction (which is better than underprediction). Second, there is not yet a proper pre-processor for calculation of the cloud stabilization height (computed as the source height). Here we set a hypothetical source height (with no physical basis), it is, for the atmospheric conditions of the cases here presented this height could be greater or smaller, depending on PBL stability, which will influence the concentration of pollutants. As an example, extra simulations were made in order to investigate the effects of reducing the source height from 150 to 50 m in C2 (not shown). The results show that the concentration at Z = 1.5 m increases remarkably up to 21 mg/m³, 300 s after release. In this situation the source height would be at Z/h = 0.1, where the eddy diffusivity coefficient is almost reaching its greatest value, while in the displayed situation (H_s = 150 m) the relative source heightis Z/h = 0.3, and the value of K_z is already smaller. Third, the wind profile formulation (even though seems to fit well on an unstable PBL) works better on the surface layer, tends to underestimate the velocities on the higher portion of the boundary layer and is also unable to represent the low level jet (e.g.Fig. 3), which might play an important role concerning mechanical turbulence at night. Fourth, the physical characteristics of the boundary layer, mostly at night, must be improved under deeper specific research and study. In the particular situation of a stable boundary layer, the analysis is quite more complex, mainly due to geographical position of the CLA site, close from both the ocean and the 45 m height cliff, leading to the formation of an internal boundary layer (IBL). This IBL is a transition zone between the boundary layers over two different surfaces in terms of roughness or heat flux, for example, which also may have different diurnal cycles due to the physical differences of the underlying surface. Added to this fact, the nocturnal boundary layer is constituted by the stable layer and the residual layer, which smoothly blend in one another (Stull 1988). However, since the rocket exhaust cloud when stabilized may reach an altitude higher than the stable boundary layer, we considered, due to all those facts, that the PBL height during night time was the top of the residual layer (Reuter et al. 2004).

CONCLUSION

An analytic model to assess dispersion of rocket exhaust clouds was presented. The results showed physical consistency regarding the dispersion and diffusion mechanisms, indicative of a promising tool for research and management. There are some further developments that may take place on the future in order to improve the model, concerning mostly the source condition: the source height, which lacks physical support (calculation of a stabilization height of the buoyant cloud), and dimensions (it was treated here as punctual, but more accurately should be considered a sphere or a cylinder).

REFERENCES

• ADMS-CERC (2012) Modelling dry deposition - ADMS 5. Cambridge Environmental Research Consultants Ltd., Cambridge, UK [accessed 01 Aug 2014]. http://www.cerc.co.uk/environmental-software/assets/data/doc_techspec/CERC_ADMS5_P17_13.pdf
  » http://www.cerc.co.uk/environmental-software/assets/data/doc_techspec/CERC_ADMS5_P17_13.pdf
• Bainy BK, Buske D, Quadros RS (2015) Analytical model for rocket effluent dispersion: sensitivity analysis. Proceedings of the 14th International Conference on Wind Engineering; Porto Alegre, Brazil
• Bianconi R, Tamponi M (1993) A mathematical model of diffusion for a steady source of short duration in a finite mixing layer. Atmos Environ 27(5):781-792. doi: 10.1016/0960-1686(93)90196-6
  » https://doi.org/10.1016/0960-1686(93)90196-6
• Bjorklund J, Dumbauld JC, Geary H (1982) User's manual for the REEDM (Rocket Exhaust Effluent Diffusion Model) compute program, NASA contractor report 3646. Huntsville, AL: George C. Marshall Space Flight Center.
• Buske D, Vilhena MT, Tirabassi T, Bodmann B (2012) Air pollution steady-state advection-diffusion equation: the general three-dimensional solution. J Environ Protect 3(9A):1124-1134. doi: 10.4236/jep.2012.329131
  » https://doi.org/10.4236/jep.2012.329131
• Chuva Project (2015) [accessed 2015 Feb 4]. http://chuvaproject.cptec.inpe.br
  » http://chuvaproject.cptec.inpe.br
• Degrazia GA, Anfonssi D, Carvalho JC, Mangia C, Tirabassi TC (2000) Turbulence parametrisation for PBL dispersion models in all stability conditions. Atmos Environ 34(21)::3575-3583. doi: 10.1016/S1352-2310(00)00116-3
  » https://doi.org/10.1016/S1352-2310(00)00116-3
• Degrazia GA, Rizza U, Mangia C, Tirabassi T (1997) Validation of a new turbulent parameterization for dispersion models in a convective conditions. Boud Layer Meteor 85(2):243-254. doi: 10.1023/A:1000474204748
  » https://doi.org/10.1023/A:1000474204748
• Gonçalves GA, Quadros R, Buske D (2013) An analytical formulation for pollutant dispersion simulation in the atmospheric boundary layer. J Environ Protect 4(8A):57-64. doi: 10.4236/jep.2013.48A1008
  » https://doi.org/10.4236/jep.2013.48A1008
• Massman WJ (1998) A review of the molecular diffusivities of H₂O, CO, CH, CO, O, SO, NH, NO, NO and NO in air, Oand N near STP. Atmos Environ 32(6):1111-1127. doi: 10.1016/S13522310(97)00391-9
  » https://doi.org/10.1016/S13522310(97)00391-9
• Moreira DM, Trindade LB, Fisch G, Moraes MR, Dorado RM, Guedes RL (2011) A multilayer model to simulate rocket exhaust clouds. J Aerosp Technol Manag 3(1):41-52. doi: 10.5028/jatm.2011.03010311
  » https://doi.org/10.5028/jatm.2011.03010311
• Moreira DM, Vilhena MT, Buske D, Tirabassi T (2009) The state-of-art of the GILTT method to simulate pollutant dispersion in the atmosphere. Atmos Res 92(1):1-17. doi: 10.1016/j.atmosres.2008.07.004
  » https://doi.org/10.1016/j.atmosres.2008.07.004
• Nyman R (2009) NASA Report: evaluation of Taurus II Static Test firing and normal launch rocket plume emissions. NASA. Final Report: environmental assessment for the expansion of the Wallops flight facility launch range. Wallops Island, VA: National Aeronautics and Space Administration [accessed 01 Aug 2014]. http://sites.wff.nasa.gov/code250/docs/expansion_ea/Appendix_G_MARS_Final_EA.pdf
  » http://sites.wff.nasa.gov/code250/docs/expansion_ea/Appendix_G_MARS_Final_EA.pdf
• Panofsky HA, Dutton JA (1984) Atmospheric turbulence. New York: John Wiley & Sons.
• Reuter ED, Fisch G, Mota GV, Cohen JC (2004) Estudo observacional da camada limite planetária marinha na região do Centro de Lançamento de Foguetes de Alcântara-MA. Rev Bras Meteorol 19(3):251-264.
• Roballo ST, Fisch G (2008) Escoamento atmosférico no Centro de Lançamento de Alcântara, CLA, Parte 1 - aspectos observacionais. Rev Bras Meteorol 23(4): doi: 10.1590/S010277862008000400010
  » https://doi.org/10.1590/S010277862008000400010
• Seinfeld JH, Pandis SN (2006) Atmospheric chemistry and physics. New Jersey: John Wiley & Sons.
• Stroud A, Secrest D (1966) Gaussian quadrature formulas. Englewood Cliffs, NJ: Prentice Hall Inc.
• Stull RB (1988) An introduction to boundary layer meteorology. Dordrecht: Kluwer Academic Publishers.
• Zhang L, Gong S, Padro J, Barrie L (2001) A size-segregated particle dry deposition scheme for an atmospheric aerosol module. Atmos Environ 35(3):549-560. doi: 10.1016/S1352-2310(00)00326-5
  » https://doi.org/10.1016/S1352-2310(00)00326-5

Publication Dates

• Publication in this collection
  Jul-Sep 2015

History

• Received
  07 July 2015
• Accepted
  21 Aug 2015