A central goal of atmospheric chemistry is to understand quantitatively how the concentrations of species depend on the controlling processes: emissions, transport, chemistry, and deposition. This dependence is expressed mathematically by the continuity equation, which provides the foundation for all atmospheric chemistry research models. In the present chaper we derive the continuity equation in its Eulerian form (fixed coordinate system) and in its Lagrangian form (moving coordinate system). We also describe the methods and approximations used to solve the continuity equation in atmospheric chemistry models. As we will see, the simple models presented in chapter 3 represent in fact drastic simplifications of the continuity equation.




    5.1.1 Derivation


We wish to calculate the number density n(X,t) of a species in a 3-dimensional frame of reference fixed to the Earth. Here X = (x,y,z) is the vector of spatial coordinates, and t is time. Consider an elemental atmospheric volume (dx, dy, dz) as shown in Figure 5-1 :


Figure 5-1 Flux through an elemental volume


The number density of the species in this elemental volume changes with time as a result of:


  • Transport in and out of the volume, characterized by the flux vector F = (Fx, Fy, Fz) which has units of molecules per cm2 per second. In the x-direction, the transport rate (molecules s-1) of molecules into the volume, across the shaded face in Figure 5-1 , is Fx(x)dydz (flux multiplied by area). Similarly, the transport rate of molecules out of the volume is Fx(x+dx)dydz. The resulting change in number density inside the volume dxdydz per unit time (molecules cm-3 s-1) is: [Fx(x)dydz - Fx(x+dx)dydz]/dxdydz = Fx(x)/dx - Fx(x+dx)/dx = - Fx/ x. In the y- and z- directions, the changes are -Fy/ y and -Fz/ z, respectively.
  • Sources inside the volume from emissions and chemical production; we denote the ensemble of these sources as P (molecules cm-3 s-1).
  • Sinks inside the volume from chemical loss and deposition; we denote the ensemble of these sinks as L (molecules cm-3 s-1).


We thus obtain the following expression for the temporal change n/t of number density inside the elemental volume:




where -- = (/x, /y, /z) is the gradient vector and --.F is the flux divergence (flux out minus flux in).


Let us now relate the flux vector F to the local wind velocity vector U = (u, v, w). The flux Fx(x) across the shaded face of Figure 5-1 represents the number of molecules nudydz crossing the face per unit time, normalized to the area dydz of the face; Fx(x) = nudydz/dydz = nu. Similar expressions apply for Fy and Fz. Replacing into (5.1) we obtain




which is the general Eulerian form of the continuity equation. This form is called Eulerian because it defines n(X,t) in a fixed frame of reference. The continuity equation is a first-order differential equation in space and time that relates the concentration field of a species in the atmosphere to its sources and sinks and to the wind field. If U, P, and L are known, then (5.2) can be integrated to yield the concentration field n(X,t).


We saw in section 4.4 how the turbulent nature of atmospheric flow leads to difficulty in characterizing the flux F = nU. This difficulty can be addressed by using the turbulent diffusion parameterization for the turbulent component FT of the time-averaged flux ( section 4.4.3 ) and solving the continuity equation only for the time-averaged concentration n(X,t). When expressed in three dimensions, equation (4.21) for the turbulent flux becomes




where K is a 3x3 matrix with zero values for the non-diagonal elements and with diagonal elements Kx, Ky, Kz representing the turbulent diffusion coefficients in each transport direction. Equation (5.2) is then written as:




where FA = n U is the mean advective component of the flux. The merit of couching the continuity equation solely in terms of time-averaged quantities is that it avoids having to resolve turbulent motions except in a statistical sense. In what follows all concentrations are time-averaged quantities and we will drop the overbar notation.


    5.1.2 Discretization


Solution of the continuity equation in 3-dimensional models of atmospheric chemistry must be done numerically. Consider the problem of calculating n(X, to+Dt) from knowledge of n(X, to). The first step involves separation of the different terms in (5.4) :







Here the "chemistry" term is intended to account for all terms contributing to P and L (these may include emissions and deposition in addition to chemistry). Each term in (5.6) is integrated independently from the others over the finite time step Dt using an advection operator A, a turbulence operator T, and a chemical operator C. The formulation of the advection operator is




