Determination of atmospheric aerosol size distribution functions from screen diffusion battery data: Mathematical aspects

J. Aerosol SCL Vol. 22. No. 3, pp. 373.388. Printed in Great Britain.

1991

0021-8502/91 s3.cQ+o.o0 Q 1991 Pcrgaman Press plc.

DETERMINATION OF ATMOSPHERIC AEROSOL SIZE DISTRIBUTION FUNCTIONS FROM SCREEN DIFFUSION BATTERY DATA: MATHEMATICAL ASPECTS V. S. BASHUROVA, K. P. KOIJTZENOGIL, A. Y. PUSEP and N. V. SHOKHIREV Institute

of Chemical (Received

Kinetics 9 October

and Combustion,

Novosibirsk,

1989; and infinalform

630090,

14 August

U.S.S.R

1990)

Abstract-This paper considers the mathematical aspects of the reconstruction of size distribution functions of fine atmospheric aerosol from screen-type diffusion battery (SDB) measurements. It is shown that the problem of reconstruction may be successfully solved using Tikhonov’s regularization method of the second order. A series of experimental data is considered, and the spectra obtained are classified. The bimodal size distribution is proved to exist in the highly dispersed fraction of atmospheric aerosol. The reliability of the solutions is analysed; it is found that, for the existing characteristics of the SDB, it is impossible to determine more than two modes. Recommendations on the planning of SDB experiments are given.

1. INTRODUCTION

Most of atmospheric processes are associated with the appearance, transformation and transport of aerosols. The dynamics of atmospheric circulation and structure development are determined by the properties of atmospheric aerosols which, in turn, depend on their microphysical characteristics, such as concentration, dispersion and chemical composition of aerosol particles (Koutzenogii, 1987). A reliable measurement of atmospheric aerosols parameters becomes particularly important as regards the problems of atmospheric pollution (Kondratiyev et al., 1983; Ivlev and Andreev, 1986; Rogers, 1976; Schreier, 1982). We are going to consider one of the above characteristics, i.e. the particle size distribution. In terms of the well-known Whitby model (1978), atmospheric aerosol consist of three modes: (i) nuclei (ultrafine) mode (particle size less than 0.1 pm); (ii) accumulation

mode (particle size 0.1-l pm); and

(iii) coarse-dispersed fraction (particle size greater than 1 pm). The three-mode distribution proposed by Whitby (1978) is related to the combined atmospheric physico-chemical processes. Thus, the nuclei mode (i.e. particle radius less than 0.1 pm) is thought to appear due to the ‘gas-to-dispersed-phase’ conversion resulting from photochemical reactions. The accumulated mode arises from the coagulation aging of the nuclei mode particles. The coarse-dispersed fraction is mostly defined by the processes of soil and water erosion. This interpretation is the simplest one concerning many complex and interrelated processes affecting aerosol size distribution. Other approaches serving as a theoretical foundation for the structure of atmospheric aerosol size distribution as a whole, or as its parts, are also available (Friedlander, 1977; Smirnov, 1977a). At the same time, detailed measurement of size distribution provides important information about the formation and space-time evolution of aerosol in real atmosphere. Thus, Davies (1974) proposed the expansion of the whole size range into a set of superimposed log-normal distributions. A number of papers (Lushnikov and Smirnov, 1975; Smirnov 1977a, b, 1987; Smirnov and Sergeev, 1980) report a theoretical basis for the power approximation of separate parts of the atmospheric aerosol size range caused by the processes of condensation, vaporization and coagulation. These indicate how the problem depends on the reliability of the measurements of the atmospheric aerosol size distribution. 373

V. S. BASHUROVA

374

ef al.

The overall atmospheric particle size range spans five orders of magnitude (from 10e3 to IO2 pm) and no common method is known to measure it. Therefore, various techniques and apparatus are applied to investigate different aerosol size sub-ranges. As studies of the laws of photochemical formation of smogs are very important, the problem of detailed measurements of the nuclei mode size distribution is very relevant. As a rule, most theoretical models lead to the one-mode (usually log-normal-approximated) size distribution (McMurry and Friedlander, 1978, 1979; McMurry, 1983). However, some numerical models (Sutugin, 1975) and experiments (Whitby, 1978; Michael, 1986; Bashurova et al., 1987; Jaenicke, 1988) indicate the possibility of multi-mode distribution in the nuclei mode. The present paper is concerned with a particular problem in the determination of the nuclei mode spectrum. The diffusion battery (DB) technique has been shown to be one of the most efficient means to study this part of the size range (Koutzenogii, 1986). The determination of the aerosol size distribution using diffusion batteries is based on the size dependence of aerosol penetration through such devices. Information about the particle size is obtained by measuring the particle flux differences through the battery. Hence, this method for measuring particle size distribution is the so-called indirect measurement. To determine the required dependence it is necessary to use special mathematical procedures (Delves and Walsh, 1974; Tikhonov and Arsenin, 1979; Vasilenko, 1979). The reliability of the information obtained depends not only on the quality of measurements but also on the method of data processing. These questions are widely discussed in the literature and different approaches have been developed to the solution of inversion problems: fitting methods (Helsper et al., 1982), regularization methods (Phillips, 1962; Tikhonov, 1963; Twomey, 1963, 1965; Crump and Seinfeld, 1982), nonlinear programming using physical constraints (Cooper and Spielman, 1976) and nonlinear iterative methods (Twomey, 1975; Cheng and Yeh, 1984; Maher and Laird, 1985). Quite opposite opinions are proposed concerning the applicability of these methods to determination of particle size spectra. Twomey (1975) prefers the iterative method, although Crump and Seinfeld (1982) consider it to be inadequate and propose a version of a smoothing algorithm. In the present paper we also use a version of the regularization method (Tikhonov and Arsenin, 1979) with a simple recipe of the choice of regularization parameter. We have successfully employed this method earlier (Pusep and Shokhirev, 1986). The most important methodical novelty is the analysis performed of the reliability of size distribution reconstruction. This analysis, based on singular value decomposition (SVD) of the integral operator does not depend on the method of reconstruction. In our opinion, this analysis should go in line with the solution of an inversion problem proper. The second section contains a successive mathematical formulation of the problem, and the third one presents the method of its solution and analysis of the DB inversion problems, which we developed earlier (Pusep and Shokhirev, 1984, 1986; Shokhirev et al., 1985). The fourth section exemplifies the atmospheric size distribution obtained by this technique. The fifth section considers the reliability of the reconstruction of different types of spectra. In Section 6 the main results are summarized and an experimental design is considered.

