# Urban Lead: Modeling Its Distribution and Effects on Children

Authors: Zhixiong Chen; Yi Ding; Andrew Getz; Bernard Lipat
Journal: Journal of Applied Mathematics
 (2017)

---

## Abstract

We model the transportation of lead from the atmosphere and from the surface of the soil simultaneously at the macroscale and mesoscale to study its health effects on children in Jersey City, NJ. We conceptualize Jersey City as an open system where lead is continuously emitted from a local smelting plant and a local power plant, deposited onto the surface soil of playgrounds, and ingested by children. The model is constructed using the diffusion-advection partial differential equation in three spatial dimensions and one temporal dimension with an initial condition and boundary conditions. The model is solved using the Crank-Nicolson numerical method at the macroscale to determine the deposition of lead from the smelting plant and the local power plant and at the mesoscale to refine the amount of lead deposition for the areas considered. We then determine the health consequences for the average child using the bioaccessibility of lead from soil to children, the bioavailability of ingested lead to the circulatory system, and the biological half-life of lead isotopes in the blood. The health effects on children from lead are directly proportional to the blood lead concentration.


---

## Body

## 1. Conceptualization of the System

On June 22, 2012, officials of the City of Jersey City announced the closure of downtown-area Mary Benson Park [1]. The reason was lead. The problem of lead in the urban environment so epitomized by Jersey City has been a long one and a challenging one. Although the second largest city in the State of New Jersey, Jersey City, is a preeminent commercial, distribution, and transportation hub for the Tristate Area [2–4], Jersey City has not been immune to the decline in manufacturing so characteristic of postindustrial, globalized America. Once home to the production of such popular brands as American Can, Emerson Radio, Lorillard, Colgate, Dixon Ticonderoga, and Chloro, Jersey City now claims only four industrial centers: Greenville Yards, Claremont Industrial Park, Montgomery Industrial Park, and Liberty Industrial Park [2–4]. It is unfortunate, however, that this reduction in manufacturing coupled with the ban on leaded gasoline has not resulted in the elimination of the threat of lead poisoning: the accumulation of lead in the soil of industrial areas has created a problem that persists.In the present study, we construct a mathematical model to determine the spatial distribution of lead in the City of Jersey City so as to determine the resultant effects on the health of children. We begin with the foundational assumptions about the system. First we assume that the City of Jersey City is an open cube, wherein lead is free to enter and exit. We assume that the two major sources of lead in the vicinity are the only sources of lead and that the effects of the other sources are negligible. To make calculations simpler, we assume that all isotopes of lead have the same atomic radius. We assume, by the insolubility of lead at reasonably balanced acidity levels, that lead will not be transported by rain or surface water, nor will it diffuse through the soil. Further, we assume that lead, because of its particulate form, will diffuse through the air, assumed to be homogeneous, and be transported by wind advection. Our initial and boundary conditions are approximately chosen to simulate data for Jersey City, keeping the simplistic 3D-City model. We focus specifically on four public parks in Jersey City.


## 2. Diffusion-Advection

