Health risks from arsenic-contaminated soil in Flin Flon–Creighton, Canada: Integrating geostatistical simulation and dose–response model

Environmental Pollution 157 (2009) 2413–2420

Contents lists available at ScienceDirect

Environmental Pollution journal homepage: www.elsevier.com/locate/envpol

Health risks from arsenic-contaminated soil in Flin Flon–Creighton, Canada: Integrating geostatistical simulation and dose–response model Hua Zhang a, Guo-he Huang a, *, Guang-ming Zeng b a b

Environmental Systems Engineering Program, Faculty of Engineering, University of Regina, Regina, Saskatchewan, Canada S4S 0A2 College of Environmental Science and Engineering, Hunan University, Changsha, Hunan 410082, PR China

Cancer risks induced by arsenic in soil were evaluated using geostatistical simulation and dose–response model.

a r t i c l e i n f o

a b s t r a c t

Article history: Received 23 January 2009 Received in revised form 9 March 2009 Accepted 10 March 2009

Elevated concentrations of arsenic were detected in surface soils adjacent to a smelting complex in northern Canada. We evaluated the cancer risks caused by exposure to arsenic in two communities through combining geostatistical simulation with demographic data and dose–response models in a framework. Distribution of arsenic was first estimated using geostatistical circulant-embedding simulation method. We then evaluated the exposures from inadvertent ingestion, inhalation and dermal contact. Risks of skin caner and three internal cancers were estimated at both grid scale and census-unit scale using parametric dose–response models. Results indicated that local residents could face non-negligible cancer risks (skin cancer and liver cancer mainly). Uncertainties of risk estimates were discussed from the aspects of arsenic concentrations, exposed population and dose–response model. Reducing uncertainties would require additional soil sampling, epidemic records as well as complementary studies on land use, demographic variation, outdoor activities and bioavailability of arsenic. Ó 2009 Elsevier Ltd. All rights reserved.

Keywords: Geostatistics Stochastic simulation Arsenic Risk assessment Cancer

1. Introduction Arsenic occurs naturally in soils. Reported background concentrations in Canada range from 4.8 to 13.6 mg/kg, and average urban and agricultural soil concentrations are mostly between 4 and 6 mg/kg (Wang and Mulligan, 2006). Arseniccontaminated soil has long been recognized as a threat to human health (Tang et al., 2007). Inorganic arsenic is the only known carcinogen for which there is adequate evidence of carcinogenic risk by both inhalation and ingestion (Kapaj et al., 2006). It has been on the first Priority Substances List 1 under the Canadian Environmental Protection Act (CEPA) since 1989. Elevated levels of arsenic have been observed in the vicinity of anthropogenic sources such as smelter and base-metal refinery facilities (Newhook et al., 2003; Schelle et al., 2008). Human health risks caused by exposure to arsenic are of increasing concerns to the regulatory agency and communities in these areas. In 2006, elevated concentrations of arsenic were detected in surface soils in the City of Flin Flon (Manitoba) and the Town of Creighton (Saskatchewan), which are adjacent to a smelter complex of the Hudson Bay Mining and Smelting Company Ltd (HBM&S).

* Corresponding author. Tel.: þ1 306 5854095; fax: þ1 306 5854855. E-mail address: [email protected] (G.-he Huang). 0269-7491/$ – see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.envpol.2009.03.014

Concentrations of arsenic at 70 of the 106 sites (66%) were above 12 mg/kg, the national human health protection guideline for residential/parkland soil (Manitoba Conservation, 2007). This base-metal mining and smelting complex has been processing sulphide ores from local mines since 1930. Although particulate emission from smelter stacks was reduced by 90% in 1995 due to a technological improvement, the eight-decade emission has resulted in elevated concentrations of heavy metals (including arsenic) in adjacent soils (Henderson et al., 1998; McMartin et al., 1999), posing a potential health threaten to the neighboring 7000 residents. Following a soil study of Manitoba Conservation in 2006, a human health risk assessment conducted by HBM&S is underway to determine whether exposures to the contaminated soil present unacceptable, long-term health risks to residents and visitors. Estimating the spatial distribution of pollutant is crucial for quantifying the risk level at different locations (Chen et al., 1998; He et al., 2008; Qin et al., 2008). The accuracy of the direct analysis of in situ data is often dubious because the field investigation is limited by time and cost and the observations contain considerable uncertainty (Lee et al., 2007a). Spatial interpolation methods such as kriging have thus been widely used for estimating environmental conditions and concentrations of pollutants at unsampled locations (Yang et al., 2004; Chen et al., 2008; Shi et al., 2008). Kriging estimates honor the known points, but do

2414

H. Zhang et al. / Environmental Pollution 157 (2009) 2413–2420

