Terrestrial ecosystem health under long‐term metal inputs: modeling and risk assessment

Abstract

Metal contamination of soils may pose long‐term risks to ecosystem health if not properly managed. Future projection of contamination trends, coupled with ecological assessment, is needed to assess such risks. This can be achieved by coupling dynamic models of soil metal accumulation and loss with risk assessment on the basis of projected metal levels. In this study, we modeled the long‐term dynamics of Cu, Zn, and Cd in agricultural topsoils of a northern Chinese catchment (Guanting reservoir) and related projected metal levels to 2060 to ecological risk. Past metal dynamics were simulated using historical metal inputs from atmospheric deposition, irrigation, fertilizers, and animal manures. Modeling future dynamics was done using scenarios of projected metal input rates. Ecological risk assessment was done using the Potentially Affected Fraction () approach to estimate the combined toxic pressure due to the three metals. Modeled labile soil metals agreed well with measurements from monitoring in 2009 following adjustment of the porewater dissolved organic concentration. Metals were predicted to be largely retained in the topsoil. Projections were sensitive to changes in imposed soil pH, organic matter, and porewater dissolved organic carbon. Modeling suggests that decreases in input rates to between 5% and 7.5% of 2009 levels are required to prevent further accumulation. Computed s suggest zinc makes the greatest contribution to ecological risk. Under the most conservative estimate of , the threshold of potential ecological risk was reached before 2060 in two of the three future input scenarios.

Key words:

• China
• dynamic modeling
• ecological risk
• metals
• Potentially Affected Fraction
• soils

Introduction

The accumulation of hazardous metals in soils poses potential risks to ecological and human health. Soils receive metal inputs due to wet and dry deposition from the atmosphere, and agricultural soils may be exposed to metal inputs due to additions of sewage sludge, fertilizers, and animal litters and manures containing metals (Nicholson et al. 2003), irrigation with contaminated waters (Singh et al. 2010), or by use of metal‐based fungicidal agents such as Bordeaux mixture (Komarek et al. 2010). Metal accumulation in agricultural soils may have adverse ecological effects (Paoletti et al. 1998, Mackie et al. 2013) or enter the food chain via accumulation in crops, posing potential risks to animal and human health (Bi et al. 2006). Projection of future trends in soil metal accumulation is particularly important in areas where intensification of agricultural practices and industrial activity is placing steadily increasing pressure on the quality of agricultural soil (e.g., Luo et al. 2009).

Since metals usually associate strongly with soil solids, repeated inputs may lead to topsoil concentrations potentially capable of adversely impacting on soil health, even if annual inputs increase concentrations by a small amount relative to the concentration already present. Present day management of metal inputs for long‐term soil sustainability may benefit from projection of future metal concentrations and consequent risks under different scenarios of input rates. Models of soil metal dynamics have been applied in recent years to project future trends in the accumulation of metals such as copper, zinc, and cadmium in agricultural soils (e.g., de Vries et al. 2004, Groenenberg et al. 2006, Chen et al. 2013).

Here, we couple a model of soil metal dynamics (IDMM; Lofts et al. 2013) with the Potentially Affected Fraction (PAF) approach to project future metal accumulation and ecological risk to soil health due to inputs of Cu, Zn, and Cd. We firstly test the capability of the model to reproduce present day soil metal concentrations based on estimated historic natural and anthropogenic inputs. Our study catchment, the Guanting reservoir (GTR), is located in northern and has been previously surveyed for these metals (Luo et al. 2007, Xu et al. 2013). We highlight the key, sensitive model parameters and discuss how their uncertainties may be managed and reduced, and present projections of possible trends in metal concentrations to 2060. The PAF approach used (Lofts et al. 2005) takes into account the effect of bioavailability on the ecological risks of the metals, and allows estimation of their combined impacts.

Materials and Methods

Catchment description

The GTR is located northwest of Beijing city (E115.43°, N40.19° to E115.97°, N40.50°) (Fig. 1). The catchment area covers about 920 km², including 98 km² of water. The average elevation of the sampling sites (see Soil sampling and analysis) was 502 m ASL. The area has a cool temperature continental monsoon climate, with a mean annual temperature between 3°C and 9°C. Average annual precipitation is between 370 and 480 mm. The primary rock types are intermediate and acidic igneous; and the soil types are fluvo‐aquic (FAS), calcareous‐cinnamon (CCS), fluvo‐cinnamon (FCS) and meadow‐wind sand (MWSS). The land use types around the GTR are woodland, orchard, farmland, and fallow, which together account for about 90% of the total area. Much of the catchment is dedicated to agricultural use, focusing mainly on maize production.

short-legendFigure 1.

Soil sampling and analysis