2. MATHEMATICAL

MODEL INVERSE

OF DB FORMULATION PROBLEM

FOR

THE

The screen diffusion battery is a tube filled with screens through which the aerosol is drawn. The number of particles left in the flow is measured after each screen. The determination of the aerosol particle size is based on the dependence of the aerosol precipitation rate (by diffusion) on the particle size. To this end a penetration coefficient is measured, i.e. the ratio between the number of particles penetrating through a given number of screens and the number of particles in the initial flow. To use these measurements for determining particle size, we must know the aerosol penetration through the screen as a

Atmospheric

aerosol

size distribution

functions

375

function of particle size. This dependence is defined by the diffusion battery design and by the experimental conditions, and may be tabulated using special calibration experiments. In practice, however, appropriate analytical dependencies are more convenient. A well-known dependence (Cheng and Yeh, 1980; Kirsch and Fuchs, 1968) has been used for monodisperse aerosol (particles of radius r) penetration through n screens of DB: P(n)=exp{

-C.n.P3(r)},

(1)

where C is the constant dependent on the flow rate and screen parameters (screen step and fibre diameter) and D(r) is the diffusion coefficient. For the TSI DB Cheng et al. (1980) used instead of equation (1) a more complicated expression, which included an interception term. We have evaluated the influence of the interception. Our estimations show, that for particle radii less than 0.1 pm, the contribution of direct interception is not more than O.l%, and the correction term for diffusion and interception is not more than 3%. An essential difference of these data from that for TSI DB is mostly conditioned by a large diameter of fibre: 72 pm in our case and 20 pm for TSI DB. It is quite clear how to use the DB for measuring the aerosol particle size. The problem is solved by measuring the exponential decay in equation (1) followed by the determination of the diffusion coefficient from the slope of the straight line In P(n), and hence, the dependence D(r) being known, the particle radius. For D(r), we have employed the simple power approximation formula suggested by Mavliev and Ankilov (1985). Later it might be assumed that: D(r)=B,/r+B,/r2,

(2)

where B, = 8.53 x 10m9, B, = 1.27 x 10m6 (r is meaasured in nm, and D in m2 s-l). Note, that Cheng et al. (1980) used the Einstein-Stokes formula with the Cunningham slip correction. The values of D(r) calculated according to these formulae differ not more than by 4.5%, the particle radii being l-100 nm. For polydispersed aerosol, the problem is more complex. This is described in terms of the particle size distribution function. If we make x = log (r), the penetration of the polydispersed aerosol through y screens is unambiguously determined by the form of the distribution function: CL K(Y, x)f(x)dx=p(y)> cGy
CD”’

(r), r = 10”.

(4)

(5)

Note that it is convenient to draw the data to a standard DB by redefinition of C in equation (5). The DB with C = 25,5 11 s2’3 m-4’3 is taken as a standard. Hence, y in equations (3) and (4) is the adjusted number of screens calculated by the formula: y=C,nlC,

(6)

where C, and n are a constant and the number of screens for a real battery, respectively. Note that the adjusted number of screens, y, is not an integer as a rule. Equation (3) may be considered as the equation whose solution, according to the measurement dependence p(y), is used to determine the particle size distribution function f(x). Equation (3) is a linear integral equation of the first kind. The problem of its solution belongs to the class of incorrect inverse problems (Vasilenko, 1979). Common to these problems is the instability of their solution to the measurement error in the right-hand part of p( y). Another feature is that thef(x) form may be determined by using both the input data [i.e. the p(y) form] and the solution procedure. In this case, dependence on the choice of the solution method exists even when the method provides the stability to the measurement error. Any f(x) function substituted into equation (3) reproducing the p(y) dependence

V. S. BASHUROVA

376

et al.

within experimental accuracy may be taken as the solution for the inverse problem. functions differ in the details whose contributions to the signal are less than the ment error. The solutions differing in such details may be obtained by various Thus, a special analysis is preferable for separating a ‘true’ part of the solution artefacts introduced by processing.

3. SOLUTION

OF DB INVERSION

All these measuremethods. from the

PROBLEM

