Application of Adaptive Grid Refinement to Plume Modeling

A. Sarma, N. Ahmad, D. Bacon, Z. Boybeyi, T. Dunn, M. Hall, and P. Lee.
Science Applications International Corporation, 1710 Goodridge Dr., McLean, Virginia, USA.
E-Mail: sarma@apo.saic.com

Abstract

Resolving the flow and pollutant concentrations at the highest possible resolution is of paramount importance in atmospheric chemistry calculations especially when dealing with chemical reactions in plumes. The Operational Multiscale Environment model with Grid Adaptivity (OMEGA) is used to explore the modeling of pollutant plumes. OMEGA is built upon an unstructured adaptive grid made up of triangular prisms. OMEGA also has an embedded Atmospheric Dispersion Model (ADM). The ADM uses a puff approach to disperse the pollutants. The particles are treated as centroids of growing puffs, with the growth determined by the ambient turbulent characteristics. OMEGA also features a particle diffusion algorithm using a Monte Carlo method with a receptor-oriented concentration calculation algorithm.

1 Introduction

Over the past few decades the field of numerical modeling has seen considerable advance. Considerable progress has also been made in the fields of algorithm development and parameterizations of physical processes. Numerical models are being used to simulate environments with spatial resolutions down to a few kilometers. Even though these resolutions may be deemed high for the purposes of weather prediction, they are inadequate for the purposes of plume transport and corresponding atmospheric chemistry. Chemical reactions are strongly dependent on the concentrations of the reactant species. To resolve the concentration peaks and variations in plumes require that the plumes be resolved at a very high resolution. Traditional numerical models that are used in meteorology use a nested grid approach to increase resolution over a small region of interest. Examples of nested grid models include the Regional Atmospheric Modeling System, RAMS (Pielke, et al.1) and Penn State/NCAR Mesoscale Model, MM5 (Grell, et al.2). This method involves solving the model equations on a small high-resolution grid in the region of interest while the entire domain may be solved at a much coarser resolution. However, the nested grid method requires a-priori knowledge of the solution so that the simulated features such as a plume or frontal activity can be contained within the high-resolution grid.

At the same time these advances were being made in the field of numerical weather prediction, researchers in computational fluid dynamics (CFD) were experimenting with adaptive grid techniques to model flow around complex shapes. These techniques marked a significant departure from the traditional structured grids and employed unstructured grids in two and three dimensions. The model described here, Operational Multiscale Environment model with Grid Adaptivity (OMEGA), merges the grid technology that evolved from this CFD research to the physical parameterizations developed for numerical weather prediction.

2 Model Description

OMEGA is a multiscale, non-hydrostatic atmospheric simulation model with an adaptive grid that permits a spatial resolution ranging from roughly 100 km to less than 1 km without the need for nested grids. The basic features of the OMEGA model are provided in Table 1. OMEGA is a fully non-hydrostatic, three-dimensional prognostic model. It is based on an adaptive, unstructured triangular prism grid that is referenced to a rotating Cartesian coordinate system. The model uses a finite volume flux based numerical advection algorithm derived from Smolarkiewicz3. OMEGA has a detailed physical model for the planetary boundary layer (PBL) with a 2.5 level Mellor and Yamada4 closure scheme. OMEGA uses a modified Kuo scheme to parameterize cumulus effects (Kuo5, Anthes6), and an extensive bulk-water microphysics package derived from Lin et al.7. OMEGA models the short-wave absorption by water vapor and long-wave emissivities of water vapor and carbon dioxide using the computationally efficient technique of Sasamori8. OMEGA uses an Optimum Interpolation analysis scheme (Daley9) to create initial and boundary conditions and supports piecewise four dimensional data assimilation using a previous forecast as the first guess for a new analysis. Finally, OMEGA contains both Eulerian (grid based) and Lagrangian (grid free) dispersion models embedded into it. The Lagrangian dispersion model can be configured in one of two ways - puff and tracer modes. In the puff mode each particle can act as the centroid of a Gaussian spheroid whose dimensions are determined by the turbulent characteristics of its environment; in tracer mode, each particle acts as a passive tracer, which feels the effects of the ambient turbulence via turbulent velocity fluctuations, which are added to the advective velocity of the particle.