The transportation of lead through the air is governed by the diffusion-advection equation:(1)∂C∂t=D∂2C∂x2+∂2C∂y2+∂2C∂z2-vx∂C∂x-vy∂C∂y-vz+vT∂C∂z,where C=Cx,y,z,t is the concentration of lead at a point x,y,z in space and at time t, D is the diffusion coefficient for lead particles in the air, vx is the velocity of the wind in the x, or eastwest, direction, vy is the velocity of the wind in the y, or northsouth, direction, vz is the velocity of the wind in the z, or vertical, direction, and vT is the terminal velocity of a falling particle of lead.The diffusion coefficient for lead in air is calculated using the mean velocity and the mean free path [5]; that is,(2)D=BkT3π2Pdair+dPb28BkTπmN,where Bk=1.380650×10-23 J/K is Boltzmann’s constant, T=298.15 K is the temperature, mN=3.44063610352×10-25 kg is the mass of a particle of lead, P=100000 kg/m3 is the pressure, dair=3.66×10-10 m is the average diameter of a particle in air, and dPb=1.4193×10-7 m is the diameter of a particle of lead [5–7].The advective coefficients, that is, velocity of the wind in thex- and y-directions, respectively, vx and vy are taken as functions of season, as in Table 1 [8], where α=π/4 and β=π/8.Table 1
Autumn
Winter
Spring
Summer
v
x
3
⋅
cos
⁡
α
4
⋅
cos
⁡
α
3
⋅
cos
⁡
β
2
⋅
cos
⁡
α
v
y
3
⋅
sin
⁡
α
4
⋅
sin
⁡
α
3
⋅
sin
⁡
β
2
⋅
sin
⁡
αThe directions of the velocities reflect the choice of origin to be the bottommost southeast corner of the cube. During the autumn, winter, and spring, the direction of the wind is southeasterly, that is, in the direction of the origin. However, during the summer, the direction of the wind is northeasterly. To maintain the strict diagonal dominance of the coefficient matrixA, which will be detailed in Section 3, the cube is reflected so that the origin becomes the bottommost northeast corner of the cube.The advective coefficient, or velocity of the wind in thez-direction, vz=dz/dt is calculated by solving the hydrostatic ordinary differential equation, which relates the vertical velocity of the wind with the ω-vertical velocity, which is the rate of change of pressure over time. The hydrostatic equation is as follows:(3)dPdz=-d·g⟹dPdtdtdz=-d·g⟹dzdt=1-d·gdPdt,where dP/dt is the ω-vertical velocity in Pa/s, d is the air density in kg/m3, and g is gravitational acceleration [9–16]. Air density d is calculated according to the ideal gas law for mixtures of ideal gases(4)d=pa287.05T+pv461.495Tpa being the partial pressure of dry air in Pa, pv being the partial pressure of water vapor also in Pa, and T being the air temperature in K and 287.05 and 461.495 being the specific gas constants for dry air and water vapor, respectively [9–16]. The partial pressure of water vapor is determined by(5)pv=RH100pvsat,where RH is the relative humidity and pvsat is the partial pressure of water vapor at saturation, which since pvsat is [17](6)pvsat=RH·6.1078×107.5T-2048.625/T-35.85solves for pv. Again, using standard atmospheric temperature and pressure, and since dP/dt=0.02 Pa/s, we have dz/dt=vz=-0.00158 m/s [16, 17].The terminal velocityvT is calculated using the definition(7)vT=2mNgdCDA,where mN=3.44063610352×10-25 kg, CD=0.47 is the dimensionless drag coefficient of a spherical particle of lead, and A=3.16423×10-14 m2/s is the cross-sectional area of a lead particle [5–7]. Thus vT=4.9267×10-6 m/s.


### 2.1. Initial Conditions

Data available on the mixing ratios of lead, measured in parts per million (ppm) in 2007 at each of the parks under consideration, was used in construction the initial conditions of the system [18]. We have taken the data of lead found in the surface of soil, which we will use to show the effects on human life. We do not consider the lead below the surface. They are summarized in Table 2.Table 2
Park
Concentration (ppm)
Hamilton Park
500
Mary Benson Park
383.34
E. F. Jones Park
840
Lincoln Park
512.5To avoid the unreasonable and distorting assumption that the concentration of lead at points in the system with unavailable data is nonexistent, the following interpolation method for the initial conditions is used:(8)C¯=∑i=1n2-pC+∑i=1pfxi,yin2⟹C=n2C¯-∑i=1pfxi,yin2-pto agree with the known concentration data and to preserve the average concentration over the surface soil.Once the data values were uncovered, two-dimensional cubic spline interpolation was used to refine the surface soil data.


### 2.2. Boundary Conditions

