Large-eddy simulation of plume dispersion under various thermally stratified boundary layers

Abstract. Contaminant gas dispersion in atmospheric boundary layer is of great concern to public health. For the accurate prediction of the dispersion problem, the present study numerically investigates the behavior of plume dispersion by taking into account the atmospheric stability which is classified into three types; neutral, stable, and convective boundary layers. We first proposed an efficient method to generate spatially-developing, thermally-stratified boundary layers and examined the usefulness of our approach by comparing to wind tunnel experimental data for various thermal boundary layers. The spreads of plume in the spanwise direction are quantitatively underestimated especially at large downwind distances from the point source, owing to the underestimation of turbulence intensities for the spanwise component; however, the dependence of the spanwise spreads to atmospheric stability is well represented in a qualitative sense. It was shown that the large-eddy simulation (LES) model provides physically reasonable results.


Introduction
Contaminant gas dispersion resulting from accidental release from industrial areas or intentional release of CBRN (chemical, biological, radiological, or nuclear) agent is of great concern to public health and social security.Plume dispersion within the atmospheric boundary layer is influenced by roughness elements, terrain, and thermal stability.In terms of thermal stability, atmospheric boundary layers are in general classified into three types; neutral boundary layer (NBL), stable boundary layer (SBL), and convective boundary layer (CBL).In an NBL turbulence is generated and maintained by wind shear, while in an SBL turbulence is not only maintained by wind shear but also constrained by negative buoyancy.In a CBL, turbulence is mainly produced by shear and/or buoyancy.
For simulating plume dispersion under various thermal conditions, there are typically two approaches: one is a wind tunnel experimental technique, and the other is a computational fluid dynamics (CFD) approach.It is well known that wind tunnel experiments are a reliable tool.There have been many studies of plume dispersion in various thermally-stratified boundary layers over the past 30 years (Fackrell and Robins, 1982;Deardorff and Willis, 1984;Fedorovich and Thäter, 2002).Ohya et al. (1997) and Ohya and Uchida (2004) also investigated turbulence structures in thermally-stratified boundary layers.On the other hand, with the rapid development of computational technology, the CFD technique also has come to be regarded as a useful tool.In particular, there are two different approaches; the Reynolds-Averaged Navier-Stokes (RANS) and LES models.Hanna et al. (2004), Demael and Carissimo (2007), and Vervecken et al. (2013) carried out RANS simulations of atmospheric dispersion and compared to observed data of mean concentrations.Sykes and Henn (1992) and Henn and Sykes (1992) performed LESs of plume dispersion in neutral or convective boundary layers and investigated the instantaneous structure of a plume and the characteristics of concentration fluctuations.
As mentioned-above, there are a relatively small number of LES studies on plume dispersion in boundary layers with a wide range of thermal stability.Applicability of LES technique to plume dispersion under different thermal stability has not been fully demonstrated.In this study, we first propose an efficient method to generate various thermally-stratified boundary layers such as neutral, stable, and unstable boundary layers, and then examine the usefulness of our approach in comparison to wind tunnel experimental data.