In May 2009, ten farmland soil samples were collected from the GTR watershed. Each sample was a composite of five sub‐samples collected from the center and at the four corners of an area of approximately 100 × 100 m². Sub‐samples were collected at a depth of 0–20 cm using a stainless steel shovel. Soil pH was determined in a 1:2.5 w/w soil water slurry using a Delta 320‐s pH meter (Mettler‐Toledo, Greifensee, Shanghai, China). Soil organic matter (SOM) was determined by titration with iron(II) sulfate after digestion with a potassium dichromate‐sulfuric acid solution. Soils were digested with a mixture of concentrated HCl–HNO₃–HF–HClO₄. Total concentrations of Cu, Cd and Zn in the digestates were determined by ICP‐MS (Agilent 7500a). A soil standard reference material (GSS‐1) was analyzed as part of the quality assurance and quality control procedures. The detection limits of Cu, Cd, and Zn were 0.02, 0.005, and 0.002 μg/g, respectively. The recovery rates of Cu, Cd, and Zn were in the ranges 95–105%, 98–103%, and 96–108%, respectively. Labile soil metal concentrations were estimated by extraction with 0.05 M EDTA (Ure et al. 1993) and measurement by Atomic Absorption Spectrophotometry (Shimadzu AA‐6300).

Model Description

Basic principles of the model

The IDMM is designed for long‐term prediction of metal dynamics within a topsoil based on a parsimonious input data set. The soil is simulated as a single well‐mixed layer. The model runs on an annual time step.

Soil metal is partitioned into labile and nonlabile pools (see Metal chemistry). The labile (or “geochemically active”) pool comprises metal dissolved in the soil porewater and metal adsorbed to the surfaces of the soil solids, whereas the nonlabile pool comprises metal in less chemically reactive forms, for example occluded into mineral lattices. Metal may be lost from the layer via vertical water movement to the soil below or as surface runoff, either dissolved in the porewater or in adsorbed or aged form in eroded soil particles (see Metal transport in the soil). Metal may also be lost from the soil due to uptake into crop biomass and subsequent removal due to harvesting (see Metal uptake by crops).

Table 1 provides the parameters used for metal partitioning, aging, and cropping.

Table 1. Parameters used in the IDMM simulations

Metal	Solid‐solution partitioning (Eq. 1)				Metal aging (Fig. 2)				Cropping^†
	α0	α1	α2	n	log kf,a	log kb,a	log kf,m	log kb,m	{M}crop (μg/g dw)
Cu	−6.37	0.64	0.87	0.57	−2.5 + 10^−3.3∙e^pHpw	−2.1 + 10^−3.5∙e^pHpw	−5.0	^‡	2.42
Zn	−4.67	0.46	0.84	0.84	−4.2 + 0.26∙pHpw	−3.2	−5.0	^‡	21.7
Cd	−5.71	0.41	0.91	0.70	−2.9 + 0.18∙pHpw	−1.6	−5.0	^‡	0.013

†f_layer = 0.275.

‡Calculated by fixing the total metal concentration at pristine steady state.

Metal chemistry

Equilibrium solid‐solution partitioning of metals is handled by a combination of two models. A Freundlich‐type expression describes the relationship between the metal adsorbed to the soil solids and the free metal ion in soil solution: 1 where {M}_ads is the adsorbed labile concentration (mol·g soil⁻¹), [M]_free is the free ion concentration in porewater (mol/dm³), pH_pw is the porewater pH, and SOM is the soil organic matter concentration (% w/w). The parameters α₀, α₁, α₂, and n may be derived by fitting to metal partitioning data. We used the parameters derived by Groenenberg et al. (2010) from 216 soils of the United Kingdom and the Netherlands, covering a porewater pH range of 3.7–8.3 and SOM range of 0.5–97.8%. The total dissolved metal in the soil solution is computed from the free ion concentration using WHAM/Model VI (Tipping 1998). The porewater pH, pH_pw, and concentrations of Na, Mg, Ca, Cl, NO₃, SO₄, and dissolved organic matter (DOM) in the porewater are specified. Concentrations of carbonate species are computed assuming equilibrium with a fixed partial pressure of CO₂ in the soil atmosphere. Activities of Al(III) and Fe(III) are computed according to Lofts et al. (2013). Ion concentrations are adjusted during model setup to balance the porewater charge and to account for the possible control of Ca concentrations by calcite (CaCO₃) (see Temporal simulation of metal dynamics).

Transfers of metal between the labile and nonlabile pools (aging and weathering) are simulated by first order kinetics (Fig. 2). The nonlabile metal is partitioned into two pools termed “weakly aged” and “strongly aged”. Metal may reversibly transfer between the labile and weakly aged pools, from the weakly aged to the strongly aged pool, and from the strongly aged pool to the labile pool. A full description of the model aging and its parameterization is provided in Appendix S1.

short-legendFigure 2.

Metal transport in the soil

The fluxes of metal into, through and out of the soil comprise labile, weakly aged, and strongly aged components. Labile metal may be transported in dissolved form or adsorbed to eroded soil, while weakly, and strongly aged metal is transported only with eroded soil. The expressions for the annual loss fluxes of each form in drainage (in mol/m²) are 2 3 4 5 where V is the annual volume of drainage (m³/m²), F_ES,drainage is the annual flux of eroded soil in drainage (kg/m³), and f_HE is the proportion of highly erodible soil in the layer. Corresponding expressions are used to calculate annual loss fluxes of metal in leaching. The parameter f_HE allows for the eroded portion of the soil to be enriched in metal relative to the bulk soil, as has been shown to occur under field conditions (Quinton and Catt 2007). In this work we set f_HE to 0.3.