not fully reflect the local variability. The inherent smoothing effect of kriging methods tends to overestimate small values and to underestimate large values. This effect is not spatially uniform because it depends on the local data configuration only. It was suggested that interpolated maps produced by kriging should not be used for situations sensitive to the presence of extreme values (Goovaerts, 1997). Therefore, rather than producing a unique map of the best unbiased prediction, stochastic simulations have recently been used to generate equally probable realizations. Each realization reproduces the same sample statistics (histogram) and spatial continuity (semivariogram), and matches the conditional data. Fluctuations between realizations provide a quantitative measure of uncertainties in both the natural phenomena and the field measurement. These realizations can serve as inputs to risk models for producing a distribution of plausible answers, enabling the decision maker to see the average, maximum and minimum likely risk (Gay and Korre, 2006). There is a variety of stochastic simulation methods in geostatistics (Goovaerts, 1997; Webster and Oliver, 2007). Due to the computational advantage, methods using the sequential simulation paradigm have become the most popular approaches for estimating the distributions of pollutants in soils (Gay and Korre, 2006; Bourennane et al., 2007; Goovaerts et al., 2008; Verstraete and Van Meirvenne, 2008) as well as in groundwater (Cinnirella et al., 2005; Lee et al., 2008). However, the performance of sequential geostatistical simulation is not obviously superior to other simulation methods (Lin et al., 2001; Lee et al., 2007b). This is possibly because of the inherent moving-neighborhood process. Because the kriging systems increase dramatically when more and more previously simulated values are involved following the sequential diagram, in practice, the conditioning values are reduced to only the nearest points around the simulated node (within a moving neighborhood). In general, the moving neighborhood entails a loss of accuracy because the screening effect of the closest values is partial; furthermore, the errors in each realization tend to spread when the simulation proceeds since previously simulated values are used in subsequent steps (Emery, 2004). Geostatistical simulation has been rarely used in human health risk assessment. Two exemptions were dealing with carcinogenic risk of ingesting arsenic in smelt fish in Taiwan (Lee et al., 2008) and unspecified health risk from contaminated land in the United Kingdom (Gay and Korre, 2006), respectively. Both of them, however, used sequential Gaussian simulation without consideration on its limitations. Furthermore, none of these studies used dose–response models to quantify the health risks. Quantitative estimates of cancer risk can be generated through either nonparametric methods (e.g. slope factor, unit risk, concentration at specific risk level) or parametric methods such as dose–response models. While the non-parametric methods are easy to use, dose– response models are more effective in reflecting the complex relationship between body characteristics (e.g. age and gender) and cancer vulnerability. Therefore, the first objective of this study is to characterize the distribution of arsenic concentrations in soil using a geostatistical circulant-embedding simulation (GCS) method. Probabilistic distributions of concentration are generated for each node of a 50 m grid covering Flin Flon and Creighton. The second objective is to assess carcinogenic risks due to exposure to soil arsenic using parametric dose–response models. Risks of skin cancer, lung cancer, bladder cancer and liver cancer are estimated at both grid scale and census-unit scale. These goals are achieved through combining the geostatistical simulation method with demographic data and dose–response models in a framework.

2. Materials and methods 2.1. Study area and field data Flin Flon and Creighton (54 460 N and 101 510 W) share the boundary between the Province of Saskatchewan and the Province of Manitoba, Canada. The climate is continental, characterized by cold winters and relatively warm summers, with an annual precipitation of 477.9 mm; wind directions and velocities are fairly well distributed but predominant to the southeast and southwest (McMartin et al., 1999). Flin Flon is immediately adjacent to the smelter with a population of 5502, and Creighton is located 2 km southwest of the smelter with 1594 residents. Uninhabited land (including major water bodies) is excluded in this study through interpreting residential area on QuickBird satellite images and modifying the boundaries of 16 census Dissemination Areas (DAs) defined by Statistics Canada. The study area (12.36 km2) is represented as a grid with 50 m spacing as shown in Fig. 1. Soil was sampled at 91 sites in Flin Flon and 13 sites in Creighton, including playgrounds, boulevards, school yards, vacant lots, and undeveloped land parcels. The majority of sites were located less that 3 km from the smelter complex. A site at the Bakers Narrows Provincial Park (17.8 km southeast of the smelter) and another at the Cranberry Portage Park (38.3 km southeast of the smelter) were selected as low impact and background sites, respectively. Three replicate samples were collected at each site during August 21–24, 2006. Chemical analysis of arsenic was performed using inductively coupled plasma mass spectrophotometry (ICP/MS) on a strong-acid digested sub-sample. Arsenic concentrations of soil samples in Flin Flon and Creighton ranged from 2.3 to 407 mg/kg; low impact and background concentrations were 6.1 and 2.5 mg/kg, respectively. Details of sampling and laboratory analysis are referred to the survey report by Manitoba Conservation (2007).

