Challenges in Using Splitting Techniques for Large-scale Environmental Modeling

Ivan Dimov(1), Tzvetan Ostromsky(1), and Zahari Zlatev(2)	
(1)Institute for Parallel Processing, Bulgarian Academy of Sciences, Acad. G Bonchev str., bl. 25-A, 1113 Sofia, Bulgaria; e-mails : ivdimov@acad.bg , ceco@parallel.bas.bg.
(2)National Institute for Environmental Research, Frederiksborgvej 399, Roskilde, Denmark; e-mail: zz@dmu.dk .  	

Abstract:	Splitting techniques are, to our knowledge, used in all operationally run large scale air pollution models with many scenarios. The modellers believe that really huge computational tasks can be made tractable on the available computers by dividing them into a sequence of smaller and simple sub-tasks. Probably, the first simple splitting procedure for partial differential equations was proposed by Bagrinovski and Godunov in 1957. Since then many different splitting schemes were proposed and studied. A significant progress in splitting analysis was done during the last 3-4 years when the Laplace transformations technique was replaced by more powerful techniques. Several splitting schemes for large scale air pollution models are analysed and tested in a set of numerical experiments with respect to the splitting error. The computational properties of splitting schemes are analysed from algorithmic point of view. Parallel properties of these splitting schemes are also discussed. 
Key words:	Air quality, weather, numerical forecast, atmospheric chemistry, atmospheric transport.

1.	Introduction
The application of splitting procedures in the treatment of large scientific and engineering problems is an excellent tool (and, very often, the only tool) by which huge computational tasks can be made tractable on the available computers. This is achieved by dividing the original problem into a sequence of simpler tasks.
A simple splitting procedure for partial differential equations (PDE) was proposed, as an example, by Bagrinovskii and Godunov (1957). This was probably the first attempt to introduce the splitting technique. Different splitting procedures have been developed and/or studied in many scientific papers, for example: Dyakonov (1962), Yanenko (1962), Marchuk (1968), Strang (1968), Penenko and Obraztsov (1976), Tikhonov and Samarski (1977), Marchuk (1980,1986), Lanser and Verwer (1999), Faragó and Havasi (2002), Csomós et al. (2003). A detailed theoretical study and analysis of some splitting procedures can be found in Marchuk (1980), Marchuk (1988), Faragó and Havasi (2002), Faragó et al. (2002). Some convergence results were recently presented in Faragó and Havasi (2003).
Splitting techniques are successfully used in large-scale air pollution modeling (see McRae (1982), Marchuk (1986), Zlatev (1995), Dimov et al.  (1999,2001), Hunsdorfer and Verwer (2003)). Results related to the use of some splitting procedures in a large air pollution model (DEM in particular) as well as the advantages and disadvantages of these procedures are discussed in this paper.
2.	General description of splitting procedures 
Consider the differential equation
 = A(c(x,t))											(1)
where t[0, T], xD and A is a differential operator.
We shall assume that operator A(c) can be represented as a sum of m operators, which are in some sense simpler than the original operator A(c):
A(c) =  A[i](c)
It should be noted here that the division of the operator A(c) to a sum of simpler operators is in general not unique.
This means that the original problem can be split into the following m sub-problems:
 = A[i](c[i](x,t))		, i=1,2…m				(2)
Moreover, the time-interval [0,T] is divided to N=T/ “small” sub-intervals, where  is called the splitting time-step. The particular splitting procedure depends on the way in which the operator A(c) is represented as a sum of simpler operators as well as on the order in which the smaller problems are treated. 
Assume that the representation of the original operator A(c) as a sum of simpler operators A[I](c) is fixed. Various types of splitting procedures can be constructed by ordering the simpler operators in different ways. Four of them, e.g. sequential splitting, symmetric splitting, weighted sequential splitting and weighted symmetric splitting are considered and analyzed in this paper.
2.1	Sequential splitting procedure
Consider, for example, the Danish Eulerian Model (DEM), described in detail in Zlatev (1995). It is formally represented by the following PDE system:
     = -  -  - +
+  +  +  +			  (3)
+ E + Q(c) − ckT , 