One of the most universal methods for solving the inversion problems is the Tikhonov regularization method (Vasilenko, 1979). The method is as follows. The initial incorrect problem in equation (3) may be written as a search for the norm minimum IIPW- jW,x)fb)dxII +min is stabilized. That is, instead of equation (7), the problem (Philips, Arsenin, 1979) IIp(y)--

jK(y,

x)f(x)dxl12

+rSCfl+min

1962; Tikhonov

(7) and (8)

is solved where S [f] is the stabilizer, and the positive value y is the regularization parameter. We have employed the second-order stabilizer: S[f]= which provides the smoothness of the reconstruction function f(x). In order to construct the algorithm of the solution we must first determine the norm of equation (7) and introduce the interval [a, b] within which the solution should lie. When measuring a set of the values pm = p(y,), the norm in equation (8) is written as IMP-

*WY, 4fb)dxl12 s0

= Cum tp,n

K,,Cfl~2~

(10)

where kCfl=

*K(y,AfWx, s (I

(11)

and u, are the weight coefficients allowing for various values of the errors in the set of measurements for pm. When seeking the solution at the points xj=a+Ax(j-l),

Ax=(b-a)/(N-l),j=l

. . . N.

(12)

Equation (8) reduces to the solution of the system of linear algebraic equations (Tikhonov and Arsenin, 1979; Vasilenko, 1979): ATA+&Q2

i= ATi, >

where AT is the transposed matrix of A. Amj = Ju,K(~rn,

Xj)Wj,

(14)

and matrix Q is the operator of the second-order finite differences, the elements of which are zero except for:

Qj,j = -2s Qj,j+ 1= Qj.j- 1= 1,

(15)

h, = & pm, and wj are the weight coefficients of the quadrature formula which substitutes the integration over x in the interval [a, b]. The possibility of applying the regularization method to solve the concrete inversion problem is highly dependent on the choice of the optimal value of the regularization

Atmospheric aerosol size distribution functions

371

parameter y, determining the solution smoothing degree. We have chosen the optimal value from the condition of the discrepancy minimum A=)~-A~lM-1’2,

(16)

where M is the number of measurements. To calculate the discrepancy the negative values of the components of the solution fj = f(xj) were substituted by zeros. This way of choosing y proved rather efficient in solving the DB inverse problem. When using the regularization method the solution depends on the interval [a, b] which includes the points of equation (12). In our procedure this interval is changed. First of all, the problem is solved within a certain initial interval. In this case the regularization parameter is taken according to the discrepancy minimum in equation (16). Then the interval [a, b] is redetermined. If, for example, f(b) does not vanish, the interval shauld be extended to the right. If f(x) vanishes near the right-hand boundary the interval should be narrowed from the right. With a new interval y is again chosen by the discrepancy minimum of equation (16). A series of such steps is sure to lead to the solution for f(x) which tends to zero (vanishes) at the boundary of the [a, b] interval. This solution is taken as optimal. To introduce correctly the weight coefficients to equation (lo), the mean errors of all the measurements must be known. When (E:) is the mean square of the measurement error in the penetration p,, it is assumed that u, = l/(&i).

(17)

In this case the mean-square discrepancy of equation (16) is of the form:

A={~~~~~“,Ap.=K,Cfl-p..

(18)

Hence, the value A>< 1 indicates whether it was possible to keep within the assumed mean range of errors. For AS 1 we believe that one of the solutions of the inverse problems was found but that it is not a unique one. Thus, the application of the procedure for solving the inversion problem in equation (3) seems to be successful if one or more solutions are possible within the error imit. It often happens that the interval in which the solution is found is beyond the boundaries of the socalled ‘reliable reconstruction range’ (RRR). This interval may usually be estimated for each concrete inversion problem, using some simple considerations. In our problem, the kernel of the integral equation, K(y, x), has (as a function x) the characteristic form of ‘smoothed steps’. A series of such functions is presented in Fig. 1. Based on their forms, we can assert that the measurement for a given number of y, screens contains information only about the range of x variations in which function K(y,, x) ‘cuts’ the distribution functionf(x). Thus, the particle size distribution in which a correct reconstruction of thef(x) form is possible

Fig. 1. The kernel of the integral equation of SDB-problem as a function x=log(r) corresponding y,,,values are denoted near each curve.

(r in nm). The

378

V. S. BASHUROVA

et 01.

may be estimated from the form of the function K(y,, x) for a concrete set of measurements of y,. Roughly speaking, particles whose sizes are beyond this range either precipitate on the first screen (small particles) or penetrate through all the screens without deposition (large particles). The above facts show that the RRR-value first of all depends on the measurement range, i.e. on the maximum and minimum number of screens in the measurement set. For this case, the measurement range is implied to be filled uniformly, as in Fig. 1. Note that, in general, the RRR magnitude is slightly dependent on the measurement error as well. A routine procedure is available that allows RRR to be estimated for any inversion problem of the type s’hown in equation (3). This procedure is based on the formalism of the singular value decomposition (SVD) of the integral operator K ( y, x) kernel (Pusep and Shokhirev, 1984). The SVD is the representation of the kernel of equation (3) in the form (Delves and Walsh, 1974) K(~,x)=Ck,z,(~)u,(x),

(19) n with the singular functions z,(y) and u,(x) forming the orthonormal (in the general case with weight functions) sets within the ranges of the y and x variations, respectively. The singular numbers may always be considered as the ordered ones k,>,k,>, Substituting

. . .>O.

m

equation (19) into (3), we derive P(Y)=cs”z”(Y). n

(21)

Thus, the functions z,(y) are the natural basis for expanding the measured function; p,(y) may always be represented as equation (21) (Pusep and Shokhirev, 1984). A second set of functions (0,) is the natural basis for representing the required function (22)

f(x)=pdL(x), n

where, for simplicity, the weight function is taken as unity. The expansion coefficients in equations (21) and (22) are interrelated by t, = s,/k,.

(23)

It is noteworthy that, unlike the expansion of equation (21), which is always valid, equation (22) holds true assuming the completeness of the basis {on}. If {u,} is incomplete (i.e. in an incomplete experiment), only that part of j(x) may be reconstructed which can be represented as equation (22), since the orthogonal to all u, (x) parts make zero contribution to the measured dependence p(y). In practice, while reconstructingf(x) from the erroneous measurements pm= p( y,) + E,, the basis {u”} is always incomplete. Indeed, the effect of the increasing noise results from equation (23) is due to the decrease of the numbers in equation (20) (Vasilenko, 1979). To avoid the predominance of the noise component over the basic one, the expansion in equation (22) should include only a few first terms. Hence, knowing the form of functions u,(x) corresponding to the maximal k,, the question of.whatf(x) may be reconstructed using the experimental data of equation (3) may be solved. An earlier paper (Pusep and Shokhirev, 1986) reported the calculation of the functions u,(x). Each has the number of zeros equal to (n- 1). Since the basic functions are liable to reflect the details of formf(x) only in the part of x variations where they oscillate, RRR is naturally defined as the region involving the zeros of functions u,(x). The region of zeros is somewhat broadened with an increasing number of basic functions, and therefore RRR is dependent on the measurement error [more terms may be included in equation (22) with a decreasing number of measurement error]. The comparison between the RRR estimate obtained by SVD and that based on the form of the functions of the K( y, x) kernel (see Fig. 1)

Atmospheric

aerosol

size distribution

functions

319

showed the latter to be fairly good. A satisfactory RRR estimate is gained assuming each measurement of y, to bear the information about thef(x) form in the range x(l) 5 x sxt2), where x(l) and x(‘) are determined from the equations: K(y,,

x”‘)=Li,

i= 1, 2.

(24)

The L, and L, values are dependent on the experimental accuracy. This dependence is, however, rather weak and, for the typical error of the measurements of penetration being no less than a few per cent, the universal values L, = 0.4 and L2 = 0.6 may be employed. Let us consider an example of the RRR determination and demonstrate the influence of the RRR value on the quality of thef(x) reconstruction. Our measurements were modelled at the points log,y,

= -2+0.5

m, m= 1 . . . 19,

(25)

remembering that the adjust number of screens y is not necessarily an integer. When determining RRR, we calculate its lower limit from K( y,, x”‘)=O.4 and the upper one from K(Y,~, xt2)) = 0.6 which yields x(r)=0 and x (2)= 2 5 A bimodal distribution (67% of particles is a log-normal distribution with rs,,= 6 nm, uB= 1.85; 33% of particles with r5,, = 28 nm, cB = 1.5) that is sure to keep within RRR was taken as an example off(x). The latter was reconstructed by the penetration set pm = p( y,)( 1 + 6,) where p( y,) are the ‘true’ values, and 6, are the values determined using the generator of random numbers that model the error distributed uniformly within the range f2%. The model measurements were completed with the normalization condition, i.e. the additional point p (0) = 1 (thus, a total number of points M=20). The present example employs a wittingly high accuracy and a great number of points, in order to elude the influence of measurement error and of insufficiency of the measurement number on the quality of reconstruction. The results of the reconstruction are given in Fig. 2a. The reconstruction quality is considered very good. Now let us see what happens when the distribution in question is beyond the RRR limits. To this end the measurements were modelled at the points log2y,=2.75+0.25m,m=1..