2.2. Geostatistical simulation Two assumptions behind the GCS method are that (1) the data is stationary over the spatial domain (in terms of mean, variance, and semivariogram) and (2) the natural variations can be modeled as a multivariate Gaussian function. While local variations are partly lost in kriging process due to the smoothing effect, GCS adds this kind of information to kriging estimation. The added variation for a specific location is a conditional distribution defined as a mean of zero and a variance equal to the kriging variance at the same location. Simulated realizations are generated by randomly sampling from these conditional distributions at all grid nodes under the spatial correlation context defined by a preliminary semivariogram model. For the whole spatial domain, each realization is a possible image of the natural reality of interest; at each location, values from a number of realizations follow a Gaussian distribution, with a mean equal to the kriging estimation and a spread given by the kriging variance. The simulation proceeds in the following steps: (1) Decluster the concentration dataset using the cell method in the case of preferential sampling. (2) Transform the clustered dataset into a standard normal distribution using Gaussian anamorphosis if the dataset is positively skewed. (3) Establish a semivariogram model using the normalized dataset. (4) Construct a covariance map on the simulation grid using the semivariogram model and perform convolution on the covariance map using fast Fourier transform (FFT). (5) Generate a random grid which is a standard Gaussian realization (with whitenoise covariance) on the same spatial domain of the simulation grid and perform convolution on this random grid using FFT. (6) Multiply the FFT results from steps (4) and (5) and perform an inverse FFT to constitute an unconditional realization. (7) Condition the unconditional realization via kriging (which is performed for the whole spatial domain in one pass). (8) Back-transform the result of conditional realization into the original scale. A new realization can be generated by repeating steps (5)–(7) with a new random grid. Specifically, the process of steps (4)–(6) is known as simulation via circulant embedding of the covariance matrix; mathematical details are referred to Kyriakidis and Yoo (2005) and Dietrich and Newsam (1997). Comparing to the sequential Gaussian simulation, the distinct feature of GCS is that the kriging configuration is fixed and includes only the original data. The conditional kriging needs to be performed only once for all the realizations to reproduce the original data, while the sequential method has to achieve one kriging per realization unless the same visiting sequence is used (Emery, 2004). Therefore, the GCS method does not suffer from forced simplification such as the moving neighborhood, and requires substantially less computation time with respect to the sequential paradigm. The soil survey in the study area has an emphasis on downtown and high-density residential areas of Flin Flon (with the average distance between two sampling sites less than 250 m). The sample histogram exhibits strongly positive skewness. Therefore both histogram clustering and normal-score transformation are necessary. Totally

H. Zhang et al. / Environmental Pollution 157 (2009) 2413–2420

2415

Fig. 1. Study area and field data (concentration unit: mg/kg).

100 realizations were obtained in the Geostatistical Analyst of ArcGIS 9.3. Each of them generated a map of arsenic concentration for the exposure evaluation.

specific for each group (Health Canada, 2004). The total daily intake of arsenic in soil is the sum of ADIoral , ADIinh , and ADIskin .

2.3. Exposure evaluation

2.4. Risk assessment

Three exposure pathways are considered: (1) inadvertent ingestion of contaminated soil, (2) inhalation of contaminated soil particles, and (3) dermal contact with contaminated soil. The average daily intake (ADI, mg/kg/day) of each pathway is calculated using the dose estimation equation recommended by Health Canada (2004) as

Increased mortality from multiple internal organ cancers (liver, kidney, lung, and bladder) and increased incidence of skin cancer were observed in populations exposed to high concentrations of arsenic in environment (Mandal and Suzuki, 2002; USEPA, 2004). The Black-foot Disease endemic area of Taiwan has historically been the principal data source for studying arsenic-induced cancer risks in North American (Kapaj et al., 2006); dose–response models were developed for skin cancer and internal cancers (Brown et al., 1989; USNRC, 1999, 2001). These parametric models are used in this study to predict the potential carcinogenic risks in a northern Canadian context. The prevalence ratio of arsenic-induced skin cancer is estimated as

ADI oral ¼

Cs  IRs  AF git  D1  D2  D3 BW  LE

(1)

ADI inh ¼

Cs  PAir  IRA  AF inh  D1  D2  D3  D4 BW  LE

(2)

Cs  SAH  SLH  AF skin  EF  D1  D2  D3 ¼ BW  LE

(3)

ADI skin

   pðc; tÞ ¼ 1  exp  q1 c þ q2 c2 ðt  mÞk Hðt  mÞ

where ADIoral , ADIinh , and ADIskin denote the average daily intakes of inadvertent ingestion, inhalation and dermal contact respectively, Cs is the concentration of arsenic in soil (mg/kg), IRs is the factor for receptor’s soil ingestion rate (kg/day), PAir is the particulate concentration in air (kg/m3), IRA is the receptor’s air intake (inhalation) rate (m3/h), SAH is the skin surface area exposed (cm2), SLH is the soil loading to exposed skin (kg/cm2 per event), EF is the exposure frequency (event/ day), AFgit , AFinh and AFskin are the unitless factors for gastrointestinal tract, inhalation absorption and dermal absorption respectively, D1 is the exposed days per week/7 days, D2 is the exposed weeks per year/52 weeks, D3 is the total year exposed to site, D4 is the exposed hours per day, BW is the body weight (kg), and LE is the life expectancy (year) for the assessment of carcinogens. The Canadian general population is divided into three age groups, i.e. children (0–11 years), teenager (12–19 years) and adult ( 20 years); values of parameters in Equations (1)–(3) are

(4)

