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 on-the-run, locally and separately for every species according to its characteristic reaction time. Unlike conventional time-scale 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 ≈105 to 106, 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, high-resolution 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 GEOS-Chem 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 NOx 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 (C3H6) 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 CO2. The time scales of the reactions vary from seconds (lifetime of OH) to minutes (lifetimes of PO2, MCO3, MO2) to hours (lifetimes of other propene intermediates) to weeks or greater (CO).
Reduced propene 
NO
x
HO
x
O
3 atmospheric chemistry mechanism
Schematic of the reduced chemical mechanism in Table 1. (a) Cycling of HOx radicals, NOx radicals, and ozone. The organic peroxy radicals in the mechanism (PO2, MCO3, MO2) perform functions similar to HO2. (b) Intermediate species in the oxidation of propene to CO2. The concentrations of all 15 species shown in this figure (excepting CO2) are computed within the chemical mechanism.
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 ith chemical species, and si 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 species-by-species 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
and
implied by Eq. 4. This is to emphasize that in more general situations, it is not possible to analytically solve for the fast reactant profiles.
However, γ and α can always be determined by using matching conditions, by requiring that both concentrations and fluxes are
continuous at the boundary between the fast and slow region. Thus, if n
A,f is the concentration of species A as computed by the full solver in the fast region, and the boundary between fast and slow occurs at x = x
b, we require that γ = n
A,f(x
b), and α = −∂x
n
A,f(x
b)/n
A,f(x
b). This method for choosing γ and α ensures that the spatially reduction conserves mass exactly.
This one-dimensional 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 one-dimensional 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).
Schematic picture of the level set curves. The boundary interface Γ is depicted by thick solid line and the two level set curves corresponding to different values of d are depicted by thin solid lines. The curves Γ⊥ are shown as dashed lines.
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 slow-reacting 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 GEOS-Chem 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 second-order 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.
One-Dimensional.
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 = 1012 molecules per cm3. We impose zero flux at the edges of the domain for all species except the pollutants C3H6 and NO, where the fluxes at the ground are JC3H6 = 1011 molecules per s cm2 and J NO = 1011 molecules per s cm2. We take the turbulent diffusion coefficient κ = 3 × 105 cm2 /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 × 104 s and ends at time t 1 = 5 × 104 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 source-sink 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 × 104 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.
Species densities versus distance at t = 4 × 104. The computation is done in 15 km physical domain resolved by 100 points. Solid line, full chemical kinetics; dots, spatially reduced chemical kinetics.
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 reaction-diffusion-advection 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 NO2) 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 PO2 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 δ.
Two-Dimensional.
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 one-dimensional simulation above, the nonzero pollutant fluxes are JC3H6 = 1011 molecules per s cm2 and JNO = 1011 molecules per s cm2.
The molar concentrations of three selected species, PRPE, MAP, and O3, with very different characteristic reaction times are shown in Fig. 4 a at time t = 4 × 104. The domains separating the fast and slow regions for these three species are shown in Fig. 4 b.
Comparison of spatially reduced and full calculations in two dimensions. (a) Concentrations of three selected species C3H6, MAP and O3 vs. (x, y). (b) Moving boundary Γ, which separates the computational domain into two subdomains D f and D r, where full and reduced models are used accordingly. (c) Comparison of these species distribution obtained by the spatial reduction algorithm (spheres and cubes) and conventional method (solid line) depicted in blue and red colors at x = 750 km and at y = 0.5625 km correspondingly, at time t = 4 × 104 for the atmosphere with linearly distributed horizontal shear layer type flow. Vertical and horizontal spatial scales are y = 0 − 15 km and x = 0 − 1,500 km (not shown in the figure) correspondingly. The computation is carried out in 1,500 km × 15 km physical domain uniformly resolved by 80 × 80 points.
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 time-dependent 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 three-dimensional 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 GEOS-Chem.
As an illustration, a three-dimensional domain splitting for the chemical PRPE (based on the PRPE density distribution produced by GEOS-Chem) 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 cross-sections 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.
Domain partitioning at the attitudes h = 0.17, 12.087, and 15.198 km.
Acknowledgments
This research was supported by the National Science Foundation Grant DMS-0417810.
Footnotes
- *To whom correspondence should be addressed. E-mail: 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
