Environ. Sci. Technol. 1987, 21, 1088-1096
(3) Hagenmaier, H.; Brunner, H.; Haag, R.; Kraft, M. (Verfahren zur Zerstiirung von Chloraromaten)Patent Pending
DP 36 23 492.3. (4) Hagenmaier, H. in Mullverbrennung und Umwelt; Thome-Kosmiensky,K. J., Ed.; EF-Verlag Tech.: Berlin, 1985. ( 5 ) Buser, H. R. Ph.D. Thesis, University of Umeh, 1978.
(6) Firestone, D. J. Agric. Food Chem. 1977, 25, 1274. (7) Buser, H. R.; Rappe, C. Anal. Chem. 1984, 56, 442.
Received for review September 4,1986. Accepted May 18,1987. This work was supported by the Ministerium fur Ernahrung, Landwirtschaft, Umwelt, und Forsten, Baden- Wurttemberg, FRG.
Historical Fluxes of Particle-Bound Pollutants from Deconvolved Sedimentary Records Erik R. Christensen” and Richard H. Goetzt Department of Civil Engineering, University of Wisconsin-Milwaukee,
rn A deconvolution method to extract the historical input records of particle-associated pollutants from their sedimentary profiles, influenced by mixing, is developed. The method operates in the time domain and requires that an equal number n of time and depth intervals be appropriately chosen. An n X n matrix E of pollutant concentrations as a function of time and depth for unit influx is generated by forward programs, and the historical influxes are calculated by inverting E and multiplying with the actual concentration profile of the pollutant. Application of the method to Pb, Zn, and Cd in several Lake Michigan sediment cores shows that pollution with these metals started around 1894 in the southern part of the lake and 15-30 years later further north. A maximum was obtained during 1954-1969 for P b and Zn and during 1939-1954 for Cd, while current levels have declined 45-35% relative to these maxima in accordance with available atmospheric loading data (Pb) and U.S.consumption figures (Zn and Cd). Introduction
Effective pollution control requires an understanding of sources, phase transformations, and pathways of the pollutants investigated. To gain such insights, it is useful to determine historical fluxes of metallic and organic pollutants to lakes or coastal zones and to correlate the fluxes with alternative historical data such as sales or consumption records. On the basis of these correlations, it is then possible to make predictions for future pollutant levels in response to management strategies to reduce the impact of selected toxic substances. Examples include the ban on DDT in the early 1970s, termination of domestic PCB production in 1977, and the federally mandated reduction of Pb levels in gasoline (1,2). Historical fluxes of pollutants may be obtained from dated sedimentary records. Common time markers are 210Pb = 22.3 yr), 13’Cs ( T l j 2= 30 yr), and 239Puand 240Pu(TI = 2.4 X lo4 and 6.6 X lo6 yr), although other radionuciides such as 32Si(TI = 276 yr), 22sTh(Tl,2 = = 1620 yri, and 14C = 5730 yr) 1.91 yr), 22eRa are used on occasion (3). In the absence of mixing, the historical influx of a particle-bound pollutant will be directly reflected in the sedimentary record, even if compaction would make it necessary to carry out a trivial mapping of the nonlinear density versus depth on a linear time scale. However, most lake and near-shore sediments are subjected to some degree of mixing of either biological or physical origin. Thus, ‘Present address: Triad Engineering, Inc., Milwaukee, WI 53222. 1088
Environ. Sci. Technol., Voi. 21, No. 11, 1987
Milwaukee, Wisconsin 5320 1
rapid changes in the influx will be attenuated, or may not even show up in the sedimentary profile, while slow variations are left intact (4-6). In order to recover the original input record, it is therefore necessary to remove the influence of mixing by deconvolution. This is the so-called inverse problem as opposed to the forward problem in which the sedimentary profile of a pollutant is determined from a given input record. Although of great practical significance, the inverse problem has until recently received very little attention. An abstract by Ruddiman et al. (7) outlines the problem, but the first fully documented successful attempt to generate an inverse solution was made by Berger et al. (8), who developed an “unmixing” equation on the basis of the simple mixing model of Berger and Heath (9). Their equation may be written as ac c, = c - maz
where c, = original tracer concentration before mixing (g/cm3), c = tracer concentration in mixed layer (g/cm3), m = thickness (cm) of mixed layer, and z = depth (cm) below the sediment-water interface. The derivative should be taken at the boundary between the mixed layer and the underlying historical layers. Equation 1 was applied to deep-sea records of the l8O isotope signal. However, this method was later criticized by Jones and Ruddiman (IO), who showed that the mixing intensity used in the deconvolution operation was overestimated. A more refined approach was subsequently applied by Schiffelbein (11), who again considered l80records in deep-sea sediments. The method is based on inversion in the frequency domain and includes an impulse response function (12),which is derived from Guinasso and Schink’s (13) mixing model. Thus, features of the record that are still primarily located in the mixing zone cannot easily be deconvolved. Also, compaction and a realistic continuous mixing function were not included. Finally, the results of a test of the method, that is to bring the ISO records of two biological species of different abundances versus depth to coincide, were not particularly convincing. A time domain deconvolution method that appears to be restricted to small time reversals At for the sedimentary record well below the sediment-water interface was investigated by Christensen et al. (14). This method may be formulated as
where c and z have the same meaning as in eq 1 and where
0013-936X/87/0921-1088$01.50/0
0 1987 American Chemical Society
Table I. Lake Michigan Core Sample Information geographic position latitude longitude north west
core SLM84-CO SLM84-D1 SLM84-E1 SLM84-F1 SLM84-H1 SLM84-I1 SLM84-J1 SLM84-K1 CLM84-LO CLM84-M0 NLM84-A0 NLM84-BO NLM84-C1 NLM84-E1 NLM84-G1 NLM84-I1
I
880
Figure 1.
I
67’
I
86’
I
850
Map of Lake Michigan with sampling locations.
u (cm/yr) is the sedimentation rate, t (yr) is the time, At
(yr) is the time reversal step, and D (cm2/yr) is the diffusion coefficient. Recently, we have developed a different time domain deconvolution method that overcomes the above difficulties and that may be used to determine historical influxes directly rather than “unmixed”sedimentary profiles. While our previous report highlighted the application of the method to hypothetical data with a low degree of mixing (14), we intend here to provide a thorough description of the method itself. In addition, we describe details of the sampling and the experimental procedures, as well as the results obtained for several sediment cores from Lake Michigan. Materials and Methods
Sampling. Sediment cores were collected in the summer of 1984 from Lake Michigan by means of RV Neeskay. A 6.7-cm inner diameter gravity corer was used for sampling. Cores undisturbed by the sampling process are obtained when the following conditions are met: (a) uppermost traces of sediment on the outside of the transparent tubular polyethylene sampler are level with the inside sediment-water interface, (b) the water column above the interface in the polyethylene tube is clear immediately after sampling, meaning that little, if any, stir-up of material occurred during the sampling process, and (c) the coring device is lowered through the water column and into the sediment with a moderate, consistent velocity. The above points were observed in the extraction of our sediment cores. In addition, videotapes of core sampling taken from the submersible Johnson-Sea-Link II in Lake Superior during the summer of 1985 demonstrate that cores may be taken with minimal physical distortion. We, therefore, conclude that errors caused by sampling artifacts are insignificant. Cores were collected from three general areas of the lake (Figure 1): southern (SLM), central (CLM), and northern
42’28’00” 42’22’40’’ 42’15’30’’ 42°13’00’’ 42’19’40’’ 42’32’00’’ 42’ 38/20” 42’46’00’’ 43’38‘20‘‘ 43’45’25’’ 44’20’00” 44’28’20” 44’31’40’’ 44’40’50’’ 44’41’40’’ 44’ 35’00”
87‘ 28/20“ 87’ 17/00’’ 87’00’00’’ ‘ 86’ 47/30’’ 86’40’35” 86’42’20‘‘ 86’ 57’00” 87’ 15’00’’ 86’46’00” 86’39’10” 86’38’20” 86’43’25’’ 86’42’20’’ 86’45’00” 86’55’00” 87’04/36’’
water depth, m
core depth, cm
88 110 100 90 72 87 160 125 111 81 235 250 263 263 216 190
51.0 59.5 36.0 31.0 66.0 73.0 84.0 84.0 60.5 55.0 90.5 84.0 78.0 73.0 28.0 65.0
Lake Michigan (NLM). Each core was given a letter designation for identification and a number to represent the individual sampler (1-3 for the legs in a triple corer and 0 for a single corer) with which the core was obtained. The number 84 is included to indicate that the samples were collected in 1984. For example, SLM84-J1 means that the core is from the J area of southern Lake Michigan and is obtained with core sampler number 1. Core depth, water depth, and coordinates of each core are listed in Table I. Experimental Procedures. The cores were sectioned aboard ship with a standard extrusion device. Briefly, the rubber stopper that was affixed to the bottom of the polyethylene tube at the time of bringing the core aboard ship was removed, and the bottom end was then subjected to a controllable water pressure in the extrusion apparatus. By gently increasing the water pressure, the entire core with the top end, i.e., the water-sediment interface, first was extruded through the upper end of the polyethylene tube. Thus, measured sections of the core were pushed out of the tube and sectioned in the process by means of a thin clean nalgene blade. Each core section was stored in a separate glass or nalgene bottle. At the end of the sampling cruise, sample containers were transferred from the ship to an onshore storage cold room kept at 4 “C. Because the extrusion was done very carefully, we feel that errors associated with this process were negligible. Porosity was determined by a gravimetric method within 48 h of sampling. Samples for radiometric and metals analysis were dried for 24 h at 103 “C and ground with a mortar and pestle. Determination of 13’Cs was carried out by a NaI(T1) detector and automated y spectroscopy (662 keV) and of 210Pbby plating of the granddaughter zlOPo on copper disks (15) and (Y spectroscopy with surface barrier detectors. The metals Pb, Zn, and Cd were measured by acid digestion (HC1/HN03) and atomic absorption spectroscopy (AAS) with Instrumentation Laboratory Model Video 12 atomic absorption spectrophotometer in the Great Lakes Research Facility. A graphite furnace method was applied for Pb and Cd, whereas flame AAS was adequate for Zn. The analytical results were checked for accuracy by analyzing National Bureau of Standards Standard Reference Material 1645, river-sediment. The measured concentrations of Pb, Zn, and Cd were generally within the certified ranges (certified value f certified standard deviation), Forward Problem. The forward problem is to determine the sediment profile s (dpm/g or g/g, i.e., gram of Environ. Sci. Technol., Vol. 21, No. 11, 1987 1089
pollutant per gram of sediment) of a tracer when the input to the sediment is given as a function of time. A solution procedure has been developed previously on the basis of the advection-diffusion equations for sediment solids and a tracer (16): at
(3)
where z (cm) is the depth below the sediment-water interface, u (cm/yr) is the sedimentation rate, t (yr) is the time, D (cm2/yr) is the diffusion coefficient, and X (yr-l) is the radioactive decay constant. The bulk density of sediment solids is defined as p = (1- 4)pa, where 4 is the porosity and pS is the density of sediment solids (2.45 g/ cm3). We shall in the following assume that the sediment characteristics are stationary so that a p / a t = & / a t = aD/at = 0. It is further assumed that the bulk density and diffusion coefficient are functions of depth: p
= P- - p1 exp(-az)
D = Doexp[-z2/(2u2)]
(5) (6)
where PO = pm- p1 is the density at the top, pmis the density at large depths, Dois the diffusion coefficient at the top, and u is an effective mixing depth. Because of the assumed steady-state sediment characteristics, the mass sedimentation rate r (g cm-2 yr-l) is a constant. The differential eq 3 and 4 are solved with proper initial and boundary conditions, using finite differences, either for steady-state (210Pb)or time-dependent (137Cs,Pb, Zn, Cd) conditions. The 137Csinventory expected from atmospheric input Ma (dpm/cm2) is given by
Me = CF, exp(-XtJ i
(7)
where Fi(dpm/cm2) is the atmospheric 137Csinflux in year i, determined by multiplying the data (17) with a factor of 1.5 ( l a ) ,A (yr-l) is the decay constant for 137Cs, and t, is the number of years the nuclide has been in the sediment. Fallout data F, for Argonne, IL, were used for the southern Lake Michigan stations; Green Bay, WI, data for the northern Lake Michigan stations; and an average of Argonne and Green Bay data for the stations in central Lake Michigan. The measured sediment inventory of 137Cs Mm(dpm/cm2) is given by
Mm = CP,S,AZL
(8)
2
where p, (g/cm3) and s, (dpm/g) are the bulk density and measured l3Vs activity, respectively, in the ith core section of depth Az,. Similarly, the measured steady-state 210Pb flux Jo to the sediment is determined from Jo = XCP,S,AZ,
(9)
k
where the symbols have the same meaning as above, except that the decay constant X and the measured activity s, refer to 210Pb. Inverse Problem. In the inverse problem, the sedimentary profile of the pollutant is given, and the input that produced this profile is to be determined. Another way of expressing this is to say that the solution s to eq 4 is known, while the time-dependent boundary condition, Le., the input flux to the sediment, is the unknown. We ap1090
Environ. Sci. Technol., Vol. 21, No. 11, 1987
proach this problem by deconvolution in the time domain. For a linear, time-invariant system, the sediment concentration of a tracer s(z,t) may be expressed in the form of a convolution integral:
where so(z,t-r) is the impulse response, f ( r ) is the influx at time r , and z and t are the actual depth and time, respectively. We apply a discrete version of eq 10:
where i and j are space and time indices, S (dpm/g or g/g) is the concentration profile in the sediment, E (g-l cm2yr) is a matrix of elementary contributions to this profile for unit influx, and F (dpm cm-2 yr-l or g cm-2 yr-l) is the historical input flux to be determined. F is calculated by inverting E and multiplying with S. The matrix inversion was carried out by Gauss elimination (19). Both eq 10 and 11 are exact. They express the tracer concentration in the sediment as a superposition of elementary contributions from all previous times, assuming that there was no influx for t < 0 or j < 1, In the discrete version, eq 11,it is assumed that the influx is constant in each of the n time intervals and that the arithmetic average concentration is considered in each of the n sediment intervals. In fact, since the physical sampling intervals are finite, i.e., usually 0.5,1.0, or 2.0 cm each, all that is known about the tracer concentrations are averages, meaning that the discrete version lends itself naturally to the solution of the inverse problem. The method can now be summarized as follows. An equal number n of space and time intervals is chosen such that each depth interval is at least as large as the physical sampling intervals and with the time intervals of such a duration T that nT covers the period of interest, i.e., at least the length of time where tracer (pollutant) inputs have occurred. Further restrictions on the depth intervals, acquired through the actual implementation of the method, will be described in the following section. Next, the matrix E is determined by forward programs as solutions for unit influx to eq 3 and 4. Finally, E is inverted and multiplied wih the concentration profile S to produce the influxes F. The correctness of the influxes F calculated by inversion from eq 11was checked by two different methods. First, a hypothetical tracer input record was applied to the cores, and the resulting sedimentary profile S was calculated by forward programs. On the basis of the calculated S profile and E matrix, the input fluxes were then regenerated by inversion from eq 11. The calculated influxes were compared with the original ones. Obviously, a good match is necessary in order for the method to work on real systems. The second test involved actual metal (Zn) data. The concentration profile S was determined from the sedimentary concentrations by subtracting appropriate preindustrial levels, The influxes F were then determined by inversion from eq 11, and these fluxes were, in turn, used to regenerate S by the forward programs. The regenerated profile was finally compared with the original tracer profile. If these profiles were identical or nearly so, this would provide strong evidence in favor of the method. Further checks on the method were made by comparing the reconstructed historical fluxes with available atmospheric loading data (Pb) or U.S. consumption data for Zn and Cd. In addition, reconstructed fluxes from areas of considerable mixing were compared with similarly obtained
Table 11. Core Parameters density parameters p-, g cm-3 pl, g cm"
core
a, cm-l
mass sedimentation rate r, g cm-2 yr-l
zloPbflux Jot dpm cm-2 yr-'
SLM84
co
D1 El
F1 H1 I1 J1 K1 CLM84 LO MO NLM84 A0 BO
c1 El G1 I1
mixing parameters Do,cm2 yr-I u, cm
0.570 0.360 0.808 0.894 0.624 0.547 0.432 0.403
0.4759 0.1898 0.4972 0.8747 0.3583 0.4009 0.2660 0.2714
0.1755 0.1237 0.2867 0.3555 0.0445 0.1594 0.0297 0.0649
0.01353 0.00970 0.00766 0.01156 0.06088 0.01405 0.01765 0.03139
0.8933 0.5213 0.1988 0.6430 2.2823 0.8458 1.1636 2.4581
20.00 50.05 50.20 0.20 0.10 2.30 0.60 120.00
10.05 0.10 0.55 3.55 0.50 0.90 1.15
0.453 0.904
0.2546 0.6822
0.2033 0.0607
0.01847 0.02086
0.7862 0.7327
50.10 0.10
0.10 0.10
0.411 0.437 0.300 0.406 0.443 0.332
0.3162 0.3247 0.2188 0.2887 0.3693 0.2546
0.0726 0.0446 0.1402 0.0707 0.1934 0.2886
0.01324 0.01586 0.01233 0.01253 0.01019 0.00465
0.9327 1.3853 0.9596 1.0550 0.5880 0.3135
12.00 18.00 1.00 19.00 0.70 50.10
0.40 0.65 0.55 0.45 0.25 0.15
TOTAL Cs 137 ACTIVITY
08
I
I
TOTAL Cs 137 ACTIVITY
EXCESS Pb 210 ACTIVITY d p m l g
POROSITY $ 07
09 1001
I
1
02
I
'
04
08 10
1 ' / ' / 1 1
20 I
'
PO
0.50
dmlo
8 0 1 0 0 200 4 0 0 6 0 0 0
10
20
30
40
2
:t ,A ! l0C
1
10 12
12
14
14 D,=O i o c r n 2 i y r
16
Do= 0 10 crn2/yr - 16
I
o=o i o crn D,=O
i o crn2iyr 16
o:355cm 1=0.06086 g/Cm21yr
a=010cm
r = O 02086 g/crn21yr
-20
%=09608 Ma
Jo=O 7327 d p m l c m 2 / y r
1
22
J0=2 2823 dprnlcm21yr
-18 -20
- 22
24
i
26
26 -
- 26
26 -
- 28
Flgure 2. Actual and calculated porosities and activities of 210Pband 137Csfor core SLM84-Hi.
fluxes from areas of relatively little mixing. 07
Results and Discussion The determination of the source function F from eq 11 requires that the sedimentation, mixing, and compaction characteristics of the sediments are well-defined so that the matrix E can be calculated accurately. The execution of the method is, therefore, a two-step process. First, tracers such as 210Pb(steady state) and 13'Cs (time dependent) are used to determine the sediment properties from the forward programs (16);next, these characteristics are used in the implementation of the method to the pollutants of interest, which are Pb, Zn, and Cd here. Determination of Sediment Characteristics. The results of the evaluation of the sediment core parameters are listed in Table 11. The largest mass sedimentation rates r and bulk densities pm are found in southern and central Lake Michigan. As expected from a particle-associated tracer such as 210Pb,there is a strong correlation between the zloPbflux Jo and the mass sedimentation rate r. Also, the variability of the 210Pbflux in the northern basin is less than in the southern basin, reflecting the less significant sediment focusing in the north (16). Examples of actual and calculated porosity and 210Pb and 137Cs profiles are given in Figures 2-5 for cores SLM84-H1, CLM84-MO, NLM84-BO, and NLM84-El. Parameters from Table I1 are indicated in these figures
08
09 1 0 0 1
0 2 , 0 4 ,, q 8 , i O
20
, PO,
,:\iQO
200 4 0 0 8 0 0
2-
46-
*
810-
6 t
121416-
18-
22 24 20
/+-i;
j
-12
14
~
26
-
0,=i8 00 crn2/yr o:065cm
1=001566 g/Cm2/yr
~ 3853 ~ d p= r n / c ~i2 / y r
IOo=1600crn21yr-16
I$=' M
o = 0 640g7 5crn
? -18 22
I
24
26 -
30
a -10
- 26 ~
I
- 26 30
Table 111. SLM84-J1 E Matrix ( p g / g ) for Unit Input Flux (fig
yr-')
time, yr
depth, cm 0-2 2-4 4-6 6-8 8-10 10-12 12-14 14-16 16-18 18-20
'
15
30
45
60
32.05 12.84 0.12 0.00 0.00 0.00 0.00 0.00 0.00 0.00
13.61 21.12 6.73 0.06 0.00 0.00 0.00 0.00 0.00 0.00
6.07 12.40 17.79 2.84 0.02 0.00 0.00 0.00
0.00 0.00
2.71 5.7315.83 11.94 1.02 0.01 0.00 0.00 . 0.00 0.00
75
90
105
120
135
150
1.21 2.57 8.77 16.46 6.33 0.31
0.54 1.15 4.12 12.48 13.26 2.70 0.08 0.00 0.00 0.00
0.24 0.51 1.86 6.90 14.42 8.29 0.95 0.02 0.00 0.00
0.11 0.23 0.83 3.31 10.42 12.92 4.09 0.28 0.00 0.00
0.05 0.10 0.37 1.51 5.89 12.74 8.95 1.63 0.07 0.00
0.02 0.05 0.17 0.68 2.91 9.18 12.09 4.87 0.54 0.02
0.00
0.00 0.00 0.00
Time (years)
361
24: 26
28-
'
~
I
I
- 24 - 26 - za
I
Depth (cm) Flgure 6. Tracer profiles at various times in core SLM84-Jl in response to the given unit input flux. The tracer profiles are used for the E matrix determination (see Table 111).
catch the tail of the tracer profile from the latest time period (150 yr). In Figure 6, depth intervals Az of 2 cm were used. Although it does not appear from the figure as if nonzero values occurred within the deepest depth interval (18-20 cm), values of the order of pg/g actually existed. Matrix element Eijis the arithmetic average of the tracer concentration (pg/g) in depth interval i, with a depth of Az(i-1) to Azi cm, and at time index j , that is, Atj years after the front of the unit input pulse entered the sediment. For example, ES5= 8.77 pg/g is the average concentration in depth interval 3 (4-6 cm) 75 years after the unit input pulse first entered the sediment. The resulting E matrix can be seen in Table 111. No row or column should contain all zeros (E would then be singular), and the last depth interval should go far enough into the tail of the curve from the last time period. Several methods for depth interval division were attempted prior to finding the described method. Variable depth increments were initially utilized. The increment size varied according to the width of the peaks for each time period and were centered on the corresponding peak to generate a matrix with the maxima along the diagonal. This method provided poor results, possibly due to the 1092
Environ. Sci. Technol., Vol. 21, No. 11, 1987
Table IV. Calculated Influxes Based upon a Generated Sediment Profile from a Hypothetical Known Influx SLM84
H1 J1 original calcd original calcd original calcd influx influx influx influx influx influx
D1
time, yr
1.0 2.0 1.0 3.0 5.0 6.0 8.0 7.0 6.0 9.0
0-15 15-30 30-45 45-60 60-75 75-90 90-105 105-120 120-135 135-150
1.00 2.01 1.02 3.00 5.06 6.01 8.05 7.03 6.03 9.07
.LO 0.0 1.0 2.0 2.0 3.0 5.0 7.0 8.0 8.0
1.00 0.00 1.00 2.00 2.00 3.00 5.00 7.00 8.00 8.01
1.0 0.0 1.0 2.0 2.0 3.0 5.0 7.0 8.0 8.0
1.00 0.00 1.00 1.98 2.04 2.95 5.05 6.97 8.00 8.01
4 c
7.0
6.5-
60-
1.0 0.0 1.0 2.0 3.0 4.0 4.0 5.0 6.0 5.5
1.00 0.00 0.99 2.02 2.94 4.09 3.89 5.09 5.95 5.56
1.0 2.0 1.0 3.0 5.0 6.0 8.0 7.0 6.0 9.0
1.00 2.02 0.96 3.12 4.87 6.20 7.94 7.04 6.07 9.01
1.0 0.0 1.0 2.0 3.0 4.0 4.0 6.0 5.0 5.5
1.00 0.00 1.00 2.00 3.00 4.00 4.01 5.99 5.01 5.49
Table V. Comparison of Actual Zn Sediment Profile with Profile Calculated by Forward Programs
depth, cm 0-2 2-4
4-6 6-8 8-10
10-12 12-14 preind. level, w m
calcd profile baaed actual concn in profile upon calcd Zn input minus preindustrial level, ppm flux, ppm SLM84-Jl NLM84-BO SLM84-J1 NLM84-BO 205.35 195.10 128.65 70.40 12.70 0.00 0.00 130.0
250.95 237.90 173.70 57.40 2.00 0.00 0.00 113.8
205.62 195.53 130.05 74.86 17.25 0.73 0.04 130.0
SLM84.Hl
.:05
-
CLM84-MO BO El original calcd original calcd original calcd influx influx influx influx influx influx
0-15 15-30 30-45 45-60 60-75 75-90 90-105 105-120 120-135 135-150
....
.... 5.5-
NLM84 time, yr
-
253.16 240.76 178.83 66.04 5.72 0.07 0.00
113.8
varying size of the depth increments. Another method also involved choosing depth intervals around the peaks for each time period, but the size of the depth increments were kept constant. This method gave satisfactory results on some cores but was not consistent. The method that was eventually devised, as described above, provides consistently good results. As previously mentioned, one test of the method is to apply a hypothetical input record, calculate the sedimentary profile, and regenerate the input record by inversion from eq 11. Results of such a test are displayed in Table IV. The errors are generally of the order of 1% , almost zero when the relative effect of mixing is mall (SLM84-H1), and up to 4% when there is a fair amount of mixing (NLM84-BO). In the other test, corrected Zn profiles were first generated from the measured one by subtracting constant preindustrial levels for each core. The inversion procedure was then applied to the corrected Zn profiles, and the calculated influxes were, in turn, used to generate Zn profiles in the sediments. As may be seen from Table V, the agreement between calculated and actual profiles is good, with less than 3% error in the upper three depth
h N '
45-
. 2
40-
m
35-
1984 1954 1924 1894 1664 YEAR
1834
t L
... -
Z 30-
k 25-
,
Figure 7. Historical inputs of Pb to Lake Michigan sediments. Full horizontal lines are calculated by deconvolutlon, and dashed horizontal lines are obtained directly from the sedimentary record. The curve (NLM84-El) shows atmospheric loading data.
intervals. The lengths of the physical sampling intervals were 1 cm for a depth of 0-10 cm and 2 cm for a depth of 12-14 cm of core SLM84-J1. Similarly, these thicknesses were 0.5 cm for a depth of 0-4 cm, 1cm for a depth of 4-10 cm, and 2 cm for a depth of 10-14 cm of core NLM84-BO. In the general application of the inverse method to reconstruct historical input fluxes, preindustrial levels of Pb, Zn, and Cd were first estimated from the deep portion of the cores where the levels have decreased to fairly constant, low values. Next, minor adjustments in these levels were made by minimizing any calculated negative input fluxes that might occur as part of the early historical record, when the anthropogenic influxes were small. Application to Lake Michigan Sediments. Results of the deconvolution for the Lake Michigan sediment cores are shown in Figures 7-9. While deconvolution of the sedimentary Pb, Zn, and Cd profiles was attempted for all of the cores listed in Table IV, only the successful attempts are shown in the figures. In other cases, significant negative input fluxes, especially in the early historical record, were indicative of errors either in the model assumptions or, more likely, in the determinations of metals. It is symptomatic that all of the Zn profiles could be deconvolved (Figure 8) whereas only three profiles out of six for Pb (Figure 7) and four out of six for Cd (Figure 9) were amenable to deconvolution. This is consistent with the fact that Zn was determined with higher precision and accuracy than it was possible to achieve for Pb and Cd. The input records, as obtained directly from the sedimentary profiles, are also shown in Figures 7-9. These inputs (dashed line) were obtained as the product of the Environ. Sci. Technol., Vol. 21, No. 11, 1987
1093
lO.Or
lO.Or
:4,0r5& 6 2.0
I
"k_
,
_-__
0.0 1984 1954 1924 1894 1864 1834
I
__---
0.0
1984 1954 1924 1894 1864 1834
YEAR
YEAR
22 0
-
20 0
-
- 2.2 - 2.0
18.0 -
- 1.8 - 1.6
1984 1954 1924 1894 1864 1834
0
5 . 9 -
W
0 7
12.0 '-
c
2 z
r
4,0k -0.0
I
-_&--
4
0.0
-
UJ 2
6.0 -
- 0.6
-
- 0.4
4.0
r.5 2.0
2
0
10.0 -
El 8.0
SLM 84-J1
6.0 c
- 0.2
20I
I
I---
U
%
- 1.4 - 1.2 - 1.0 5 - 0.8 v1
YEAR
r
.
'
0.01
I
1
-_
-7 I Y
1 0.0
Figure 8. Historical inputs of Zn to Lake Michlgan sediments. Full horizontal lines are calculated by deconvolutlon, and dashed horizontal lines are obtained directly from the sedimentary record. The curve (SLM84-H1) shows U.S. consumption data.
mass sedimentation rate r (g cm-2 yr-l) and the average metal concentration s (kg/g),where the time intervals are given by the cumulative sediment mass (g/cm2) divided by r (g cm-2 yr-l). I t is clear that deconvolution has the effect of sharpening peaks, as for example in core SLM84-H1, and shifting the onset time of pollution toward more recent times. The latter effect is of the order of one to two time intervals, Le., 15-30 years. Furthermore, when the effect of mixing is sufficiently strong, the deconvolved records may show a peak where none were evident from the unprocessed records. Examples are P b in core NLM84-E1 (Figure 7) and Cd in core NLM84-BO (Figure 9). Superimposed on these records are atmospheric lead loading data for the Great Lakes Region and U.S.consumption data for zinc and cadmium. The P b data were obtained by combining the source data of Edgington and Robbins (21) with the data for lead in gasoline from 1942 to 1983 given by Trefry et al. (2). The latter data were normalized to corresponding values from ref 21 and were used as gasoline-derived P b data from 1942 to 1983. Loading data for P b from coal burning (21) were assumed to have remained constant at 600 tons/yr since 1974. U.S. consumption data for Zn and Cd were taken from ref 22. 1094
Envlron. Sci. Technol., Vol. 21,
No. 11, 1987
The source functions have been added to one of the reconstructed records for which they appeared to provide a good match, that is, NLM84-E1 for lead (Figure 7), SLM84-H1 for zinc (Figure 8), and CLM84-MO for cadmium (Figure 9). The reconstructed input records are fairly consistent throughout the lake for Pb and Cd (Figures 7 and 9). Lead shows a maximum during 1954-1969, and Cd has one during 1939-1954. The Pb records are in accordance with atmospheric loading data for the Great Lakes Region, and the Cd records approximately follow U.S. consumption data. In contrast to Pb and Cd, the input records for Zn show some local variation throughout the lake (Figure 8). Cores in the southwestern part (CLM84-MO, SLM84-J1, and SLM84-H1)have a maximum during 1954-1969, while the northern ones (NLM84-BO and NLM84-El) exhibit no distinct maximum but rather a constant plateau during the last -45 years. The reason for this pattern is not known. However, since the three metals Pb, Zn, and Cd enter the lake primarily through the atmospheric route (23, 24), it would appear that the sources and atmospheric transport paths of Zn are more local in nature than those of Pb and Cd. The input record of Zn at core CLM84-MO
NLM84.BO
008
00 1984 1954
0
1
0
y
1924 1894 1864 1834 YEAR
1 10
1984 1954 1924 1894 1884 1834 YEAR
._ L_
1984 1954 1924 1894 1864 1834 YEAR
Figure 9. Historical inputs of Cd to Lake Michigan sediments. Full horlzontal lines are calculated by deconvolution, and dashed horizontal lines are obtained directly from the sedimentary record. The curve (CLM84-MO) shows US. consumption data.
forms not only an approximate average for Lake Michigan but also follows the U S . consumption data well, thus supporting the validity of the reconstructed historical records. Future research will focus on improvements in the time resolution of deconvolved records from the present 15 years to 2-3 years and on the development of error measures for deconvolved records. Conclusions
A deconvolution method to recover the historical input records of particle-bound pollutants from their sedimentary profiles, influenced by mixing, has been developed. The method operates in the time domain and is based on the principle of superposition for a linear system. A major advantage of this method over frequency domain methods is that features of the profile that are still predominantly in the mixing zone easily can be deconvolved. By contrast, frequency domain methods require that the sedimentary profile has left the mixing zone and has been immobilized in the historical layers. The developed method has been tested for internal consistency with hypothetical input records and for accuracy by means of actual Zn data. The method relies on accurately calculated pollutant concentrations at n time intervals and n depth intervals from a unit input flux. These forward calculations are carried out by a previously published method based on solution of the advectiondiffusion equations for sediment solids and a tracer (pollutant) by means of finite differences. Specific conclusions regarding the method itself and its application to Pb, Zn, and Cd in Lake Michigan sediment cores are as listed: (1)The n time intervals of the same length T should cover the time period of interest, that is, the time when
anthropogenic inputs have occurred. We have here chosen n = 10 and T = 15 years. The n depth intervals should be of equal depth Az and should, furthermore, be chosen such that they are at least as deep as any physical sampling interval and such that the lower boundary n A z of the deepest interval catches the front tail of the sediment profile from the unit influx of tracer (pollutant), the front of which has spent the maximum time nT in the sediment. Typical values of Az were 2 (SLM84-J1, NLM84-BO), 3 (SLM84-H1), and 1.5 cm (SLM84-MO, NLM84-El). (2) The effects of deconvolution on sedimentary records of pollutants are to sharpen peaks (SLM84-H1),to restore peaks lost because of mixing as, for example, for Pb in core NLM84-E1 and Cd in core NLM84-BO8 and to shift the onset times of pollution toward more recent times. From the reconstructed input records it appears that pollution with Pb, Zn, and Cd started around 1894 in the southern end of Lake Michigan and 15-30 years later in the central and northern parts. The records for Pb and Cd are consistent throughout the lake and show a maximum during 1954-1969 for P b and during 1939-1954 for Cd. The pattern of Zn pollution is more local in nature with a maximum during 1954-1969 in the southeastern part of the lake and a constant high plateau during the last 30 years in northern Lake Michigan. The deconvolved records for P b are in remarkable agreement with atmospheric P b loading data for the Great Lakes Region. Similarly, there is a good correlation between the reconstructed historical records for Cd and the average one throughout the lake for Zn on one hand and the U.S. consumption data for these metals on the other. These facts support the validity of the approach and point to the usefulness of conducting similar deconvolution studies in other aquatic systems where the sedimentary profiles have been distorted by mixing. Acknowledgments
We thank Mark Hermanson for assistance with chemical analytical procedures and Donald Nelson for cooperation in 13'Cs counting. Prasanta Bhunia helped in the initial phases of the solution procedure for the inverse problem. Registry No. Pb, 7439-92-1; Zn, 7440-66-6; Cd, 7440-43-9. Literature Cited (1) Eisenreich, S. J.; Metzer, N. A.; Urban, N. R.; Robbins, J. A. Environ. Sci. Technol. 1986, 20, 171-174. (2) Trefry, J. H.; Metz, S.; Trocine, R. P.; Nelsen, T. A. Science (Washington, D.C.)1985, 230, 439-441. (3) Krishnaswami, S.; Lal, D. In Lakes: Chemistry, Geology, Physics; Lerman, A., Ed.; Springer-Verlag: New York, 1977; Chapter 6, pp 153-177. Goreau, T. J. Nature (London) 1980,287, 620-622. Schiffelbein, P. Nature (London) 1984, 31I , 651-653. Miller, K. M.; Heit, M. Limnol. Oceanogr., in press. Ruddiman, W. F.; Dicus, R. L.; Glover, L. K. Geological Society of America, Meeting Abstract 1976; Geological Society of America: Boulder, CO, 1976; p 1079. (8) Berger, W. H.; Johnson, R. F.; Killingley, J. S. Nature (London), 1977,269, 661-663. (9) Berger, W. H.; Heath, G. R. J. Mar. Res. 1968, 26(2), 134-143. (10) Jones, G. A.; Ruddiman, W. F. Quat. Res. (N.Y.) 1982,17, 148-172. (11) Schiffelbein, P. Mar. Geol. 1985, 64, 313-336. (12) Officer, C. B.; Lynch, D. R. Mar. Geol. 1983, 52, 59-74. (13) Guinasso, N. L., Jr.; Schink, D. R. J. Geophys. Res. 1975, 80(21), 3032-3043. (14) Christensen, E. R.; Bhunia, P. K.; Hermanson, M. H. (4) (5) (6) (7)
Proceedings of the 5th International Conference on Heavy Metals in the Environment;CEP Consultants: Edinburgh, UK, 1985; Vol. 1, pp 116-118. Environ. Sci. Technol., Vol. 21, No. 11, 1987
1095
Environ. Sci. Technol. 1907, 21, 1096-1 102
MacKenzie,A. B.; Scott,R. D. Analyst (London) 1979,104, 1151-1158. Christensen, E. R.; Bhunia, P. K. J . Geophys. Res. C: Oceans 1986, 91 (C7), 8559-8571. Health Saf. Lab. Environ. Q. (U.S. Energy Res. Dev. Adm.) 1977, HASL-329. Health Saf. Lab. Enuiron. Q. (U.S. Energy Res. Dev. Adm.) 1972, HASL-258. Kreyszig, E. Advanced Engineering Mathematics; Wiley: New York, 1983. Krezoski, J. R.; Robbins, J. A. J . Geophys. Res. C: Oceans 1985, 90(C6),11999-12006. Edgington, D. N.; Robbins, J. A. Environ. Sei. Technol. 1976, 10(3),266-274.
(22) U.S. Bureau of Mines Minerals Yearbook;U.S. Department of the Interior: Washington, DC, 1984. (23) Eisenreich, S. J. Water, Air, Soil Pollut. 1980,13, 287-301. (24) Christensen, E. R.; Chien, N. K. Enuiron. Sci. Technol. 1981, 15(5),553-558.
Received for review November 14,1986. Accepted July 13, 1987. Although the information described in this paper has been funded by the U.S. Environmental Protection Agency under Assistance Agreement R810419 to E.R.C., it has not been subjected to the Agency’s required peer and administrative review and, therefore, does not necessarily reflect the views of the Agency, and no official endorsement should be inferred.
A Microscale System for Estimation of Model Parameters for Fixed-Bed Adsorbers Walter J. Weber, Jr.,* and Chang Keun Wang Environmental and Water Resources Engineering Program, The University of Michigan, Ann Arbor, Michigan 48 109
W
A technique for prediction of the performance charac-
teristics of fixed-bed adsorbers using parameters estimated from microscale experiments (micro diameter-depth adsorption systems, MIDDAS) is described. This technique employs specially designed “microcolumn” adsorbers containing small adsorbent particles to develop information for scale-up and design from relatively straightforward experimental measurements and associated model simulations.
Introduction The design of fixed-bed adsorbers generally requires information regarding the breakthrough characteristics of specific “target” compounds in systems containing several adsorbable species. Such information can be gained empirically from pilot-plant investigations, but tests of this nature are typically lengthy and expensive. Mathematical models that can utilize more fundamental data to predict adsorber performance for alternative design configurations and different conditions of operation provide useful adjuncts to pilot-scale studies. At a minimum, they can be used to optimize investigations involving pilot-scale fixed-bed reactor (FBR) systems and to facilitate interpretation of resulting data. Predictive FBR models require input of adsorption equilibrium and mass transport parameters characteristic of the particular thermodynamic and kinetic properties of the system(s) of interest (1). A variety of experimental techniques and numerical correlations for evaluation or estimation of such parameters are available, but many of these involve unacceptable levels of difficulty, imprecision, and uncertainty. The micro diameter-depth adsorption system (MIDDAS) technique described here was developed in response to the need for a relatively simple yet accurate means for determination of the equilibrium and mass transport parameters required for FBR adsorber modeling. The technique employs “deep-bed adsorber” microcolumns (MIDDAS-DBA’s)for equilibrium parameter determinations and ”short-bed adsorber” microcolumns (MIDDASSBAs) for mass transport or rate parameter evaluations. Both the MIDDAS-DBA and MIDDAS-SBA microcolumn systems employ adsorbent particles of significantly smaller diameter than typically used in full-scale FBR systems. 1096
Environ. Scl. Technol., Vol. 21, No. 11, 1987
The use of “micro” particle diameter and bed depth dimensions promotes high rates of adsorption and early bed exhaustion in MIDDAS systems, thus facilitating rapid parameter measurements while still affording sufficiently representative system configuration and flow conditions to ensure accurate parameter characterization. A scale-up procedure has been developed to utilize parameters estimated from the MIDDAS technique to predict the breakthrough properties of target compounds in FBR adsorbers comprised of larger diameter columns, greater bed depths, larger adsorbent particles, and adsorbents of polydisperse particle size.
MIDDAS Development The MIDDAS technique employs a modified version of the “high-pressure minicolumn” (HPMC) method described by Rosene et al. (2) in combination with the microcolumn or short-bed adsorber (SBA) technique of Weber and Liu (3, 4 ) to determine capacity and mass transport parameters, respectively. The SBA is defined as an adsorber of sufficiently shallow depth to ensure “immediate” breakthrough of the compound(s) of interest. Mass transfer thus occurs over the entire bed depth after one theoretical “bed volume” or unit detention time. Weber and Liu ( 3 , 4 )employed sensitivity analyses of SBA breakthrough curves with respect to the extraparticle “film”-transfer coefficient (123 and intraparticle “surface diffusion” coefficient (D,) associated with a widely used two-resistance FBR adsorber model (1) to demonstrate that (1)the initial stage of an incipient breakthrough event in an SBA is dominated by film transfer, regardless of the eventually predominant mass transport control, and (2) the response of an SBA to changes in the parameters kf and D,is more sensitive than is that of a deep FBR adsorber. The authors verified the SBA approach for single-solute systems by successfully predicting FBR adsorber breakthrough curves using models calibrated with parameters obtained from SBA measurements and corresponding model simulations. Liang and Weber (5) and Smith et al. (6)subsequently extended the SBA procedure to estimation of mass-transfer coefficients for mixed-solute systems. The general validity of the SBA approach to rate parameter estimation has since been confirmed by other investigators (7, 8).
0013-936X/87/0921-1096$01.50/0
0 1987 American Chemical Society