. . 19,

(a) 10

0.5

s

0

u 2 2

1.0

.E 5 DL x a

0.5

:fi

(b)

0 1.0

05

0

Particle Fig. 2. The dependence measurement interval;

of solution (b) 6=2%,

radius (nm)

quality on accuracy and measurement interval. (a) 6 = 2%. a full reduced measurement interval; (c) 6= IO%, full measurement interval.

(26)

V. S. BASHUR~VA

380

et al.

with the same ‘noise’ level. In this case the lower RRR limit x’“=O.9, and hencef(x), fails to keep within RRR. It is readily observed (see Fig. 2b) that the reconstruction quality of the nuclei mode is much worse than in the first case. Thus, to solve in terms of a given set of measurements, it is necessary to estimate RRR. The reliability of the information may be estimated comparing RRR and the interval including the reconstructed distribution. The diffusion battery resolution is most likely to depend on the experimental accuracy. A detailed analysis of the influence of the experimental accuracy on the resolution is also possible, based on the SVD formalism. Assuming a set of singular functions {on} to be the solution basis, the resolution is believed to be determined by the period of singular function oscillations. As mentioned above, with increasing experimental accuracy the number of functions u, that may be left in the expansion in equation (22) increases, and these with large numbers have more oscillations. Hence, expansion (22) acquires additional functions with a smaller period that is sure to increase the resolution. Note that the increase in number n causes just a slight broadening of the region involving the zeros of V, (x). Consequently, the decrease of the experimental error has practically no effect on RRR and leads to the increase of the resolution. Let us demonstrate the effect of the experimental accuracy on the reconstruction quality by using the above bimodal distribution. Figure 2c depicts the f(x) distribution reconstructed in terms of the model data, with errors 6, distributed uniformly within the range f 10%. The comparison between this solution and that derived with error f2% (see Fig. 2a) shows that the increase in the measurement error reduces substantially the reconstruction quality. The model experiments with the error proportional to square root of p(y,), yield an analogous conclusion. The previous example demonstrates that the bimodal distributions may be reconstructed by solving the inversion problem. It is hardly possible, however, in our opinion, to distinguish between the three modes. We have performed a numerical experiment by using the distribution being the sum of three log-normal ones with rSOequal to 20 nm, 90 nm and 250 nm and with the same gB = 1.4. The contribution of each distribution was, respectively, 56%, 30% and 14%. Note that the choice of the contributions was determined by the fact that the ‘small’ modes are always measured against the background of the ‘big’ ones. Therefore, they are difficult to reconstruct. Hence, the decrease in the fraction of big particles should improve the reconstruction quality. The measurements were modelled at points [see equation (25)] without the ‘noise’ added [note that the calcualtion error for p( y,) is no more than 0.1 %] and were supplied with the normalization condition p(O) = 1. The reconstruction results are given in Fig. 3. To our mind, this example demonstrates the impossibility of distinguishing between the three modes. Indeed, the three-mode distribution in reconstruction is substituted by a bimodal one. In this case the difference in the penetrations calculated from ‘true’ and reconstructed distributions is no more than l%, which is undoubtedly less than the usual measurement error.