where x = (x,y,z) is the position vector in the 3-D space domain, c = (c1, c2..., cq)T are the concentrations of different pollutants, E = (E1, E2..., Eq)T are the emissions of different pollutants, Q = (Q1, Q2..., Qq)T is the chemical reactions function and k = (k11+k21 , k12+k22 ,..., k1q+k2q)T are the dry/wet deposition coefficients.
(3) was split in 5 separate submodels in the older versions of DEM, when the pseudospectral algorithm was applied in the numerical treatment of the horizontal transport (advection) and the horizontal diffusion (see for more details Zlatev (1995)). 
In the new splitting procedure used in DEM finite elements are used in the discretization of both the horizontal advection and the horizontal diffusion. Therefore, there is no need to split these two processes and they are combined in one (advection-diffusion) submodel. Furthermore, the deposition process is added to the chemical part. In this way, by implementing the new splitting procedure, the number of the sub-models is reduced from five to three. These three sub-models are given below:
 =  -+,								  (4)
 = --++,		(5)
	= 	E + Q(c[3]) − c[3] kT											(6)
These three systems are solved successively at every time-step (t = tn, n=1,2, ...,N), with t0=0 and  tn - tn–1 = . Assume that some approximation of the concentration vector c is known for t = tn . Then it is used as an initial value for solving (4) at time-step  n+1. The solution c[1] of (4) at t = tn+1 is used as initial value when solving (5) at time-step n+1 (i.e., when c[2] is computed at t = tn+1). Similarly, the solution c[2] of (5) at t = tn+1 is used as initial value when solving (6) at time-step n+1 (i.e., when c[3] is computed at t = tn+1). Finaly, the solution c[3] of the last sub-model at  t = tn+1 = tn+  is accepted as an approximation of the solution c of the original system at time-step  n+1. This means that the splitting time-step  n+1  is completed and everything is ready to carry out the same kind of computations at splitting time-step  n+2.
The simple splitting procedure, described above and used by default in the current version of DEM, is called sequential splitting. It is very often useful to describe how the splitting procedure is used by  (i) giving the order, in which the operators  A[i]  are applied, and (ii) indicating the splitting time-step , which is actually used. For the particular case discussed in this sub-section (i.e., the sequential splitting procedure) such a description is given by the following sequence:
(A[1]) (A[2]) ...... (A[m-1]) (A[m]) .									(7)
2.2	Symmetric splitting procedure
Assume that the first m−1 sub-problems in (2) are solved by using a splitting time-step /2. Assume that the last of the simpler sub-problems, the m-th sub-problem, is after that treated by using a splitting time-step . Finally, assume that the computations at the splitting time-step under consideration are completed by performing (again, but in reverse order) the computations with the first m-1 simpler problems by using a splitting time-step /2. It is easy to express these computations by showing just the sequence of operators, as in (7) for the sequential splitting. This is done in (8) for the symmetric splitting (note that only the index numbers of the operators and the splitting time-steps are given; in other words all c[I], i=1,2,..., m, are omitted in the formula given below):
(A[1])/2 (A[2])/2 ...... (A[m-1])/2 (A[m]) (A[m-1])/2 ...... (A[2])/2 (A[1])/2 . 			(8)
2.3	Weighted sequential splitting procedure
Consider again (2). Assume that the computations at the splitting time-step under consideration are carried out by applying the following two procedures (c[I], i=1,2,..., m, are again, as in the previous section, omitted):
(A[1]) (A[2]) ...... (A[m-1]) (A[m])                    		          (9)
and
 (A[m]) (A[m-1]) ...... (A[2]) (A[1])  .								 (10)
It is seen from (9) and (10) that only the order in which the simpler sub-problems are treated is different. In the first of these two formulae the computations are performed by taking sequentially the operators (the index numbers of the operators being increasing). The splitting procedure in the second formula (10) is also sequential, but now the index numbers of the operators are decreasing.
Denote by cforward and cbackward the results obtained when (9) and (10) are applied at a given splitting time-step. Then an approximation:
 c =  cforward + (1−) cbackward        							(11)
of the solution of the original problem (1) at the splitting time-step under consideration is calculated by using some weighting parameter  and one can proceed by the calculations for the next splitting time-step.
The splitting procedures based on the performance of the calculations at each splitting time-step by using the formulae (9), (10) and (11) are called weighted sequential splitting procedures.
The definition of weighted sequential splitting procedures given in this sub-section can be slightly generalized. One can apply ordering with upper indices {r1, r2,..., rm-1, rm} instead of ordering with upper indices {m, m-1,…, 2, 1}, where indices {r1, r2,... , rm-1, rm} are any permutation of {1, 2,..., m-1, m}. Then the procedure defined by (15) and (16) is a special case of the generalized weighted splitting procedures. Probably, it is most natural to use in practice that special case. 

