10.5. Eulerian models

10.5.1. Solving atmospheric transport equations

The main idea of any Eulerian models is to solve numerically the atmospheric transport equation (equation 10.9). The atmospheric transport equation is mathematically a second order differential equation, and its solution with the appropriate initial and boundary conditions provides the spatiotemporal evolution of the concentration, i.e., c = c(t,x). There are several numerical methods to solve PDE, one of the most powerful and common method is the so-called “method of lines”. The method consists of two steps: (i) spatial discretization and (ii) the temporal integration of the derived ordinary differential equations (ODEs). Spatial discretization of PDE (the atmospheric transport equation) is performed on a mesh (grid). This reduces the PDE to a system of ordinary differential equations (ODEs) in one independent variable, time. The system of ODEs can then be solved as an initial value problem, and a variety of powerful methods and software tools exist for this purpose.

Let’s consider a simple and 1D form of the atmospheric transport equation

.

(10.23)

This equation describes a simplified situation, where the spread of a pollutant occurs by advection (wind) and turbulent diffusion. We will show the method of lines technique on this example.

Spatial discretization:

There are different discretization techniques (e.g., finite difference method (FDM), finite volume method (FVM), finite element method (FEM)). The basic idea of the spatial discretization is approximating the first and the second order spatial derivatives in (10.23). The simplest one is the finite difference method, and here the derivatives will be approximated by the Taylor’s theorem. Let’s divide the 1D space equidistantly into N grid point (Figure 10.8). We assume that c(x,t) is at least twice continuously differentiable, and examine its first order Taylor expansion

Spatial discretization.

Figure 10.8: Spatial discretization in 1D.

,

(10.24)

.

(10.25)

It can be seen that the first order derivatives can be easily derived that from (10.24) and (10.25) if we assume , thus

(forward difference),

(10.26)

(backward difference).

(10.27)

Those two approximations have truncation error of . It is easy to construct other schemes that have smaller truncation error. Let’s subtract equation (10.25) from equation (10.24), thus we get the central difference scheme

(central difference).

(10.28)

This scheme has truncation error .

For approximation of second order derivatives, let’s add equation (10.26) to equation (10.27). It can be seen that the first order derivatives will be eliminated, and we construct a central difference scheme for second order derivatives

(central difference).

(10.29)

Here the truncation error is .

Temporal integration:

After spatial discretization we obtain a set of ODEs (each ODE describes the variation of concentration in time in a given grid point). The next step is to solve these ODEs (using appropriate initial conditions) with a numerical integrator. There are three numerical integrator schemes: explicit (forward, equation 10.30), implicit (backward, equation 10.31), and Crank–Nicolson (semi-implicit, equation 10.32) method

,

(10.30)

,

(10.31)

,

(10.32)

where and are the concentration at t and tt, respectively. Δt is the time step (integration step).

10.5.2. Operator splitting

We have seen that the basis of nearly all deterministic air pollution models is a system of partial differential equations, the right-hand side of which contains several terms. These terms describe the different physical/chemical subprocesses of the air pollution transport, namely, advection, diffusion, chemical reactions, deposition and emission. (The definition of a sub-process is not unique. For example, we can call advection one sub-process, but horizontal and vertical advection can be viewed as two different sub-processes as well.) The several different terms on the right-hand side make the model problem rather complicated, especially in complex Eulerian models, where all known sub-processes are taken into account. The high complexity of this system requires the application of numerical solution techniques.

The numerical integration of the transport equations is a rather difficult computational task, especially in large-scale and global Eulerian models, where the number of grid-points can range from a few thousand to a few hundred thousand, and the number of chemical species is typically between 20 and 200. If some readily available solver is applied directly to the system, we cannot obtain a sufficiently accurate numerical solution within reasonable computational time. Therefore, procedures that allow us to lead the solution of the system back to the solution of simpler systems for which sufficiently accurate as well as efficient methods exist, are of great importance.

The problems describing the sub-processes have been thoroughly investigated, and appropriate solution techniques have been developed for their numerical solution. I.e., there exist efficient and accurate methods for the diffusion equation, for the advection equation etc., but not for a system which describe all the sub-processes at the same time. This gives the idea of operator splitting: let us divide the right-hand side of the original system into a few simple terms, and solve the corresponding systems – which are connected to each other through the initial conditions – successively in each time step of the numerical integration. This means that we replace the original model with one in which the different sub-processes take place successively in time. This de-coupling procedure allows us to solve a few simpler systems instead of the original, complicated one.

Operator splitting has been successfully applied in the numerical solution of the Navier-Stokes equations. Splitting is also popular in numerical models of ocean circulation, which typically include multiple time-scale motions. The governing equations of the oceans are split into sub-systems that model the fast and slow motions separately. Further fields where operator splitting is successfully applied include cloud physics, biomathematics, astrophysics, and medical science.

As an example of operator splitting from a complex transport-chemistry model, we introduce a splitting method used in the Danish Eulerian Model (DEM splitting) (Dimov et al., 2008). In this long-range air pollution model the subsystems describe the horizontal advection (1), the horizontal diffusion (2), the chemical reactions to which also emissions are added (3), the deposition (4) and the vertical exchange (5):

,

(10.33)

,

(10.34)

,

(10.35)

,

(10.36)

,

(10.37)

for i = 1,…, m, where m is the number of chemical species. In the above splitting procedure, the original system of PDEs has been split into five simpler systems, which are to be solved in each time step one after the other in the following way. Assume that some approximation to the concentration vector (c1,…, cm) at the beginning of the time step has been found. System (equation 10.33) is solved by using this vector as a starting vector. The obtained solution will serve as the initial vector in the treatment of system (equation 10.34), and so on. The solution of the fifth system is accepted as an approximation to the concentration vector at the end of the time step.

The DEM splitting is based on a separation of the different physical processes of the air pollution transport and the separation of the vertical and horizontal directions for advection and diffusion. An alternative of the DEM splitting is the so-called physical splitting, where the sub-operators belong to the five basic air pollution processes (i.e., advection, diffusion, chemistry, deposition and emission), and there is no directional separation.

Operator splitting has the great advantage that the sub-systems are easier to solve numerically than the original system. We can exploit the special properties of the different sub-systems and apply a suitable method for each. In some situations one of the sub-problems may even be solved analytically (semi-analytical method). On the other hand, this de-coupling procedure may cause inaccuracy in the solution, since we have replaced the original model with one in which the sub-processes take place one after the other. This error is called splitting error. Obviously, the success of operator splitting is primarily determined by the effect of the splitting error. Therefore, further splitting schemes with reduced splitting error have been constructed, such as the Marchuk–Strang splitting or the weighted splitting.