## Abstract

The mean flow characteristics of a turbulent round jet and a plume are simulated in a high-resolution regional atmospheric model. The model used for the study is the Advanced Regional Prediction System (ARPS). It is observed that the predicted flow characteristics are very sensitive to the turbulence closure scheme used. Among the three closure options in ARPS, namely, the constant eddy viscosity scheme, the Smagorinsky diagnostic closure, and the 1.5-order total kinetic energy scheme, the constant eddy viscosity scheme predicts the jet characteristics in better agreement with experiments. For the plume, all three schemes are unsatisfactory. It is shown that a modification of the constant eddy viscosity scheme incorporating a length-scale variation as suggested by theory predicts plume characteristics in good agreement with experiments. The simulations are carried out at one fixed grid-box size and flow inlet conditions; extending the present simulation results to other cases is discussed.

## 1. Introduction

Round jets and plumes are flows resulting from a source of momentum and buoyancy, respectively, in a large body of constant density fluid. These flows are generally turbulent, expand conically in the axial direction, and appear to originate from a point source (Fig. 1). Jets and plumes are geometrically simple flows, amenable to experimental investigation and theoretical analysis. They provide the basis for modeling a variety of practical and natural flows, including combustion, waste disposal, cooling towers, and cumulus clouds (Fisher et al. 1979; Chen and Rodi 1980; Arakawa and Schubert 1974). For these reasons, jets and plumes have been investigated extensively (Turner 1973; Hinze 1975; Fisher et al. 1979; Chen and Rodi 1980; Malin and Spalding 1984; Cho and Chung 1992).

Jets and plumes are deceptively simple flows and their numerical simulation is a challenging task (e.g., Cho and Chung 1992). The majority of the studies are based on the *K*–ɛ turbulence closure model or its variants (e.g., Chen and Chen 1979; Cho and Chung 1992). The standard *K–*ɛ model has not been found to be very successful in predicting jet and plume characteristics (Cho and Chung 1992). Modifications to the *K*–ɛ model to account for (a) the effect of vortex stretching (Pope 1978), (b) the effect of normal stress (Hanjalic and Launder 1980), and (c) the effect of entrainment (*K*–ɛ–*γ* model;Cho and Chung 1992) have enabled the characteristics of the round jet to be accurately predicted. The *K*–ɛ–*γ* model (Cho and Chung 1992) is found to be superior to other two in the case of a jet, but all are unsatisfactory for a plume. The *K*–ɛ–*t*′^{2}–*γ* model, which accounts for the entrainment and the buoyancy transport by the fluctuating quantities, predicts the plume characteristics in good agreement with experiments (Dewan et al. 1997). Thus, accurate prediction of both jet and plume characteristics by means of numerical methods is now possible by employing closure schemes specific to these flows.

So far, the emphasis of the modeling efforts has been on predicting jet and plume flows under laboratory-type conditions. A couple of issues are to be addressed before such models can be extended to atmospheric flows. In the atmosphere, the flow conditions are much more complex compared to that assumed in the successful simulations so far. For example, the latent heat release, due to the condensation of water vapor, and the stable stratification of the surrounding environment are very important in the case of cumulus clouds, and a realistic model for a cumulus cloud needs to take these factors into account. The numerical models that are successful in jet and plume simulations assume a constant density ambient fluid and there is no phase change of water. Therefore, modifying the successful jet and plume models mentioned earlier for application to atmospheric problems is not a trivial task. Second, the turbulence closure schemes in these models are involved and computationally expensive, and it may not be practical to include them in atmospheric models.

Now, limited area (regional) atmospheric models are available that directly solve the Reynolds-averaged equations governing the motion of the air in the atmosphere. These models contain important physical processes relevant to the atmosphere, and a variety of problems can be attempted. Therefore, simulation of jet and plume flows in a limited area model seems feasible and attractive. The advantage of this is that we can test the effect of phase change of water, ambient stratification, and ambient wind field (shear) on jet and plume behavior. The authors are not aware of any previous work on the simulation of round jets and plumes in a regional atmospheric model. Therefore, our first objective was to see if the basic characteristics of these flows can be simulated in a regional atmospheric model without adding new equations or making changes to the model equations. If the effort is not successful, we will then explore minimum modifications that give satisfactory results.

