Environ. Sci. Technol. 1991, 25, 474-481
(15) Farkas, L.; Lewin, M.; Bloch, R. J. Am. Chem. SOC.1949, 71, 1988-1991. (16) Peintler, G.; Nagypal, I.; Epstein, I. R. J . Phys. Chem., in press. (17) Gordon, G. J . Pure Appl. Chem. 1989, 61, 873-878. (18) Werdehoff, K. S.; Singer, P. C. J.-Am. Water Works Assoc. 1987, 79, 107-113. (19) Singer, P. C.; O’Neil, W. K. J.-Am. Water Works Assoc. 1987, 79, 75-76. (20) Schmitz,G.; Rooze, H. Can. J . Chem. 1981,59,1177-1187. (21) Schmitz,G.; Rooze, H. Can. J . Chem. 1984,62,2231-2234. (22) Cady, G. H. In Inorganic Syntheses; Moeller, T., Ed.; McGraw-Hill Book Co., Inc.: New York, 1957; Vol. V, pp 156-165. (23) Tang, T.-F.; Gordon, G. Anal. Chem. 1990,52, 1430-1433. (24) Suzuki, K.; Gordon, G. Anal. Chem. 1978,50,1596-1597. (25) Ikeda, Y.; Tang, T.-F.; Gordon, G. Anal. Chem. 1984,56, 71-73. (26) Gordon, G.; Cooper, W. J.;Rice, R. G.; Pacey, G. E. J.-Am. Water Works Assoc. 1988,80, 94-108. (27) Gordon, G.; Ikeda, Y. J.-Am. Water Works Assoc. 1984, 76, 98-101.
(28) Silverman, R. A.; Gordon, G. Anal. Chem. 1974, 46, 178. (29) Silverman, R. A.; Gordon, G. J . Chem. Educ. 1973, 50, 654-655. (30) Moore, J. W.; Pearson, R. G. Kinetics and Mechanism, 3rd ed; John Wiley & Sons: New York, 1981; p 21. (31) Taube, H.; Dodgen, H. J . Am. Chem. SOC.1949, 71, 3330-3336. (32) Lister, M. W. Can. J. Chem. 1956, 34, 465-478. (33) Lister, M. W. Can. J . Chem. 1956, 34, 479-488. (34) Lister, M. W.; Petterson, R. C. Can. J . Chem. 1962, 40, 729-733. (35) Halperin, J.; Taube, H. J. Am. Chem. SOC.1952,74,375-380. (36) Fabian, I.; Gordon, G., unpublished results, Miami University, Oxford, OH, 1990. (37) Lister, M. W. Can. J . Chem. 1952, 30, 879-889. (38) Yokoyama, T.; Takayasu, 0. Kogyo Kagaku Zasshi 1967, 70, 1619-1624. (39) Tasaka, A.; Akimura, A.; Taguchi, T.; Mimoto, A. Kagaku Kaishi 1989, 2, 174-180.
Received for review March 6,1990. Revised manuscript received October 15, 1990. Accepted October 19, 1990.
Geostatistical Modeling of Dissolved Oxygen Distribution in Estuarine Systems Donato Posat and Mario E. Rossi”8t IRMA, CNR, Universiti degli Studi di Bari, 70125 Bari, Italy, and Fluor Daniel Inc., 10 Twin Dolphin Drive, Redwood City, California 94065
The probabilistic random function model is used in geostatistics to infer the attribute under study a t unsampled locations. In many earth sciences applications, sparse information and incomplete sampling are not uncommon. This problem sometimes can be approached with an appropriate definition of the stationary random function model being used. In spatial statistics this decision is usually influenced by practical factors, such as availability and distribution of the samples. This paper describes a procedure that trades the third dimension for a larger number of samples in the two-dimensional space. This allows for a better characterization of the distribution of dissolved oxygen within the Mar Piccolo of Taranto, Italy. Introduction
To accurately assess the risk and plan the remediation of contaminated estuarine systems, it is necessary to model the spatial distribution of the relevant physical and chemical variables within that system. More often than not, incomplete knowledge of all the processes involved prevents the analyst from using deterministic models. In addition, there are usually too few and too sparse samples to aid in the modeling process. In such cases, a probabilistic model that takes into account the uncertainty associated with the estimates is needed. Geostatistics is the name given to a toolbox of statistical techniques that deal with data correlated in space. This relatively new branch of statistics has gained widespread acceptance within the mining, petroleum, and, lately, hydrogeology fields. Geostatistical techniques are based on the probabilistic concept of a random function model ( I ) , built under the following assumptions: (i) the data are statistically hodegli Studi di Bari. *Fluor Daniel, Inc. t Universitd
474
Environ. Sci. Technol., Vol. 25, No. 3, 1991
mogeneous (stationary), and (ii) the statistics inferred from the samples are representative of the global population. The spatial dependence between two variables at different locations within the area of interest is inferred and modeled from the sample data. A correlation function (either a variogram, covariance, or correlogram) is expressed as a function of the distance between samples. “Inference” of the correlation function implies obtaining from the samples a few experimental points, which are then modeled with an analytical equation. This mathematical model thus measures how “close” two samples are in space and uses this information to interpolate values a t unknown locations by use of a generalized least-squares algorithm ( 2 ) called kriging (ref 3, pp 303-4431, In practice, and as is usually the case in any data-based methodology, geostatistical analysis is adversely affected by lack of information. In particular, obtaining an experimental variogram that describes the spatial variability of the data can be extremely difficult if there are not enough samples available. This is so because each experimental variogram value requires a minimum number of sample pairs for the variogram function to be stable and statistically significant (ref 3, pp 194). The problems associated with inference of statistical moments of a random function model have been discussed by several authors (see refs 4-6 for environmental applications). Any analysis that requires statistical inference will suffer from sparse information. The study presented in this paper deals with dissolved oxygen (DO) measurements, taken a t 16 stations on the Mar Piccolo of Taranto, Italy, on March 19, 1986, over an area of approximately 3000 m by 3800 m. The basin is made up of two bays, which communicate a t the Punta Penna and Punta Pizzone (Figure 1). The first bay, which is the focus of this study, is connected with the Mar Grande and the Ionian Sea by two channels: Canale Porta
0013-936X/91/0925-0474$02.50/0
0 1991 American Chemical Society
r
1Z.5 2o
r
15 12,s
._
f
2
10
Second Bay
7,s
s
Navigabile
~~~
Prob [Z(x)
PuntauntaPenna Pizzone
2*5
0 -20S
-5
Mar Grande -5
0
5
10
15
20 25 Eosting
30
35
40
45
Figure 1. The Mar Piccolo of Taranto.
Taranto
fc/ I+
+ + + + + + + + + + + + + + + +
+
L
2*2
-0.2 Navigabile
Mar Grande -5 -5 - 2 * 3 Oe4 3.1
5.8
The set of samples z(x) spread over A is called a regionalized variable (ref 1 and pp 27-29 in ref 3). Each regionalized variable is interpreted as a realization of a random variable Z(x),with stationary univariate cumulative density function (cdf)
8Q5 1102 13.9 16.6 19+3 Eos+ing
Figure 2. First bay (Primo Seno) of the Mar Piccolo.
Napoli and Canale Navigabile. The bay has an ellipsoidal shape and is fairly shallow, with a major axis of -4 km, a minor axis of -3 km, and a maximum depth of -13 m. Though the depth a t each station varies, there are some 20-25 DO measurements per station, with samples down to 10-m depth (7). This project was prompted by concern about ecological conditions wtihin the bay. Parenzan (8)published a detailed study on Mar Piccolo, after a mussel farm located in the area of the Canale Navigabile (Figure 2) was contaminated. Soon after, the Istituto Tassalografico-C.N.R. of Taranto and the Istituto per Richerche di Matematica Applicata-C.N.R. of Bari were involved in the research. The objective is to understand the dynamics of the dominant physical processes within Mar Piccolo. This will be used to characterize the degree of contamination within the basin. Other variables, besides DO, were sampled as well: salinity, temperature, pH, and current velocity. The case study presented describes the customized statistical techniques and methodology applied to compensate for data shortage. The existence of a larger number of samples in the vertical direction, the uncertainty associated with the exact coordinates of the station locations, and the high continuity of the variable under study have all been exploited to obtain more information in the x-y plane. It is assumed here that it is sufficient to describe the DO field within Mar Piccolo in several twodimensional planes, located at different depths within the bay. A full three-dimensional study is not possible with the information available. Stationarity Decision
Consider the DO value z(x),x E A , where x is the location coordinate vector and A represents the domain of interest.
z] = F ( z )
The set of all the random variables Z(x) over A is called the random function [Z(x),xE A] or, more simply, Z(x). The upper case Z(x) is reserved throughout this paper for both the random variable and the random function, while the lower case z(x) is used to designate the regionalized variable, i.e., the sample value a t location x. The random function Z(x) is said to be strictly stationary if its multivariate distribution [Z(x,),Z(x2),..., Z(x,)] is invariant under translation. Note that this is a property of the model, not of the real samples. Stationarity is not an intrinsic property of the data, but rather a fundamental statement about the model. Geostatistics is based on the idea that the random variables (RVs) Z(x) are dependent. Thus, a spatial correlation model can be inferred. This model is the basic tool for all studies in spatial statistics, and in particular, it allows for the implementation of interpolation-extrapolation techniques (9). The probabilistic random function model requires that significant number of samples be available to infer the characteristics of the variable a t unsampled locations x. The stationary area within which the random function model is applied has to be chosen such that the data are statistically homogeneous and representative of the area A . Despite its rigorous mathematical definition, the stationary decision is almost always influenced by practical considerations, most notably lack of information. In layman's terms, and from a practical standpoint, stationarity can be seen as a decision that allows for statistical inference. The random function Z(x) is said to be stationary of order 2 if its mean E[Z(x)] = m exists and is independent of the location x, and its covariance C(h) exists and depends only on the separation vector h. C(h) is used to describe the average spatial continuity between two RVs Z(x) and Z(x + h) and is defined as
C(h) = E[Z(x + h)Z(x)]- m2 for all x E A
(1)
The covariance is a measure of spatial continuity. The larger the covariance, the less different (on average) the pair of values z(x),z(x h). The corresponding measure of spatial variability is the variogram, defined as
+
2y(h) = E [[Z(X)- Z(X + h)I2]
(2)
C(h)and r(h) are moments of the bivariate distribution, which characterizes the joint variability of Z(x) and Z(x + h). Under the second-order stationarity assumption, they are related by r(h) = C(0) - C(h)
(3)
where C(0) is the covariance value a t h = 0. Model for Dissolved Oxygen
In practice, both the covariance and variogram must be estimated from available samples. A sufficient number of samples is needed to estimate r(h) as an average. 1
N(h)
T(h) = - C [ ~ ( x+j h) - ~(xj)]' 2N(h) i = l Environ. Sci. Technol., Vol. 25, No. 3, 1991
(4) 475
8.65
5
14.0
c
5
*P: 5
8.8
7.75
6.85
++++++++++
6.4
7.3
*
3,9Q*. %
*
5
0
1
2
3
5
4
6
7
8
9
10
Depth
Figure 3. Vertical profiles found at four stations of the Mar Piccolo. The station with almost constant and low DO values is atypical and corresponds to a station located near submarine springs.
where N(h) is the number of pairs separated a vector distance h. An empirical rule states that -30 pairs is adequate in order to obtain a reliable estimate (ref 3, p 194). Any two samples in space are seldom or never exactly distant a vector h; thus, angle and distance tolerances are defined. These tolerances are dependent on the amount of experimental information available within A : the less samples within A , the larger the tolerances needed to estimate y(h). The experimental function obtained for different h is then modeled with a closed analytical and admissible function (10, 11). An example of the flexibility a user has in defining a statistically homogeneous sample set to deal with is presented hereafter. The 16 stations sampled for DO in Mar Piccolo are not sufficient for variography analysis in the horizontal plane, as was shown in Rossi and Posa (6). The area under study has been gridded with a 200 m by 200 m mesh. Inside 16 of these cells a sampling station was established, but the exact ( x , y ) coordinates of the stations are not known. It is known, however, that the station is located inside the cell, and arbitrarily, each station has been located at the center of its cell (7). The uncertainty involved in the location of the stations allows for a solution to the inference problem. The third dimension, present in all 16 stations, can be “traded in” for additional information in the horizontal plane. Four different layers were defined: layer 1, between 0- and 2-m depth; layer 2, between 2.1- and 4-m depth; layer 3, between 4.1- and 6-m depth; and layer 4, between 6.1- and 8-m depth. A n y sample falling within these layers was relocated to t h e horizontal reference p l a n e , with elevation equal to the mean d e p t h of each layer, and uniformly distributed over its 200 m by 200 m cell. After
the relocation procedure was completed, the number of samples available at each level increased to -880-90 values. Remarks. 1. The procedure described assumes that the short-distance variability of DO is isotropic, Le., is similar in horizontal and vertical directions. This isotropy assumption is commonly made in geostatistical applications, since there is usually little information available on horizontal short-scale variability. 2. The relocation procedure is acceptable if the variable is highly continuous, Le., if it varies smoothly in space. Because there is high intermixing of DO within the bay due to currents and other natural phenomena, the distribution of DO is very continuous. The coefficient of variation (CV) (defined as a i m ) is a good indication of the 476
*
5
*5
++++++
5.95 5.5
t I-
7.3
Environ. Sci. Technol., Vol. 25,
No. 3, 1991
l
0.5
0.5
2.6
c
4
*5
5 c
1
** ,I
4.6
I
6.7
I
8.7
1
I
I
I
10.8 12.8 14.9 16.9 19.0
Figure 4. Posting map of the relocated data. Note the high degree of clustering.
intrinsic variability of the attribute. Here CV = 0.09 (9%), which is very low if compared to the CV of most other variables studied in geostatistics, where CVs of the order of 50-100% are common. 3. Figure 3 shows a group of vertical profiles for 4 of the 16 stations within Mar Piccolo. It is possible to choose the thickness of the layer according to the relative variation of the variable within that depth. Also, layers need not be of the same thickness. 4. The layer thickness determines the number of samples whose coordinates are transformed. The thicker the layer, the more samples are available for data analysis. On the other hand, the thinner the layer, the smaller the approximation introduced. 5 . Since the exact location of the station is not known within its 200 m by 200 m cell, the associated uncertainty minimizes the consequence of the sample’s relocation within the cell. 6. The kriging algorithm can be applied now on a two-dimensional grid, at four different horizontal planes, a t depths of 1, 3, 5, and 7 m. 7. This procedure yields a “transformed” set of samples that is spatially clustered. In practice, clustered samples introduce biases, since they are not representative of the entire area. In particular, inference of moments of the random function model becomes complicated. Clustering poses more serious problems when analyzing variables with a high coefficient of variation. Clustering. Detailed discussion of the effect of spatial aggregation of samples in geostatistical studies (12) is beyond the scope of this paper. However, some relevant remarks can be made. The degree of clustering observed after the relocation procedure described is largely dependent on the following: the thickness of the layer chosen, the overall number of stations in the study area, and the size of the cell where the samples are redistributed. Figure 4 gives a posting map of relocated stations in layer 1 (mean depth 1 m) and shows the severity of clustering. The experimental variogram function is known to be greatly affected by clustering. There exist other estimators for spatial correlation, more robust with respect to clustering; see Parker and Srivastava (13). Refer to the following section and Figures 7 and 8 for a description of an improved estimator. Exploratory Data Analysis and Structural Analysis
All the figures presented hereafter correspond to layer 1. A similar study was carried out for each of the other
three layers defined.
47.9
N
dm min qo.25
94 9.489 0.975 0.104 6.200 9.400
q0.w
9.700
q0.75
9*900 10.700
m o2
max
'Class width = 0.600
0.600 2.100 3.600 5.100 6.600 8.100 9.600 11.100 12.600
DO Figure 5. Histogram and univariate statistics of the data corresponding to layer 1.
::::
J"i
11.9 10.0 8.1
*
6.2
*
* * *
1 .o
*
* * *
*
*
4.3
*
+
0
0.0
1010
5.0
0.5
2.8
5.1
7.4
9.8
12.1
14.4
16.7
19.0
Figure 6. Contour map of the initial information available in layer 1. Contour interval 0.4 mg/L.
The initial st,ep of any case study is an exploratory analysis of the information available. Traditional statistical tools are used in this step to describe the data: histograms, univariate statistics, normal probability plots, contour plots, scatter plots, etc. Figure 5 presents the univariate statistics and histogram of the relocated data in layer 1, and Figure 6 gives the corresponding contour map using a traditional inverse square distance interpolator. Such map allows for a fast qualitative appreciation of the spatial distribution of the data. The trend model adopted later for kriging is partly based on this preliminary information. Note in Figure 5 the outlier low 110 values, also evident in the north area of the contour map of Figure 6. These correspond to a station close to submarine freshwater springs, as described in Parenzan (8). Figure 7 presents the experimental variograms in two directions of layer 1. The clustering effect is quite evident in the destructured variogram. The directions shown are north-south (NS) and east-west (EW). An improved variogram estimator, called "nonergodic" variogram, is known to be more resistant to clustering (14). Figure 8 shows the nonergodic variograms for the same layer, and for the same two directions NS and EW. A structure is now evident. These experimental variograms are then modeled and used to inject the correlation model
15:O
20:o
Lag D i s t a n c e
I I I I I I I I I I I I
Figure 7. Experimental variograms for two main directions of anisotropy of layer 1, E-W and N-S. Note the lack of structure. 3.0
c)
3
****e
*****
E-W N-S
Non-ergodic vario Non-ergodic vario
*
2.0
E E o
c3
4
*
1 .o
*
,++
*
0.0
0.0
10.0
20.0
30.0
Lag D i s t a n c e Flgure 8. Nonergodic experimental variograms for the two main directions of layer 1. The solid line corresponds to the model retained. Note that the variogram in the direction of the trend (N-S) is more erratic. Also, compare with Figure 7.
into the interpolation algorithm. A detailed discussion on variogram modeling can be found, for example, in Journel and Huijbregts (ref 3, pp 148-235) and Isaaks and Srivastava (ref 11, pp 369-399). The anisotropies observed in the variograms are conEnviron. Sci. Technol., Vol. 25, No. 3, 1991 477
16.5 15.5 14.5 13.5 12.5 11.5 10.5 ~
9.5
S I
.-
5
8.5
3
1.5
i
z 6.5 5.5 4.5 3.5 2.5 1.5 0.5 0.5
' 5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 12.5 13.5 14.5 15.5 : 6 . 5 17.5 18.5
tasting Figure 9. Contour map of the DO model, layer 1. Contour interval 0.2 mg/L.
sidered negligible, and four slightly different spherical isotropic variogram models were adopted: layer 1: r l ( h ) = C1 + 0.85SphU,=,,,(h) C, = 0.3 layer 2: rz(h)= C2 + 0.90Sph,,=,,o(h)
Cz = 0.01
layer 3: ~ ~ (= hC3)+ 0.95Sph,3=,,o(h) C3 = 0.02 layer 4: -y4(h)= C4 + 0.82Sph,,=,,0(h)
known coefficients. R(x) is a stochastic stationary residual, with zero mean: E[R(x)] = 0. The trend model adopted for the distribution of DO is linear:
m(x)= a.
+ alx + a,y
The kriging estimate is written
C, = 0.08
(5)
with
Sph,(h) = 1.5h/a - 0 . 5 ( h / ~ ) if~ h I a The function Sph,(h) = 1 if h > a. The Ci,i = 1, ..., 4 are the nugget constants corresponding to the white-noise components, and the ai,i = 1, ..., 4 are the ranges of the models, Le., the distance beyond which spatial correlation vanishes. Note that the structural model of the deepest layer has the largest range of the four models. Kriging with a Trend Model
where m(x) is a deterministic trend of known form: L
CaJi(x) 1=0
Here the fi's are specified functions and the ai's are un478
strained normal equations: n
After the variograms are modeled, a kriging with a trend model algorithm [commonly known as universal kriging (UK)] is implemented for each layer. The kriging algorithm is briefly described hereafter. For a more detailed presentation, refer, for example, to Journel and Huijbregts (ref 3, pp 313-320). The UK algorithm assumes a two-component model Z(x) = m ( x ) + R(x)
m(x) = E[Z(x)l =
It is evident from the eq 5 that kriging is an interpolation technique that uses a weighted linear combination of the available samples. This type of estimator is common in interpolation algorithms. Kriging differs, however, in the way the weights A, are determined. A Lagrange formalism is used to minimize the error - Z(x)] under the unbiasedness convariance Var [Z*(x) dition E[Z*(x) - Z(x) = 01, resulting in a system of con-
Environ. Sci. Technol., Vol. 25, No. 3, 1991
c X p Y d X , - Xp)
-
Po
-
PlX,
-
PZYa =
p=1
Y R b - x,) ( a = I , ..., n )
n
CA,
= 1
p=1
n
cx,x,
=x
B=1 n
CAaYp = Y
p=1
where -yR represents the residual (detrended) variogram (ref 3, pp 315-316). The @!'sare the three Lagrangian
16.5 15.5 14.5 13.5 12.5 11.5 10.5
o,
9.5
.-r
0.5
* i
0 Z
7.5 6.5 5.5 4.5 3.5 2.5 1.5 0 5
0.5
1.5
2.5
3.5
4.5
5.5
6.5
7.5
8.5
9.5 1 0 . 5 1 1 . 5 1 2 . 5 1 3 . 5 1 4 . 5 1 5 . 5 16.5 1 7 . 5 1 8 . 5
Easting Figure 10. Contour map of the DO model, layer 2 . Contour interval 0.2 mglL. 16.5 15.5 14.5 13.5 12.5 11.5 10.5
o,
.-c
9.5 8.5
u i
0 7.5
Z
6.5 5.5 4.5 3.5 2.5 1.5
0.5
1.5
2.5
3.5
4.5
5.5
6.5
7.5
8.5
9.5 10.5 11.5 12.5 13.5 14.5 15.5 16.5 17.5 18.5
Easting Figure 11. Contour map of the DO model, laver 3. Contour interval 0.2 ma/L. I I I U I L I ~ I I b~ Icull
espuiiuirig
LU
LIW
Lmee unuiaseaness con-
straints introduced, and x and y are the coordinates of the samples, chosen as the functions fL's.
Bystem 0 is soivea for eacn point or the grid being interpolated, using neighboring samples, found according to some predefined search strategy. Fnvirnn Sei Taehnnl
Vnl 3 5 Nn ?
19QI
A70
16.5 15.5
14.5 13.5 12.5 11.5 10.5 9.5
c .-I
8.5
u i
0 7.5
Z 6.5 5.5 4.5 3.5 2.5 1.5 0.5 0.5
1.5
2.5
3.5
4.5
5.5
6.5
8.5
7.5
9.5
10.5 11.5 12.5 ‘3.5 14.5 15.5 16.5 17.5 18.5
Easting Figure 12. Contour map of the DO model, layer 4. Contour interval 0.2 mg/L.
Results and Conclusions
The contour maps of the DO kriging estimates for each layer are presented in Figures 9-12. The spring area (8) in the north side of the basin is quite marked, as well as a previously unnoticed NS channel of low values in the west side of the basin. It is interesting to note that this two-dimensional interpolation provides a more “rugged image of the variable, due to the smaller ranges of the 2-D variograms, as compared to the 3-D model described in Posa (7). The trend, although always present, is more evident in the bottom layer. This last remark agrees with the prevailing hypothesis that the trend in DO values is mostly caused by the submarine springs, with tidal currents having only minor influence. These tidal currents are likely to produce the NS channel of low DO values in the west side of the basin, particularly for the shallow layers. This is expected, since deeper water currents are depleted in DO (15). Another feature of Figures 9-12 is a trough of very low DO values, near the center of the basin (whose approximate coordinates are 3c = 5 and y = 7). This “hole” in the DO surface, although located within the channel mentioned above, cannot be explained yet by any known natural phenomena. It is very marked in the upper two layers, there is only a hint of it in Figure 11 (layer 3), and it is missing in Figure 12 (layer 4). Interestingly, the port and city of Taranto are located nearby (recall Figures 1 and 2). This feature calls for further investigation. The case study introduces a relocation procedure that enhances the inference of correlation of the data. The three-dimensional information available is traded for better inference in the horizontal plane. This, coupled with a correlation estimator that is resistant to the clustering introduced, allows for construction of a two-dimensional dissolved oxygen model, using universal kriging (or kriging with a trend model), a t selected depths. This relocation 480
Environ. Sci. Technol., Vol. 25, No. 3, 1991
procedure proves useful in cases where there is very limited information in the horizontal plane, the variable under study is highly continuous, and the exact locations of the stations are not known. Acknowledgments
The Dipartimento di Matemdtica della Universiti degli Studi of Bari, IRMA-CNR of Bari, Prof. V. Capasso of IRMA-CNR, and Prof. A. G. Journel of the Department of Applied Earth Sciences, Stanford University, are gratefully acknowledged for supporting the collaboration of the authors. Istituto Tassalografico-CNR of Taranto, Italy, is acknowledged for providing the data set. Thanks are also due to both reviewers. Their comments and suggestions significantly improved the original manuscript. Literature Cited (1) Matheron, G. The Theory of Regionalized Variables and Its Applications; Les Cahiers du Centre de Morphologie Mathgmatique, Fasc. 5 ; CG: Fontainebleau, France, 1971. (2) Goldberger, G. Best Linear Unbiased Prediction in the Generalized Linear Regression Model. J. Am. Stat. Assoc. 1962,57, 369-375. (3) Journel, A. G.; Huijbregts, Ch. Mining Geostatistics; Academic Press: New York, 1978. (4) Matgrn, B. Spatial Variation;Lecture Notes in Statistics 36; Springer Verlag: New York, 1960. ( 5 ) Myers, D. E. To be or not to be ... stationary. Math. Geol. 1989, 21, 347-362. (6) Rossi, M.; Posa, D. 3-D Mapping of Dissolved Oxygen in Mar Piccolo. J. Enuiron. Geol. Water Sci. 1990,16, 209-219. ( 7 ) Posa, D. Variography Analysis for Anisotropies Identification of Dissolved Oxygen Measurements. Bollettino di Oceanologia TeBrica ed Applicata 1988, 6, 4. (8) Parenzan, P. I1 Mar Piccolo di Taranto. Report for the Camera di Comercio, Bari, Italy, 1984. (9) Ripley, B. D. Spatial Stetistics; J. Wiley and Sons: New York, 1981.
Environ. Sci. Technol. 1991, 25, 481-488
Christakos, G. On the Problem of Permissible Covariances and Variogram Models. Water Resour. Res. 1984, 20, 251-265. Isaaks, E. H.; Srivastava, R. M. A n Zntroduction t o Geostatistics; Oxford University Press: New York, 1989. Rossi, M. Impact of Spatial Clustering on Geostatistical Analysis. M.Sc. Thesis, Dept. of Applied Earth Sciences, Stanford University, Stanford, CA, 1988. Parker, H.; Srivastava, R. M. Relative Variograms and Robust Spatial Continuity Measures; Third International Congress of' Geostatistics Avignon, France, 1988; Vol. 2.
(14) Isaaks, E. H.; Srivastava, R. M. Spatial Continuity Measures for Probabilistic and Deterministic Geostatistics. Math. Geol. 1988, 20, 313-341.
(15) Broecker, W. S. Chemical Oceanography; HBJ Inc.: New York, 1974.
Received for review November 9, 1989. Revised manuscript received May 2 1 , 1990. Accepted October 23, 1990. Financial support was provided to D.P. by a grant from the McGee Fund, School of Earth Sciences, Stanford University.
Metal( II)Ion Binding onto Chelating Exchangers with Nitrogen Donor Atoms: Some New Observations and Related Implications Arup K. Sengupta," Yuewei Zhu, and Diane Hauze
Fritz Engineering Laboratory, Environmental Engineering Program, Lehigh University, Bethlehem, Pennsylvania 18015 Chelating exchangers with only nitrogen donor atoms exhibit some unusual properties in relation to metal ion uptake from the aqueous phase. Specifically, for these chelating exchangers, (a) the metal ion uptake increases with an increase in competing Ca2+and Na+ ion concentrations; (b) metal removal capacities remain almost unaffected a t pH as low as 1.0; (c) regeneration or desorption of metal ions with ammonia is very efficient while acid regeneration is ineffective; and (d) both metal cations and other anions can be removed simultaneously from the aqueous phase. The foregoing behaviors are not observed a t all for other commonly used chelating exchangers with carboxylate, iminodiacetate, or aminophosphonate functionalities. A molar-exchange model, where the individual nitrogen donor atoms bind metal ions independently on molar basis, can explain the anomalous characteristics of these specialty chelating exchangers. Because of their unusual properties, these chelating exchangers with only nitrogen donor atoms offer new application opportunities for heavy-metals removal unattainable through other exchangers. ~~
Introduction Chelating Exchangers: General Background and Metal Ion Selectivity. Chelating exchangers are, in general, coordinating copolymers with covalently bound side chains, which contain one or multiple donor atoms (Lewis bases) that can form coordinate bonds with most of the toxic metal ions (Lewis acids). Due to coordination-type interactions, all such chelating exchangers offer extremely high selectivity toward commonly encountered toxic M(I1) cations, namely, Cu2+,Pb2+,Ni2+,Cd2+,and Zn2+,over competing alkaline (Na+,K+) and alkaline-earth (Ca2+,Mg2+)metal cations. Depending on the number of donor atoms present in a pendant ligand of the polymer, the repeating functional groups are often referred to as mono-, bi-, or polydentate. In the recent past, a great deal of interest has been observed in relation to the applicability of these chelating polymers for removal, separation, and purification of metal ions from heavy-metal-contaminated water, wastewater, and solid wastes (2-6). State of the art reviews on synthesis and properties of a great number of chelating polymers including the recently developed ones have been provided by Hudson (7) and Warshawsky (8). Table I provides functionalities and other salient information for five different chelating exchangers used in this study, of which the first three (iminodiacetate, car0013-936X/91/0925-0481$02.50/0
boxylate, aminophosphonate) are quite popular and traditionally used in many industrial applications (3,5).The last two chelating exchangers (Dow Chemical, Midland, MI) have only nitrogen donor atoms and are weakly basic; on protonation, they develop fixed positive charges. For convenience, the exchanger with three N atoms and the one with two N atoms per functional group in Table I will be referred to as DOW 3N and DOW 2N, respectively. Figure 1provides experimentally determined Cu(II)/Ca and Ni(II)/Ca separation factors a t pH 4.0 for three different chelating exchangers, and the high selectivities of the metal ions can be readily noted ( 4 ) . Separation factor is a dimensionless measure of relative selectivity between two competing ions and, in this case, equal to the ratio of distribution coefficient of the metal(I1) concentration between the exchanger phase and the aqueous phase to that of calcium ion and is given as follows: separation factor
[m]
aMICa =
[RM][ ~ a 2 + 1 [M2+][RCa]
(1)
where and [M2+] correspond to exchanger- and aqueous-phase concentrations of M(II), respectively. It is true that Coulombic and hydrophobic interactions are present in such chelating exchange processes, but their roles are relatively insignificant compared to Lewis acidbase interaction in contributing toward such high metal ion selectivities. Another generic similarity among all chelating exchangers lies in the fact that these functionalities are weak acid or weak base, and therefore, all of them exhibit high affinity toward hydrogen ions. Differences in metal ion selectivities among different chelating exchangers in Figure 1 are closely related to the metalligand stability constants of different functionalities; this aspect is, however, outside the central objective of this article. Metal Ion Binding onto Chelating Polymers and Premises of the Study. In a chelating polymer, the functional groups with the donor atoms are rigidly (covalently) bound to repeating monomers (like styrene) which, again, are fixed as part of a three-dimensional network cross-linked through divinylbenzene. As a result, the ligands in the polymer phase experience considerable strains to orient themselves spatially around the receptor metal ions. This strain may not allow the individual functionalities in the polymer phase to reproduce their aqueousphase metal-ligand configurations. Extensive experimental results are available for more widely used chelating
0 1991 American Chemical Society
Environ. Sci. Technol., Vol. 25, No. 3, 1991 481