The concentrations at the two sources of lead from 2007 to the present were available from the Environmental Protection Agency [19] and are summarized in Table 3.Table 3
Year
Source concentration
2007
2008
2009
2010
2011
2012
Kearny point (ppm)
2325
2194
1515
680
2259.75
2259.75
Hudson generating (ppm)
210
581
407
55
313.25
313.25Consistent with the assumption of openness, Neumann boundary conditions were imposed at the six faces of the cube:(9)∂C∂xx=0=0∂C∂yy=0=0∂C∂zz=0=0,∂C∂xx=L=0∂C∂yy=L=0∂C∂zz=L=0,where L is the distance to the opposite face of the cube.The boundary conditions are incorporated into the numerical scheme at each of the boundaries. The emissions from a local smelting factory and a local power plant are taken into account.


## 2.1. Initial Conditions

Data available on the mixing ratios of lead, measured in parts per million (ppm) in 2007 at each of the parks under consideration, was used in construction the initial conditions of the system [18]. We have taken the data of lead found in the surface of soil, which we will use to show the effects on human life. We do not consider the lead below the surface. They are summarized in Table 2.Table 2
Park
Concentration (ppm)
Hamilton Park
500
Mary Benson Park
383.34
E. F. Jones Park
840
Lincoln Park
512.5To avoid the unreasonable and distorting assumption that the concentration of lead at points in the system with unavailable data is nonexistent, the following interpolation method for the initial conditions is used:(8)C¯=∑i=1n2-pC+∑i=1pfxi,yin2⟹C=n2C¯-∑i=1pfxi,yin2-pto agree with the known concentration data and to preserve the average concentration over the surface soil.Once the data values were uncovered, two-dimensional cubic spline interpolation was used to refine the surface soil data.


## 2.2. Boundary Conditions

The concentrations at the two sources of lead from 2007 to the present were available from the Environmental Protection Agency [19] and are summarized in Table 3.Table 3
Year
Source concentration
2007
2008
2009
2010
2011
2012
Kearny point (ppm)
2325
2194
1515
680
2259.75
2259.75
Hudson generating (ppm)
210
581
407
55
313.25
313.25Consistent with the assumption of openness, Neumann boundary conditions were imposed at the six faces of the cube:(9)∂C∂xx=0=0∂C∂yy=0=0∂C∂zz=0=0,∂C∂xx=L=0∂C∂yy=L=0∂C∂zz=L=0,where L is the distance to the opposite face of the cube.The boundary conditions are incorporated into the numerical scheme at each of the boundaries. The emissions from a local smelting factory and a local power plant are taken into account.


## 3. Numerical Method