2.4	Weighted symmetric splitting procedure
Assume that the calculations at each splitting time step are carried out by using the following two formulae:
(A[1])/2 (A[2])/2 ...... (A[m-1])/2 (A[m]) (A[m-1])/2 ...... (A[2])/2 (A[1])/2         (12)
 and
(A[m])/2 (A[m-1])/2 ...... (A[2])/2 (A[1]) (A[2])/2 ...... (A[m-1])/2 (A[m])/2 .       (13)
It is seen that the same principles as in the previous section are used. However, symmetric splitting procedures are used both in the forward mode and in the backward mode; in fact, the formula (12) is precisely the same as (8).
Denote again (as in the previous subsection) by cforward  and cbackward the results obtained when (12) and (13) are applied at a given splitting time-step. Then, as in the previous section, the approximation:
c =  cforward + (1-) cbackward           				       (14)
of the solution of the original problem (1) at the splitting time-step under consideration is calculated by using some weighting parameter  and one can proceed by the calculations for the next splitting time-step.
3.	Advantages and disadvantages of the splitting procedures
The replacement of the original problem (1) by a series of simpler problems (2) has several advantages:
•	It is in general much easier to find optimal (or, at least, good) numerical methods for each of the simpler systems in (2) than for the big system arising after the direct discretization of operator A(c(x,t)) in the right-hand side of (1).
•	It is easy to change numerical algorithms and physical mechanisms used in the different sub-problems arising as a results of the implementation of the splitting procedure chosen. Thus, improvements of the code can easily be achieved.
•	The application of splitting procedure facilitates the preparation of efficient codes for different parallel computers. In many cases, the splitting procedure is leading directly to parallel tasks. Consider, for example, the first two sub-models defined by (4) and (5) in the default DEM splitting procedure. It is easily seen that each of these sub-models consists of  q  parallel tasks (where  q  is the number of chemical species studied by the model under consideration).
4.	Comparison of tHE splitting procedures
In this section we are trying to compare the above splitting procedures and to give some recommendations about the choice of a suitable splitting procedure according to two important criteria: the computational cost and the accuracy of the results.
4.1	Comparing the computational costs
•	The sequential splitting procedures have the cheapest computational cost (among the four types of splitting procedures) when the number of the simple operators m is fixed.
•	If a symmetric splitting procedure is to be selected, then it is most efficient to put the most expensive operator in the middle of (8). In this way the total computational cost of the symmetric splitting procedure is reduced, but this cost will always be higher than the cost of the corresponding sequential procedure.
•	The computational cost of a weighted sequential splitting procedure is higher, by a factor approximately equal to two, than that of the corresponding sequential procedure. If two processors are available when the weighted sequential procedure is used and only one processor is applied with the sequential splitting procedure, then the computer time needed will be nearly the same for the two splitting procedures. This is due to the fact that the two tasks (9) and (10) can be carried out in parallel. However, if the code based on a sequential splitting procedure is parallelized and if the fact that two processors are available is utilized when the sequential splitting procedure is run, then this is no more true. It will be shown in Section 5 that an efficient parallel code can be developed in the case when the sequential splitting procedure is used. Of course, a similar conclusion can be drawn if 2p processors are available when the weighted sequential splitting procedure is used, while the number of processors available when the sequential splitting procedure is used is only p. If the parallel codes for the two splitting procedures are scalable, then approximately the same computing time will be used for the two splitting procedures.
Similar conclusions (as those for the relationship between weighted sequential splitting procedures and sequential splitting procedures) are valid for the relationship between weighted symmetric splitting procedures and symmetric splitting procedures. Some numerical results in support of these conclusions are given in Table 3 in the next section (Section 5).
4.2	Comparing the accuracy properties
•	All sequential splitting procedures are of first order of accuracy, i.e., of order O(). This means that the errors due to these splitting procedures are proportional to the splitting time-step . Thus, if the time-step  is sufficiently small and if the error due to the use of the splitting procedure chosen is dominant, then one should expect the errors caused by the splitting procedure to be reduced by a factor of two when the splitting time-step is reduced from   to /2.
•	Symmetric splitting procedures of second order of accuracy, O(2), can be constructed. This means that if the splitting time-step   is sufficiently small and if the error due to the use of the splitting procedure chosen is dominant, then one should expect the errors caused by the splitting procedure to be reduced by a factor of four when the splitting time-step is reduced from  to /2.
•	Weighted sequential and symmetric splitting procedures of order higher than two can be constructed.
The discussion in this sub-section shows that if the accuracy is the most important issue, then symmetric splitting procedures and even weighted (either sequential or symmetric) splitting procedures are to be selected.
Combining the conclusions in this sub-section with those in the previous sub-section makes it clear that the selection of a splitting procedure is a trade off between two requirements (the requirement for reducing the computing time and the requirement for achieving sufficiently accurate results). Unfortunately, these two requirements are working in opposite directions. This is why in general it is very difficult to find a choice which satisfies sufficiently well both requirements.
If some assumptions are satisfied, then the errors due to the splitting procedures can vanish. Such results were proved by Lanser and Verwer (1999) by using the concept of "L-commutativity of two operators". This property, which is rather important from a practical point of view, will be studied in further works. 
5.	real data experiments with DEM – comparison with observations
One-year experiments with real data were carried out by using DEM and the results were compared with the observations, obtained from a number of measurement stations scattered throughout Europe. The results for two major pollutants, NO2 and SO2, are presented in Table 1 and Table 2, respectively.
 The scalability and the performance of DEM on a Sunfire 6800 parallel machine for the above one-year experiments are shown in Table 3. The chemistry-deposition stage is perfectly parallelized as it is independently executed for each grid point. Moreover, it shows superlinear speed-up because of the cache utilization of the processors.