where p(c, t) is the fraction of people with skin cancer, c is the arsenic concentration (mg/kg), t denotes age (year), and H denotes the Heaviside function (i.e. Hðt  mÞ ¼ 0 for t < m and Hðt  mÞ ¼ 1 for t  m); the parameters q1, q2, k, m are nonnegative. The incidence rate of arsenic-induced internal cancer is estimated as   hðc; tÞ ¼ k q1 c þ q2 c2 ðt  mÞk1 Hðt  mÞ

(5)

where h(c, t) denotes an incidence rate and the variables and parameters are defined as in (4). The parameter values of q1, q2, k, m are specific for each cancer and are referred to Yu et al. (2003). Both p(c, t) and h(c, t) are calculated for each gender at each age, and summarized for the entire population based on its demographic structure: pðcÞ ¼

 X X fij pij ðc; tÞ ; i

j

hðcÞ ¼

 X X fij hij ðc; tÞ i

j

(6)

2416

H. Zhang et al. / Environmental Pollution 157 (2009) 2413–2420

Fig. 2. Semivariogram graph (a) and surface (b).

where fij is the fraction of population of gender i at age j to the total population, and pij(c, t) and hij(c, t) denote the risks associated with this sub-population defined by Equations (4) and (5). Population density is assumed to be uniform within a DA and is calculated by dividing the number of residents by the inhabited area. Demographic information is extracted from the database of 2006 Census. We that assume the total cancer risk induced by arsenic in soil is represented by the sum of risks of skin cancer, liver cancer, lung cancer, and bladder cancer. The total cancer risk is firstly determined at each grid node, and point estimates are then averaged for each DA to generate regional risk level that could be readily used by regulatory agency and community.

3. Results The experimental semivariogram was fitted through the regionalization of an anisotropic Gaussian model (major range ¼ 7067.42 m, minor range ¼ 3851.23 m, direction ¼ 352.15 , partial sill ¼ 0.63), an isotropic spherical model (range ¼ 1033.21 m, sill ¼ 0.48) and a nugget effect of 0.34, as shown in Fig. 2a. The lag distance was 250 m. An elliptic searching neighborhood with four sectors was used. The minimum and maximum numbers of neighboring points considered in each sector were 2 and 4 respectively. The goodness of fit was checked by cross-validation showing a satisfactory result. The standardized mean of prediction errors approached zero and the average standard error was close to the root-mean-squared prediction error. Zonal anisotropic features were clearly revealed in the semivariogram surface, and the maximum continuity was located at the E–W direction (Fig. 2b). The ratio of nugget to sill was 23.45% indicating strong spatial dependence in arsenic distribution according to Wang et al. (2001).

The long range of Gaussian model indicated that the regional spatial autocorrelation prevailed. Information of arsenic concentrations contained in the set of 100 realizations is summarized in Fig. 3. The maps of the mean and standard deviation reflected the uncertainty of predictions based on the 104 sample sites. While local spatial patterns were different among simulated maps, some special features were more certain as seen on most of the simulated maps. Arsenic showed anomalously high concentrations close to the smelter, decreasing with distance until the background value. The most distinct area of high concentrations was the V-shape zone connecting the two communities, which coincided with the prevailing wind direction. Another high-concentration area was located between the smelter complex and the Ross lake, probably because it was the nearest residential region to the stack. The uncertainty about predicted arsenic concentrations, as measured by the differences among realizations, was also the largest on these hotspots. Averaged concentrations varied from 5.18 to 157.93 mg/kg, and were above the national guideline (12 mg/kg) for 71.80% of the study area. The total cancer risk based on the 100 concentration maps is summarized in Fig. 4. Estimated risks were compared with a 1-in100 000 incremental risk which is the risk level of negligible impact for Canadian contaminated site exposures (Health Canada, 2004). The regional pattern generally coincided with the distribution trend of arsenic concentrations, because risk increased with arsenic concentrations in the dose–response models. At the DA scale, the V-shape hotspot area shrank dramatically, because the related DA had a relatively large coverage and concentrations varied

Fig. 3. The mean (a) and the standard deviation (b) of the distribution of arsenic concentrations from 100 realizations of GCS (unit: mg/kg).

Fig. 4. Summary statistics for the distribution of risk levels at the nodes of a 50-m grid (left) or at the scale of Dissemination Area (right): mean (a, b), deviation (c, d), and probability that the risk level exceeds 1/100 000 (e, f). Units for mean and deviation are 1/100 000.

2418

H. Zhang et al. / Environmental Pollution 157 (2009) 2413–2420