Our study shows that it is possible to predict the mean characteristics of both turbulent jets and plumes using a limited area atmospheric model. The materials in this paper are organized as follows. Section 2 contains a brief description of the time mean characteristics of turbulent jets and plumes. The details of the regional atmospheric model relevant to the present study are given in section 3 followed by the results in section 4. Section 5 analyzes why the model results are very sensitive to the choice of the turbulence closure scheme, section 6 contains a brief discussion, and section 7 contains the conclusions.

## 2. Round jet and plume: Theory

Here we consider the steady-state behavior of a vertically pointed turbulent round jet and plume, and follow the treatments given in Turner (1973), Fisher et al. (1979), and Shabbir and George (1994). The initial volume flux is zero in the case of an ideal jet and plume;however, the jets and plumes produced in practice have finite source size and initial volume flux. Observations show that, after some distance from the source (typically 8–10 source diameters, called the zone of flow establishment), the flow is turbulent, expands conically in the horizontal, and appears as if originating from a point source (Fig. 1). The instantaneous flow field is three-dimensional, highly complex, and varies from one instant to another (e.g., Papantoniou and List 1989; Venkatakrishnan et al. 1999). However, when averaged over a period that is much larger than the timescale given by *l*/*U* (the eddy turnover time), where *l* and *U* are characteristic length and velocity scales in the flow, the resulting mean flow field is steady and two-dimensional (e.g., Papantoniou and List 1989; Bhat and Narasimha 1996). The mean flow is symmetric about a vertical axis passing through the center of the orifice, referred to as the axis and also as the centerline. In the following, the distance along the axis is denoted by *z* and that in the horizontal by *r.* Beyond the zone of flow establishment, flow has the tendency to become self-similar; that is, when normalized by the appropriate velocity and length scales in the flow, the shape of mean horizontal profiles become independent of the axial location (Hinze 1975). Normally, the time mean velocity at the centerline (*W*_{c}) is taken as the velocity scale, and a length characterizing the jet–plume width (denoted by *b* here) is taken as the length scale. Dimensional arguments together with experimental observations suggest the following forms for the mean flow variables, which are known as similarity solutions (Fisher et al. 1979; Hussein et al. 1994; Shabbir and George 1994): jet,

plume,

for both jet and plume,

In the above equations, *d* and *W*_{o} are, respectively, the nozzle diameter and the vertical velocity at the inlet (Fig. 1); *W* is the vertical velocity; *b* is the radius at which *W*/*W*_{c} = 1/*e*; *A* and *β* are constants; *B* is the buoyancy flux; *T* is temperature; *g* is acceleration due to gravity; and *z*_{o} is the virtual origin, that is, the distance above/below the orifice where the flow appears to originate (Fig. 1). Here, *B* = (*π*/4)*gα**d*^{2}*W*_{o}*θ*_{o}, where *α* is the coefficient of thermal expansion and *θ*_{o} is the perturbation temperature at the inlet. The subscripts *j* and *p* refer to the jet and plume, respectively. Various experimental studies agree on the form of the equations above; however, the values of the constants are not unique (Fisher et al. 1979). Experimentally measured spread rate for the jet (*β*_{j}) varies in the range 0.1–0.13, and that of the constant *A*_{j} between 5 and 7 (Fisher et al. 1979; Hussein et al. 1994). The spread rate of the plume (*β*_{p}) is comparable to that for the jet (Turner 1986). The constant *A*_{p} varies from 3.9 to 4.7 (Fisher et al. 1979; Shabbir and George 1994; Dai et al. 1995).

## 3. Numerical model

Our aim in these numerical simulations is to see if the similarity solutions of round turbulent jets and plumes can be produced in a limited area atmospheric model. Jets and plumes are axisymmetric flows. Though a cylindrical coordinate system is the natural coordinate system for them, we find successful jet and plume simulations carried out in a rectangular box using the equations in a Cartesian system (Cho and Chung 1992). Therefore, the choice of the model coordinate system is not a crucial issue as long as appropriate boundary conditions are applied.