2.1 The OMEGA Grid Structure

A unique feature of OMEGA is its unstructured grid. Unstructured grid cells in the horizontal dimension can increase local resolution to better capture topography or the important physical features of the atmospheric circulation and cloud dynamics. In many fields of application there is recognition that these methods are more efficient and accurate than the structured logical grid approach used in more traditional codes (Baum and Löhner10, Schnack, et al.11). To date, however, unstructured grids and grid adaptivity have not been used in the atmospheric science community (Skamarock and Klemp12). OMEGA represents the first attempt to use this CFD technique for atmospheric simulation.

OMEGA is based on a triangular prism mesh that is unstructured in the horizontal dimension and structured in the vertical (Figure 1). The rationale for this mesh is the physical reality that the atmosphere is highly variable horizontally, but always stratified vertically. In addition, the structured vertical grid enables the use of a tri-diagonal solver to perform implicit solution of both vertical advection and vertical diffusion relaxing the limitation on the timestep.

2.2 The OMEGA Grid Adaptivity

Since the accurate solution of any complex computational problem depends on fine spatial discretization of the computational domain, the accurate representation of multiscale events in numerical models has long been a principal issue in computational fluid dynamics. For example, one typically desires to capture not only the development and evolution of small-scale features but also their interaction with and influence upon the larger scale flow. This is a particularly important requirement in atmospheric models, because numerous events such as fronts, clouds, and plumes are not only relatively localized with respect to their environment, but are also forced on scales larger than their own. Because practical limitations in computer size and speed prohibit the use of uniformly high spatial resolution appropriate for the smallest scales of interest, numerous techniques have been developed to deal with multiscale flows.

Grid nesting techniques involve the sequential placement of multiple finer scale meshes in desired regions of the domain so as to provide increased spatial resolution locally. A principal limitation of grid nesting technique is that one must know a priori, and for the duration of the calculation, which regions of the domain will require high spatial resolution. Also, the propagating waves discontinuously change their speeds upon passing from one mesh to the next and reflect off the boundaries of each nest. This has been a topic of great concern during the past two decades.

OMEGA co-ordinate system

Figure 1: OMEGA coordinate system and the vertical grid structure.

The flexibility of unstructured grids and their ability to adapt to transient physical phenomena are features that make unstructured grid algorithms for partial differential equations so powerful. Grid adaptivity improves the fidelity of the numerical schemes, eliminating some of the errors that plague existing models. The improvement comes from the ability to adapt the grid structure in the vicinity of rapidly changing flow features in the atmosphere.

Using unstructured grids also eliminates the limitations of nested grid techniques. The main advantage of unstructured grids is the ease with which dynamic or solution adaptation can be implemented. There is no longer a need for extensive user interaction for creating topologies of complex terrain features. The whole procedure can be fully automated. The unstructured grid also has a smooth and continuous transition from coarse to fine regions within the whole domain, thus eliminating wave propagation errors, which occur due to interactions with the different nest boundaries in structured grids.

To the best of our knowledge, OMEGA is the only operational atmospheric flow model which fully exploits the advantages and flexibility of unstructured grids. It can adapt its grid both statically and dynamically to different criteria such as fronts, clouds, hurricanes, plumes, etc. For real-time flow predictions, the capability of grid adaptivity, given the computational constraints, becomes important. This capability is also crucial in responding to emergency scenarios such as release of hazardous materials. OMEGA with its grid adaptation capability has a unique advantage over other atmospheric flow models in providing accurate solutions quickly in an operational setting. The grid adaptivity in the OMEGA model takes place in two different ways: (1) static grid adaptivity, and (2) dynamic grid adaptivity. The static adaptivity refers to the method of adapting the grid to static features such as terrain and land-water features, while the dynamic grid adaptivity refines the grid in response to time-dependent features such as plume or frontal movement.

The grid can be made to adapt to several features. The adaptation criteria will vary with the fields that are used as a guide to adaptation. As different fields may have different levels of significance, adaptation is performed in OMEGA by using a single cost function. This function is built by combining the weights of the various fields. Wherever the cost function exceeds a set threshold, vertices are added, and wherever the cost function is less than a threshold vertices are deleted. These steps are followed by vertex reconnection and relaxation steps. In the case of dynamic adaptation, all the model fields are interpolated from the old grid to the new grid using a pseudo-Laplacian technique. In this study we explore the utility of dynamic adaptation in plume simulation and its implications on plume chemistry.