NO2 concentrations in air for 1998:   DEM results in comparison with measured values	
Splitting             method:	Jan	Feb	Mar	Apr	May	June	July	Aug	Sept	Oct	Nov	Dec	year	
	NO2  [mol/l]  (averaged over 35 places with measurement stations)	
1. Sequential	3.20	4.96	1.83	1.76	2.05	2.33	1.86	1.99	2.42	2.38	3.97	4.10	2.740	
2. Symmetric	2.35	3.83	1.38	1.39	1.66	1.95	1.56	1.81	2.11	1.85	3.04	3.10	2.169	
3. Weighted seq.	3.21	4.97	1.84	1.76	2.05	2.34	1.86	1.99	2.42	2.39	3.97	4.11	2.742	
Measured value	2.79	3.24	2.18	1.98	1.81	1.58	1.57	1.70	1.84	1.96	3.10	3.33	2.258	
	DISCREPANCY   [%]	
1. Sequential	13
	35
	-19
	-12
	12
	32
	15
	15	24	18
	22
	19
	17
	
2. Symmetric	-19
	15
	-58
	-42
	-9
	19
	-1
	6
	13
	-6
	-2
	-8
	-4
	
3. Weighted seq.	13
	35
	-19
	-12
	12
	32
	16
	15
	24
	18
	22
	19
	17
	
	CORRELATION FACTOR	
1. Sequential	0.67
	0.61
	0.69
	0.59
	0.64
	0.56
	0.50
	0.59
	0.65
	0.57
	0.65
	0.70
	0.692
	
2. Symmetric	0.62
	0.59
	0.75
	0.49
	0.66
	0.47
	0.49
	0.57
	0.56
	0.53
	0.70
	0.71
	0.679
	
3. Weighted seq.	0.66
	0.61
	0.69
	0.59
	0.64
	0.56
	0.50
	0.59
	0.65
	0.57
	0.65
	0.69
	0.690
	
Table 1: The concentrations of NO2 in the air for 1998, calculated by DEM by using three different splitting schemes. The results are compared with the values, obtained from 35 permanent measurement stations in Europe at the corresponding grid points. The (96 x 96) grid is used in the discretization of the spatial domain.
SO2 concentrations in air for 1998:   DEM results in comparison with measured values	
Splitting             method:	Jan	Feb	Mar	Apr	May	June	July	Aug	Sept	Oct	Nov	Dec	year	
	SO2  [mol/l]  (averaged over 35 places with measurement stations)	
1. Sequential	2.70
	2.02
	1.49
	1.06
	1.13
	0.72
	0.63
	0.73
	0.85
	0.68
	1.59
	2.89
	1.375
	
2. Symmetric	2.49
	1.91
	1.41
	1.01
	1.12
	0.72
	0.62
	0.75
	0.86
	0.67
	1.58
	2.64
	1.317
	