4. PROCESSING

OF EXPERIMENTAL

DATA

We have used experimental data obtained by Koutzenogii et al. (1986) in Bulgaria. Note that particle concentration in these experiments was in the range of 103-lo5 cmd3. When processing the experimental data, all the P,,, points were assumed to have the same relative error, i.e. the weight coefficients u, a pm, The quality of experimental data fitting was estimated from the value of the mean-square relative discrepancy 6, thus ii2 = & 1 (AP,,,/P,)~. m

(27)

Each set consisted of six points and was completed with the normalization condition p(0) = 1 (i.e. the total number of points M =7). The total number of particles at the DB entry was measured at the beginning and at the end of each experiment. The experiments were analysed for which the change in the number

Atmospheric

aerosol

size

distribution

functions

381

15 r

0

1

10

Particle Fig.

3. Reconstruction

of

three-mode

distribution: constructed

radius

100

t nm)

solid function.

line-initial

function;

dashed

line-re-

of entering flows was no more than 13%. We have processed dozens of sets. The size distribution functions obtained for atmospheric aerosols may be divided into four types. The first one involves the unimodal symmetric distribution, exemplified in Fig. 4. Such distributions are successfully fitted by the log-normal function. Indeed, in the above example, 6 = 5.5% for the solution derived by the regularization method and 6 = 7.2% for the log-normal distribution. The integral characteristics of these distributions are rather close: the mean particle radius (r) =42.0 nm for the solution obtained in terms of the regularization method and (r)=41.6 nm for the log-normal one, the mean-mass size (r3)li3 is 59.7 nm and 56.9 nm, respectively. The second type involves the unimodal asymmetric distribution. Using a number of experimental data sets, we have obtained several types of asymmetric distributions (Figs 5

Particle Fig.

4. One-modal

symmetric distribution: line-log-normal

solid distribution;

radius

trim)

line-solution rsc = 35.6

by regularization nm, uI = 1.X.

method;

dashed

V. S. BASHUROVAet

382

0

al.

I 10

1

I 100

Particle radius (nm) Fig. 5. Distribution with straightened left wing, 6 = 3.6%.

1.0 -

0

1

10

100

Particle radius (nm) Fig. 6. Distribution with a more sloping right wing, 6=8.4%.

1.0 -

0

1

10

100

Particle radius Inm) Fig. 7. Distribution

with shoulder to the left: solid line-solution by regularization 6 = 8.7%; dashed line-expansion into two modes, 6 =9.09%.

method,

and 6) that are difficult to fit by the log-normal ones (in comparison with the previous distributions). Table 1 lists the discrepancy values for the solution in terms of the regularization method and for the log-normal fittings. The third type includes the ‘shouldered’ distributions. The distributions were obtained with the ‘shoulders’ in both regions (Figs 7 and 8). In the case of the pronounced ‘shoulder’ (see Fig. 7), a bimodal distribution may be assumed, and one may try to distinguish between these modes.

Atmospheric aerosol size distribution functions Table 1. Values of discrepancy (per cent) for log-normal regularization solutions

383

fitting and

(a) Log-normal fitting Y.W

Fig. 5

10

12.4 1.1

30

50 70 115

160 5

Fig. 6

Fig. I

2.9 6.2 -32.3 -25.4 -51.1 -61.3 34.03

-20.9 -19.4 -31.0 -34.2 21.0

12.8 2.4 -26.6

-5.7 12.8 -18.4 14.2

Fig. 8 30.4 39.0 26.2 32.5 -15.0 1.6 25.1