The numerical code selected here is the Advanced Regional Prediction System (ARPS) developed by the Center for Analysis and Prediction of Storms, Norman, Oklahoma. ARPS is a nonhydrostatic atmospheric prediction model based on time-averaged (but unsteady) compressible Navier–Stokes equations. Prognostic equations for *u, υ, w, θ*′, *p*′, and *q*_{ψ} are solved, which are respectively, the *x, y,* and *z* (vertical) components of the Cartesian velocity, the perturbation potential temperature, perturbation pressure, and six categories of water substance. The state variables (e.g., potential temperature and pressure) are defined as the sums of base-state variables and the deviation from the base state. The base state is assumed to be horizontally homogeneous, time invariant, and hydrostatically balanced. The governing equations are represented in a curvilinear coordinate system, which is orthogonal in the horizontal. Through coordinate transformation, these equations are solved in rectangular computational space. A comprehensive description of the model can be found in Xue et al. (1995).

In the present simulations, the bottom surface is flat and hence the coordinate system reduces to the standard Cartesian system. There is no Coriolis force or water vapor. The equations actually solved using the ARPS code are the following:

where Div* = ∂*u**/∂*x* + ∂*υ**/∂*y* + ∂*w**/∂*z, **ρ* is the mean density, *u** = *ρ**u, υ** = *ρ**υ, w** = *ρ**w, α* is the divergence damping coefficient introduced to attenuate the acoustic waves, *D* is the subgrid-scale turbulence and is expressed in terms of the Reynolds stress tensor *τ*_{ij} for the velocity as *D*_{i} = ∂*τ*_{ij}/∂*x*_{j}, *D*_{θ} = ∂*H*_{i}/∂*x*_{i} (*H* being the heat flux), and *S*_{θ} is the heat source. In addition, *B* is the buoyancy force and *c*_{s} is the acoustic wave speed.

The nonhydrostatic formulation used in ARPS supports acoustic waves, which are meteorologically unimportant but limit the time step size whenever an explicit time integration scheme is used (as is the case here). ARPS employs a mode-splitting time integration technique presented by Klemp and Wilhelmson (1978). This technique divides the big integration time step into a number of computationally small time steps to update the acoustically active terms while all other terms are updated every big time step. The big time step integration uses the leapfrog time differencing scheme.

Subgrid-scale turbulence is important for the flows being simulated here. In the ARPS model, turbulence stress *τ*_{ij} is expressed as

where *k*_{m} is the turbulent mixing coefficient and *D*_{ij} is the velocity deformation tensor (∂*u*_{i}/∂*x*_{j}). The closure schemes differ depending on how *k*_{m} is calculated. ARPS version 4.2.3a, on which the present study is based, has options for three turbulence closure schemes:(a) constant eddy viscosity, (b) Smagorinsky diagnostic closure, and (c) 1.5-order total kinetic energy (TKE) scheme. Descriptions below are very brief, relevant to an isentropic base state, and more details can be found in Xue et al. (1995).

Constant eddy viscosity model (CEV): As the name suggests, the mixing coefficient is constant; that is,

where *k*_{c} is a constant.

Smagorinsky scheme (SK): The mixing coefficient is given by

where *k*_{s} is an empirical constant assigned a value of 0.21 in the model, Δ is a measure of the grid size, and |Def|^{2} is related to the square of deformation tensor *D*_{ij} and is given by

1.5-order TKE closure (TKE): In this closure scheme, the mixing coefficient is given by

where *k*_{t} is a constant set equal to 0.1 in ARPS, *E* is the turbulent kinetic energy, and *l* is the length scale proportional to the grid spacing in the version of ARPS used for this study. When TKE is used for turbulence closure, an additional prognostic equation for the turbulent kinetic energy is solved.

### a. Boundary conditions