3. Weighted seq.	2.72
	2.04
	1.51
	1.07
	1.14
	0.72
	0.64
	0.74
	0.88
	0.69
	1.60
	2.91
	1.387
	
Measured value	1.36
	1.41
	1.11
	0.80
	0.73
	0.51
	0.49
	0.56
	0.60
	0.58
	1.32
	1.41
	0.906
	
	DISCREPANCY   [%]	
1. Sequential	50	30	26	25	35	29	22	23	29	16	17	51	34	
2. Symmetric	46	26	22	21	35	30	20	25	30	15	17	47	31	
3. Weighted seq.	50	31	26	25	36	30	23	24	29	17	18	52	34	
	CORRELATION FACTOR	
1. Sequential	0.83
	0.61
	0.64
	0.60
	0.71
	0.72
	0.82
	0.63
	0.82
	0.75
	0.29
	0.39
	0.712
	
2. Symmetric	0.83
	0.69
	0.64
	0.64
	0.72
	0.70
	0.83
	0.71
	0.81
	0.81
	0.45
	0.52
	0.775
	
3. Weighted seq.	0.83
	0.61
	0.64
	0.60
	0.71
	0.72
	0.83
	0.63
	0.82
	0.75
	0.28
	0.39
	0.712
	
Table 2: The concentrations of SO2 in the air for 1998, calculated by DEM by using three different splitting schemes. The results are compared with the values, obtained from 35 permanent measurement stations in Europe at the corresponding grid points. The (96 x 96) grid is used in the discretization of the spatial domain.
	Time [sec.]  (Speed-up)  of 2-D version of DEM on Sunfire 6800	
Splitting method:	1 proc.	2 proc.	4 proc.	8 proc.	16 proc.	
1. Sequential	37893	18570  (2.04)  
	8406  (4.51)  	4588  (8.26)  	2362  (16.0) 	
2. Symmetric	45864	23116  (1.98)  
	11780  (3.89)  	5684  (8.07)  	2942  (15.6) 	
3. Weighted seq.	68746	33465  (2.05)  
	16227  (4.24)  	7700  (8.93)  	4127  (16.7) 	
Table 3: Time in seconds and speed-up (in brackets) of the three splitting schemes for a one-year experiment with the 2-D version of DEM on a (96 x 96 x 1) grid. 
The superlinear speed-up of the chemistry-deposition stage has stronger influence on the total speed-up in the sequential and weighted sequential schemes, where this stage takes more than 70% of the total time. In the symmetric scheme it is put in the middle and is surrounded by 2 advection-diffusion stages with a half time-step, so its relative weight decreases as well as the total speed-up.
6.	experiments on test problems with known solution
Three important cases are studied by these experiments:
•	The behaviour of the results in the case where only splitting errors appear in the computations.
•	The behaviour of the results in the case where the splitting errors are appearing together with errors due to the spatial discretization and the errors due to the time integration.
•	The behaviour of the results in the case where there are no splitting errors (i.e., all pairs of the involved operators are L-commuting). 
6.1	Information about the presentation of the results
Let  ci,j,n  ≈  c(xi,yj,tn) ; i{l,2,...,Nx}, j{l,2,...,Ny}, n{l,2,...,N}. The error ERRn at time-step n is evaluated by:
,      where , .
The global error ERRMAX (which is the maximal error over the whole time-integration interval) is evaluated by .
The two-dimensional case is treated above. It is clear that values of ERRMAX for three-dimensional problems can be calculated in a quite similar manner. The values of ERRMAX, rounded to the third digit, will be presented in Tables 4-7.
6.2	Only splitting errors appear in the computations
Consider the following problem:
 ,      t  (0,1] ,     c(x,y,0) = x+y ,            (15) 
with exact solution
c(x,y,t) = x+y+t  .  												(16)
Assume that a sequential splitting procedure is applied in the solution of (15) by using the following two simpler operators:
, 											(17)
.                                         			(18)
Linear finite elements are used in the semi-discretization of the first sub-problem (17). The resulting system of ODE’s is treated by a third-order predictor-corrector method. The second sub-problem (18) is solved exactly. Since the exact solution (16) of (15) is a linear function in  x, y  and  t,  this means that the numerical methods are not causing errors. Thus, the only error is caused by the splitting procedure.
Two important conclusions can be drawn from the results given in Tables 4 and 5:
•	As mentioned above, there are no errors caused by the numerical algorithms. Therefore, all errors (excepting the rounding errors, which are negligible, because the computations were carried out in double precision) are caused by the use of the splitting procedure. This implies that refining the grids should have no effects on the accuracy of the approximate solution. This statement is fully confirmed by the results shown in Table 4.
•	Since the sequential splitting procedure used is of first order, one should expect the errors to depend linearly on the splitting time-step. The results shown in Table 5 indicate that this is true (reducing the splitting time-step by a factor of ten is leading to a reduction of the errors of the numerical solution by a factor approximately equal to ten).