3 Description of the Simulations

To study the effects of dynamic adaptation on plume structure, three simulations were performed - 1) a baseline simulation with no adaptation in which the grid was not changed during the simulation, 2) one with adaptation in which the grid was allowed to adapt around particle locations, and 3) one in which the grid resolution was held uniformly at a high resolution. To separate out the effects of terrain, the domain was chosen with a flat surface. A simple steady-state flow was imposed in the domain with a wind speed of about 5 ms-1. Figure 2 shows the initial grid and the initial wind conditions.
 
 

Coarse grid Initial velocity field
Hi-resolution grid Adapted grid at +12 hr

Figure 2: (Top left) Grid used for the baseline simulation., (top right) initial velocity distribution, (bottom left) grid used for the high-resolution simulation and (bottom right) the grid at 12 hours for the simulation with grid adaptation.
 

To simulate the plume a continuous source was provided in one cell (Figure 2a). The source strength was arbitrarily set equal to 1.0 unit s-1; the same value was used for all the simulations. The initial grid has a resolution of about 25 km over most of the domain. To better define the source, the grid resolution was increased to roughly 2 km in its vicinity. Figures 2c and 2d show the high-resolution grid as well as the grid that had adapted to the solution after 12 hours. Note the high-resolution region along the center resulting from the adaptation. The adaptation was performed based on the location of the Lagrangian particles that were released at the centroid of the cell containing the Eulerian source. These particles were advected by the Atmospheric Dispersion Model built into OMEGA. As the flow was uniform and steady, the particles formed a very narrow plume along the centerline of the domain. A region of influence was specified around each particle location, within which the model attempts to increase the resolution to specified levels. After the grid refinement, all the model variables are interpolated to the new grid geometry.

4 Simulation Results

Figure 3 shows the concentration distribution after 12 hours of transport for the baseline simulation, the simulation with grid adaptivity, and the high-resolution simulation. The baseline simulation lets the plume advect in a grid without any dynamic adaptation. Hence once the plume moves out of the high-resolution source region, it diffuses very rapidly due to the coarse resolution. The simulation in which the grid was allowed to adapt to the solution was better able to capture the higher concentrations in the plume as well as the concentration gradients. These are further evident in the vertical sections, shown in the right-hand panels of Figure 3, taken along the centerline of the plume. Figure 4 shows a comparison of the plume concentrations along the centerline. These are taken along a line that extended along the center of the plume at an altitude of about 100m. Note that the baseline simulation shows concentrations of only about 20% when compared to the high-resolution case. The adaptation case is able to resolve up to 60%.

The ability to resolve plume concentrations has a direct impact on atmospheric calculations. The concentrations of the product species are directly proportional to the reactants. With a multitude of reactions and with species being produced and used at the same time, the end result on atmospheric composition can be significant. To explore this premise, we introduced a simple single-reaction chemistry module to OMEGA. This reaction converted SO2 to H2SO4 by the OH radical (Lurmann et al.13) and is given by the equation

An iterative scheme with OH and HO2 radical concentration convergence as termination criteria was used to calculate pseudo-steady state concentrations of the gas phase species. In this experiment the source was set to release 21 tons of SO2 per day, which corresponds to the emission from a 1000MW coal-fired power plant (Burns and McDonnell Engr. Co.14). The total amount of sulfate produced within the computational domain as a function of time is shown for the three cases in Figure 5. The simulation with no adaptation produce only 60% of the total sulfate when compared with the high-resolution run while the adaptation was able to increase this to 85%.

The principal reason for adaptive grid is the efficient use of computer resources by providing high resolutions only where they are needed. In this case the savings in computer time was significant. The high-resolution case took 12 times as much CPU time as the baseline case. The adaptive grid case able to achieve a significantly improved solution with only a factor of 2.5 increase in CPU time when compared to the baseline case (a factor of 5.3 lower than the high-resolution case.)
 