Numerical model
The model used here is LOcal-scale High-resolution atmospheric DIspersion Model using LES (LOHDIM-LES) developed by Japan Atomic Energy Agency (Nakayama et al., 2011(Nakayama et al., , 2012(Nakayama et al., , 2013)).The basic equations are the filtered continuity equation and the Navier-Stokes equation as follows; where u i , t, p, ρ, ν, τ ij , g, β, θ, θ 0 , δ ij , and are the wind velocity, time, pressure, density, kinematic viscosity, subgrid-scale (SGS) Reynolds stress, gravitational acceleration, thermal expansion coefficient, temperature, reference temperature, Kronecker delta, and grid-filter width, respectively.The subscript i stands for coordinates.ν SGS and C s are the eddy viscosity coefficient and the model constant of the flow field, respectively.Upper bars () denote application of the spatial filter.The subgrid-scale turbulent effect is represented by the standard Smagorinsky (1963) model with the constant value of 0.1.It should be noted that the application of the static Smagorinsky model requires some caution.For example, Moin and Kim (1982) mentioned that it is difficult to capture near-wall turbulent behaviors from LESs of turbulent channel flows using the conventional Smagorinsky model.Basu et al. (2008) conducted LESs of a diurnally varying atmospheric boundary layer and mentioned that the static type of SGS model doesn't reproduce a clear low-level jet in the nighttime although no significant differences between the static and dynamic type of models are found in the daytime.On the other hand, we have performed LESs of turbulent flows over a two-dimensional hill (Nakayama and Nagai, 2010), in building arrays with different obstacle densities (Nakayama et al., 2013), and in an actual urban area (Nakayama et al., 2014), and showed that the standard Smagorinsky model can produce reasonable results in comparison to wind tunnel experimental data.Therefore, the standard Smagorinsky model is used in our LES model because of its wide applicability to various environmental flows.
The filtered temperature and concentration transport equations are as follows; where c is concentration.Pr SGS , Sc SGS , α, and κ are the turbulent Prandtl and Schmidt numbers, and molecular diffusivity for temperature and concentration fields, respectively.The subgrid-scale scalar flux is also parameterized by an eddy viscosity model.Although the turbulent Prandtl and Schmidt numbers also should be dynamically determined depending on the flow type, both numbers are set to a constant value 0.71.The coupling algorithm of the velocity and pressure fields is based on the marker-and-cell method with the secondorder Adams-Bashforth scheme for time integration.The Poisson equation is solved by the successive over-relaxation method.For the spatial discretization in the basic equations, a second-order accurate central difference scheme is used.However, for the advection term of the concentration transport equation, cubic interpolated pseudo-particle (Takewaki et al., 1985) is used.

Wind-tunnel experimental data
In this study, we compare the simulated results with the experimental data of Fackrell and Robins (1982) for a NBL case, Ohya et al. (1997) for a SBL case, and Ohya and Uchida (2004) for a CBL with a capping inversion case.The Reynolds numbers Re defined as U ∞ δ/ν are 320 000 for the NBL case, 110 000 for the SBL case and 19 400 for the CBL case.U ∞ and δ are a free-stream velocity and boundary layer thickness, respectively.The bulk Richardson number Ri δ defined as g δ / 0 U 2 ∞ is 0.12 and −0.45 for the SBL and CBL cases, respectively.Here, 0 , and are average absolute temperature in boundary layer, and temperature difference between temperature of ambient air ∞ and surface temperature s .The bulk Richardson number based on the boundary layer thickness is generally defined using the potential temperature.However, within a scale of wind tunnel experiments which have a height of a few meters, the atmospheric pressure can be considered to be constant in the vertical direction, which indicates that the absolute temperature equals the potential temperature.Therefore, the Richardson number is defined using the absolute temperature here.
In the experiment of Fackrell and Robins (1982), a neutral boundary layer flow was generated by spires and a plume was released from the point source located at the downwind position of 6.7δ and elevated at the height of 0.19δ.In the experiment of the SBL case by Ohya et al. (1997), first, air flow with a temperature profile was imposed and a turbulent boundary layer was generated by an obstacle.Then, a stable boundary layer was produced by cooling the floor at a large downwind distance from the entrance of the test section.In their experiment, a downwind distance of 36.2δ was required to obtain a fully developed stable boundary layer.In the CBL with a capping inversion experiment by Ohya and Uchida (2004), first air flow with a temperature profile was imposed at the entrance of the test section.The convective boundary layer was produced by the obstacle and the entirely heated floor.In their experiment, a downwind distance of 24.1δ was required to obtain a fully developed convective boundary layer flow.

Computational settings
A driver region is set to generate a basic turbulent boundary layer flow by a recycling technique of Kataoka and Mizuno (2002) at the upstream part of model domain for each case.In this turbulent inflow technique, only fluctuating components are extracted at a recycling station and recycled back to the inlet boundary.
At the inlet boundary, mean wind velocity profile of the power law 0.14 for the NBL and CBL cases, and 0.25 for the SBL case is imposed.At the recycle station, a nearly target temperature profile is imposed for the SBL and CBL cases, and a thermally-stratified boundary layer is spatially developed.For the velocity field, the Sommerfeld radiation condition is applied at the exit.At the top, a free-slip condition is imposed for streamwise and spanwise velocity components, and vertical velocity component is set to be zero.At the side, a periodic condition is imposed.At the ground surface, a no-slip condition is imposed for each velocity component.For temperature and concentration fields, zero-gradient is imposed at all boundaries.However, for a temperature field, θ s = 0 and θ s = ∞ is set at the ground surface for the SBL and CBL cases, respectively.θ s is instantaneous surface temperature.The time step interval t U ∞ /δ is about 0.003 ( t: time step).The maximum Courant-Friedrich-Levy number is about 0.1.The length of the simulation run to calculate the time averaged values T U ∞ /δ (T : averaging time) is 200.The length of the simulation run before releasing a plume is 300.
The sizes of the computational models for the NBL and SBL cases are 26δ × 2δ × 2δ in the streamwise, spanwise and vertical directions.The number of grid points is 650 × 100 × 90 for both cases.The model size and the number of grid points for the CBL case is 26δ × 5δ × 3δ and 650 × 250 × 96 in the streamwise, spanwise and vertical directions, respectively.Because large-scale convective WT (Fackrell,1982)_u' WT (Fackrell,1982)_v' WT (Fackrell,1982) motions are formed in CBL, a model domain size should be large enough to capture them.Therefore, more number of grid points is set than that of the NBL and SBL cases.The streamwise and spanwise grid spacing is uniform and the vertical grid spacing is stretched from 0.002δ to 0.09δ for each case.Re is set to 12 000 for each case and Ri δ is set to the same value as that of the experiments for the SBL and CBL cases.The release point of a tracer gas is set at a distance of 11.0δ downstream from the inlet boundary and elevated with a height of z/δ = 0.19 for each case.

Turbulence statistics under various thermal conditions
Figure 1 compares the LES results of mean wind velocity (U ) and turbulence intensities (u , v , w ) with the wind tunnel experimental data of Fackrell and Robins (1982) for the NBL case.The mean wind velocity profile of the LES is consistent with the experimental data.Although the turbulence intensities of the LES are a little underestimated, the shape of each vertical profile is similar to those of the experimental data.
Figure 2 compares the LES results with the experimental data of Ohya et al. (1997) for the SBL case.They conducted wind tunnel experiments of SBL flows with a wide range of Ri δ and showed that SBLs show two different distribution patterns of vertical profiles of turbulence statistics based on Ri δ ; a weak stability flow for Ri δ < 0.25 and a strong stability flow for Ri δ > 0.25.Based on this classification, the present LES case for the SBL flow corresponds to a weak type.The mean wind velocity profile of the LES is consistent with the experimental data.The turbulence intensity especially for the vertical component of the LES is much smaller than the experimental data of Ri δ = 0.12.It should be noted that the ratio of the vertical component turbulence intensity to the free-stream velocity of the stable boundary layer flow by Ohya et al. (1997) is comparable to that of the neutral z/δ U/U ∞ WT (Ohya, 1997)  boundary layer flow by Fackrell and Robins (1982).However, the shape of the vertical profile corresponds to that of the experiments of a weak stability type Ri δ < 0.25.The vertical profile of the mean temperature is in good agreement with experimental data of Ohya et al. (1997).The r.m.s.(root mean square) temperature and vertical heat flux are entirely overestimated and underestimated in comparison to the experimental data of Ri δ = 0.12.As described in Sect.2, it is pointed out that the conventional Smagorinsky model cannot accurately capture near-wall turbulent behaviors especially for stable boundary layers (Basu et al., 2008).These differences are due to the use of the static Smagorinsky model.However, the shape of the vertical profile corresponds well with that of the experiments for a weak type.The values of the ratio u * /U ∞ and θ * / in this LES are 0.021 and 0.026, while those are 0.026 and 0.025 in the experiments of Ohya et al. (1997).Here, u * and θ * are friction velocity and friction temperature, respectively.Those values are found to be comparable to those of the experiments.However, in order to examine a turbulent Reynolds number effect in detail, turbulence statistics profiles in inner scaling have to be investigated in comparison to wind tunnel experimental data.
The LES results for the CBL case are shown in Fig. 3. U m is mean wind velocity in the range where the values are constant in the middle part of the CBL.w * is convective velocity scale defined as (g Q s δ/ 0 ) 1/3 .Here, Q s is the maximum value of vertical heat flux w θ near a ground surface.The mean wind velocity and turbulence intensities are in good agreement with the experimental data of Ohya and Uchida (2004).The characteristics such as a constant profile of mean temperature in the main portion of the CBL and the patterns of vertical distributions of the r.m.s.temperature and vertical heat flux are similar to those of the experiments.Deardorff (1972) showed that turbulent eddies have tendencies to become elongated by the mean wind shear and form rolls for a weakly unstable boundary layers −z i /L < 4.5 z/z i w' 2 /w * 2 WT (Ohya, 2004)_Rib=-0.74WT (Ohya, 2004) WT (Ohya, 2004)_Rib=-0.74WT (Ohya, 2004)_Rib=-0.45LES_Rib=-0.45WT (Ohya, 2004)_Rib=-0.74WT (Ohya, 2004) WT (Ohya, 2004)_Rib=-0.74WT (Ohya, 2004)_Rib=-0.45LES_Rib=-0.45 (corresponding to u * /w * > 0.45) while the ordering effect of the wind shear gives way to the random orientation of convective structures for strongly unstable boundary layers −z i /L > 10 (corresponding to u * /w * < 0.34).z i and L are the inversion height and the Obukhov length, respectively.Moeng and Sullivan (1994) performed LESs of planetary boundary layers including extremes of shear and buoyancy forcing and intermediate cases, and implied that CBL flows have shear-and buoyancy-driven flows with u * /w * around 0.65.In this study, u * /w * is 0.46 which slightly exceeds the limit shown by Deardorff.From their results, the CBL flow in our LES model is considered to be driven by not only buoyancy flows but also shear flows.
From these results, for the NBL case, the boundary layer is found to be reasonably simulated.For the SBL case, the vertical component turbulence intensity and vertical heat flux are underestimated, and the r.m.s.temperature is overestimated.However, the shape of those profiles is similar to the experimental data for a weak type.For the CBL case, the turbulence characteristics are generally similar to the experimental data.These indicate that the LES model provides physically reasonable results depending on atmospheric stability.

Plume spreads under various thermal conditions
The most common stability classification scheme is the Pasquill-Gifford (P-G) (Turner, 1970), which defines six stability classes namely A (highly unstable), B (moderately unstable), C (slightly unstable), D (neutral), E (moderately stable), and F (extremely stable).A comparison of the plume spreads with Pasquill-Gifford curves and those of the wind tunnel experimental data (Fackrell and Robins, 1982) where σ y and σ z are spanwise and vertical plume spreads, respectively.y c and z c are the peak locations of the spanwise and vertical distributions of mean concentration at each downstream position, respectively.According to a classification for the P-G category by Snyder (1981), the bulk Richardson number values evaluated based on differences between 2 m and 10 m height of -0.02, -0.01, 0, 0.004 and 0.05 correspond to the categories B, C, D, E, and F. In this study, the bulk Richardson number values are estimated −0.018 and 0.025 for the CBL and SBL cases, which correspond to the categories B and between E and F, respectively.The σ y and σ z of plumes in the experimental data of Fackrell and Robins (1982) are found to be distributed around the category D. On the other hand, those of the LES are overestimated especially near the point source for each case.Michioka et al. (2003) examined the sensitivity of the plume spreads to grid resolution for a point source which has 1.0 and 10 times the real diameters of the point source in comparison to wind tunnel experimental data.They showed that the plume spreads especially near the point source are largely overestimated in the coarse grid resolution, while those in the fine grid resolution are consistent with the experimental data.The point source is a tube with a diameter of d s /δ = 0.007 in the experiments of Fackrell and Robins (1982), while the mean size is d s /δ = 0.016 in the LES (d s : a tube diameter of a point source).Therefore, it is considered that the overestimation especially near the point source is due to numerical diffusion by coarse grid resolution for the point source.The differences at large downwind distances from the point source are due to the underestimation of the turbulence intensities.These tendencies are observed for the SBL and CBL cases.The increases of the vertical spreads in the CBL case are constrained at a large downwind distance because of the capping effect of inversion.However, the distribution patterns of the spanwise and vertical plume spreads in response to the difference in atmospheric stability are well reproduced.This indicates that our LES model provides physically reasonable results.

Conclusions
We performed LESs of plume dispersion under various thermally-stratified boundary layers and examined the usefulness of our approach by comparing the simulated results with wind tunnel experimental data.For the NBL case, the turbulent boundary layer is reasonably produced.For the SBL case, the turbulence statistics such as the vertical component turbulence intensity, r.m.s.temperature, and vertical heat flux especially in the lower part of the boundary layer are quantitatively different from the experimental data.It is pointed out that the Smagorinsky model cannot reproduce near-wall turbulent behaviors by Moin and Kim (1982) and Basu et al. (2008).These differences are due to the use of the static Smagorinsky model.However, the shape of the distribution patterns corresponding to that of experimental data for a weak stability flow is obtained.This indicates that a weaktype SBL flow is successfully produced.For the CBL case, the turbulence characteristics are generally reproduced well.Focusing on dispersion fields, the spanwise and vertical spreads of plumes are overestimated near the point source for each case.This is due to numerical diffusion by coarse grid resolution for the point source.At a large downwind distance from the point source, the spanwise plume spreads are underestimated for each case.This is due to the underestimation of the turbulence intensity for the spanwise component.However, the plume spreads of the LES are distributed nearly along the P-G curves depending on atmospheric stability.
From these results, the plume dispersion patterns are successfully simulated depending on atmospheric stability.It can be concluded that our LES model provides physically reasonable results.

Figure 1 .
Figure 1.Vertical profiles of (a) mean wind velocity and (b) turbulence intensi 25 wind component of a NBL flow.26

Figure 1 .
Figure 1.Vertical profiles of (a) mean wind velocity and (b) turbulence intensities for each wind component of a NBL flow.

Figure 2 .
Figure 2. Vertical profiles of (a) mean wind velocity, (b) turbulence intensity for streamwise component, (c) turbulence intensity for vertical component, (d) mean temperature, (e) r.m.s.temperature, and (f) vertical heat flux of a SBL flow.

Figure 3 .
Figure 3. Vertical profiles of (a) mean wind velocity, (b) turbulence intensity for streamwise 24 component, (c) turbulence intensity for vertical component, (d) mean temperature, (e) r.m.s. 25 temperature, and (f) vertical heat flux of a CBL with a capping inversion flow.26

Figure 3 .
Figure 3. Vertical profiles of (a) mean wind velocity, (b) turbulence intensity for streamwise component, (c) turbulence intensity for vertical component, (d) mean temperature, (e) r.m.s.temperature, and (f) vertical heat flux of a CBL with a capping inversion flow.

Figure 4 .
Figure 4. Streamwise variation of (a) spanwise and (b) vertical plume spreads under NBL, 24 SBL, and CBL with a capping inversion flows.The characters A-F indicate the thermal 25 stability classification scheme of the Pasquill-Gifford.26

Figure 4 .
Figure 4. Streamwise variation of (a) spanwise and (b) vertical plume spreads under NBL, SBL, and CBL with a capping inversion flows.The characters A-F indicate the thermal stability classification scheme of the Pasquill-Gifford. σ