Grid
	Without splitting
	With splitting
	
8 x 8
	3.34E-14
	1.70E-4
	
16 x 16
	3.79E-14
	1.80E-4
	
32 x 32
	3.99E-14
	1.80E-4
	
96 x 96
	7.21E-14
	1.97E-4
	
288 x 288
	4.57E-14
	1.99E-4
	
480 x 480
	4.28E-14
	1.99E-4
	
Table 4: The values of ERRMAX obtained when (15) is solved without splitting and with splitting (the number of time-steps is 1500 for all grids, but the values of ERRMAX are rather insensitive to variations of the number of the time-steps when no splitting is used). The values of ERRMAX are close to the machine precision when no splitting is used. The values of ERRMAX are of order Ο(10-4) for all grids when the splitting procedure is applied.

Steps
	Step
	Without splitting
	With splitting
	
10
	 0.1
	4.35E-14
	2.42E-2
	
100
	 0.01
	4.36E-15
	2.39E-3
	
1000
	 0.001
	5.76E-14
	2.39E-4
	
10000
	 0.0001
	7.24E-13
	2.39E-5
	
Table 5: The values of ERRMAX obtained when (15) is solved by using different numbers of time-steps without splitting and with splitting. The 16 x 16 grid is used in the discretization of the spatial derivatives. The experiment shows that the values of ERRMAX do not depend on the number of time-steps when no splitting is used. In the case where splitting is applied, the values of ERRMAX are reduced when the number of steps is increased, because, as mentioned above, Δt = .
6.3	Splitting errors versus other numerical errors
The combined effect of splitting errors and other errors, caused by the numerical algorithms applied in the discretization of the spatial derivatives and by the numerical methods applied in the treatment of the systems of ODE’s arising after the semi-discretization, is  illustrated in this sub-section.
Consider the following three-dimensional problem:
 + y(cos z - sin z) + x(sin z + cos z) ,   t[0,1] ,        (19)
with an initial condition
c(x,y,z,0) = ez(x2+y2). 											(20)
The analytical solution of this problem is
c(x,y,z,t) = ez+t(x2 + y2) + t[ y(cos z − sin z) + x(sin z + cos z) ].          (21)
A sequential splitting procedure is used in the solution of this problem. This procedure is based on the sequential treatment of the following sub-problems:
 ,    								        (22)
  y(cos z − sin z) + x(sin z + cos z) ,	     				     (23)
 .     									      (24)
Results obtained with this splitting procedure for   = Δt  for different resolution grids and different values of   Δt  are given in Table 6.
The results in Table 6 indicate that the error due to the spatial discretization is very quickly becoming dominant (except the results for the finest grid). Further reductions of the time-step  Δt  do not lead to reduction of the total error, because the errors caused by the spatial discretization become dominant.

Time-step =Δt
	Grid size	
	16x16x16
	32x32x32
	64x64x64
	128x128x128
	
8.00E-3
	1.58E-2
	1.92E-2
	unstable
	unstable
	
4.00E-3
	7.99E-3
	9.03E-3
	unstable
	unstable
	
2.00E-3
	4.03E-3
	4.55E-3
	5.07E-3
	unstable
	
l.00E-3
	2.22E-3
	2.33E-3
	2.53E-3
	unstable
	
5.00E-4
	1.32E-3
	1.22E-3
	1.28E-3
	1.35E-3
	
2.5E-4
	l.OOE-3
	6.73E-4
	6.58E-4
	  6.74E-4 •
	
1.25E-4
	1.02E-3
	4.02E-4
	3.47E-4
	3.43E-4
	
6.25E-5
	1.03E-3
	2.96E-4
	1.91E-4
	1.77E-4
	
3.12E-5
	1.03E-3
	2.95E-4
	1.14E-4
	9.32E-5
	