Metal uptake by crops

The IDMM simulates crop uptake by assuming a constant metal concentration, {M}_crop (mol/g fresh weight), in harvested crop material. The flux of metal out of the soil due to crop uptake and harvesting is then 6 where F_crop is the annual loss flux due to cropping in mol/m², Y_crop is the annual crop yield in tons fresh weight/ha, and f_layer is a term giving the proportion of the total crop metal sourced from the soil layer.

Temporal simulation of metal dynamics

The model is run from a starting year in which metal within the soil is assumed to be in steady state, representing “pristine” conditions prior to significant anthropogenic metal inputs. At steady state, input and output fluxes of metal are assumed to balance, that is, 7

The input flux of metal comprises the natural deposition flux F_deposition, which is assumed to be entirely labile, and a replenishment flux F_{replenishment} equal to the adsorbed and aged metal losses due to erosion. Erosion is assumed not to cause a net loss of soil mass; instead it is assumed that erosional losses of soil are balanced by an addition of the same mass of soil, having the same properties and containing solid‐phase metal at the same concentrations at that soil material lost. Accounting for soil due to erosion would cause the depth of the soil layer to decrease, so this soil addition effectively shifts the lower boundary of the soil layer deeper into the profile to maintain a constant soil depth. This balancing of erosional soil losses means that the flux balance at steady state can be rewritten in terms of dissolved metal fluxes, and a steady‐state concentration of dissolved metal in the porewater can be computed: 8 9

The steady‐state adsorbed labile concentration is then backcalculated, and the steady state weakly and strongly aged metal pools computed assuming equilibrium among all three soil pools. Soil porewater chemistry setup accounts for the possibility of calcite control of porewater pH. Initial porewater speciation is computed using fixed concentrations of Na, Mg, Cl, SO₄, and DOC, and fixed pH and pCO₂, and adjusting the Ca concentration to achieve charge balance. If the system is oversaturated with respect to calcite, CaCO_{3 (s)}, then the Ca²⁺ activity is fixed according to the pCO₂ and calcite solubility product, and the porewater is re‐speciated adjusting pH to achieve charge balance. A standard solubility product of 10^−8.48 and an enthalpy change of −8.0 kJ/mol are used for CaCO_{3 (s)} precipitation.

The model is then run in dynamic mode for the desired simulation period. The metal mass balance is given by 10 where ΔM_layer is the annual change in the metal pool within a layer and F_input is the total annual metal input from all sources. The replenishment flux is assumed to be constant through the simulation period at the value computed at steady state. All the flux terms have units of mol/m².

Model Application

In this study we apply the IDMM to simulate Cu, Zn, and Cd dynamics in GTR topsoils (0–20 cm depth) with maize cultivation as the land use. We primarily aimed to predict generic patterns of metal dynamics across the catchment. Simulations were done using individual soil compositions in order to assess cross‐catchment uncertainty in the “average” metal dynamics. For consideration of model sensitivity and future projections, including risks, simulations were done using a single soil composition calculated using the arithmetic mean pH and geometric mean % soil organic carbon (SOC) were used for modeling. We term these catchment‐mean simulations.

The starting year for simulation was determined by studying the trends in estimated anthropogenic metal input rates (see Time series of metal inputs) to establish a point at which anthropogenic inputs are likely to be negligible relative to natural inputs. On the basis of the estimated anthropogenic inputs, we used 1900 as the starting year and assumed that anthropogenic inputs were negligible to 1950. Although there will certainly have been pre‐1950 agricultural activity within the catchment, which will have influenced soil physicochemical and hydrological properties, for simplicity, and given the lack of available information on pre‐1950 land use, we have assumed that the soils were in a steady state with respect to natural metal inputs in 1900. We have assumed that physicochemical and hydrological variables from 1900 onwards are constant at their present day values, with the exception of the SOM content.

The purpose of this work, in addition to simulating present day metal concentrations, is to make projections of future metal accumulation under different scenarios of future inputs. Therefore, we have run simulations to the year 2060. The scenarios of metal inputs from 2009 to 2060 are described in Metal Input Scenarios: Future inputs.

Time series of metal inputs

We assumed that metal inputs to the soils derive from four possible sources: atmospheric deposition (both natural and anthropogenic), impurities in applied fertilizers, application of manures containing metal derived from animal feed additives, and irrigation using surface waters.