To solve the diffusion-advection PDE, the Crank-Nicolson numerical scheme is employed, and we used MATLAB to assist us. It is the average of the forward finite difference method and of the backward finite difference method with the index shifted to agree with the forward difference method. Essentially, the scheme uses the relation, as set out by the two discretizations of the diffusion-advection PDE, between the values of the concentration at the six surrounding points in 3-dimesional space at one time step and the values of the concentration at the same six surrounding points in the next time step, to approximate the value of the concentration at the central point over time. The resulting finite difference scheme is desirable over either of its constituents because of its unconditional stability [20].Then×n×n cube is reshaped into an n3×1 column vector by iterating first over x, then over y, and finally over z. The following matrix equation results:(10)A·Ci,j,kl+1=B·Ci,j,kl⟹Ci,j,kl+1=A-1B·Ci,j,kl,where l is the time step, i is the space step in the x-direction, j is the space step in the y-direction, k is the space step in the z-direction, C is the n3×1 column vector of concentration values, and A and B are n3×n3 square matrix of coefficients [21]. The invertibility of A, a strictly diagonally dominant matrix, is guaranteed by the Levy-Desplanques Theorem [22] and the reflection described in Section 2.The boundary conditions are incorporated into the matrix using the centered difference method at each boundary, that is, at each of the six faces of the cube. For example,(11)∂C∂xx=0≈C0,j,kl-C2,j,kl2h=0⟹C0,j,kl=C2,j,kl,∂C∂xx=L≈Cn-1,j,kl-Cn,j,kl2h=0⟹Cn,j,kl=Cn-1,j,kl,where the space step is h=1300 m at the macroscale and h=100 m at the mesoscale. Let ai1n3 be the sequence of row of the column vector C. Given the set-up of the cube, the x=0 boundary occurs where ai mod n=1, x=L occurs, where ai mod n=0, y=0 occurs, where 1≤ai mod n2≤n, y=L occurs, where n2-n≤ai mod n2≤0, z=0 occurs, where 1≤ai mod n3≤n2, and z=L occurs where n3-n2≤ai mod n3≤0.The matrix equation is programmed into MATLAB and iterated with a time step consistent with the stability requirement [8](12)D·lh2≤12from autumn 2007 to autumn of 2012, that is, twenty-four seasons.The results are depicted in Figures1 and 2. Although the model takes into consideration the seasonal variations in the directions and the velocities of the wind, there is, as expected, a southeasterly tendency for lead transportation across the surface soil. In fact, even though this lead is blown across the surface soil, the increase in the concentration of lead, particularly at the parks, indicates that the rate of lead deposition from the two sources of lead is greater than that of lead transportation out of the parks. As in Figure 1, if we calculate in macroscale, we can only see the lead deposition of Lincoln Park clearly; after we used mesoscale as shown in Figure 2, it is clear to see the lead deposition of E. F. Jones Park, Mary Benson Park, and Hamilton Park.Figure 1Figure 2Further, a sensitivity analysis was conducted on the macroscale numerical solution. A ten percent (±10%) change in each parameter (e.g., source concentrations, initial concentrations, vertical velocity, and angle of advection) yielded only slight percent differences in the final concentrations of lead. It was found that the model was relatively insensitive to most parameters. However, it should be noted that the model was moderately sensitive to a change in the angle of advection. This is depicted in Figure3. The absolute differences are so small, that is, on the order of 10-74, that the significance of a moderately large percent change (at the very maximum about -70% and only in the upper atmosphere 13 km above the surface) should not be conflated.Figure 3Of note, the model reflects the fact that in 2010 Hamilton Park underwent renovations that included removal and replacement of the surface soil by setting the concentration at that point to 100 ppm, which the authors thought reasonable in the absence of data to reflect a significant decrease in the concentration of lead but so as not to distort the lead deficiency and therefore the rate of lead deposition over time. The continued increase of the concentration of lead at Hamilton Park after 2010 at a rate commensurate with the rate of increase before 2010 implies, through something of a natural experiment, that remediation—the removal and replacement of surface soil—is not an effective long-term strategy because it will only be a matter of years before the concentration of lead at, say, Hamilton Park, returns to elevated and, perhaps, prohibited levels.Figure4 demonstrates the comparison of the data generated by the numerical solution of the model with data published by the City of Jersey City [1], and these are seasonal averages. Here we use our numerical solutions to compare with sample 2 (410 mg/Kg) and sample 4 (563 mg/Kg) from [1]. “Sample 2 and sample 4 from [1] were taken by an outside geotechnical and environmental materials testing consultant and were taken in four distinct locations within the playground area” [1]. The numerical solution of the model herein developed agrees with the data available from the aforementioned closure of Mary Benson Park in June 2012.Figure 4


## 4. Health Effects