1.56E-5
	1.03E-3
	2.94E-4
	8.11E-5
	5.15E-5
	
7.81E-6
	1.03E-3
	2.94E-4
	8.09E-5
	3.07E-5
	
3.91E-6
	1.03E-3
	2.94E-4
	8.08E-5
	2.15E-5
	
Table 6: Values of ERRMAX obtained when the problem defined by (19) and (20) is solved with the splitting procedure (22) – (24) by using different numbers of time-steps and different spatial grids.
The splitting time-step and the discretization time-step used in the integration of the system of ODE’s arising after the semi-discretization were, as stated above, equal when the results shown in Table 6 were calculated (i.e.,  = Δt). This does not allow us to distinguish between errors due to the use of different splitting steps    and errors caused by use of different discretization steps  Δt  in the solution of the system of ODE’s. That is why some experiments, in which  is not equal to Δt, are also needed and were carried out. More precisely, the number of splitting steps was kept constant and the number of time-steps per splitting step was varied. Some results for ERRMAX in dependence with  and Δt, obtained in such experiments, are given in Table 7 for the (161616) grid.

Δt		
	5.0E−3
	4.0E−3
	2.0E−3
	1.0E−3
	5.0E−4
	2.5E−4
	
8.00E-3
	1.59E-2
	-
	-
	-
	-
	-
	
4.00E-3
	1.59E-2
	7.99E-3
	-
	-
	-
	-
	
2.00E-3
	1.59E-2
	7.99E-3
	4.03E-3
	-
	-
	-
	
l.00E-3
	1.59E-2
	7.99E-3
	4.03E-3
	2.22E-3
	-
	-
	
5.00E-4
	1.59E-2
	7.99E-3
	4.03E-3
	2.22E-3
	1.32E-3
	-
	
2.50E-4
	1.59E-2
	7.99E-3
	4.03E-3
	2.22E-3
	1.33E-3
	1.00E-3
	
1.25E-4
	1.59E-2
	7.99E-3
	4.03E-3
	2.22E-3
	1.33E-3
	1.01E-3
	
6.25E-5
	1.59E-2
	7.99E-3
	4.03E-3
	2.22E-3
	1.32E-3
	1.01E-3
	
3.12E-5
	1.59E-2
	7.99E-3
	4.03E-3
	2.22E-3
	1.32E-3
	1.00E-3
	
1.56E-5
	1.59E-2
	7.99E-3
	4.03E-3
	2.22E-3
	1.32E-3
	1.01E-3
	
7.81E-6
	1.59E-2
	7.99E-3
	4.03E-3
	2.22E-3
	1.32E-3
	1.01E-3
	
3.91E-6
	1.59E-2
	7.99E-3
	4.03E-3
	2.22E-3
	1.32E-3
	1.01E-3
	
Table 7: Values of ERRMAX obtained when the problem defined by (19) and (20) is solved by using the splitting procedure (22) – (24) with different numbers of discretization time-steps Δt per splitting time-step .
It is seen from the results presented in Table 7 that if splitting time-step  is fixed, then it does not matter how many discretization time-steps of size  Δt  will be carried out per splitting time-step. In other words, the results shown in Table 6 indicate that the errors due to the splitting (with time-step ) are dominant in this example.
7.	Conclusions and plans for future research
>From the real data experiments with DEM by using different splitting schemes and comparing the results with data from the measurement stations (see Section 5) the following conclusions can be drawn:
•	In general, the results, obtained with the weighted sequential splitting, are quite similar to those of the sequential splitting. This could be an indication for a small commutator of the two operators.
•	For some pollutants (NO2, SO2) the symmetric splitting scheme (due to Marchuk and Strang) gives results, which are in most cases closer to the measurements than those, obtained with the other two splitting methods.  More experiments are needed in order to investigate the consistency of such behavior.
•	Although the weighted and the symmetric splitting schemes are of higher order of accuracy with respect to the splitting error, this cannot be clearly seen by the experimental results.  This is probably because of the influence of other sources of error as well as the uncertainty of the measurements.
The examples given in the previous section show that the splitting procedure in general causes splitting errors. Therefore, it is worthwhile to carry out work in two important directions:
•	attempt to avoid the use of splitting procedures when this is possible,
•	finding reliable and robust methods for evaluating the splitting errors.
Both problems are still open when general large-scale air pollution problems are to be handled. Success in the solution of any of these two problems is a very desirable, but also a very challenging task. It should be noted that even a partial solution of any of the two problems will lead to a very considerable progress in the numerical treatment of large-scale air pollution models and many other large-scale computational problems arising in science and engineering.
Better understanding of the interference of errors due to the splitting procedures and errors that are caused by other reasons (numerical methods, uncertainties of the input data, etc.) is highly desirable. Research efforts in this direction are needed.