Robust data on metal inputs to the catchment soils over time are absent, necessitating the generation of plausible estimates of inputs. To do this, we adopted the general assumption that inputs correlated with (1) industrial production, in the case of anthropogenically derived atmospheric deposition; (2) fertilizer use and metal content, in the case of fertilizer inputs; (3) manure application rates and manure metal concentrations; (4) irrigation volumes and irrigation water metal concentrations. To derive trends in industrial production and fertilizer use, we used data available in the national and provincial Statistical Yearbooks. Yearbook figures are typically available both for China as a whole and for individual provinces, although the timescales of data availability are typically shorter for provinces than for the whole country. Data were not available for land areas smaller than province level. The GTR catchment has sections within both the Hebei and Beijing provinces. We chose Hebei as the source for province‐level data, as we considered the trends in industrial and agricultural activity to be more plausibly representative than those for Beijing municipality, which is dominated by the urban area of Beijing and thus less representative of the desired trends in industrial and agricultural activities.

Details of the estimation of time trends in metal inputs are provided in Appendix S1.

Crop yield

Dry weight yields of maize were available at the national scale for 1949–2009 (National Bureau of Statistics of China, 1981–2010) and for Hebei province from 1985 to 2009 (Bureau of Statistics of Hebei Province, 1996–2010). National and Hebei yields for 1985–2009 were significantly correlated, thus national trends were used to estimate Hebei yields for 1949–2009 by linear regression (Appendix S1: Fig. S7).

Copper, zinc, and cadmium concentrations in maize cobs, estimated from Wang et al. (2008) and Zhang et al. (2010) are shown in Table 1, along with the value of f_layer used.

Soil hydrology and erosion

Annual water fluxes from the topsoil, as vertical drainage and surface runoff (m) were estimated using the Pesticide Root Zone Model (PRZM; Carousel et al. 2003). Information on the PRZM simulation, and a comprehensive set of input parameters used, are given in Appendix S1. The simulated mean water balance is shown in Fig. 3.

short-legendFigure 3.

Surface erosion losses were computed for a square field of 0.5 ha (5000 m³) area and three land slopes: 0%, 1%, and 5%, to allow modeling of IDMM sensitivity to erosion rate. Computed mean annual surface erosion losses for these slopes were 0.000, 0.027, and 0.193 kg/m², respectively. These figures correspond to concentrations of eroded soil in runoff of 0.0, 1.2, and 8.4 g/dm³, respectively. In running the IDMM, the default land slope used was 0%.

Metal Input Scenarios

Metal sources

We constrained the number of possible metal input scenarios by using the land use information provided by Zhang et al. (2013) for the Guanting catchment. Fertilizer use was assumed to be a mixture of N, P, and K fertilizers and manure inputs were assumed to be of either chicken or cattle manure. We simulated two scenarios of input, both comprising inputs due to atmospheric deposition, irrigation, and fertilizer, together with either chicken manure (scenario S1) or cattle manure (scenario S2).

Future inputs

For scenario S1, we simulated soil response to possible future scenarios of metal input, from 2010 to 2060. We derived three trends in metal inputs:

1. inputs constant to 2060 at the input rate estimated for 2009;

2. inputs increasing linearly from 2009 to 2060, with an annual increase equal to the average annual increase from 2000 to 2009;

3. inputs increasing according to the expression 11 which produces a rising trend in emissions, but a gradually decreasing rate until a threshold is reached in 2060. The input in 2060 (F₂₀₆₀) is fixed to the mean of the input values for that year under trends (1) and (2).

An example of the projected future trends under these scenarios is shown in Appendix S1: Fig. S8.

These trends were applied to create three future scenarios: A, B, and C. In all the scenarios inputs due to irrigation were projected using trend (1). Other inputs were projected in scenarios A, B, and C using trend numbers (1), (2), and (3), respectively.

Additionally, we computed the metal input rates from 2010 to 2060 that would be required to prevent any further accumulation of metals beyond the concentrations modeled in 2009 (the “standstill load” approach; de Vries et al. 2002). Computations were done for each soil by iteratively adjusting the input rates from 2010 to 2060 until the change in predicted total metal between 2010 and 2060 was minimized.

Soil physicochemistry

The observed soil pH values were converted to porewater pHs using the expression of de Vries et al. (2004). Observed SOC was assumed to comprise 50% carbon for conversion to SOM.

We considered the influence of changes in SOM over time due to agricultural practice. Zhang et al. (2013) presented information on SOC in GTR soils (0–20 cm) in 1980 and 2005 under varying land uses. In maize‐growing land with irrigation, mean SOC was 6.78 g/kg in 1980 and 9.30 g/kg in 2005, a mean increase of 0.10 g/kg/a. This gives a net organic matter addition for the 0–20 cm layer of 0.054 kg/m²/a. Zhang et al. (2013) associate the observed SOC increases with changes in agricultural practice in the catchment around 1980, entailing increased use of manure for fertilization and decreased removal of waste crop biomass from the soil. To take these changes in SOC content into account in modeling, we fixed the 2009 SOM for each soil to its observed value and computed annual SOM for each year to be simulated, assuming it to be constant from 1900 and 1979 and to increase by 0.054 kg OM/m²/a from 1980 onwards. For the “catchment average” soil composition the same calculation was done, fixing the SOM in 2009 to the geometric mean of the observed values.