significantly within it. Dividing this DA into several sub-DAs would better represent the variation of risk levels if demographic data is available at sub-DA scale. For decision making, the uncertainty of cancer risk was supplemented by a measure of the likelihood that the criterion of Health Canada is exceeded. This likelihood was derived at both grid and DA scales by counting the proportion of 100 simulated risks that were higher than 105. The probability map indicated a large probability of exceeding the threshold at areas immediately adjacent or downwind to the smelter complex. Among the nine DAs with an average simulated risk above 105, the probability of exceedance ranged from 0.23 to 1. This information is critical for the risk management, especially for delineating hazardous area or identifying locations for additional sampling. These probabilities could be further interpreted by comparing to two threshold p1 and p2 defining tolerable risks for classifying a site as safe (probability of exceedance < p1) or hazardous (probability of exceedance > p2), while probabilities in the range [p1, p2] indicate sites whose classification is uncertain (Goovaerts et al., 2008). Assuming threshold values as 0.2 and 0.8, 49.72% of the study area would be classified as safe at the grid scale, while 18.68% as dangerous; at the regional scale, two DAs could be in dangerous situation and seven could be safe. The average risk level of each community is also presented for each cancer (Table 1). For male residents in both communities, arsenic-induced cancer risks were dominated by skin cancer, while female residents were mainly threatened by skin cancer and liver cancer. Estimated health effects were generally more serious in Creighton than in Flin Flon; this is possibly because that Creighton is much smaller than Flin Flon in terms of inhabited area and the high-concentration zone accounts for a relatively large portion of total area in Creighton. Special attentions of regulatory agency and community would be needed for areas where the estimated risks are over 10 times higher than the criteria of Health Canada. 4. Discussion The proposed integration of geostatistical simulation and dose– response models produces point estimates of cancer risks under consideration of spatial uncertainty. The results are also impacted by other uncertainties in the information and assumptions underlying the health issues of interest. These uncertainties cannot be incorporated into the present assessment because of limited availability of required data and models. They are discussed in this section from three aspects: (1) the present and future concentrations of arsenic in soil, (2) the characteristics of target population in the health risk assessment, and (3) the dose–response models for estimating cancer risks from arsenic in soil. 4.1. Arsenic concentrations in soil Geostatistical simulation generates concentration maps with a balance between the estimation of point concentration and the

description of regional variation. Because the simulation is based on the preliminary analysis of sample data, it is crucial to judge whether the soil samples are representative of the real situation. In this study, the answer largely relies on two issues. The first is the spatial configuration of sampling sites. An ideal scheme would be sampling uniformly and densely, which is always constrained by resources availability and research objectives. The majority of samples sites were located in the City of Flin Flon because the survey was conducted by Manitoba Conservation; only 13 of 104 sites (12.5%) lied in the extent of Creighton which accounted for 26.5% of the total study area. Sampling work was also focused on the densely inhabited blocks, while few attentions were given to the wide suburban area. Such a preferential sampling strategy reflected the primary concern on human health effects, but would not be the best plan for capturing the regional trend of arsenic concentrations. Although this impact could be partly corrected through using declustering techniques, the variation of arsenic concentrations in sparsely sampled areas would still not be captured effectively. The second issue is the uncertainty in the measurements at individual sites. Even within a short range around a site, the transportation and degradation processes of arsenic compounds could varies largely due to natural variations and/or human activities. This could be exaggerated by operational errors in fieldwork and laboratory analysis. As found during the survey (Manitoba Conservation, 2007), the difference among three replicate samples was relatively high at a number of sites. The uncertainties in the conditional dataset will propagate through the simulation process and decrease the accuracy of estimations. A much greater uncertainty lies in the distribution of arsenic concentrations over time as discussed by Yu et al. (2003). Since the exposure duration in cancer risk assessment is the lifetime, one important assumption of present study is that the local population will expose to the same concentrations as those at present. However, the future concentration could be largely different from the present level. For the study area, decreased concentrations could occur at the regional level when emission from the smelter complex is further reduced using improved technologies; at the local level, arsenic concentrations in soil could vary dramatically due to land use activities (e.g. developing bare land into residential lot) or remediation processes (e.g. restoring the natural vegetation cover). Considering this uncertainty in the risk assessment would need information regarding potential land cover changes, remediation measures and emission control plans, which are hardly available at present. 4.2. Exposure to arsenic in soil The extent of inhabited area is a crucial factor for estimating cancer risk at the regional scale. The accuracy of satellite images interpretation could be improved if rigorous ground verification was available. Furthermore, because DA is the smallest census unit available for the study area, we assume the distribution of

Table 1 Estimated cancer risks in Creighton and Flin Flona (unit: 1/100 000). Cancer

Flin Flon

Creighton

Male

Female

Skin Lung Bladder Liver

1.5 3.7 2.2 1.5

4.2 2.7 3.3 5.9

Total

1.5 (1.0)

a

(1.0) (2.5)  102 (1.2)  105 (0.9)  105

Standard deviations shown in parenthesis.

(2.8)  101 (1.4)  105 (1.6)  105 (3.9)  101

1.0 (0.7)

Both

Male

Female

1.9 3.7 5.5 5.9

2.5 6.3 7.9 4.8

6.9 8.5 1.0 9.7

(1.3) (2.5)  102 (2.8)  105 (3.9)  101

2.5 (1.7)

(1.8) (4.6)  102 (4.5)  105 (2.8)  105

2.5 (1.8)

(5.0)  101 (4.8)  105 (0.6)  104 (7.1)  101

1.7 (1.2)

Both 3.1 6.3 1.8 9.7

(2.3) (4.6)  102 (1.0)  104 (7.1)  101

4.2 (3.0)