Acknowledgement: This research was supported in part by grant I-901/99 from the Bulgarian NSF, grant ICA1-CT-2000-70016  Center of Excellence and by NATO grant NATO ARW  project ”Impact of future climate changes on pollution levels in Europe”.
8.	references
Bagrinovskii, K. A. and Godunov, S. K., 1957, Difference schemes for multidimensional problems. Dokl. Akad. Nauk USSR, 115, pp. 431-433.
Csomós, P., Faragó, I. and Havasi, Á., 2003, Weighted sequential splittings and their analysis. Comput. Math. Appl. (to appear).
Dimov, I., Faragó, I., Havasi, Á. and Zlatev, Z., 2001, L-commutativity of the operators in splitting methods for air pollution models. Annalles Univ. Sci. Budapest, Sec. Math., 44, pp. 127-148.
Dimov, I., Faragó, I. and Zlatev, Z., 1999, Commutativity of the operators in splitting methods for air pollution models, Report 04/99, Central Laboratory for Parallel Processing, Bulgarian Acad. Sci..
Dyakonov, E. G., 1962, Difference schemes with splitting operators for multidimensional stationary problems. Journal of Computational Mathematics and Mathematical Physics, 2(4), pp. 57-79 (in Russian).
Faragó, I. and Havasi, Á., 2002, The mathematical background of operator splitting and the effect of non-commutativity. In: Large Scale Scientific Computations III (Margenov, S., Yalamov P. and Wasniewski, J. eds.), Springer, Berlin, pp. 264-271.
Faragó, I. and Havasi, Á., 2003, On the convergence and local splitting error of different splitting schemes. Progress in Computational Fluid Dynamics, (to appear).
Faragó, I., Havasi, Á. and Korotov, S., 2002, Mathematical analysis of the weighted splitting. Reports of the Department of Mathematical Information Technology, University of Jyväskylä, No. B/12.
Hunsdorfer, W. and Verwer, J. G., 2003, Numerical solution of time-dependent advection-diffusion-reaction equations. Springer, Berlin.
Lanser, D. and Verwer, J. G., 1999, Analysis of operators splitting in advection-diffusion-reaction problems in air pollution modelling. Journal of Computational and Applied Mathematics, 111, pp. 201-216.
Marchuk, G. I., 1968, Some application of splitting-up methods to the solution of mathematical physics problems. Applik. Mat., 13(2).
Marchuk, G. I., 1980, Methods of computational mathematics. Nauka, Moscow, 1980.
Marchuk, G. I., 1982, Mathematical Modeling for the Problem of the Environment. Nauka, Moscow (in Russian).
Marchuk, G. I., 1986, Mathematical modelling for the problem of the environment. Studies in Mathematics and Applications, No. 16, North-Holland, Amsterdam.
Marchuk, G. I., 1994, Applications of Adjoint Equations to Problems of Global Change and Environmental Protection. In: Proceedings of HERMIS'94, Vol. 1, (E.A. Lipitakis, ed.), Hellenic Mathematical Society, Athens, pp. 1-20.
McRae, G. J., Goodin, W. R. and Seinfeld, J. H., 1982, Numerical solution of the atmospheric diffusion equations for chemically reacting flows. Journal of Computational Physics, 45, pp. 1-42.
Penenko, V. V. and Obraztsov, N. N., 1976, A variational initialization model for the fields of meteorological elements. Soviet Meteorol. Hydrol., 11, pp. 1-11.
Strang, G., 1968, On the construction and comparison of difference schemes. SIAM J. Numer. Anal., 5, pp. 505-517.
Tikhonov, A. N. and Samarski, A. A., 1977, Equations of the Mathematical Physics. Nauka, Moscow.
Yanenko, N. N., 1962, On convergence of the splitting method for heat equation with variable coefficients. Journal of Computational Mathematics and Mathematical Physics, 2(5), pp. 933-937 (in Russian).
Zlatev, Z., 1995, Computer treatment of large air pollution models, Kluwer.