The “default” dissolved organic carbon (DOC) concentration in the soil porewater was fixed initially to 1.0 mg C/dm³, based on trial modeling which showed good agreement between observed and modeled present day metal concentrations. The porewater dissolved organic matter was modeled as 65% chemically active (as fulvic acid) and 35% inert.

Model sensitivity and uncertainty

Sensitivity computations were done with three parameters (erosion rate, DOC, pH, and soil organic matter content) varied. The erosion rates were set to the values predicted by PRZM for varying land slopes (see Soil hydrology and erosion). Soil pH was varied to two SDs below the mean value. The pristine SOM content was adjusted to high and low values, two SDs above and below the geometric mean soil SOM measured in 2009, with correction to account for temporal changes in SOM.

The DOC concentration was varied to 0.5, 2.0, 10.0, and 50.0 mg C/dm³. This wide range was chosen in order to investigate the influence of this parameter not only within the range of the default value selected from initial model testing, but also within ranges previously measured in field studies.

The influence of hydrological variability on the predictions was investigated by running the model with annually varying drainage and runoff as predicted by PRZM. A 20‐yr cycle of hydrology, simulated by PRZM based on modeling from 1975 to 1994, was used and predictions compared with those obtained assuming the constant, mean hydrology as given in Fig. 3.

Uncertainty in model predictions was evaluated by comparing the distribution of soil‐specific predictions against the distribution of measured labile metal concentrations. Statistical analysis was done using MINITAB (Minitab 2010).

Ecological risk assessment

We assessed the potential risks to soil organisms of metal accumulation in the Guanting soils by computing multisubstance Potentially Affected Fractions (msPAFs) for Cu, Zn, and Cd. The Potentially Affected Fraction (PAF) for a single substance is calculated from the species sensitivity distribution (SSD) for that substance and is an estimate of the proportion of the species within the ecosystem adversely affected by that concentration of the substance (Fig. 4), that is, a measure of toxic pressure on the system (Van Straalen 2001). Here, we have used the mechanistic approach of Lofts et al. (2005) to calculate single substance PAFs for the metals. This approach calculates PAFs as a function of the labile metal concentration above the natural background and accounts for metal bioavailability on the basis of the soil porewater pH and organic matter content. It also allows for the calculation of a PAF at a defined level of statistical confidence, after Aldenberg and Jaworska (2000).

short-legendFigure 4.

Potentially Affected Fractions were computed for the IDMM projections of metal concentrations from scenario S1 and future input scenarios A, B, and C, using catchment‐mean predicted metal concentrations. We assumed the computed labile metal concentrations at steady state to be the natural background concentrations, thus by definition PAF = 0 for all the metals at steady state. The PAFs were calculated individually for the metals at three confidence levels (50%, 15.9%, and 84.1%), where a given confidence level represents the probability that the true PAF is smaller than the value calculated. Single substance PAFs were combined to give msPAFs at each confidence level, assuming independent combined toxic action of the metals: 12

Results

Soil physicochemistry

The geometric mean organic carbon content of the soils was 0.92% w/w, giving a mean organic matter content of 1.84% w/w assuming organic matter to be 50% carbon. The mean measured pH (water extraction) in the study area was 7.65, giving a geometric mean pH_pw of 7.72. Individual soil properties and metal concentrations are given in Table 2. Mean total metal concentrations were in the order Zn > Cu > Cd, following the general crustal abundance of these metals. Geometric mean EDTA‐extractable concentrations were in the order Cu > Zn > Cd. Extractable metal as a proportion of the total followed the trend Cd > Cu > Zn. Extraction with 0.05 M EDTA has been previously shown (Marzouk et al. 2013a) to be a reasonable method for extracting the labile or “geochemically active” pool of soil metal. In the slightly alkaline soils of the GTR catchment, the proportions of Cu, Zn, and Cd extractable by EDTA were broadly similar to proportions of isotopically labile forms of these metals previously measured in agricultural and mining‐impacted soils of similar pH (e.g., Smolders et al. 2012, Marzouk et al. 2013b).

Table 2. Properties and metal concentrations (total and EDTA‐extractable Cu, Zn, and Cd) of field soils (0–20 cm) in the Guanting catchment

Sample	pHH2O	pHpw^†	TOC (%)	Bulk density (g/cm^3)	Metal concentrations (μg/g d.w. soil)
					Copper		Zinc		Cadmium
					Total	EDTA‐extractable	Total	EDTA‐extractable	Total	EDTA‐extractable
GTR‐01	7.63	7.70	1.60	1.24	21.4	2.67	77.0	2.41	0.32	0.086
GTR‐07	7.79	7.87	0.88	1.42	15.6	2.89	63.8	1.38	0.27	0.080
GTR‐09	7.97	8.05	0.59	1.54	11.5	3.01	59.4	1.59	0.27	0.082
GTR‐33	7.70	7.77	1.04	1.37	12.8	2.03	48.3	1.93	0.31	0.074
GTR‐37	7.54	7.60	0.49	1.59	9.12	1.49	54.6	0.94	0.29	0.068
GTR‐41	7.66	7.73	0.80	1.45	12.3	2.73	50.9	1.82	0.30	0.096
GTR‐42	7.40	7.46	0.83	1.44	14.8	3.26	53.6	3.13	0.61	0.19
GTR‐58	7.56	7.62	1.12	1.35	19.7	2.90	57.9	1.56	0.32	0.093
GTR‐60	7.65	7.72	1.38	1.29	16.2	2.49	58.1	2.76	0.40	0.074
GTR‐61	7.64	7.71	1.00	1.38	14.8	2.01	55.3	2.29	0.30	0.072