Lead’s neurotoxic effects have been widely studied [14, 23–32]. They result mainly from lead’s similarity to and consequent competition with calcium, which is vital for the effective delivery of action potentials along axons. Replacement of calcium ions in the brain with lead results in retardation of this process and can have severe neurological effects.The following integral is used to calculate the cumulative exposure, that is, blood lead level (BLL), to lead at a pointx1,y1,z1, for a time period t1,t2:(13)BLL=Bacc·Bav·∫t1t2Cx1,y1,z1,t12t/30dt,where Bacc=0.154897 is the bioaccessibility factor and Bav=0.108 is the bioavailability factor [30, 32]. Bioaccessibility is the measure of the percentage of lead deposited onto the surface soil accessible for absorption into a biological system, and bioavailability is the percentage of lead absorbed into a biological system available for uptake into the circulatory system [33].The health effects of lead considered in this paper are those directly proportional to the blood lead level (BLL), with the proportionality summarized in Table4 [8, 15–17, 19, 23, 27, 33–36].Table 4
Health effects
Proportionality
Fine motor composite score
−4.6
IQ
−0.25
Mortality
0.0035
Average lifetime earning
0.005975Table5 depicts the cumulative exposure and resultant health effects under the assumption of an hour of exposure per day after the course of the simulation.Table 5
Summary of health effects
Park
Blood lead level (μg/dL)
Fine motor decrease
IQ decrease
Excess relative risk (%)
Earnings decrease (%)
Downtown
6.0428
0.1314
1.5107
2.115
0.9026
Hamilton
1.6447
0.0358
0.4112
0.5756
0.2457
Lincoln
7.0903
0.1541
1.7726
2.4816
1.0591
Mary Benson
7.136
0.1551
1.784
2.4976
1.0659
E. F. Jones
12.6348
0.2747
3.1587
4.4222
1.8873


## 5. Conclusions

This model is tuned to the conditions presented in the City of Jersey City, but it is not unreasonable to generalize these results to similar postindustrial urban centers. Still, future research may wish to refine the simplifying assumption regarding the bioaccessibility of lead by modeling saltation of granular particles near the surface of the soil to more precisely quantify the probability of uptake. A more accurate bioavailability factor may be achieved by modeling the competition between lead and calcium, for example, for neural receptors. The initial conditions are derived from samples taken on location and so, while necessarily accurate, are limited as any samples are limited by questions of representativeness and quality of equipment. The interpolation thereof, however, while preserving the invariant of average concentration, assumes uniform distribution. Further research may choose to assume that concentrations of lead are higher near highways and factories than they are away from such. Still, the diffusive and advective processes, while possibly weak in the short-term, may have resulted in a redistribution of lead in, say, the four decades since the ban of leaded gasoline.Nevertheless, the qualitative predictions of the solution to the model are no less worrisome. Of local interest, our model predicts especially elevated levels of lead concentration at many highly trafficked and popular parks, such as E. F. Jones Park. Public health officials may wish to take action to protect all who benefit from public recreation areas, especially children, by collecting and testing soil samples to validate experimentally the results predicted and to take the actions necessary to safeguard the public from adverse health effects. These actions may include public notices about the potential threats of lead ingestion and inhalation, covering exposed topsoil, prohibition of access to the parks, and soil remediation. Although none of the aforementioned are likely to be effective long-term solutions, they may be effective in the short run. Cost-benefit analyses may be done to determine which options are best.However, the results of the model’s predictions for Hamilton Park, while a single data point, reveal that even the most drastic of the short-term solutions, that is, soil remediation, is ineffective in the long run because more lead is being deposited by the emissions than is transported away by advection. Long-term solutions may include inspections of the local refineries, smelting factors, power plants, and so forth, to ensure compliance with EPA emissions regulations, reduction in lead emissions, and elimination of lead emissions altogether. There may be a steady-state solution with nonzero lead emissions, but, considering the current concentrations of lead in local parks, the health consequences may outweigh any potential economic benefit.How prescient it was of Dr. Clair C. Patterson, who wrote: “Sometime in the near future it probably will be shown that the older urban areas of the United States have been rendered more or less uninhabitable by the millions of tons of poisonous industrial lead residues that have accumulated in cities during the past century.” [37]



---

*Source: 6107430-2017-06-11.xml*