H. Zhang et al. / Environmental Pollution 157 (2009) 2413–2420

population is uniform within each DA; thus the risk estimate for a DA is the arithmetic average of all point estimates within the inhabited area in this DA. However, we recognize that the actual population density is not uniform and high-density residential areas should be given more weights in regional estimation. This reality could be better reflected only if information at finer scale is available; for example, Gay and Korre (2006) used digital Ordnance Survey maps and EDINA Digimap Land-Line Plus data to identify residential buildings and to allocate community population to each building. We assume the entire population of Creighton and Flin Flon to be exposed to the arsenic in soil. More people will be threatened by the arsenic-induced cancer risks in the future if the population grows. In addition, since cancer risk increases exponentially with age in dose–response models, the proportion of senior citizens is decisive to the risk estimation of the entire population. However, uncertainties exist in both the quantity and the age structure of future population. The population of the study area has decreased since 2001, but both Saskatchewan and Manitoba have seen strong population growth since 2007 due to gains from immigration. It is thus unsure whether the exposed population will continue to go down, and to what extent this vagueness can affect the estimation of future demographic structure. Another uncertainty lies in the mobility of population. It is practical to allocate each resident to a DA for census purpose, but it is not rational to assume his/her exposure also occurs in this DA exclusively. Specifying the daily footprints for every resident would be overcomplicated, but establishing a regional pattern of population mobility should help reduce the uncertainties in exposure evaluation. This is especially valuable in a northern Canadian context where outdoor activities are largely constrained due to the inclement weather. Exposure equations developed by Health Canada are based on the Canadian average conditions, and parameter values need to be modified to match the site conditions in the study area. Intakes of arsenic from all pathways contribute to the health risk in an accumulative way. Besides the soil-related pathways considered here, other possible routes include inhalation of contaminant vapors, drinking of contaminated water and ingestion of contaminated food (Mandal and Suzuki, 2002; Dube´ et al., 2004; Lambert and Lane, 2004; Zhu et al., 2008). Furthermore, a number of factors unrelated to arsenic also contribute to the development of cancers (e.g. smoking). From the perspectives of communities and government agencies, the estimated risks from this study would be more significant if they are compared to risks resulted from other exposure pathways and carcinogens. 4.3. Dose–response models The selected dose–response models were originally developed for cancer risks from drinking arsenic-contaminated groundwater based on data from Taiwan. Uncertainties of dose–response models are distinct from and independent of uncertainties in the geographical distribution of arsenic or in the exposure of people to the arsenic (Yu et al., 2003). Two questions could be raised. Are these models applicable to the Flin Flon–Creighton population? Are they applicable to arsenic in soil? The first question concerns the difference between the previously studied Asia population and the Canadian population studied here. The risk level generated by Equation (4) is a prevalence ratio of skin cancer based on the assumption of zero background prevalence, and the incidence rate by Equation (5) is an additional rate to the non-zero background value (Yu et al., 2003). These assumptions are supported by the evidence from detailed epidemiological survey in Taiwan, but such a survey has not been reported in Flin Flon–Creighton. Another assumption is that the

2419

susceptibility of people in Flin Flon–Creighton should be similar to that of people in Taiwan. This could be questioned by differences in factors such as nutrition conditions, genetic features, and age structure, but related information is not available. Uncertainty also lies in the exposure duration. In the Taiwan case people had ingested arsenic from groundwater wells for at least 60 years, while people in Creighton and Flin Flon has exposed to the arsenic for eight decades. However, because harsh weather keeps the residents indoors during winter, the actual exposure duration in this area should be shorter that that in Taiwan given the same age. The second question addresses the transferability of existed dose–response models. This study is the first attempt of using parametric dose–response models for health risks induced by soilbound arsenic. We assume that the carcinogenic effect of arsenic in groundwater will be the same as that of soil-derived arsenic if the amounts of arsenic in both systems are also the same. The arsenic concentrations of soil samples (mg/kg) are thus converted into the equivalent arsenic concentrations in groundwater (mg/l) based on daily intakes of soil and groundwater. The daily consumptions of drinking water are adopted from the national data of Health Canada (2004), and the daily intakes of soil are determined through the exposure equations. The estimated arsenic concentrations in this study vary from 0 to 567 mg/g, equivalent to a range of 0–7.41 mg/l of arsenic in drinking water for a Canadian adult; it is within the range with which the dose–response models were developed and applied (Brown et al., 1989; Yu et al., 2003). We recognize that uncertainty exists here because of the bioavailability of arsenic, which is the fraction of arsenic freely available to be absorbed into a body. It was reported that the bioavailability of arsenic in soil is less that of the soluble arsenic (Dube´ et al., 2004); however, as discussed by Valberg et al. (1997), results indicating less bioavailability of soilbound arsenic were only derived from animal studies and tests of statistical significance were not provided. Furthermore, it is not clear that to what extent those reported differences in arsenic bioavailability could apply to the northern regions where extreme natural conditions prevail. Bioavailability of arsenic is affected by many factors that control the fate and behaviour of arsenic in environment (e.g. pH and the reduction/oxidation potential), while related soil and biological studies are not reported yet for cold regions like Flin Flon–Creighton. Therefore, before parametric dose–response functions specific to soil-bound arsenic are available, we would suggest the best use of what is available. 5. Conclusions This study integrates geostatistical simulation and dose– response models to assess the cancer risk induced by soil-derived arsenic in two northern Canadian communities. High-concentration zones of arsenic were found at areas immediately adjacent or downwind to the smelter. Risks of skin cancer and three internal cancers are not only provided at grid scale and census-unit scale, but also presented as the probability of exceeding national criterion. The results showed that the spatial pattern of risk levels generally coincided with the spatial distribution of arsenic, and indicated that non-negligible cancer risks could exist in the two communities due to the exposure to high concentrations of arsenic in residential soil. Uncertainties associated with the risk estimates were discussed from the aspects of arsenic concentrations, exposed population and dose–response model. Reducing these uncertainties would require additional soil sampling, epidemical records as well as complementary studies on land use, demographic variation, outdoor activities and bioavailability of arsenic. The proposed approach could be a viable option to the human health risk assessment not only for Flin Flon and Creighton, but also for other communities threatened by arsenic-contaminated soil.