† Calculated from the measured pH_H2O: pH_pw = 1.0462pH_H2O − 0.2847 (de Vries et al. 2005).

Metal inputs

Fig. 5 shows the time trends in estimated metal inputs. For all sources except irrigation there is a continuous increase in estimated inputs, commensurate with the pattern of industrial expansion and agricultural intensification in China since 1949. The atmospheric deposition increases steadily from the earliest available estimate (1953) and the rate of increase itself accelerates in the early 2000s. The fertilizer‐associated input rates remain relatively flat until 1960, when they also begin to rise steadily. Irrigation inputs show a maximum in the 1980s.

short-legendFigure 5.

Fig. 6 shows the magnitudes of estimated metal inputs from each source in 2009. Atmospheric deposition is the largest single input for all the metals. The estimated deposition fluxes for 2009 are 136, 944, and 2.08 μmol/m², which are broadly comparable with the values estimated by Luo et al. (2009), for China as a whole, of 170, 989, and 3.56 μmol/m². The Guanting inputs are, respectively, 156, 700, and 104 times greater than those estimated for geogenic deposition. For Cu and Zn, chicken manure application is estimated to be the second largest input, but is still relatively small, being approximately one‐fifth of atmospheric deposition for both metals. The smallest estimated input is consistently NPK fertilizer, at 2.5–5% of atmospheric deposition.

short-legendFigure 6.

Model simulations

Simulations with scenario S1, including uncertainty on the predicted labile metal, are shown in Fig. 7. The model predicted an increase in both total and labile metal concentrations in response to anthropogenic inputs. In all cases the observed (present day) geometric mean labile metal was within a single SD of the predicted geometric mean. Mann–Whitney testing indicated no significant differences between observed and modeled labile metal (df = 18; Cu: W = 112.0, P = 0.623; Zn: W = 107.0, P = 0.910; Cd: W = 116.0, P = 0.427). Total Cu, Zn, and Cd concentrations in 2009 were predicted to have increased by 0.62, 4.3, and 0.018 mg/kg, respectively. These represent increases of 4.5%, 8.0%, and 5.6%, respectively, over the pristine total metal concentrations. Mann–Whitney testing showed no significant differences between the observed present day concentrations and the estimated pristine total concentrations (df = 18; Cu: W = 98.0, P = 0.623; Zn: W = 88.0, P = 0.212; Cd: W = 91.0, P = 0.308). This indicates that when the measured cross‐catchment variability in total metal is accounted for, the predicted changes in total metal are not statistically significant. The labile metal concentrations were predicted to increase by 0.34, 0.95, and 0.009 mg/kg, respectively, representing increases of 17%, 118%, and 15% over the pristine mean labile metal concentrations. Mann–Whitney testing indicated that the increases in labile Cu and Cd were not significant when cross‐catchment variability was taken into account (df = 18; Cu: W = 95.0, P = 0.473; Cd: W = 96.0, P = 0.521), whereas the increases in labile Zn were significant (df = 18, W = 58.0, P = 0.0004).

short-legendFigure 7.

Comparison of the increases in total and labile metal concentrations shows that the model is predicting extensive transformation of the input metal into aged forms (Fig. 8). The proportions of Cu and Cd predicted to have entered the aged pool are comparable to those remaining in the labile pool; in the case of Zn, the majority of metal retained in the soil was predicted to be nonlabile. Total losses due to drainage and cropping were below 10% for all the metals. Cropping was predicted to be a more important loss process than drainage for Cu and Zn.

short-legendFigure 8.

All the future scenarios (Fig. 9) projected further metal accumulation, including the “best case” scenario of no future increase in inputs. This emphasizes the high metal retention characteristics of the GTR soils. Mean “standstill” metal input rates for 2010–2060 were 12.9, 74.1, and 0.124 μmol/m², respectively. These are, respectively, 7.3%, 6.3%, and 5.0% of the estimated inputs for 2009.

short-legendFigure 9.

Model sensitivity

The sensitivity of present day catchment‐mean labile metal concentrations to model variables is shown in Table 3. Reducing the soil pH reduced the predicted labile Cu and Cd concentrations somewhat, while predicting a slightly higher labile Zn concentration. Adjusting the soil OM had the largest effect on the predicted soil concentrations of all the metals. Increasing soil OM increased the predicted labile Cu, Zn, and Cd concentrations while reducing the soil OM had the opposite effect. The model response to adjustments in soil properties was largely due to the effect on the predicted pristine, steady‐state concentrations—effects on prediction accumulation rates were small. Varying the erosion rate had a negligible influence on the predicted concentrations, which is expected given the simulated replenishment of eroded soil to maintain a constant layer depth.