Coarse: Concentration field vertical x-section
Adaptive: Concentration field vertical x-section
Hi-res: Concentration field vertical x-section

Figure 3: Plume concentration distribution for (top row) no-adaptation/low-resolution, (center row) adaptation and (bottom row) no-adaptation/high-resolution cases. The left-hand column shows a horizontal slice at an altitude of about 100m (AGL). The right-hand column shows a cross section along the plume centerline. The vertical dimension has been exaggerated and represents a total height of 250m. The contours represent values of 0.1 (outermost contour), 0.5, 1, 5, 10, 50, 100, 500, and 1000 (arbitrary units).
 

5 Conclusions

In this paper we have the application to atmospheric chemistry, in particular to plume modeling, of a new atmospheric modeling system which is built upon the concept of unstructured adaptive grids. This methodology allows for the efficient use of computing resources by increasing the spatial resolution of the numerical model only in areas of the domain where an increase in resolution is warranted. The resolution adjustment is done during the simulation as the solution evolves and is performed in an automatic manner. The ability to resolve pollutant plumes at high-resolution increases the fidelity of any atmospheric calculations that may be performed in the plume.
 

Centerline Concentration

Figure 4: Plume concentration along the plume centerline for the baseline, high-resolution and adaptation cases.
 
 

SO2

Figure 5: Total sulfate formed within the domain for the baseline, high-resolution and grid adaptive cases.
 

References

Pielke, R. A., Cotton, W. R., Walko, R. L., Tremback, C. J., Lyons, W. A., Grasso, L. D., Nicholls, M. E., Wesley, D. A., Lee, T. J., & Copeland, J. H., A Comprehensive Meteorological Modeling System - RAMS. Meteorology and Atmospheric Physics, 49, pp. 69-91, 1992.

Grell, G. A., Dudhia, J., & Stauffer, D. R., A Description of the Fifth-Generation Penn State/NCAR Mesoscale Model (MM5). NCAR/TN-398+IA, National Center for Atmospheric Research, Boulder, Colorado, USA, 1994.

Smolarkiewicz, P. K., A fully multidimensional positive definite advection transport algorithm with small implicit diffusion. Journal of Computational Physics, 54, 325, 1984.

Mellor, G. L. & Yamada, T., A hierarchy of turbulence closure models for planetary boundary layers. Journal of Atmospheric Science, 31, pp. 1791-1806, 1974.

Kuo, H. L., On formation and intensification of tropical cyclones through latent heat release by cumulus convection. Journal of Atmospheric Science, 22, 40-63, 1965

Anthes, R. A., A cumulus parameterization scheme utilizing a one-dimensional cloud model. Monthly Weather Review, 105, 270-286, 1977.

Lin, Y-L., Farley, R. D., & Orville, H. D, Bulk parameterization of the snow field in a cloud model. Journal of Climate and Applied Meteorology, 22, pp. 1065-1092, 1983.

Sasamori, T., A linear harmonic analysis of atmospheric motion with radiative dissipation. Journal of the Meteorological Society of Japan, 50, pp. 505-518, 1972.

Daley, R., Atmospheric Data Analysis. Cambridge University Press, 457 pp., 1991.

Baum, J. D., & Lohner, R., Numerical simulation of shock-elevated box interaction using an adaptive finite element shock capturing scheme. Proceedings of the 27th Aerospace Science Meeting, AAIA-89-0653, 1989.

Schnack, D. D., Lottati, I., Mikic, Z., & Satyanarayana, P., MHD simulation on an unstructured, adaptive mesh (abstract), EOS, Transactions of the American Geophysical Union, 74, SH11A-16 (fall meeting), 1993.

Skamarock, W. C. & Klemp, J. B., Adaptive grid refinement for two-dimensional and three-dimensional nonhydrostatic atmospheric flow. Monthly Weather Review, 121, 788-804, 1992

Lurmann, F. W., Lloyd, A. C., & Atkinson, R., A Chemical Mechanism for use in long range transport / acid deposition computer modeling, Journal of Geophysical Research, 91, pp 10905 - 10936.

Burns and McDonnell Engr. Co., Report on the environmental analysis for the Laramie River Station Steam Generating facility of the Missouri Basin power project. Report No. 74-012-4-001, 1975.