Fig. 9 -12.8 8.7 -3.6 32.7 0.9 -32.4 18.4

Fig. 10 16.4 -4.6 -16.3 -1.5 -27.6 -6.8 14.0

(b) Regularization solution Ylll 0 10

Fig. 5

Fig. 6

-0.1

30 50 70 115

0.4 3.6 - 8.4 2.3 -0.3

160 6

5.2

0.9 -5.1 10.7 -16.8 1.5 -2.1 0.9 8.4

1.3

Fig. I -0.7 4.0 2.6 -21.4 5.8 4.0 -0.3 8.7

Fig. 8 0.6 2.8 1.2 -2.2 10.2 - 19.0 6.5 8.6

Fig. 9 9.6 16.4 -5.4 -16.6 19.2 6.7 -11.2 13.2

-

Fig. 10 -1.1 4.8 -3.2 -9.2 10.9 - 16.0 1.6 8.9

(c) Bimodal parameters Fig. 7 rg

nm

6.17 2.14 17% 34.67

6l1)s

41

r!$ o12) s

Fig. 9 23.17 nm

4.42 nm 1.76 31%

1.64 89%

nm

1.78 83% 9.09%

1.0

Fig. 10

167.88 nm 1.33 11%

67.61 nm 1.86 63% 10.45%

15.00%

-

10

Particle

100

radius

(nm)

Fig. 8. Distribution with shoulder to the righ, Jbtained by regularization method, 6 = 8.6%.

The fourth type involves the bimodal distribution with separate modes, the two types of which are shown in Figs 9 and 10. The corresponding discrepancy values for the distributions of the third and fourth types are listed in Table 1.

V. S. BASHUROVA

384

Particle Fig.

9. Bimodal

distribution 6= 13.2%;

Fig.

10. Bimodal

distribution 6= 8.9%;

radius

with separate modes: dashed line-expansion

Particle with dashed

separate line-xpansion

et cl.

trim)

solid into

radius modes:

solid into

line-solution two modes,

by regularization 6= 15%.

method,

by regularization 6 = 10.45%.

method,

trim) line-solution two modes,

5. DISCUSSION In inversion problems of the type described by equation (3), the reconstruction precision is determined by two factors: the experimental accuracy and the width of measurement interval (in the DB problem this interval is the range of the variations of the adjust screen number). In all the sets of experimental data, the number of screens varied from 10 up to 160 (see Table l), the screen parameters coinciding withthe standard ones. RRR was estimated in terms of the methods presented earlier to be in the range from 14 to 300 nm. However, the physical model describing DB is known to fail for a particle radius exceeding 100 nm. The upper RRR limit cannot exceed this physical threshold. As mentioned above, for the case of symmetric unimodel distributions, the solutions derived by fitting and regularization methods are in fair agreement (see Section 4). On comparison, it should be taken into account that, unlike the regularization method, the reconstruction ones assuming a priori a definite j(x) form may also be valid if the true distribution is far beyond the RRR limits. For example, Mavliev et al. (1984) recommend the measurements in the range 0.001 < pm < 0.1. For this case RRR covers only the right-hand

Atmospheric

wing of the distribution

aerosol

size

distribution

function. If the distribution

functions

385

form coincides with that assumed a

priori (log-normal), the experimental evidence is sufficient for determining two parameters, rsO and erg.When there is no such coincidence, the fitting method may either lead to error or

make it impossible to keep within the error limit during fitting. Consider now the question of the reliability of the bimodal solutions. The above examples indicate (see Figs 7 and 8) that the ‘small’ mode with a small weight is always to the left from RRR. Thus, the question arises of how reliable the conclusion is about the bimodality of such a distribution. To answer this question we have changed the extent of the solution smoothness to transform the bimodal into unimodal distribution (Fig. 11). The solution whose mean-square discrepancy (Fig. 4) is minimal has two separate modes (Fig. 1la). The other two, the smoother ones (Fig. 1 lb and c), have a less pronounced distribution bimodality that, in our opinion, is noticeable in the most smooth solution (Fig. 1 lc). The dependence of the mean-square discrepancy on the regularization parameter always has a minimum (Vasilenko, 1979). This dependence for the solution analysed is given in Fig. 12. The minimum exists because at small y the solution oscillates, and at large y it becomes so smoothed that the true distribution form is concealed. Therefore, the solution with the minimum 6 is always taken as a true one. From this viewpoint the conclusion on the bimodality of the obtained distribution seems fairly reliable. Note, however, that an attempt to fit the measurements by the log-normal distribution leads to a 1.5-fold increase in the mean-square discrepancy (see Fig. 1 lc). Of interest is an essential difference between the highly smoothed solution and the optimal log-normal distribution. Nevertheless, when the regularization method yields a unimodal solution, such a noticeable difference from the optimal log-normal distributions is absent (see Fig. 4). The reliability of the bimodal solution obtained by the regularization method is verified by the model data. The bimodal solution was never derived (with the level of the ‘model noise’ comparable to the real experimental errors) for unimodal true distribution. The closeness of the discrepancies 6 in Fig. 1 la and b shows, however, that the width of the ‘small mode cannot actually be determined. Note that other bimodal solutions are available in

1.0

0.5

0 .-:, z s

10

E T

0.5

p_ > 2

(b)

0