Table 3. Catchment‐mean simulated labile Cu, Zn, and Cd concentrations in 2009, based on different metal input and parameter scenarios

Simulation code	Metal inputs included in simulation^†					Soil physicochemistry^‡	Predicted labile concentration in 2009 (μg/g)
	AD	IRR	NPK	CH	CT		Cu	Zn	Cd
Pristine							1.97	0.81	0.059
S1	×	×	×	×			2.31	1.76	0.067
S2	×	×	×		×		2.28	1.74	0.068
S3	×	×	×	×		[DOC]pw = 50 mg/dm^3	0.28	1.06	0.011
S4	×	×	×	×		[DOC]pw = 10 mg/dm^3	0.49	1.51	0.036
S5	×	×	×	×		[DOC]pw = 2.0 mg/dm^3	1.42	1.64	0.064
S6	×	×	×	×		[DOC]pw = 0.5 mg/dm^3	3.85	1.77	0.069
S7	×	×	×	×		pHpw = 7.40	2.10	1.82	0.058
S8	×	×	×	×		Soil OM = 3.18% w/w	4.77	2.72	0.146
S9	×	×	×	×		Soil OM = 0.29% w/w	0.90	1.19	0.025

† Natural atmospheric deposition always included. AD, anthropogenic atmospheric deposition; IRR, irrigation; NPK, nitrogenous, phosphate, and potash fertilizer; COMP, compound fertilizer; CH, chicken manure; CT, cattle manure.

‡ Unless stated, physicochemical parameters were set to their default values (see main text for details).

Varying the porewater [DOC] in the range 0.5–50 mg C·dm⁻³ caused the predicted labile Cu to vary by a factor of 13.8, that is, over an order of magnitude. Smaller variations in labile Zn and Cd were seen, with factors of variation being 1.7 and 6.3, respectively. These differences in present day predictions were due to predicted differences in the computed steady‐state labile metal. Using a higher [DOC] concentration gave lower predicted labile metal concentrations. The relatively large variability for Cu results from the relatively strong association of this metal with DOC.

Simulations done using variable hydrology had a negligible influence on predicted concentrations. An example is given in Appendix S1: Fig. S9.

Ecological risk assessment

Trends in msPAF for input scenario S1 and future scenarios A, B, and C are shown in Fig. 10. The figure shows the msPAF of 0.05, which is typically considered to be the threshold above which potential risks to the ecosystem are indicated. In none of the future scenarios is the msPAF at a 50% confidence level (msPAF₅₀) predicted to exceed this threshold, although in scenarios B and C the upper confidence level on the msPAF reaches 0.05 before 2060.

short-legendFigure 10.

Fig. 11 shows the msPAFs predicted in 2060, along with the single substance PAFs for Cu, Zn, and Cd for the same year. The PAFs for Cd are small, being <10⁻⁴ for all three future scenarios. The Zn PAFs are the highest, being between five and eight times higher than the Cu PAFs and three to four orders of magnitude higher than the Cd PAFs.

short-legendFigure 11.

Discussion

Estimation of soil metal inputs suggested that atmospheric deposition is the main source of Cu, Zn and Cd to the GTR catchment. In the case of atmospheric deposition, this estimation was critically dependent on the availability of deposition data within the catchment. For broader scale application, larger scale measurements or modeled estimates of atmospheric metal deposition are essential.

The agreement obtained between the measured EDTA‐extractable metal concentrations, and the predicted labile metal concentrations, was most dependent upon the porewater DOC concentration used in the simulations. Using a DOC concentration of 1.0 mg·C·dm⁻³ produced good agreement between observations and predictions. This concentration is lower than those previously reported in situ in the porewaters of agricultural soils (e.g., Keller et al. 2002, Ostermann et al. 2015). Agreement of predictions with observations is clearly reliant on the prediction of steady‐state pristine labile metal, which is strongly dependent upon the porewater DOC, particularly for Cu. There is currently insufficient knowledge on the variables controlling in situ porewater DOC concentrations in agricultural soils. This is a priority for future research and monitoring, since it currently limits more extensive application of the model for risk assessment. Broader monitoring and modeling studies on a wider range of well‐characterized systems are needed to assess the degree to which the assumption of pristine steady state is generally reasonable.

On the basis of the modeling, the Guanting catchment soils are not currently to be at ecological risk. The model predicts that accumulation has raised the total soil metal concentrations by <10% in approximately 60 yr, suggesting that most of the soil metal is still natural in origin. The model predicts greater proportional changes in the labile pools, which are thus more closely related to risks to soil health than the total metal (e.g., Lofts et al. 2004). Despite this, on a catchment‐wide basis only the labile zinc concentration was predicted to have increased significantly to 2009. Continued metal accumulation, which is predicted to occur even under the “best case” scenario of no further increase in metal inputs, would if allowed to continue eventually result in risk being predicted, although exceedance of the PAF₅₀ is not predicted to occur up to 2060.

