New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Spatial reduction algorithm for atmospheric chemical transport models

Communicated by James R. Rice, Harvard University, Cambridge, MA, June 19, 2007 (received for review January 17, 2007)
Abstract
Numerical modeling of global atmospheric chemical dynamics presents an enormous challenge, associated with simulating hundreds of chemical species with time scales varying from milliseconds to years. Here we present an algorithm that provides a significant reduction in computational cost. Because most of the fast reactants and their quickly decomposing reaction products are localized near emission sources, we use a series of reduced chemical models of decreasing complexity with increasing distance from the source. The algorithm diagnoses the chemical dynamics ontherun, locally and separately for every species according to its characteristic reaction time. Unlike conventional timescale separation methods, the spatial reduction algorithm speeds up not only the chemical solver but also advection–diffusion integration. Through several examples we demonstrate that the algorithm can reduce computational cost by at least an order of magnitude for typical atmospheric chemical kinetic mechanisms.
The modeling of tropospheric oxidants (ozone and the hydroxyl radical OH) is of central importance for addressing issues of air quality, aerosol formation and evolution, acid formation, and global budgets of greenhouse gases such as methane (1). Although the chemical factors controlling tropospheric ozone and OH are fairly well established, the computational challenge of modeling concentrations in relation to changes in emissions or climate is enormous. Current photochemical mechanisms describing oxidant chemistry include hundreds of coupled chemical species reacting on time scales ranging from milliseconds to many years (2).
Typical chemical transport models (CTM) combine chemical reactions with advection by a meteorologically predicted flow velocity. The resulting system of equations is extremely stiff, nonlinear and involves a large number of chemically interacting species. The difficulty of solving these equations imposes severe limitations in the spatial resolution of the CTMs, as well as the timescales that can be simulated. This difficulty represents a major obstacle for progress in atmospheric chemistry research. Even with a simplified chemical mechanism of ≈100 species, present computational resources limit the number of CTM grid points to ≈10^{5} to 10^{6}, corresponding to a horizontal resolution ≈100 km and a vertical resolution ≈1 km (3). This is not adequate for resolving the vertical stratification of the troposphere (4), the dynamics of the planetary boundary layer (5), convective processes (6), or sharp horizontal gradients across fronts and coastlines (7). The push for fast, highresolution CTMs will be exacerbated over the next decade by the need to interpret satellite observations of tropospheric ozone and related species (8).
In this paper, we present a numerical algorithm that addresses two major aspects of this problem: the large number of interacting chemical species and extremely large stiffness of the chemical model. Our algorithm is inspired by the observation that atmospheric emissions have most of their fast reactants and their quickly decomposing reaction products localized near the emitter, typically near the ground. Far from the emitter, the fast reactants do not play a significant role. There is thus the opportunity to spatially reduce the chemical mechanism, where a full chemical mechanism is used near emission sources and reduced mechanisms are used far from sources.
To give an idea of the scale of the computational savings that could result from such an algorithm, consider the chemical mechanism of the GEOSChem CTM (9), extensively used for tropospheric chemistry applications (10): ≈65 of the 120 species describe nonmethane hydrocarbons (NMHC) and their oxidation products. These NMHCs are emitted by anthropogenic and natural processes at the surface of continents. Most have short lifetimes. The most important NMHC is biogenic isoprene, which has a lifetime of <1 h against oxidation and whose oxidation products account for 25 species in the mechanism. Although proper simulation of isoprene and other NMHC chemistry is critical near the ground (11), in the remote troposphere, concentrations of NMHCs and their reactive oxidation products are in general very low and a simpler chemical mechanism including 15 species can be sufficient (12).
The idea of reducing chemical mechanisms is quite old. Traditional methods address spatially uniform situations, and exploit the large range of time scales between different chemical reactions, by splitting the species ensemble into subsets that evolve over similar time scales (13–19). This splitting speeds up computation by reducing the size and the stiffness of the individual Jacobian matrices to be inverted.
The spatial reduction mechanism proposed here gradually and continuously discards the equations for fast reacting species with distance from a pollution source. We track a series of chemical boundary layers (CBL), which separate regions where the concentration of a given species is computed with its equation of motion from those where the concentration is determined by extrapolation from the source. The number of domains (multiconnected in general) can be as large as the number of chemical species. The variable CBLs represent the region of the atmosphere affected by fresh emissions; they include not only the meteorological boundary layer but also possibly pollution plumes injected to high altitudes.
The Chemical Mechanism
Emissions of NMHCs and NO_{x} are the principal drivers of complexity in tropospheric oxidant chemistry. The former are emitted by vegetation, combustion and industrial processes whereas the latter are emitted by combustion, soils and lighting. Because the oxidation schemes of all NMHCs follow similar schematics (20), we use propene (C_{3}H_{6}) as representative NMHC for a simplified representation of the chemistry, and include only the most important reactions in its oxidation sequence. This greatly decreases the number of species in the chemical mechanism while retaining the essence of the mechanism complexity.
Table 1 lists the reactions in the reduced mechanism, and Fig. 1 illustrates the major features of the reaction scheme. Oxidation of propene by OH produces a suite of carbon compounds leading eventually to CO_{2}. The time scales of the reactions vary from seconds (lifetime of OH) to minutes (lifetimes of PO_{2}, MCO_{3}, MO_{2}) to hours (lifetimes of other propene intermediates) to weeks or greater (CO).
The evolution of n _{i}(x, t), the number density of the ith chemical species, couples chemical reactions with advection and diffusion through the continuity equation where u is the wind velocity, κ is the turbulent diffusivity, ω̇_{i} is the production rate of the i^{th} chemical species, and s_{i} describes local emissions and nonchemical sinks. The local net production rate R(i) ≡ ω̇_{i} + s _{i}(x, t) of ith species can be separated into the local gross production rate R ^{+}(i) and the local loss R ^{−}(i) terms, i.e., R(i) = R ^{+}(i) − R ^{−}(i). Emissions are in general near the ground but can extend to all altitudes (aircraft, lightning, volcanic emissions, etc.).
Spatial Reduction Algorithm
Introduction.
The spatial reduction algorithm partitions the computational domain for every chemical species into a “fast region” where we calculate the species concentration by solving the full system of equations, and a “slow region” where we calculate the species concentration by extrapolating from the “fast region.” Because the analysis is done on a speciesbyspecies basis, the spatial reduction is gradual with fast equations being eliminated closer to the source.
We explain our algorithm in the context of a very simple example. Consider a three species reaction system, comprised of chemicals A, B, and C. The reactant A decays quickly through A + B → C. The concentration of B is assumed fixed, and there is a (slowly changing) flux f(t) of A at the ground. We assume that A is the fastest reactant, so that τ_{C},τ_{B} ≫ τ_{A}, where τ_{B}, τ_{C}, and τ_{A} are the characteristic decay times for the elements B, C, and A respectively. Additionally, we assume that the only transport mechanism is 1−D diffusion.
Under these conditions, the concentration of A is localized near the ground. The dynamics of A and C are described by
where κ is the diffusion coefficient, and the boundary conditions are −κ∂_{x}
n
_{A}_{x = 0} = f(t). As long as f(t) varies more slowly in time than τ_{A}, the solution to Eq. 2
is approximately
This equation for n
_{A}(x, t) can then be used in Eq. 3
, so that we now have a single equation for n
_{C}(x, t) which is valid far from the source. This represents a (simple) spatial reduction because we now need to solve only one equation far from the ground. This equation is
Eq. 5
uses γ and α as free constants, not necessarily given by
This onedimensional example illustrates the basic principles of spatial reduction: in the fast region, we solve the full chemical mechanism, whereas in the slow region we extrapolate concentrations of the fast reactants from the fast region. In the present example the full reaction system is comprised of two coupled partial differential equations (PDE), whereas the reduced system involves solving only one PDE. The computational cost of solving a system of PDEs using an implicit solver increases at least quadratically with the number of dependent variables. In the present simple example, the CPU time reduction will thus be a factor of four. We note that, in practice, the CPU savings should be greater because by eliminating the fast variables we are eliminating the primary source of stiffness.
Domain Partitioning.
How do we choose the borderline between the fast and slow region? Let us denote the fast and slow regions as D _{f} and D _{s}, respectively, with the regions separated by a moving boundary Γ. In the above example, Γ corresponds to the point x _{b}. It is intuitively clear that in the fast region D _{f} the species density n should either be sufficiently large, or have the potential to become large in a short period since its production rate is large. Our domain partitioning criterion thus relies on the local magnitudes of n and the local net R(i) or gross production R ^{+}(i) rates. Because R(i) ≤ R ^{+}(i), the criteria based on the local gross production rate R ^{+}(i) is more conservative than the one based on the local net production R(i). One definition of D _{f} uses an absolute threshold for n(i) and R(i): for each chemical species n(i) we require that the chemical concentration n(i) > n(i)_{0} and the reaction rate R(i) > R(i)_{0} or R ^{+}(i) > R(i)_{0}, where n(i)_{0} and R(i)_{0} are prescribed thresholds. Alternatively, we can define D _{f} using relative thresholds. Here, we first find the local maxima of n(i) and R(i) (R ^{+}(i)), and then define D _{f} by requiring that n(i) and R(i) are larger than some fraction δ of their values at the local maximum. The computations presented below use a relative threshold.
Matching Conditions and Extrapolation Algorithm.
We discussed above, how we match the solution between the fast domain and the slow domain by imposing continuity of chemical concentrations and fluxes. In the onedimensional example these conditions fixed the parameters α and γ, which allowed extrapolation of the fast species into the slow region through exponential decay. There are two additional complications in multidimensions.
First, our argument for the exponential decay of the fast reactants into the slow region was in the context of a one dimensional model with eddy viscosity. One might question whether this would still hold in multidimensions with advection. In fact, the exponential decay of the fast reactant is a generic feature of a chemical source in multidimensions, with transport by either diffusion or advection. To see this, consider an isolated chemical source of strength q located at the origin. The air is moving with a local velocity u along the x axis and the value of turbulent diffusion coefficient is κ. Assume that the chemical decays with characteristic reaction time τ. The concentration n satisfies the continuity equation which has the exact solution Here, φ is an angle between r and x. Because exponential decay dominates over geometric decay 1/r far from the source, the concentration decays exponentially. This derivation assumes that the chemical source is steady; in practice this means that as long as the source varies on a time scale much longer than the chemical decay time τ, this derivation holds.
The second major complication in multidimensions is that now the matching conditions must be satisfied on the surface separating the fast region and the slow region. Therefore, we need to find a consistent method for extrapolating information on this surface into the slow region through the expected exponential decay. Stated mathematically, given a chemical species denoted by g defined on Γ, we need to extrapolate g so that it decays exponentially with distance from Γ, with the decay rate determined by the rate of change of the g normal to the interface. Note that α_{n} is not constant, but varies over the interface Γ.
The problem of efficiently extrapolating information defined at an interface to a full domain has been previously addressed extensively in the context of the level set method (21). We first define a function d, which measures distance from Γ in the region we wish to extrapolate. As a distance metric, we require that ∇d = 1; additionally it must vanish on the interface (d_{Γ} = 0). Then we define the set of curves Γ that includes the boundary interface Γ as well as other equally distanced from Γ curve i.e., {Γd = const} as well as the set Γ^{⊥}, which are locally perpendicular to Γ (see Fig. 2). To extrapolate α_{n} from Γ to the entire region, we require that α_{n} is constant along each curve of Γ^{⊥}, and hence that on these curves g decays exponentially. Because the normal vector n = ∇d to Γ is tangent to Γ^{⊥}, these requirements are equivalent to the equations Eq. 9 states that the gradient of the local exponential decay rate is orthogonal to the local normal vector n, so that α_{n} is constant along each curve of Γ^{⊥}. Eq. 10 requires that g exponentially decays with rate α_{n} along Γ^{⊥}. Eqs. 9 and 10 can be solved simultaneously and explicitly by marching outwards from the separating interface Γ into the domain D _{s} using the Fast Marching Method (22).
Computational Savings.
How much will the spatial reduction algorithm reduce computation time in the full problem? To estimate the savings, let us assume that there are N reacting species with N _{s} slowreacting species, that the volume of the computational domain is V and that of the fast reacting domain is V _{f}. An implicit solver needs matrix inversion, which causes computational cost increase quadratically or even faster with the number of species. Assuming a constant mesh spacing, the reduced system has a fraction of the number of equations of the original full system. For the GEOSChem mechanism, we expect both V _{f}/V and N _{s}/N to be of order one tenth, hence implying an order of magnitude expected speed up in the computation time. In making this estimate we have used the fact that the spatial reduction does not substantially increase the CPU time; i.e., the computational cost of Eqs. 9 and 10 integration with described above explicit method is negligible (which is indeed true), and have also neglected additional benefits related to stiffness reduction.
Summary of Algorithm.
To summarize, the numerical algorithm consists of three stages to implemented at each timestep. First, we construct a domain partitioning of every individual species based on the current species density n(x) and source S(x) distributions. This domain partitioning implies that there are different numbers of equations (and variables) to be integrated at each grid point. In the second stage, we integrate the fast equations at each grid point to advance the solution to the next time step. Finally in the third stage, we extrapolate the fast species concentration into the slow domains. Note that the domain partitioning is recomputed at every time step, so that the distribution of fast and slow species changes dynamically throughout the simulation.
In the test simulations that follow, we advance Eq. 6 using secondorder Strang operator splitting for the time integration, separating the advection–diffusion part of the calculation from the reaction part. The advection diffusion calculation uses second order spatial central differences and the linearized trapezoidal method for time stepping; the reaction part of the problem is solved by using a Bulirsch–Stoer (23) scheme.
Numerical Experiments
We now proceed to demonstrate the implementation of the algorithm in both one and two spatial dimensions.
OneDimensional.
We first consider the performance of the algorithm in an example with one spatial dimension. To demonstrate the algorithm's ability to locate and adapt to pollution sources, we choose an example where the location of the pollution source changes in time. The challenge here is that the algorithm must automatically detect the regions where a new pollutant plume occurs.
We consider an isothermal atmosphere (T = 298 K), and have the chemical transport take place in the altitude range 0 < x < L, where L = 15 km. The computational domain is uniformly resolved with 100 points. The initial concentration of all species is taken to be 0 except n _{O3} = 10^{12} molecules per cm^{3}. We impose zero flux at the edges of the domain for all species except the pollutants C_{3}H_{6} and NO, where the fluxes at the ground are J_{C3}H_{6} = 10^{11} molecules per s cm^{2} and J _{NO} = 10^{11} molecules per s cm^{2}. We take the turbulent diffusion coefficient κ = 3 × 10^{5} cm^{2} /s.
In addition to the stationary source, we also implement a simple model of a thunderstorm partway through the simulation. Our thunderstorm rapidly transports pollutants from the lower to the upper troposphere. We assume that the thunderstorm starts at the time t _{0} = 3 × 10^{4} s and ends at time t _{1} = 5 × 10^{4} s. The thunderstorm is modeled by a pollutant sink that scavenges chemicals in the lower troposphere (defined to be 0 ≤ x ≤ x _{0} = 1.5 km) and instantly transports them in a “chimney” upwards by a distance x _{1} = 9 km. The pollutant sink of the ith species is taken to be s − _{i}(x, t) = a _{0} n _{i}(x, t) for 0≤ x ≤ x _{0}, and the pollutant source is s _{+} ^{i}(x, t) = s _{−} ^{i}(x − x _{1}, t). The sourcesink pair simulates rapid transport upward by x _{1}. We choose a _{0} = 2.5 × 10^{−5} s^{−1} so that during the thunderstorm half of the density of the original plume is transported upward.
Fig. 3 shows the test results for this model at a time t = 4 × 10^{4} s, when ≈25% of pollutants are transferred from the lower to the upper plume. The solid lines represent a simulation of the full chemical kinetics and the solid dots denote the spatially reduced model, with the relative tolerance parameter δ = 10^{−2}. The gray regions in Fig. 3 denote the automatically generated domain partitioning. The test results demonstrate that the algorithm detects the changes in the source distribution, properly rearranges the domain splitting and maintains accuracy throughout the simulation.
It is worth noting that in the simulation all of the chemical concentrations are nonzero at every point in space, even in the slow regions where their governing reactiondiffusionadvection equations are not calculated explicitly. This is a subtle but important point for maintaining accuracy. For example, it is well known that there are chemicals such as nitrogen oxides (e.g., NO and NO_{2}) in the reaction scheme that cannot be excluded from the chemical mechanism, even if their concentrations are very low.
What is the error in the spatial reduction? If we define the error as the norm of the discrepancy between the spatially reduced model n _{i} ^{0}(x) and the full model n _{i}(x) the largest error is obtained for PO_{2} density distribution e _{PO2} = 3.1 × 10^{−2}. Numerical experiments demonstrate that this error increases by a factor of five if the extrapolation algorithm in the slow region is not used, and the fast reactants are required to vanish in slow regions. Additionally, numerical experiments demonstrate that the error depends linearly on the thresholds δ.
TwoDimensional.
Now we turn to a two dimensional simulation, where the extrapolation between fast and slow regions is not so simple. Our test case assumes that the horizontal and vertical components of the velocity field (u(x, y, t), v(x, y, t)) form a shear layer, u(x, y) = u _{0} + (u _{1} − u _{0})y/y _{max} and v(x, y) = 0, where x and y are the horizontal and vertical coordinates and u _{0} = 10 m/s, u _{1} = 30 m/s. The simulation is carried out in a domain of dimensions 1,500 km × 15 km, with the computational domain uniformly resolved by 80 × 80 points.
We assume the pollution source is two grid spacings wide, occurring between 485 km ≤ x ≤522 km. Like the onedimensional simulation above, the nonzero pollutant fluxes are J_{C3}H_{6} = 10^{11} molecules per s cm^{2} and J_{NO} = 10^{11} molecules per s cm^{2}.
The molar concentrations of three selected species, PRPE, MAP, and O_{3}, with very different characteristic reaction times are shown in Fig. 4 a at time t = 4 × 10^{4}. The domains separating the fast and slow regions for these three species are shown in Fig. 4 b.
Fig. 4 c shows a quantitative comparison of the conventional technique (solid lines) and the spatial reduction algorithm (dots). The blue line shows a one dimensional slice at x = 750 km, whereas the red line depicts a slice at y = 0.5625 km. When δ = 10^{−2} the largest error is e _{PO2} = 3.7 × 10^{−2}.
Conclusions
The spatial reduction algorithm described herein allows great computational savings in the accurate simulation of chemical dynamics in the atmosphere. The algorithm relies on the fact that the concentration of fast reacting species decays rapidly with distance from the pollutant source. Therefore, for each chemical species, we partition the computational domain into two regions, where the full and reduced models are used accordingly. We have demonstrated the versatility of the method by applying it to a reduced mechanism for tropospheric oxidant chemistry. The algorithm works robustly and automatically in one and two spatial dimensions, with timedependent or independent pollution sources.
Our current goal is to implement the algorithm into actual atmospheric solvers, while keeping the details simple enough that it can be easily implemented by the research community as a whole. The most difficult part of the implementation is the extrapolation between fast and slow regions. Although computationally efficient fast marching algorithms have been constructed for two and threedimensional problems (24–27), the technical difficulties associated with complex geometries shape of the moving separating interface could provide prohibitive for widespread implementation. For this reason, we have instituted a major simplification of our algorithm, in which we require that the domain boundaries separating the fast and the slow regions have simple geometries. In our discussion above, we required the shapes of the domain boundaries to correspond to constant relative or absolute thresholds; however, there is no fundamental reason to do so. Increasing the size of the fast domain slightly can result in enormous simplification of the domain shape with only a slight increase in computational cost. The simplest possible shapes are rectangular patches made to be coincident with the existing grids used in current CTMs. We have adopted the fast marching algorithm for this geometry, and are currently working on adaptating and implementating the algorithm for rectangular shapes in GEOSChem.
As an illustration, a threedimensional domain splitting for the chemical PRPE (based on the PRPE density distribution produced by GEOSChem) is shown in Fig. 5. The computational mesh covering the fast subdomain is depicted by sets of dots at different horizontal layers. In the code, the earth's troposphere is uniformly resolved by 91 × 144 points in horizontal x − y planes and by 20 layers of increasing thickness in vertical z direction. Fig. 5 depicts grid crosssections at z = 1, 14, and 17, corresponding to the attitudes h = 0.17, 12.087, and 15.198 km. The multiconnected fast subdomain is well captured although the rectangular patches; for example, Fig. 5 a shows that the fast region essentially coincides with the locations of the continents, where the pollutant concentration is the highest.
Acknowledgments
This research was supported by the National Science Foundation Grant DMS0417810.
Footnotes
 *To whom correspondence should be addressed. Email: brenner{at}seas.harvard.edu