2420

H. Zhang et al. / Environmental Pollution 157 (2009) 2413–2420

Acknowledgments We would like to thank Manitoba Conservation for providing the data of arsenic concentrations and site locations, and Dr. Antoni Magri of Cornell University for the discussion on geostatistical simulation. This work was funded by the Major State Basic Research Development Program of MOST (2005CB724200 and 2006CB403307), the Special Research Grant for University Doctoral Programs (20070027029), the Canadian Water Network under the Networks of Centers of Excellence (NCE), and the Natural Science and Engineering Research Council of Canada. The authors are grateful to the reviewers for their constructive criticisms that contribute to the improvement of this article.

References Bourennane, H., King, D., Couturier, A., Nicoullaud, B., Mary, B., Richard, G., 2007. Uncertainty assessment of soil water content spatial patterns using geostatistical simulations: an empirical comparison of a simulation accounting for single attribute and a simulation accounting for secondary information. Ecological Modelling 205, 323–335. Brown, K.G., Boyle, K.E., Chen, C.W., Gibb, H.J., 1989. A dose–response analysis of skin-cancer from inorganic arsenic in drinking-water. Risk Analysis 9, 519–528. Chen, A., Huang, G.H., Chakma, A., 1998. Integrated environmental risk assessment for petroleum–contaminated sites – a North American case study. Water Science and Technology 38, 131–138. Chen, T., Liu, X.M., Zhu, M.Z., Zhao, K.L., Wu, J.J., Xu, J.M., Huang, P.M., 2008. Identification of trace element sources and associated risk assessment in vegetable soils of the urban–rural transitional area of Hangzhou, China. Environmental Pollution 151, 67–78. Cinnirella, S., Buttafuoco, G., Pirrone, N., 2005. Stochastic analysis to assess the spatial distribution of groundwater nitrate concentrations in the Po catchment (Italy). Environmental Pollution 133, 569–580. Dietrich, C.R., Newsam, G.N., 1997. Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix. SIAM Journal on Scientific Computing 18, 1088–1107. Dube´, E.M., Boyce, C.P., Beck, B.D., Lewandowski, T., Schettler, S., 2004. Assessment of potential human health risks from arsenic in CCA-treated wood. Human and Ecological Risk Assessment: An International Journal 10, 1019–1067. Emery, X., 2004. Testing the correctness of the sequential algorithm for simulating Gaussian random fields. Stochastic Environmental Research and Risk Assessment 18, 401–413. Gay, J.R., Korre, A., 2006. A spatially-evaluated methodology for assessing risk to a population from contaminated land. Environmental Pollution 142, 227–234. Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation, first ed. Oxford University Press, New York. Goovaerts, P., Trinh, H.T., Demond, A., Franzblau, A., Garabrant, D., Gillespie, B., Lepkowski, J., Adriaens, P., 2008. Geostatistical modeling of the spatial distribution of soil dioxins in the vicinity of an incinerator. 1. Theory and application to Midland, Michigan. Environmental Science & Technology 42, 3648–3654. He, L., Huang, G.H., Zeng, G.M., Lu, H.W., 2008. An integrated simulation, inference, and optimization method for identifying groundwater remediation strategies at petroleum-contaminated aquifers in western Canada. Water Research 42, 2629–2639. Health Canada, 2004. Federal Contaminated Site Risk Assessment. In: Canada Part I: Guidance on Human Health Preliminary Quantitative Risk Assessment (PQRA). Health Canada. Henderson, P.J., McMartin, I., Hall, G.E., Percival, J.B., Walker, D.A., 1998. The chemical and physical characteristics of heavy metals in humus and till in the vicinity of the base metal smelter at Flin Flon, Manitoba, Canada. Environmental Geology 34, 39–58.