The numerical domain is a rectangular box having *nx, ny,* and *nz* number of grid points along the *x, y,* and *z* (vertical) directions, respectively. The base state is an isentropic atmosphere at rest. The flow enters from the bottom through a small circular area (nozzle) having diameter *d* and centered at (*nx*/2, *ny*/2). In the present simulations, *nx* = *ny* = 40 and *nz* = 81. Both *nx* and *ny* are even, so that the jet–plume axis is a “*w*” point in order to maintain the symmetry about the flow axis (the dependent variables are discretized on the Arakawa C grid). The boundary conditions are applied as follows (Fig. 2).

Bottom (Fig. 2a): for *r* ⩽ *d*/2,

and for *r* > *d*/2, free-slip rigid wall, that is (Xue et al. 1995),

The velocity and temperature are uniform within the nozzle area, a distribution known as the tophat profile (e.g., Turner 1973, p. 172). Since the computational grid boxes are rectangular and the nozzle is circular, the velocity and temperature perturbations for the grid boxes partially intercepting the nozzle area are obtained by multiplying *W*_{o} and *θ*_{o} by the fraction of the grid area that the nozzle covers. Also, in order to avoid any sudden impulses and the associated numerical problems, *W*_{o} and *θ*_{o} are increased from zero to their maximum value in time *τ* (Fig. 2b). The value of *τ* is chosen to be much smaller than the total time of integration so that the effects of the initial transients do not influence in the final flow characteristics.

Top: zero normal gradient boundary condition, that is (Xue et al. 1995),

Lateral sides: we assume that the radial distance of the lateral side *r*_{l} is large compared to the jet and plume width (results presented in the next section support this). From Eq. (6), it is seen that *W* rapidly approaches zero as *r*/*b* increases. The horizontal velocities are also small;however, they cannot be made identically equal to zero as it amounts to zero lateral entrainment. Theoretical arguments show that at large values of *r*/*b,* horizontal velocity is directed toward the centerline and inversely proportional to *r* (Schneider 1981). Therefore, the following boundary condition is applied at the lateral boundaries (Fig. 2c):

The *K* in Eq. (14) is the laterally entrained mass at each vertical level. It is calculated from the net mass entering a circle located inside but just touching the last but one outermost grid boxes in the lateral direction.

## 4. Results

The grid sizes Δ*x,* Δ*y,* and Δ*z,* along the *x, y,* and *z* directions, respectively, are 1 m each, and *d* = 2 m. Compared to the laboratory scale, the grid-box size is rather large; however, for an atmospheric model, this can be regarded as ultrahigh resolution. Extending the results to other grid-box sizes is straightforward and is discussed in section 6. The big time step is 0.04 s. The acoustic wave speed is taken to be a constant at 180 m s^{−1} (results are found to not be sensitive to the exact value chosen as long as it is very large compared to the flow velocities) and the small time step of integration is set at 0.002 s. In jet simulations, *W*_{o} = 5 m s^{−1}, *θ*_{o} = 0, and *τ* = 10 s. In plume simulations, *W*_{o} = 3.3 m s^{−1}, *θ*_{o} = 10 K, and *τ* = 10 s. The equations were integrated for 360 s starting from *t* = 0. The solution at *t* = 360 s indicated that the flow had attained a steady behavior in the sense that the difference between the fields at the end of 360 s and their average calculated for the past 30 s of integration are almost indistinguishable. For time *t* ≫ *τ,* the forcing terms (i.e., momentum source in the case of the jet, and in addition, the buoyancy force in the case of plume) become constant. Since the model equations are averaged Navier–Stokes equations, in the absence of a time-dependent forcing, we expect the solution to attain a steady behavior. The numerical results confirm this.

### a. Jet