(cl

1.0

05

-f\ 0

Fig. 11. Solution transformation (b) y=O.O5, 6=9.7%; (c) y=O.15,

100

10

Particle

radius

(nm)

on change of regularization parameter y. (a) y=O.Ol; 6=8.9%; 6= 10.2%. Dashed line-log-normal distribution; rs,, = 35.6 nm, a,=2.44; 6=14%.

386

V. S. BASHUROVA

et al.

Fig. 12. Dependence of mean-square discrepancy b on regularization

parameter y,

which the ‘large’ mode is observed to the right from the physical threshold (see Fig. 9). Only the additional model experiments, to our mind, can verify the bimodality of these solutions. For the ‘shouldered’ and bimodal distributions we have divided the functions obtained into the sum of two log-normal ones to bring the distribution to the common language for describing the function of aerosol size distributions. For this discrimination we have distinguished a main modef”’ (x) (in this mode the maximum off(x) is reached), found its width ch” and maximum position r Lj . Having these parameters, we have found a preliminary weight of the mode 4i. Assuming that fb)=41f”‘b)

+ 42.P’ (XL

and having known rrh” and I&, we have determined forft”(x) its parameters CJ~’ and r\2. Then, employing the method of the least squares, the weights q1 and q2 were calculated (in line with that the value qi was determined more precisely). In the case where bimodality is not expressed obviously, the simple minimization algorithm suggested in the paper by Helsper et al. (1982) may be effective. Figure 11 demonstrates that an additional mode may change substantially with y, the values of 6 being practically the same. This causes a large error in the values of the additional mode parameters since the accuracy of their determination cannot exceed the accuracy of the additional mode reconstruction in the initial solution. The trimodal solution cannot be obtained by processing the experimental data using the regularization method. This is confirmed by both the numerical experiment (see Fig. 3) and the processing of a great body of data sets for atmospheric aerosols. 6. CONCLUSIbNS

The analysis outlined in this paper testifies that the regularization method allows the reconstruction of aerosol particle size distribution based on DB measurements. In this case RRR may be determined, and the reliability of the solution derived may be studied. In particular, our examples verify the existence of two modes in the submicron fraction of atmospheric aerosols. Generally, our results are in fair agreement with those described in the paper by Michael (1986). When using the procedure for solving the inverse problem, it is necessary to take account of the dependence of the solution quality on the [a, b] interval in which the solution is

Atmospheric aerosol size distribution functions

387

seeked. The procedure for choosing the optimal [a, b] interval, given in Section 3, is convenient to realize by taking the [a, b] interval coinciding with RRR as the initial one. On solution, the interval is corrected and the final variant does not usually coincide with RRR. Note that the model numerical experiments confirm that the solution with the minimal discrepancy [equation (16)] is possible when the [a, b] interval is approximately equal to that involving the true distribution function. For this case three situations are possible: RRR lies within [a, b]; RRR and [a, b] overlap; finally, [a, b] lies within RRR. The solution obtained may be confirmed by comparing [a, b] and RRR. Of importance in our method is the way of taking into account the information about the positiveness of the distribution function. Probably, the true minimum of the function [equation (8)] cannot be determined by our method with taking into account the condition f(x)>O. In addition, a ‘rough’ substitution of negative components in the solution offj by zeros results in a dramatic dependence of discrepancy [equation (16)] on the [a, b] interval and the regularization parameter which hampers the search for the optimal parameters. At the same time, application of the method to model problems confirms that it provides successful solutions within the error limit. The question of choosing the design parameters of a practical DB may be solved by analysing the mathematical model (see Sections 2 and 3). In particular, to reconstruct the distributions close in form to the log-normal ones (or the sum of log-normal ones) the sets of screens y, distributed uniformly over the logarithmic scale should be used. In this case, the functions K(y,, z) form a series of the uniformly distributed ‘sloping steps’ (see Fig. 1). Therefore, the quality of thef(x) reconstruction is believed to be the same within RRR. Another problem of planning the measurements is the choice of the screen parameters for a real aerosol. If the range of aerosol particle sizes is known, the screen parameters should be chosen so that RRR coincided with this range. Note that, in practice, the kernel [equation (4)] of the integral equation may be modified (i.e. the g(x) form may be changed in the exponent) to take into account other mechanisms of the deposition of large particles. This would afford the measurements of the distribution function within a wider size range. However, when allowing for interception and inertial impaction deposition the minimum emerges on the curve g(x), i.e. the inverse dependence x(g) is ambiguous. The reconstruction problem becomes ambiguous if the distribution required is situated in the region of the g(x) minimum. The problem of an optimum choice of g(x) has been considered in paper by Cheng and Yeh, (1984). Finally, the main results are summarized as follows: (a) a method has been developed for reconstructing distribution obtained using a diffusion battery (DB);

the highly-disperse

(b) it has been shown that, for the available measurement determination of more than two modes is impossible; (c) the specifications reliable results; and

have been defined for planning

aerosol size

range and accuracy, the

the DB experiments

to obtain

(d) the application of this technique has allowed us to classify the size range of the highlydisperse atmospheric aerosols. Acknowledgements-The fruitful discussions.

authors are grateful to A. N. Ankilov for experimental data and to R. A. Mavliev for

REFERENCES Bashurova, V. S., Koutzenogii, K. P. and Pusep, A. Yu. (1987) Aerosols and their application in national economy. hoc. V AI/-Union Conf, Vol. 1, p. 26. Yurmala (in Russian). Cheng, Y. S. and Yeh, H. C. (1980) J. Aerosol Sci. 11, 313. Cheng, Y. S. and Yeh, H. C. (1984) Am. ind. Hyg. Assoc. J. 45, 556. Cheng, Y. S., Keating, J. A. and Kanapilly, G. M. (1980) J. Aerosol Sci. 11, 549. Cooper, D. W. and Spielman, L. A. (1976) Atmos. Enuir. 10, 723. Crump, J. G. and Seinfeld, J. H. (1982) Aerosol Sci. Technol. 1, 15. AS 22:3-K

388

V. S. BASHUROVA

et al.

Davies, C. N. (1974) J. Aerosol Sci. 5, 593. Delves, L. M. and Walsh, J. (Editors) (1974) Numerical Solution oflnteyral Equations. Claredon Press, Oxford. Friedlander, S. K. (1977) Smoke, Dust and Haze. John Wiley, New York. Helsper, C., Fissan. H., Kapadia, A. and Lui, 9. Y. H. (1982) Aerosol Sci. Technol. 1, 135. Ivlev, L. S. and Andreev, S. D. (1986) Optical Properties of Atmospheric Aerosols. Leningrad (in Russian). Jaenicke, R. (1988) Aerosol physics and chemistry. In Numerical Data and Functional Relationships in Science and Technology (Edited by Fischer, G.), Vol. 4, pp. 391457. Springer-Verlag, Berlin. Kirsch, A. A. and Fuchs, N. A. (1968) Kolloid Zh. 30, 1144. Kondratiyev, K. Ya., Moskalenko, N. I. and Pozdnyakov, D. V. (1983) Afmospheric Aerosols. Leningrad (in Russian). Koutzenogii. K. P. (1986) Transformation and distant transfer of gaseous and aerosol admixtures in atmosphere and on creation of pollution models. Proc. All-Union Con/ p. 129. Vilnius (in Russian). Koutzenogii, K. P. (1987) Methods for determining aerosol size distribution and concentration. Review N. 4398. Moscow (in Russian). Koutzenogii, K. P., Ankilov, A. N. and Veltchev, K. (1986) Transformation and distant transfer of gaseous and aerosol admixtures in atmosphere and on creation of pollution models. Proc. All-Union Con/ p. 139. Vilnius (in Russian). Lushnikov, A. A. and Smirnov, V. I. (1975) Izu. AN SSSR. Fiz. Atm. Okeana 11, 139. Maher, E. F. and Laird, N. M. (1985) J. Aerosol Sci. 16, 557. Mavliev, R. A., Ankilov, A. N., Baklanov, A. M., Gorbunov, 9. Z., Kakutkina, N. A., Koutzenogii, K. P., Paschenko, S. E. and Makarov, V. I. (1984) Kolloid Zh. 46, 1136. Mavliev, R. A. and Ankilov, A. N. (1985) Kolloid Zh. 47, 523. McMurry, P. H. (1983) J. Colloid Interjace Sci. 95, 72. McMurry, P. H. and Friedlander, S. K. (1978) J. Colloid Interface Sci. 64, 248. McMurry, P. H. and Friedlander, S. K. (1979) J. atmos. Sci. 13, 1635. Michael, C. G. (1986) J. atmos. oceanic Technol. 3, 230. Phillips. D. L. (1962) J. Assoc. comput. Mach. 9, 84. Pusep, A. Yu, and Shokhirev, N. V. (1984) Opt. Spektr. 57, 792 (in Russian). Pusep, A. Yu. and Shokhirev, N. V. (1986) Kolloid Zh. 48, 108. Rogers, R. R. (1976) A Short Course in Cloud Physics. Oxford. Schreier, D. R. (Editor) (1982) Heterogeneous Atmospheric Chemistry. Washington. Shokhirev, N. V., Rapatsky, L. A. and Raitsimring, A. M. (1985) Chem. Phys. 105, 117. Smirnov, V. I. (1977a) Izu. AN SSSR, Fiz. Atm. Okeana 13, 274. Smirnov, V. I. (1977b) Formation of size distribution of atmospheric aerosol particles, clouds and rain drops. Author’s abstract of PhD. thesis. Leningrad (in Russian). Smimov, V. I. (1987) Meteorology and climatology (Results of science and technique VINITI AN SSSR). 15, 1 (in Russian). Smirnov, V. I. and Sergeev, 9. N. (1980) Izu. AN SSSR. Fiz. Atm. Okeana 16, 1310. Sutugin, A. G. (1975) Formation of aerosols at volumetric condensation of a fast type. PhD. Thesis NIIFKhI im. Karpova, Moscow (in Russian). Tikhonov, A. N. (1963) Dokl. Akad. Nauk SSSR. 151, 501; 153,49. Tikhonov, A. N. and Arsenin, V. Ya. (1979) Methods o/ Solution of Incorrecr Problems. Nauka, Moscow (in Russian). Twomey, S. (1963) J. ACM 10, 97. Twomey, S. (1965) J. Franklin Inst. 279, 95. Twomey, S. (1975) J. Camp. Phys. 18, 188. Vasilenko, G. I. (1979) Theory ojSigna/ Recovery. Soviet Radio, Moscow (in Russian). Whitby, K. T. (1978) Atmos. Enuir. 12, 135.