with analogous formulations for T and C. The three operators are applied in sequence to obtain n(X, to+Dt):



This approach is called operator splitting. Separating the operators over large time steps makes the calculation tractable but is based on the assumption that advection, turbulence, and chemistry operate independently of each other over the time step. The time step must be kept small enough for this assumption to be reasonable.


To carry out the integration involved in each operator, one can discretize the spatial domain over a grid ( Figure 5-2 ):


Figure 5-2 Spatial discretization of the continuity equation (only two dimensions are shown). Dots indicate gridpoints at which the concentrations are calculated, and lines indicate gridbox boundaries at which the transport fluxes are calculated.


In this manner the continuous function n(X,t) is replaced by the discrete function n(i,j,k, t) where i, j, k are the indices of the grid elements in the x, y, z directions respectively. The advection operator can then be expressed algebraically, the simplest algorithm being:




There are much better, more stable (but more complicated) algorithms to calculate advection over a discrete grid. The turbulence operator is applied on the discrete grid in the same way as the advection operator; concentration gradients at the gridbox boundaries can be calculated by linear interpolation between adjacent gridpoints or with higher-order methods.


The chemistry operator integrates a first-order ordinary differential equation (ODE) describing the chemical evolution of the species at each gridpoint (i,j,k) over the time step [to, to + Dt]. In general, the chemical production and loss of a species p depend on the ensemble {nq} of local concentrations of other species q = 1,...m that produce or react with species p. We may therefore have to solve a system of m coupled ODE's, one for each of the interacting species:




There are various approaches for integrating numerically such a system.


In conclusion, we have arrived through a series of approximations at a numerical solution of the continuity equation. The quality of the solution depends on the grid resolution, the time step, and the algorithms used for the transport and chemistry operators. At its extreme, discretization of the Eulerian form of the continuity equation leads to the simple box models described in chapter 3.




An alternative expression of the continuity equation for a species in the atmosphere can be derived relative to a frame reference moving with the local flow; this is called the Lagrangian approach. Consider a fluid element at location Xo at time to. We wish to know where this element will be located at a later time t. We define a transition probability density Q(Xo, to|X, t) such that the probability that the fluid element will have moved to within a volume (dx, dy, dz) centered at location X at time t ( Figure 5-3 ) is Q(Xo, to|X, t)dxdydz.


Figure 5-3 Motion of fluid element from (Xo, to) to (X, t), including a source P(X',t') along the trajectory


By conservation of total air mass,




where the integration is over the entire atmosphere. Consider an atmospheric species with a source field P(X, t) (molecules cm-3 s-1), The species has a concentration field n(Xo, to) at time to. The concentration field n(X, t) is given by




Equation (5.12) is a general Lagrangian form of the continuity equation. It can be readily modified to include first-order loss terms.


In the Lagrangian form of the continuity equation, transport is described not by the wind velocity U but by the transition probability density Q. Since Q accounts for the stochastic nature of turbulence, there is no need for parameterization of the turbulent transport flux; this is a significant advantage. On the other hand, whereas U is an observable quantity, Q is more difficult to evaluate. Also, the Lagrangian form is not amenable to representation of non-linear chemical processes, because the chemical rates are then affected by overlap between trajectories of the different fluid elements over the time step [to, t]; one cannot assume that the trajectories evolve independently and then sum over all trajectories, as is done in equation (5.12) .


A frequent application of the Lagrangian form of the continuity equation is to the dilution of a pollution plume emanating from a point source. If the turbulent component of U has a Gaussian frequency distribution, then it can be shown that Q also has a Gaussian form. In this case it is possible, with some additional assumptions, to achieve an analytical solution to (5.12) for the Gaussian spreading of the pollution plume. The simple puff models presented in chapter 3 are applications of the Lagrangian formulation of the continuity equation using trivial forms of the transition probability density Q.


Further reading:


Jacobson, M.Z., Fundamentals of Atmospheric Modeling, Cambridge University Press, Cambridge, 1998. Discretization of the continuity equation, transport and chemistry operators for 3-dimensional atmospheric models.


Seinfeld, J.H., and S.N. Pandis, Atmospheric Chemistry and Physics, Wiley, 1998. Eulerian and Lagrangian forms of the continuity equation, analytical Gaussian plume solutions to the Lagrangian form.