Figure 3a shows the variation of jet width *b* with the axial distance *z,* both normalized by the nozzle diameter *d.* Predictions from Eq. (1) are shown for comparison (open circle), with the error bars indicating the experimental spread (e.g., Fisher et al. 1979). As mentioned in section 2, *z*/*d* < 10 is the zone of flow establishment and the jet behavior is not self-similar. Therefore, it is relevant to consider the region *z*/*d* > 10 in Fig. 3a to see if the model produces the self-similar jet. The Smagorinsky scheme predicts a linear growth of *b* with *z,* but underpredicts the spread rate by an order of magnitude. The prediction from the TKE closure scheme is in good agreement with the experiments for *z*/*d* up to 15. (The value of *k*_{t} is set equal to 0.25 in the present study.) It is seen from Fig. 3a that the TKE scheme fails to maintain the spread rate as *z* increases beyond 15d. The constant eddy viscosity model (*k*_{c} = 0.12) predicts the axial variation of jet width well. Further computations have shown that the experimentally observed range of the jet spread rate can be reproduced in the model by varying *k*_{c} around 0.12 (the dependence of the spread rate on *k*_{c} is given in the next section).

Figure 3b shows the axial variation of *W*^{−1}_{c}. The predictions using the three turbulence closure schemes are shown by the solid lines. The line with open circles is Eq. (2) with *A*_{j} = 5.4 and *z*_{0}/*d* = −2.5. For reasons already stated, we will consider the model behavior for *z*/*d* > 10. The Smagorinsky scheme overpredicts *W*_{c} (i.e., lower values of *W*_{o}/*W*_{c}) throughout. The TKE scheme does well for *z*/*d* < 15 but not beyond. The constant eddy viscosity scheme performs better than the other two schemes, and the predicted *W*_{c} is in good agreement with the experimental data. Figure 3c shows the horizontal profile of the vertical velocity along the radial direction at *z*/*d* = 30. The continuous line is the Gaussian distribution [Eq. (6)]. In the central parts of the jet, the predicted horizontal profiles are nearly Gaussian; however, the velocity decays at a slower rate toward the edges. The slower velocity decay toward the edges in numerical simulations of jet and plume flows has been a problem even in models that employ complex turbulence closure schemes and that otherwise had performed well (e.g., Dewan et al. 1997). To the knowledge of the authors, only the *k*–ɛ–*t*′^{2}–*γ* model predicts the shape of the horizontal profiles of jets and plumes in good agreement with the experiments (Dewan et al. 1997).

### b. Plume

The axial variation of predicted plume width is shown in Fig. 4a along with corresponding experimentally observed range. For the Smagorinsky closure scheme, the rate of growth of *b* with *z* (*db*/*dz* ∼ 0.01) is too small compared to the experimental value (∼0.12). The predictions using the TKE scheme are (relatively) better for *z* < 20 *d,* but depart from observations for larger *z.* The CEV scheme predicts a nonlinear variation of *b* with *z.* Therefore, all three closure schemes are unsatisfactory. Several model runs were carried out to explore if the departure of the predictions from the standard plume behavior is due to deficiencies in the specification of the boundary conditions, insufficient number of grids, and/or the presence of acoustic waves. The differences among the solutions were minor in the sense that none of these led to the reproduction of the correct basic trend. Since the model results are very sensitive to the turbulence closure scheme employed, modifying the closure scheme seemed reasonable. Based on a theoretical analysis (details given in the next section), the mixing coefficient was prescribed as

subjected to the condition that *k*_{m} ≥ *k*_{min}. Parameter *k*_{min} was chosen to be comparable to that used for the jet simulation. Here the turbulence closure scheme with *k*_{m} given by Eq. (15) is called the modified eddy viscosity scheme (MEV) for the convenience of reference. In the MEV case shown in Fig. 4a, *k*_{1} = 0.03, *k*_{2} = 0.01, and *k*_{min} = 0.12. The predicted plume width variation is in good agreement with the experimental data.