Author contributions: Y.R., M.P.B., and D.J.J. designed research; Y.R., M.P.B., and D.J.J. performed research; Y.R., M.P.B., and D.J.J. analyzed data; and M.P.B. and D.J.J. wrote the paper.

The authors declare no conflict of interest.
 Abbreviation:
 CTM,
 chemical transport model.
 © 2007 by The National Academy of Sciences of the USA
References
 ↵

↵
 Atkinson R ,
 Baulch DL ,
 Cox RA ,
 Crowley JN ,
 Hampson RF ,
 Hynes RG ,
 Jenkin ME ,
 Rossi MJ ,
 Troe J

↵
 Brasseur GP ,
 Orlando JJ ,
 Tyndall G
 ↵

↵
 Chatfield RB ,
 Brost RA
 ↵
 ↵
 ↵

↵
 Holmes CD ,
 Jacob DJ ,
 Yang X ,
 ↵
 ↵
 ↵
 ↵

↵
 Gong WM ,
 Cho HR
 ↵

↵
 Sandu A ,
 Verwer J ,
 Loon MV ,
 Carmichael GR ,
 Potra FA ,
 Dabdub D ,
 Seinfeld JH
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Press WH ,
 Teukolsky SA ,
 Vetterling WT ,
 Flannery BP
 ↵

↵
 Malladi R ,
 Sethian JA
 ↵
 ↵
 ↵
Citation Manager Formats
Sign up for Article Alerts
Jump to section
You May Also be Interested in
More Articles of This Classification
Physical Sciences
Related Content
 No related articles found.
Cited by...
 No citing articles found.