The role of aging in removing metal from the labile pool is of general significance for metal dynamics, bioavailability, and risks to soil health. In the previous application of the IDMM (Lofts et al. 2013), aging reactions were not represented, and weathering inputs were simulated as constant fluxes from an infinite mineral pool. The current version simulates both rapid aging as demonstrated in laboratory studies (Crout et al. 2006, Ma et al. 2006) and temporal changes in the total metal pool. The model is based around the hypothesis that in a well‐developed soil, long‐term aging, and weathering fluxes are in steady state prior to anthropogenic inputs. This is a parsimonious formulation, providing a straightforward means for the model to simulate temporal changes in the total metal pools. Time series data on soil metal concentrations and solid‐phase speciation in soils, ideally in the context of both increasing and decreasing inputs over time, would be required to further test the hypothesis of pristine steady state among the metal pools.

The model predicts that metal removal by cropping is of comparable significance to losses in leaching. This is in agreement with other modeling studies (e.g., de Vries et al. 2004, Groenenberg et al. 2006) and emphasizes the need to consider plant uptake of metal, both for budgeting purposes and in the case of crops to calculate concentrations in edible parts for human health risk assessment. Studies exist (e.g., Hough et al. 2003) showing influences of soil chemistry and metal concentration on metal uptake by crops; however, given the relatively limited scope of such studies, we took the simplest possible approach to calculating uptake by assuming a constant concentration of metal in cropped material. This is an oversimplification, and further research is needed to provide more generally applicable expressions relating crop metal contents to soil chemistry and metal concentrations.

Combining the IDMM outputs with the PAF approach provides a powerful tool for estimation of future risks to soil health due to long‐term metal inputs. The adaptation of the approach by Lofts et al. (2005) is particularly amenable to linking with the IDMM since it computes risk on the basis of the “added” labile metal concentration, that is, the increment above the natural background. This can be readily computed from the IDMM outputs assuming the “pristine” steady‐state metal concentration to be the natural labile background. Furthermore, the PAF approach allows the computation of combined risk due to multiple metals, it provides a more thorough assessment of risk than do approaches that consider only single substances, and provides information on the relative risks due to each substance so allowing prioritization.

A long‐term goal for IDMM development is to simulate spatial patterns of metal dynamics and ecological/human health risks across different agricultural soil types. This would allow future projections of metal accumulation to be evaluated against land use practices. There is a general lack of spatially explicit, land use‐specific information relevant to estimating metal inputs, but some progress has been made in estimating generic inputs and outputs for general farm types across Europe (e.g., Eckel et al. 2008). An inventory of metal inputs to Chinese agricultural soils has been computed (Luo et al. 2009) but considers only nationally averaged values. Given the apparent dominance of metal inputs due to atmospheric deposition, knowledge of larger‐scale spatial patterns in deposition of metals to Chinese soils is the most pressing current need. Better knowledge of the spatial variability in other metal sources is also needed, particularly as the relevant importance of sources may also vary spatially. Generating spatially resolved knowledge on fertilizer and manure‐related inputs will be more challenging. A more flexible approach, for example by using a probabilistic consideration of inputs, may prove valuable.

Conclusions

1. Estimates of historic metal inputs to Guanting catchment soils indicate a steady increase in inputs from 1950 to the present day. Estimation strongly suggests that atmospheric deposition has been and remains the dominant source of anthropogenic metal inputs to the catchments soils.

2. The Intermediate Dynamic Model for Metals predicted present day mean topsoil (0–20 cm) concentrations of labile (“geochemically active”) Cu, Zn, and Cd well in the GTR soils, using plausible estimates of unknown parameters.

3. Modeling estimates that anthropogenic metal inputs from 1950 onwards have raised Cu, Zn, and Cd concentrations in the soils by <10% over the natural concentrations, however, the model predicts continued accumulation of metal even if input rates are maintained at present day (2009) levels.

4. The input rates required to prevent further accumulation of Cu, Zn, and Cd are all over one order of magnitude lower than the present day rates.

5. Modeling predicts that aging reactions of metals are important in limiting their accumulation in labile, bioavailable forms.

6. Predictions are most sensitive to the porewater DOC concentration. Improved knowledge of in situ DOC fluxes from agricultural topsoils is a research priority.

7. Projections of combined toxic pressure using the PAF approach suggest that potential risks to soil health in the Guanting catchment may arise on a decadal timescale if current input trends continue. Of the three metal studied, Zn makes the largest contribution to risk.

8. The coupled IDMM‐PAF approach has considerable potential for evaluation of spatial and temporal patterns in metal dynamics and ecological risks in agricultural soils, given suitable driving data.

Acknowledgments

This study was supported by the International Scientific Cooperation Program with Grant No. 2012DFA91150, the National Natural Science Foundation of China under Grant No. 414201040045 and No. 41371488, the Key Project of the Chinese Academy of Sciences under Grant No. KZZD‐EW‐TZ‐12, and the UK Natural Environment Research Council.