It is seen from Eq. (4) that the quantity *B*/*W*^{3}_{c} (denoted by *ϕ* here) varies linearly with *z.* As it is easy to judge the departure from a linear trend, *ϕ* is plotted against *z* in Fig. 4b. Since the buoyancy *B* is kept constant in the simulations, the axial variation of *ϕ* basically reflects the axial variation of *W*_{c}. The experimentally observed variation of *ϕ* is shown by the open circles. The prediction using the Smagorinsky scheme is very poor. Both the TKE and constant eddy viscosity schemes predict *ϕ* in reasonable agreement with experiments for *z* < 15*d* but not beyond. The prediction using the MEV scheme are in good agreement with experimental data except near the top boundary where the model tends to slightly overpredict. It may be possible to improve the results further by fine-tuning the constants *k*_{1} and *k*_{2}.

## 5. Theoretical analysis

Simple geometry and the tendency of the flow in jets and plumes toward a similarity behavior enables us to simplify the Navier–Stokes equations governing these flows to obtain ordinary differential equations describing the axial variation of width *b* and the centerline velocity *W*_{c}. While *W*_{c} dependence on *z* is different for the jets and plumes, that of *b* is linear for both with comparable spread rates (section 2). Therefore, the equation for the spread rate is considered in the following. The details of the derivation are given in the appendix where it is shown that the spread rate of these flows is approximately given by

In the above equation, *τ*_{12} is the turbulent shear stress, modeled in the ARPS code using the eddy viscosity formulation; that is,

The velocity gradient in Eq. (17) can be expressed making use of the axisymmetric nature of jet and plume flows. Along the centerline *U* = 0, therefore, its derivative with respect to *z* vanishes. Since *W* is an even function of radial distance *r,* it can be expressed in the neighborhood of *r* = 0 as

neglecting higher-order terms. For a Gaussian profile, *a*_{1} is unity and this value is assigned to *a*_{1} henceforth. Now let us consider how the three closure schemes of ARPS model the turbulent stress term in Eq. (17) and predict *db*/*dz.*

### a. The constant eddy viscosity model

For this closure scheme, *k*_{m} = *k*_{c} = constant, and *τ*_{12} = −2*k*_{c}*W*_{c}*r*/*b*^{2}. Substituting this into Eq. (16) we get

In the case of the jet, product *bW*_{c} is independent of *z* [Eq. (A7)]. Since *k*_{c} is a constant, the right-hand side of Eq. (20) is constant and by tuning *k*_{c} any required spread rate for the jet can be obtained. For the plume, the product *bW*_{c} keeps increasing with *z* because of the action of the buoyancy force [Eq. (A7)] and therefore *db*/*dz* decreases with *z.* Thus, we expect the spread rate to decrease with *z,* which is in agreement with the model results shown in Fig. 4a.

### b. Smagorinsky scheme

For this closure scheme, *k*_{m} = *k*^{2}_{s}Δ^{2}|∂*W*/∂*r*|,

and we get

Therefore, the main term responsible for maintaining the spread rate becomes zero in the Smagorinsky scheme. It may be noted here that the model-predicted spread rate is small (∼0.01) but not identically zero (Figs. 3a and 4a). This results from the terms that have been neglected in deriving Eq. (16).

### c. TKE closure scheme

where *l* is the length scale taken to be proportional to the grid size Δ and *E* is the turbulent kinetic energy [=(*u*^{2} + *υ*^{2} + *w*^{2})/2]. In the similarity region, we expect *E* to be proportional to *W*^{2}_{c}. Hence, *db*/*dz* can be expressed as

where the constant *c*_{1} is introduced to account for the proportionality constants relating *l* and *E*^{1/2} to Δ and *W*_{c}, respectively. In Eq. (23), all except *b* are independent of *z*; since *b* increases with *z,* the spread rate decreases.

Thus, the conclusions from Eq. (16) are consistent with model results and provide an insight into why the results are so sensitive to the turbulence closure schemes used in the ARPS model.

Equation (16) was used for improving the plume simulations shown in the previous section. Our main objective was to look for a minimum of modification that will improve the model predictions. Since the constant eddy viscosity scheme works best for the jet, here its modification for predicting the plume flow is considered. The product *bW*_{c} in Eq. (20) increases as *z*^{2/3} for the plume. Therefore, if *k*_{c} is made proportional to *z*^{2/3} rather than being kept fixed, then the model predictions could improve. In MEV, the mixing coefficient given by Eq. (15) has been used, and for this scheme, the equation for the spread rate becomes

Here *k*_{1} is chosen to be much smaller than the expected plume spread rate *β*_{p}. As *z* increases, *db*/*dz* becomes independent of *z.* As already shown in the previous section, predictions made with MEV are in good agreement with experimental observations.

## 6. Discussion