Kapaj, S., Peterson, H., Liber, K., Bhattacharya, P., 2006. Human health effects from chronic arsenic poisoning – A review. Journal of Environmental Science and Health Part A – Toxic/Hazardous Substances & Environmental Engineering 41, 2399–2428. Kyriakidis, P.C., Yoo, E.H., 2005. Geostatistical prediction and simulation of point values from areal data. Geographical Analysis 37, 124–151. Lambert, T.W., Lane, S., 2004. Lead, arsenic, and polycyclic aromatic hydrocarbons in soil and house dust in the communities surrounding the Sydney, Nova Scotia, tar ponds. Environmental Health Perspectives 112, 35–41. Lee, J.J., Jang, C.S., Liang, C.P., Liu, C.W., 2008. Assessing carcinogenic risks associated with ingesting arsenic in farmed smeltfiish (Ayu, Plecoglossus altirelis) in aseniasis-endemic area of Taiwan. Science of the Total Environment 403, 68–79. Lee, J.J., Jang, C.S., Wang, S.W., Liu, C.W., 2007a. Evaluation of potential health risk of arsenic-affected groundwater using indicator kriging and dose response model. Science of the Total Environment 384, 151–162. Lee, S.Y., Carle, S.F., Fogg, G.E., 2007b. Geologic heterogeneity and a comparison of two geostatistical models: sequential Gaussian and transition probability-based geostatistical simulation. Advances in Water Resources 30, 1914–1932. Lin, Y.P., Chang, T.K., Teng, T.P., 2001. Characterization of soil lead by comparing sequential Gaussian simulation, simulated annealing simulation and kriging methods. Environmental Geology 41, 189–199. Mandal, B.K., Suzuki, K.T., 2002. Arsenic round the world: a review. Talanta 58, 201–235. Manitoba Conservation, 2007. Concentrations of Metals and Other Elements in Surface Soils of Flin Flon, MB and Creighton, SK, 2006. Manitoba Conservation, Winnipeg, MB. McMartin, I., Henderson, P.J., Nielsen, E., 1999. Impact of a base metal smelter on the geochemistry of soils of the Flin Flon region, Manitoba and Saskatchewan. Canadian Journal of Earth Sciences 36, 141–160. Newhook, R., Hirtle, H., Byrne, K., Meek, M.E., 2003. Releases from copper smelters and refineries and zinc plants in Canada: human health exposure and risk characterization. Science of the Total Environment 301, 23–41. Qin, X.S., Huang, G.H., Li, Y.P., 2008. Risk management of BTEX contamination in ground water – an integrated fuzzy approach. Ground Water 46, 755–767. Schelle, E., Rawlins, B.G., Lark, R.M., Webster, R., Staton, I., McLeod, C.W., 2008. Mapping aerial metal deposition in metropolitan areas from tree bark: a case study in Sheffield, England. Environmental Pollution 155, 164–173. Shi, G.T., Chen, Z.L., Xu, S.Y., Zhang, J., Wang, L., Bi, C.J., Teng, J.Y., 2008. Potentially toxic metal contamination of urban soils and roadside dust in Shanghai, China. Environmental Pollution 156, 251–260. Tang, X.Y., Zhu, Y.G., Shan, X.Q., McLaren, R., Duan, J., 2007. The ageing effect on the bioaccessibility and fractionation of arsenic in soils from China. Chemosphere 66, 1183–1190. USEPA, 2004. IRIS Substance File for Arsenic, Inorganic. USEPA. USNRC, 1999. Arsenic in Drinking Water. In: Council, U.N.R. National Academy Press, Washington, D.C. USNRC, 2001. Arsenic in Drinking Water: 2001 Update. In: Council, U.N.R. National Academy Press, Washington, D.C. Valberg, P.A., Beck, B.D., Bowers, T.S., Keating, J.L., Bergstrom, P.D., Boardman, P.D., 1997. Issues in setting health-based cleanup levels for arsenic in soil. Regulatory Toxicology and Pharmacology 26, 219–229. Verstraete, S., Van Meirvenne, M., 2008. A multi-stage sampling strategy for the delineation of soil pollution in a contaminated brownfield. Environmental Pollution 154, 184–191. Wang, J., Fu, B.J., Qiu, Y., Chen, L.D., Wang, Z., 2001. Geostatistical analysis of soil moisture variability on Da Nangou catchment of the loess plateau, China. Environmental Geology 41, 113–120. Wang, S., Mulligan, C.N., 2006. Occurrence of arsenic contamination in Canada: sources, behavior and distribution. Science of the Total Environment 366, 701–721. Webster, R., Oliver, M.A., 2007. Geostatistics for Environmental Scientists, second ed. Wiley, Chichester. Yang, J.S., Wang, Y.Q., August, P.V., 2004. Estimation of land surface temperature using spatial interpolation and satellite–derived surface emissivity. Journal of Environmental Informatics 4, 40–47. Yu, W.H., Harvey, C.M., Harvey, C.F., 2003. Arsenic in groundwater in Bangladesh: a geostatistical and epidemiological framework for evaluating health effects and potential remedies. Water Resources Research 39. Zhu, Y.G., Williams, P.N., Meharg, A.A., 2008. Exposure to inorganic arsenic from rice: a global health issue? Environmental Pollution 154, 169–171.