Jets and plumes are flows dominated by turbulence, and turbulence controls the development of the mean flow characteristics. Therefore, in the numerical models for their simulation, proper parameterization of the subgrid-scale turbulence is very critical. In the ARPS model, the simplest closure scheme, namely, the constant eddy viscosity model, is observed to be superior to the Smagorinsky and TKE schemes in predicting jet and plume mean flow characteristics. The Smagorinsky scheme is not suitable for jet and plume simulations. In the case of jets, the predications using the constant eddy viscosity scheme are in good agreement with experimentally observed spread rate and the decay of the axial component of velocity (in the axial direction). Theoretical analysis showed that this is more of a coincidence (than due to any inherent merit of this scheme over the TKE scheme) owing to the special nature of the jet flow where, *bW*_{c}, the product of length and velocity scales, remains constant. When *bW*_{c} does not remain constant, for example, as in the case of a plume, its performance is not satisfactory. It may be noted here that many highly complex closure schemes are also not fully satisfactory in predicting plume flows (Dewan et al. 1997). Therefore, with the simplicity that it offers, the CEV (MEV) scheme is a good compromise for jet (plume) simulations in atmospheric models. Further, over a narrow range of *z*/*d,* the constant *k*_{c} in the CEV closure scheme can be tuned to obtain the desired spread rate of both jets and plumes [using Eq. (20)].

In the present simulations, grid-box size and the jet–plume inlet conditions were held fixed. The values assigned to the constants such as *k*_{c}, *k*_{1}, and *k*_{2} are perhaps optimum for the specific cases of jets and plumes. It is desirable to generalize these results to other grid-box resolutions and inlet conditions. The equations given in section 5 are useful for this purpose. Here we consider extending the constant eddy viscosity closure scheme for the jet and the MEV scheme for the plume as these schemes worked best for jets and plumes, respectively. From Eqs. (20) and (25), we see that the spread rate does not explicitly depend on the grid resolution. The product *bW*_{c} depends on the inlet conditions. Therefore, the values of the constants [*k*_{c} in Eq. (20) and *k*_{1} and *k*_{2} in Eq. (25)] need to be changed depending on the inlet conditions, namely, the nozzle diameter, velocity, and buoyancy.

For the jet, the product *bW*_{c} is related to the inlet conditions through the conservation of momentum; that is,

The parameter *α*_{j} is introduced to account for the shape of the horizontal profile of *W.* For a Gaussian profile, *α*_{j} = 2. We have seen in section 4 that the numerically simulated velocity profile decays at a slower rate compared to the Gaussian distribution, and the present simulations suggest a value of *α*_{j} = 1.82. Substituting (26) into (20) we get

Similarly, Eq. (24) can be used for guessing the value of *k*_{2}. Combining Eqs. (3), (4), and (25) we get (for larger *z*)

where *α*_{p} is introduced to account for the departure of the numerically simulated horizontal profiles from the Gaussian distribution. Our results suggest a value of 2 for *α*_{p}. The value of *k*_{1} is about three times *k*_{2}. Equations (27) and (28) can be used for obtaining the first-guess values of *k*_{c} and *k*_{2}, respectively. Some tuning around these values may be necessary to get a required spread rate.

## 7. Conclusions

It is shown that turbulent jet and plume mean flow characteristics, namely, the spread rate and velocity decay rates, can be simulated in the ARPS model. The results are found to be very sensitive to the turbulence closure scheme employed. The jet flow characteristics could be simulated using the constant eddy viscosity scheme without requiring any change in the model equations. Among the three closure options in the model, the Smagorinsky closure scheme is not suited for simulating the jet and plume flows. For simulating the plume characteristics, the closure schemes in ARPS require modifications. The theoretical analysis carried out here provided insight into why model results are sensitive to the choice of the closure scheme, and this knowledge is used for improving the plume predictions. Extending the present results to other inlet conditions and grid sizes is suggested.

## Acknowledgments

We thank Dr. Kelvin K. Droegemeier and Dr. Ming Xue, Center for Analysis and Prediction of Storms, University of Oklahoma, for providing the ARPS code and for their quick response to many questions regarding the use of the code. We thank Prof. Louiz M. Lourenco for the computational facilities and Prof. T. N. Krishnamurti for many useful discussions. The first author thanks the FSU Foundation and Don Fuqua Chairfund of The Florida State University for supporting his one-year visit to The Florida State University. We also thank two referees for their comments and suggestions for improvement.

## REFERENCES

### APPENDIX

#### Spread Rate of an Axisymmetric Jet and Plume

Here we derive an expression for the jet and plume spread rates from the approximate forms of the equations governing the flow. The model solutions had reached a steady-state behavior (i.e., time variations of velocity, pressure, etc., are negligible) at the time of the results shown here. Round jets and plumes are axisymmetric flows. Therefore, steady, time-averaged Navier–Stokes equations (Reynolds equations) in cylindrical coordinates have been considered here. Using the incompressibility condition, the momentum equations can be expressed at large Reynolds numbers as (Hussein et al. 1994)

In the above equations, the mean and the fluctuating quantities are denoted by the uppercase and lowercase letters, respectively. The velocity components along the axial, radial, and azimuthal directions are *W, U,* and *V,* respectively. Due to the axisymmtric nature of the flow, *V* = 0 but not the root-mean-square of its fluctuations. Other symbols have the same meaning as in section 3. Equations (A1) and (A2) can be combined to eliminate the pressure term and the resulting equation can be integrated across the entire jet or plume to yield an integral equation that can be simplified further. Detailed steps are given in Hussein et al. (1994). Based on the results of Hussein et al. (1994), Eqs. (A1) and (A2) are combined to yield

The magnitudes of *w*^{2}, *u*^{2}, and *υ*^{2} are comparable and *w*^{2}/*W*^{2}, etc., are about 5%–7% in jets and plumes (Wygnanski and Fiedler 1969; Papanicolaou and List 1988), and therefore their net contribution in Eq. (A3) can be ignored and we get

In the similarity region of the flow (typically for *z*/*d* > 20), the normalized profiles of velocity and temperature are functions of *η* = *r*/*b* only, and we can express *W* and *θ* as

where Φ_{1} = ∫^{∞}_{0} *f*^{2}_{1}*η **dη* and Φ_{2} = ∫^{∞}_{0} *f*_{2}*η **dη. *Equation (A7) is the statement of the conservation of momentum along the axial direction for the flow. For the jet, *θ*_{c} = 0, and therefore the product *bW*_{c} is independent of *z,* whereas for the plume, it increases with *z.* The left-hand side contains the derivatives of *b* and *W*_{c} and it is convenient to eliminate one of them. For this, we consider the vertical momentum equation along the centerline. Expanding Eq. (A1),

The term inside the parentheses in (A8) is identically zero for an incompressible fluid. The contributions of *w*^{2} (which is typically 6%–7% of *W*^{2}_{c}) and pressure gradient [which is not larger than that of *w*^{2}; see Hussein et al. (1994)] terms can be ignored. Then, along the centerline, (A8) becomes

Now, −*uw* is the Reynolds shear stress and is denoted by *τ*_{12} for the convenience of reference. Equations (A7) and (A9) are combined to obtain an expression for *db*/*dz.* The result is

where *σ* = Φ_{2}/Φ_{1}. It is to be noted from the above equation that the spread rate depends on buoyancy (first term on the right-hand side) and the gradient of the turbulent shear stress. For the jet, *θ*_{c} = 0, and the first term on the right-hand side of Eq. (A10) vanishes. In the case of the plume, from Eqs. (3)–(5), and experimentally observed range of the constants that appear therein, it can be shown that the first term is independent of *z* and is less than 0.01, that is, not more than 10% of the observed spread rate of around 0.1. Therefore, the spread rate is mainly controlled by the second term, and for a plume as well as a jet, we can express the rate of change of *b* with *z* by Eq. (16).

## Footnotes

* Permanent affiliation: Centre for Atmospheric and Oceanic Sciences, Indian Institute of Science, Bangalore, India.

*Corresponding author address:* Dr. A. Krothapalli, Dept. of Mechanical Engineering, The Florida State University, 2525 Pottsdamer Street, Tallahassee, FL 32310.

Email: kroth@fmrl.fsu.edu