Modelling of long-term Zn, Cu, Cd and Pb dynamics from soils fertilised with organic amendments

Soil contamination by trace elements (TEs) is a major concern for sustainable land management. A potential source of excessive inputs of TEs into agricultural soils are organic amendments. Here, we used dynamic simulations carried out with the Intermediate Dynamic Model for Metals (IDMM) to describe the observed trends of topsoil Zn (zinc), Cu (copper), Pb (lead) and Cd (cadmium) concentrations in a long-term (>60-year) crop trial in Switzerland, where soil plots have been treated with different organic amendments (farmyard manure, sewage sludge and compost).

The observed ethylenediaminetetraacetic acid disodium salt (EDTA)-extractable concentrations ranged between 2.6 and 27.1 mg kg⁻¹ for Zn, 4.9 and 29.0 mg kg⁻¹ for Cu, 6.1–26.2 mg kg⁻¹ for Pb, and 0.08 and 0.66 mg kg⁻¹ for Cd. Metal input rates were initially estimated based on literature data. An additional, calibrated metal flux, tentatively attributed to mineral weathering, was necessary to fit the observed data. Dissolved organic carbon fluxes were estimated using a soil organic carbon model. The model adequately reproduced the EDTA-extractable (labile) concentrations when input rates were optimised and soil lateral mixing was invoked to account for the edge effect of mechanically ploughing the trial plots. The global average root mean square error (RMSE) was 2.7, and the average bias (overestimation) was −1.66, −2.18, −4.34 and −0.05 mg kg⁻¹ for Zn, Cu, Pb and Cd, respectively. The calibrated model was used to project the long-term metal trends in field conditions (without soil lateral mixing), under stable climate and management practices, with soil organic carbon estimated by modelling and assumed trends in soil pH. Labile metal concentrations to 2100 were largely projected to remain near constant or to decline, except for some metals in plots receiving compost. Ecotoxicological thresholds (critical limits) were predicted to be exceeded presently under sewage sludge inputs and to remain so until 2100. Ecological risks were largely not indicated in the other plots, although some minor exceedances of critical limits were projected to occur for Zn before 2100.

This study advances our understanding of TEs' long-term dynamics in agricultural fields, paving the way to quantitative applications of modelling at field scales.

1 Introduction

Trace elements (TEs) are naturally present in soils due to mineral weathering and biogeochemical cycles. Several TEs, such as zinc (Zn), copper (Cu) and nickel (Ni), play important roles in biochemical processes and are essential for living organisms at low concentrations, though they can become toxic to biota at high concentrations; therefore, their presence in soil can be tolerable in a relatively narrow range of values (Adriano, 2005). Wan et al. (2020) and Wang et al. (2015) reported that the concentrations providing ecological safety in China ranged from 38.3 to 263.3 mg kg⁻¹ for Zn and from 13.1 to 51.9 mg kg⁻¹ for Cu when shifting from acidic soils to alkaline non-calcareous soils. In contrast, other TEs, such as lead (Pb) and cadmium (Cd), which are not physiologically active, may be toxic to living organisms at low concentrations, and their accumulation in soil is of particular concern. For example, background cadmium levels in world topsoils range from 0.01 to 2.7 mg kg⁻¹, and in Europe the mean cadmium concentration in cultivated soils is 0.5 mg kg⁻¹ (EFSA, 2009). Excessive uptake of trace elements by crop plants and enrichment in edible parts can pose significant risks to human health by entering into the food chain (McGrath et al., 2015).

Accumulation of TEs in cultivated soils is widespread and is mainly caused by application of low-grade agrochemicals, organic fertilisers and sewage sludge (Toth et al., 2016). In an EU-wide survey, Ballabio et al. (2018) reported that agricultural soils have among the highest potential to become enriched in Cu compared with other land uses, and that land cover and management are better predictors of soil Cu concentrations than natural soil formation factors.

Organic amendments are considered more sustainable than inorganic mineral fertilisers (Diacono and Montemurro, 2010); in fact, current industrial processes for N-fertiliser production are based on fossil fuels, and P fertilisers are manufactured from phosphate rocks which are naturally limited (Roberts, 2014). However, the application of organic amendments, such as farmyard manure, compost and digestates of bio-wastes, can also introduce TEs into agricultural soils. In the EU, typical levels of trace elements in cattle manure are 63–175 mg kg⁻¹ for Zn, 15–75 mg kg⁻¹ for Cu, 1.4–4.3 mg kg⁻¹ for Pb and 0.1–0.4 mg kg⁻¹ for Cd (NEBRA, 2015), but pig and poultry manure can be much more enriched in trace metals. Mean concentrations of Zn, Cu, Pb and Cd in green waste compost in Germany were reported to be 168, 33, 61 and 0.7 mg kg⁻¹, respectively (NEBRA, 2015). Application of sewage sludge into agricultural soils can be even more problematic as sewage sludge can contain trace metals up to 30 times of the amount of their concentrations in soil (Hudcova et al., 2019).

Once in the soil, multiple factors may control TE speciation, solubility and mobility (Gu and Evans, 2008), such as soil pH, soil and dissolved organic matter (SOM and DOM) contents and the quantity and chemical composition of clay minerals and metal (oxy)hydroxides. Furthermore, speciation can influence the TE toxicological hazard, particularly to organisms that are directly exposed to soils, such as plants and earthworms. Models, such as the biotic ligand model (Paquin et al., 2002), postulate that metal toxicity is related to the uptake of specific metal species in competition with other solution ions, rather than to total dissolved metals. Repeated applications of organic amendments can lead to the accumulation of TEs in agricultural soils, particularly through direct reactions with the soil solids (adsorption), formation of precipitates, or physical occlusion within the organo-mineral aggregates (fixation). This can lead to TE concentrations exceeding environmental legislation thresholds.

In the context of long-term TE accumulation due to regular application of organic amendments or other additions, predicting the long-term speciation and dynamics of TEs is useful to support decisions on ecosystem management and human health protection. Dynamic models, if reliable, are essential for this purpose. Models for TE dynamics exist at a number of levels of complexity, from those with a mechanistic approach requiring highly detailed input information and calibration (Bonten et al., 2011), to simple mass balance approaches generally applicable at large scales but relatively unsuitable for understanding and unravelling complex metal dynamics (Six and Smolders, 2014). Empirical models have been used to simulate the dynamics and uptake of TEs at specific agricultural sites based on site-specific calibration, (Bergkvist and Jarvis, 2004; Ingwersen and Streck, 2006), but such models can lack a reliable generalisation of the parameters to different climatic conditions and hydrological and soil physio-chemical properties.

Among models for determining TE dynamics in soil, the Intermediate Dynamic Model for Metals (IDMM; Lofts et al., 2013; Xu et al., 2016) is an example of a semi-empirical dynamic model which allows general application, given a reasonably parsimonious set of input data. It is intended for long-term application from decades to centuries. The IDMM describes metal dynamics from a past year in which metal inputs can be assumed to be uninfluenced by anthropogenic activities. Processes influencing metal dynamics, including solid solution partitioning, fixation into soil solid phases and leaching, are described in a semi-empirical manner. Hydrological variables (e.g. annual volume of soil drainage) are specified as a time series. Similarly, key properties influencing metal dynamics, such as the pH of the soil solution, dissolved organic carbon (DOC) flux, SOM content and soil erosion rate, may be fixed to single values or varied annually.

The objective of this study was to assess the capability of the IDMM to reproduce metal dynamics at a well-characterised location receiving a range of organic amendments. This application of the model is at a smaller and more detailed scale than previous evaluations (Lofts et al., 2013; Xu et al., 2016). We simulated the dynamics of Zn, Cu, Pb and Cd in the topsoil of a long-term (>60-year) agricultural trial in Switzerland, which comprises a series of plots receiving either farmyard manure, sewage sludge, green waste compost, or no amendment. After optimising the model parameters and driving the data for the site, we then made projections of metal dynamics into the future, under the same agronomic practices, in order to assess their sustainability in terms of environmental risk.

2 Materials and methods

2.1 The study site

The Zurich Organic Fertilization Experiment (ZOFE) is a long-term agricultural plot trial that was started in 1949 by the Swiss Federal Agricultural Research Institute (Agroscope) at Zurich-Reckenholz, Switzerland, to compare different fertilisation schemes in the following 8-year crop rotation: (1) winter wheat/intercrop, (2) maize, (3) potato, (4) winter wheat/intercrop, (5) maize, (6) summer barley, (7) clover grass ley and (8) clover grass ley (Oberholzer et al., 2014). Ploughing has been carried out to a depth of at least 20 cm, from north to south and vice versa, alternating the direction of adjacent passes (Fig. 1). The site is located at 420 m a.s.l. (above sea level), the mean annual precipitation is 1054 mm, and the mean annual temperature is 9.4 ^∘C. The soil is a carbonate-free, loamy (14 % clay) Luvisol (IUSS Working Group WRB, 2006), with a soil organic carbon (SOC) content of 1.43 % $w/w$ and a pH (H₂O) value of 6.5 prior to the experiment. The field trial consists of 12 treatments replicated in five blocks, of dimensions 7 m × 5 m, in a systematic design. The same cultivation and plant protection regimes have been applied to all the treatments. In the present study, we investigated the following four treatments: control (NIL no. 1), with no fertilisation and no amendment, farmyard manure (FYM no. 2), with an application of 5 t OM ha⁻¹ every second year, sewage sludge (SS no. 3), with an application of 2.5 t OM ha⁻¹ every year, and green waste compost (COM no. 4), with application of 2.5 t OM ha⁻¹ every year (Fig. 1).

Figure 1 The experimental design of ZOFE, with 12 treatments replicated in five blocks. Only the four treatments highlighted in grey were investigated.

2.2 Trace element time series

The NIL, FYM, SS and COM soils (top 20 cm; sieved to 2 mm) were sampled from the Agroscope ZOFE soil archive and analysed for total and ethylenediaminetetraacetic acid disodium salt (EDTA)-extractable concentrations of Zn, Cu, Pb and Cd. Analysed soils were sampled from the years 1972, 1979, 1982, 1991, 1995, 2000, 2003, 2007 and 2011. Before 2011, samples from the five replicate plots per treatment had been bulked, so variability among replicate plots could not be assessed. To determine the total and EDTA-extractable soil TE concentrations, respectively, soils were digested in aqua regia and extracted using EDTA and analysed by inductively coupled plasma optical emission spectrometry (ICP-OES; sequential PerkinElmer Optima 2000 DV). The EDTA-extractable pools were obtained with the extraction protocol described by Quevauviller (1998). The total metal concentrations were compared with the one-point-in-time measurements from the same plots carried out with inductively coupled plasma mass spectrometry (ICP-MS) from an independent laboratory, so that interferences of As with Cd in the readings (McBride, 2011) were ruled out (as shown in Fig. S1 of the Supplement for Cd). Quality control of the ICP-OES was done every 10 readings on the calibration curve by measuring the TE concentrations in the blank samples and in the standard sample at concentration of 1 ppm (parts per million). The limits of quantification of total and EDTA-extractable concentrations for each TE are given in Table S1 (see the Supplement).

Samples of farmyard manure from 2011 and 2014, sewage sludge from 2008 and 2012 and compost from 2011, 2013 and 2014 were also analysed for total and EDTA-extractable concentrations of Zn, Cu, Pb and Cd, as described above.

2.3 The Intermediate Dynamic Model for Metals with lateral mixing

A detailed description of the Intermediate Dynamic Model for Metals (IDMM) is given in the Supplement; here, we present an overview of the most relevant concepts. The IDMM simulates annual concentrations of total and geochemically active metals within a topsoil and metal fluxes from the soil due to porewater leaching and crop uptake. It distinguishes between a pool of geochemically active (labile) TEs, comprising dissolved and adsorbed forms, and a nonlabile (aged) pool that comprises chemically less reactive and occluded forms (Fig. 2a). The labile metal pool is partitioned into dissolved and adsorbed forms, assuming chemical equilibrium. A Freundlich-type isotherm (Groenenberg et al., 2010) describes the relationship between free and adsorbed TE ions, and the relationship between free TE ions and TEs complexes in the porewater is computed using the Windermere Humic Aqueous Model (WHAM) model VI (Tipping, 1998). Transformations between labile and aged pools follow first-order kinetics. A total of two pools of aged metal are defined and termed “weakly aged” and “strongly aged”. The weakly aged pool has relatively rapid, reversible transformation kinetics with the labile pool. Weakly aged metal may transfer into the strongly aged pool, and strongly aged metal may transfer back into the labile pool. Both of these transformations have relatively slow kinetics compared to the transformations occurring between the labile and weakly aged pools.

Figure 2 The IDMM model with lateral mixing. (a) The description of the model structure. (b) The description of the soil lateral mixing method.

The IDMM is driven by annual TE input rates, which we assume to be entirely in labile form. Simulations start from a past year (1750 in this case) in which all metal inputs are assumed to be natural and where the soil is in steady state, i.e. metal input and output fluxes balance (Tipping, 1998). Metals can be removed from the soil due to leaching of the dissolved form in drainage water, erosion of the soil, and uptake into the harvestable parts of crops. In this study, erosion was neglected in the consideration of the site's geomorphological characteristics (i.e. the site has a negligible slope). Since the soil samples were relevant to the homogenised ploughing depth, the soil was modelled as a single, well-mixed layer of 20 cm.

The work of McGrath (1987) and McGrath and Cegarra (1992) on plot experiments at the Rothamsted Research (UK) showed that the lateral movement of soil among adjacent plots, due to regular ploughing, exerted a significant influence on the temporal trends of metal concentrations within the plough layer. Therefore, we extended the IDMM to enable the lateral transfer of metal across the plots to be simulated. Figure 2b illustrates the lateral mixing model set-up. The plots are subdivided into a number of strips within which the soil composition is physicochemically homogeneous. Each year, a proportion of the topsoil, and its associated TEs, in each strip is exchanged with the adjacent strip. Mean soil TE concentrations in each plot are computed as the sum of the TE mass in each strip divided by the sum of the soil mass in each strip. The number of strips per plot, and the amount of soil being exchanged, can be varied. Figure 2b shows an example set-up with five strips of 1 m width and an exchange width of 0.2 m, equal to the plough depth. For simulations, we ordered the plots as NIL-FYM-SS-COM. We fixed the exchange width to this value for all simulations. A sensitivity analysis was carried out on the number of strips per plot in order to understand the impact of the choice of strip number on the predicted TE trends.

2.4 Metal input rate estimation

Besides TE inputs from the organic amendments, TE inputs to the plots were assumed to comprise geogenic and anthropogenic deposition from the atmosphere and mineral weathering from the coarse (>2 mm) fraction of the soil. Estimation and optimisation of the input rates, based on available measurements and data from the literature, are described in detail in the Supplement. For the sewage sludge amendment only, two approaches were used to estimate the metal inputs. The first approach, termed the “Swiss trend”, was derived from the literature trends specific to Switzerland; considering the unsatisfactory simulation results, the metal input rates were adjusted to the ZOFE plots, and this approach was termed the “idealised trend” (more details in the Supplement).

2.5 Model parameterisation

The IDMM requires the time series inputs of annual soil porewater pH, SOM and DOC fluxes in porewater and annual drainage volume in order to simulate soil TE concentrations. For computing TE losses by crop uptake and removal, annual crop yields and TE concentrations in the harvestable parts of crops are required.

Soil porewater pH trends for 1949–2014 were obtained by, firstly, converting measurements in aqueous extracts to porewater pH values, according to the formula provided by de Vries et al. (2008; see Fig. S7 in the Supplement). The trends in pH were then smoothed to reduce the effect of noise on the simulations; a 5-year period, four-pass averaging approach was used for smoothing. For future projections, we ran simulations assuming either (i) no change in pH from 2014 to 2100 or (ii) an exponential decline in pH from 2014 onwards to reach a constant pH that is 0.5 units below the 2014 value.

Prior to 1949, SOM was assumed to be constant and set to the mean of the first measurement in each plot. Fitting of SOM observations before 2014 and projections to 2100 was done by applying the two-pool SOM model ICBM (introductory carbon balance model; Andrén et al., 2004). Briefly, after fitting the NIL treatment with the coefficients obtained by Menichetti et al. (2016) for ZOFE, the humification coefficient was calibrated to fit the other amended treatments. Each plot was divided into five strips, and the soil lateral mixing approach, described above, was incorporated into the modelling. The resulting trends are shown in Fig. S7. The SOC was converted to SOM by assuming the latter to be 50 % C by mass.

No porewater dissolved organic carbon (DOC) concentration or flux data were available for ZOFE. The total carbon respiration flux from 1949 onwards, obtained from the ICBM modelling, was assumed to be proportional to the DOC loss flux. The annual DOC flux prior to the start of the experiment was estimated, and DOC fluxes during the experiment were then estimated by scaling. As the site was grassland prior to the experiment, we collated published data on DOC fluxes from improved grasslands (Fu et al., 2019; Kindler et al., 2011; Buckingham et al., 2008) and derived a median flux of 8.1 g m⁻² per year for scaling. The computed flux trends were divided by the annual drainage to obtain DOC concentration trends (Fig. 3). The annual drainage was assumed to be constant and was calculated by averaging the difference between rainfall measurements at local stations and evapotranspiration estimated with a locally calibrated Primault equation.

Figure 3 Estimated DOC concentration time trends in ZOFE topsoil (20 cm) for the treatments NIL, FYM, SS and COM.

Crop metal removal was assumed to be a function of crop biomass (Fig. S7), as crop metal concentrations were assumed not to vary. Crop yields have been measured in ZOFE yearly, based on the harvest from a subplot in each plot. Shoot biomass was estimated by scaling the crop yields linearly, according to the functions provided by Bolinder et al. (2007). The Zn, Cu and Cd concentrations in winter wheat grains and shoots were measured at harvest in 2014 and 2015, and the average values were taken to represent the respective metal contents in the grains and shoots of wheat and barley over the entire simulation period. The TE concentrations for Pb, and for the other crops in the 8-year rotation, were estimated from previous reports (de Vries et al., 2008; EFSA, 2009, 2010; SAEFL, 2003; SCAN, 2003a, b). Table S8 reports the assumed TE concentrations for all crops.

2.6 Analysis of the soil and organic amendment FTIR and XRD

NIL, FYM, SS and COM amended soil samples from 1972 and 2011 and SS samples from 2013 were analysed by Fourier transform infrared spectroscopy (FTIR) to detect any change in the time of the soil organic fraction between the treatments. The diffuse reflectance infrared Fourier transform (DRIFT) spectra were obtained using a fast-scanning Spectrum GX (PerkinElmer, Monza, Italy) Fourier transform infrared spectrometer (FTIR) in the mid-infrared spectral range (4000 to 450 cm⁻¹). The spectrometer was equipped with a Peltier-cooled deuterated triglycine sulfate (DTGS) detector and an extended-range KBr beam splitter. Soil samples of 50 mg were placed in a stainless steel sample cup, located in a PerkinElmer diffuse reflectance accessory and scanned for 60 s. A silicon carbide (SiC) reference disk was used as the background sample (PerkinElmer). The most noticeable peaks were attributed according to D'Acqui et al. (2015) and Niemeyer et al. (1992), as reported in the Supplement. The same soil samples were also analysed by X-ray diffraction (XRD) to investigate the soil mineral fraction, as also described in the Supplement.

2.7 Projections of the long-term effects of organic amendment applications

After evaluating and calibrating the IDMM for the period 1949–2014, the model was run in predictive mode from 1750–2100 to derive projections of the long-term influence of organic amendment application on topsoil TE concentrations. For the period 2015–2100, the following modelling assumptions were made: (i) stable annual temperature and topsoil drainage, (ii) constant rates of TE input rates via all sources, at the 2014 rates, including the idealised trend inputs via sewage sludge addition and (iii) constant crop yields at the 2014 values. Future SOM, DOC flux and pH trends were generated as previously described.

The projected labile TE concentrations were compared against ecotoxicological critical limits calculated for each metal according to the methodology of Lofts et al. (2004). This method assumes that the free metal ion concentration is the most appropriate indicator of toxicity, combined with a protective effect of the soil porewater pH. The resulting critical limit functions, expressed as labile metal concentrations and functions of soil pH and SOM, aim to be protective of 95 % of soil species. The critical limit functions are reported in the Supplement.

2.8 Statistical analysis

All statistical analyses were carried out in R (version 3.5.0). The Mann–Kendall test (package “Kendall”) was used to assess monotonic trends in the TE time series. Increasing trends (Kendall's τ statistic >0) and decreasing trends (Kendall's τ statistic <0) were considered significant when the two-sided p value was less than 0.05. The bias function and the “dplyr” package were used to calculate the bias and the root mean squared error, respectively, of the simulated labile concentrations versus the observed data.

3 Results and discussion

3.1 TE measurements in soil and organic amendments

Despite the continuous application of organic amendments, total TE concentrations (Fig. 4) showed no significant (P<0.05) accumulation patterns over time in the topsoil according to the Mann–Kendall trend test, except for the Pb concentration in the NIL treatment topsoil, which increased significantly due to atmospheric deposition. While the total concentrations were basically stable, the trends were decreasing for Zn in the SS treatment and for Cu in all the treatments, with Cu displaying large decreases, from 60–102 mg kg⁻¹ in 1972 to 30–57 mg kg⁻¹ in 1995. The concentrations measured in 1972 were clearly elevated when compared to Ballabio et al. (2018), in which an average total Cu concentration of ca. 17 mg kg⁻¹ was reported from more than 21 000 topsoils of EU countries. The observed higher Cu concentration in 1972 could be ascribed to past anthropogenic contamination, such as the application of Cu-based fungicides. The reasons for the extensive loss in total Cu are unknown, as the rate of loss (5–8 mg kg⁻¹ in 35 years) was larger than expected by leaching, which is of the order of 0.2–0.3 mg kg⁻¹ yr⁻¹ (Vulkan et al., 2000). For comparison, the temporal trend of the total P concentration in the topsoil (Fig. 4) exhibited lower temporal variability than the TE concentrations. Phosphorus showed significant accumulation over time in the FYM and SS treatments but not in the NIL and COM treatments.

Figure 4 Total concentration time trends of Zn, Cu, Pb, Cd and P in ZOFE topsoil (20 cm) for the treatments NIL, FYM, SS and COM.

Increases in the EDTA-extractable concentrations (Table 1) were significant for Zn and Cd in the NIL and FYM treatments and for Cu in the NIL treatment (P<0.05). Conversely, in the SS treatment, all the metals showed steady declines in their EDTA-extractable concentrations over the measurement period. Observed metal lability, calculated as the ratio of the EDTA-extractable concentration and total concentration in the same year, is reported in Fig. S8.

Table 1 Soil EDTA-extractable concentrations of Zn, Cu, Pb and Cd in ZOFE topsoil (0–20 cm) for the treatments NIL, FYM, SS and COM.

^* Mean of five replicates. NA: not available.

The total TE concentrations of the organic amendments are reported in Table 2. The farmyard manure and compost samples had comparable levels of total TE concentrations, with the green waste compost presenting higher Pb enrichments. The sewage sludge had higher total TE concentrations than the compost and farmyard manure, explaining why the magnitude of all total TE concentrations ranked in the order NIL < COM = FYM < SS in 2011. The analysed sewage sludge showed also the highest variability in the TE lability (Fig. S8), with Zn lability varying from 0.39 to 0.15 and Cu from 0.48 to 0.22 in the samples from 2008 and 2012. The lower lability of Cu and Pb in the organic amendments than in the topsoil could be ascribed to their stronger affinity for organic matter (McBride et al., 1999).

Table 2 Total and EDTA-extractable concentrations of Zn, Cu, Pb, Cd and P relative to the total dry matter of the organic amendment samples from available years.

^* Lower than or close to the limit of detection.

3.2 Simulations of the metal trends

The IDMM model was run from 1750 (pristine conditions) to simulate the observed metal concentrations data up to 2014. The observed and simulated total concentrations are shown for reference in the Supplement (Fig. S9). Here we focus, firstly, on the predicted labile metal concentrations and briefly discuss the modelling of the total metal concentrations.

From a mechanistic point of view, the IDMM does not directly simulate the total soil metal concentration but explicitly considers separate labile, weakly aged and strongly aged forms, each with their own dynamics. The IDMM sets the initial labile concentration, assuming a steady-state balance of natural input and output metal fluxes, including kinetic transfers to and from the nonlabile pool. In a dynamic simulation, changes in the labile concentration are then driven by changes in inputs and the soil parameters influencing solid-solution partitioning and rate of net ageing (pH, soil organic matter and DOC flux). The initial aged concentrations are set by a combination of equilibrium with the initial labile concentration coupled with the adjustment of one kinetic constant (see the Supplement for details). Changes in the ageing concentration are then driven by ageing transformations of the labile form. Transformations between the labile and weakly aged forms are sufficiently rapid that the weakly aged form responds rapidly to changes in the labile concentration, but the slower transformations involving the strongly aged pool result in a highly delayed response to changes in the labile pool. An example of predicted trends in the three pools is shown in Fig. S10.

Using five strips of the dimension 1 m × 7 m to simulate each plot when considering lateral mixing, the IDMM gave the predictions shown in Fig. 5 when using the Swiss trend and idealised trend inputs, both with and without lateral mixing.

Figure 5 Observed EDTA-extractable concentrations (bullets), simulated labile concentrations with mixing (solid line) and without mixing (dashed line) with Swiss trend inputs (top) and idealised trend inputs (bottom) for Zn, Cu, Pb, Cd in ZOFE topsoil (0–20 cm) for the treatments NIL, FYM, SS and COM.

Fitting of apparent model weathering rates was done for each scenario. Only small differences were seen among the fitted rates, so here we focus only on the rates obtained when assuming lateral mixing and using the idealised trend. Fitted apparent weathering rates were 9.28, 1.62, 0.0798 and 0.914 mg m⁻² per year for Cu, Zn, Cd and Pb, respectively. These are generally higher than the natural atmospheric deposition rates of 0.0699, 0.0752, 0.0109 and 0.522 mg m⁻² per year . Clearly, the contribution of the apparent weathering rate is not negligible relative to other natural inputs. Estimates of topsoil metal weathering inputs are uncommon in the literature, but Imseng et al. (2018, 2019) published estimates of Zn and Cd weathering at Swiss grassland locations. Generally, the weathering rates they estimate (0.001–0.5 mg m⁻² per year for Zn; <0.01–0.39 mg m⁻² per year for Cd) are lower than our fitted Zn rate but comparable for Cd. It is reasonable to assume that our fitted rates are apparent weathering rates only that are subject to a number of uncertainties, particularly the lack of knowledge regarding anthropogenic metal additions to the field prior to setting up the ZOFE experiment. This is particularly clear for copper, for which the apparent rate looks unrealistically high and is likely to be compensating for unquantified past inputs as per what is known about the site history. Nonetheless, the fitting exercise is useful in allowing us to both calibrate model predictions to observations and to emphasise where the key gaps in knowledge exist that should be tackled to allow for more comprehensive and plausible future applications of the IDMM.

The Swiss trend of metal inputs predicted neither the magnitude nor the trends in labile metal concentrations in the SS plot with or without lateral mixing enabled, with the partial exception of Cd (Fig. 5). The downward trend in labile metal concentrations in SS cannot be explained due to metal removal in leaching and/or crop uptake using the key variables (porewater pH and DOC concentration; metal contents of harvestable crops). Simulations were therefore repeated with lateral mixing and optimisation of the inputs to SS (the idealised trend). The optimised inputs (total inputs from 1949 to 2014) were factors of 3.2, 2.3, 1.5 and 6.6 times higher than the literature-derived inputs. With lateral mixing invoked, the model was reasonably successful in describing the downward observed trends in labile metal in the SS plot, although there was a tendency to overestimate post–1990 observed concentrations, particularly for Pb. To illustrate the influence of lateral mixing on predictions, we also ran a simulation using the idealised trend inputs without lateral mixing. The results suggest that lateral mixing due to ploughing has been a key determinant of the observed metal concentrations in SS. In the absence of lateral mixing, higher labile metal concentrations are consistently predicted from 1949 onwards (Fig. 5), and the observations are consistently overestimated. For example, in 2011, the observed labile metal concentrations were 17.6, 18.5, 0.36 and 12.8 mg kg⁻¹ for Cu, Zn, Cd and Pb, respectively. The model predictions, when invoking lateral mixing, were 21, 20, 0.46 and 25.2 mg kg⁻¹, respectively, which are in reasonable agreement, with the exception of Pb. In contrast, the predicted concentrations, when lateral mixing was not invoked, were 40.6, 42.2, 0.87 and 49.7 mg kg⁻¹ respectively, which are all approximately double the concentrations predicted in the presence of lateral mixing. These results are consistent with the lateral distribution of metals in ploughed experimental plots at Rothamsted Experimental Station (McGrath, 1987). This influence of ploughing on metal distribution is, however, likely to be confined to the specific management conditions of experimental plots, such as ZOFE, and not to be relevant for modelling commercial agricultural systems at field and higher scales.

Observed labile metal concentrations in NIL, FYM and COM are consistently lower than those in SS due to the lower rates of input. The IDMM generally reproduces the magnitude and trends in concentrations reasonably well in these plots but tends to slightly underestimate concentrations, particularly in the FYM and COM plots. In the absence of lateral mixing, there is a negligible distinction between the Swiss trend and idealised trend for these plots. Invoking lateral mixing increases the predicted labile concentrations somewhat due to the predicted net transfer of metal from the SS plot to the other plots. This results in an overestimation of the observed labile metal in FYM and COM. The effect of lateral mixing on the predictions for NIL is smaller, since it is not immediately adjacent to SS. There are marginal improvements in the predictions, particularly for Cu and to an extent Zn, when lateral mixing is invoked and the idealised trend SS input trend is applied.

Model performance with soil lateral mixing under both the Swiss trend and the idealised trend is reported in Table 3. The model adequately reproduced the EDTA-extractable concentrations when the input rates were optimised (idealised trend), resulting in a global average RMSE of 2.7 and an average bias (overestimation) of −1.66, −2.18, −4.34 and −0.05 mg kg⁻¹ for Zn, Cu, Pb and Cd, respectively.

Total metal simulations are shown in Fig. S9. Generally, the model predictions are similar in trend to the prediction of labile metal, with increases after 1949 due to ZOFE inputs, and in SS, with a decline from the late 1970s onwards as a result of lateral mixing. Trends in total Zn are well reproduced, although there is a bias towards low predictions in SS and marginally high predictions in COM. For Pb, particularly in FYM, SS and COM, temporal trends in observed concentrations are unclear. For example, there is no clear trend towards decreasing concentrations in SS, as predicted by modelling. Generally, concentrations in FYM, SS and COM are overestimated, while the trend in NIL is well reproduced. Observed cadmium concentrations also show appreciable noise; a clear (declining) trend is only seen in SS. Modelled cadmium concentrations are almost consistently biased to be high, although the trend in SS is reproduced. Inspection of the fit for the NIL plot suggests that the initial total Cd may be overestimated as a result of the fitting approach, as the mean of the first three observations is greater than the subsequent measurements. Total copper exhibits a distinctive declining trend in all the plots, as noted previously, but the reason for these observed trends is unclear. The rates of decline cannot be explained by any of the metal loss processes in the IDMM. The concentrations in SS, while showing the same declining trend as the other plots, stabilise at higher concentrations, reflecting the higher copper input.

3.3 Lateral mixing sensitivity analysis

We investigated the sensitivity of the predictions with lateral mixing to the number of homogeneous soil strips used per plot (Sect. 2.3) by running simulations using the idealised trend and either two or 10 soil strips per plot, maintaining a margin width of 0.2 m. Figure 6 shows the simulated labile concentrations of Zn, Cu, Pb and Cd in 2014 across the transect. With two strips per plot, the labile concentrations in SS were, on average, 28 % lower than with five strips when considering all the metals; the concentrations were comparable in FYM and COM but were higher in NIL by, on average, 81 %. With 10 strips per plot, the labile concentrations in SS were, on average, 40 % higher than with five strips but lower in NIL, FYM and COM. Therefore, increasing the number of strips per plot reduced the predicted redistribution of the TEs from the SS to the adjacent plots. The choice of the number of strips per plot clearly has an effect on both the determination of the idealised trend inputs and the simulation performance. The pragmatic choice of five strips for modelling provides an example of the magnitude of the effect of ploughing on the redistribution of metals across the treatments, but in general, the presence of the lateral mixing effect does somewhat limit the usefulness of the plot data for model evaluation as it requires additional modelling and parameterisation not necessary for true field application.

Table 3 Model performance in terms of the root mean square error (RMSE) and bias of simulated labile concentrations vs. observed EDTA-extractable concentrations with lateral mixing under the Swiss trend and the idealised trend.

Figure 6 Simulated labile concentrations of Zn, Cu, Pb and Cd in 2014 across a transect comprising the plots NIL, FYM, SS and COM in series with five and 10 strips per plot with idealised trend inputs.

3.4 Soil spectroscopy analysis and long-term effects of organic amendment applications

The IDMM model assumes that the TEs present in the organic amendments are fully labile, so they are added to the labile pool when they are introduced in the soil. While metal lability observations in the organic amendment samples (Fig. S8) show that the metals are not entirely labile, they also indicate high variability in the metal lability in the SS amendment samples from 2008 and 2012. Additionally, Fig. S9 shows that Pb lability in the SS-amended soils decreased over time, while, for the other metals in the SS treatment, the lability remained constant despite a significant pH decrease. At least for Pb, a consistent decrease in the incoming metal lability would improve the simulations of the labile metal concentrations in the SS treatment, including under the idealised trend inputs and soil lateral mixing. To support this hypothesis, the soil mineral and organic fractions of soil samples from 1972 and 2011 were investigated to detect eventual changes between treatments and over time caused by the application of the organic amendments.

The XRD diffractograms from 1972 and 2011 (and 2013 for SS only) did not reveal noticeable differences among the soil samples (see Fig. S11 in the Supplement). This would indicate that the long-term application of the organic amendments, including the sewage sludge, did not introduce any exogenous minerals, such as clay minerals and Fe-(oxy) hydroxides, capable of modifying the TE lability in the soil, with particular reference to the decrease in the SS plots.

The FTIR spectra from 1972 and 2011 were also not suggestive of differing organic fraction composition between treatments and over time, except for a varying peak in the sewage sludge samples. To focus on the time change of the soil organic fraction composition, the FTIR spectra from the 2011 and 2013 SS samples were subtracted from the 1972 spectra, thus resulting in differential spectra. As shown in Fig. 7, the differential spectra from 1972–2011 evidenced a peak at 1040 cm⁻¹, associated with the functional group of polysaccharide-like compounds, that was not present in the differential spectra from 1972 to 2013. This means that polysaccharide-like compounds, which are reported to have high affinity for TEs (Geddie and Sutherland, 1993; Veglio et al., 1997), were present in the soil in 1972 and 2013, but not in 2011.

In conclusion, the scarce observations available are more suggestive of organic fraction variability in the sewage sludge amendments applied to ZOFE than to any consistent (in time) metal lability trend. This organic matter composition change could have an impact on the lability of the incoming TEs and, hence, on simulation of labile TE dynamics, but more research is needed on the lability of metal in organic amendments and the response of that lability following addition to soil.

3.5 Long-term effects of organic amendment applications

Figure 8 presents the projected labile metal concentrations from 1940 to 2100 under the model conditions described in Sect. 2.7. In NIL, FYM and SS, labile metal concentrations are projected to be steady or declining, regardless of the pH trend. In COM, steady increases in Cu, Zn and Pb are projected, but there is a small decline in Cd. The projected increases in COM are likely to be due, at least in part, to the smaller decline in soil OM in this plot, relative to the other plots. Soil–porewater partitioning in the IDMM is a function of soil OM, and so a temporal decline in organic matter will lower soil solution partition coefficients and increase metal leaching. The proportional declines in labile metal concentration in NIL, FYM and SS are consistently of the order of Cd > Zn > Pb–Cu. This reflects the higher tendency of Zn and Cd to be lost from the soils, due to their lower binding affinities for the soil OM.

Figure 7 Differential FTIR spectra of the SS samples. SS 2011 spectra subtracted from the 1972 spectra (red line), and SS 2013 spectra subtracted from the 1972 spectra (blue line).

Figure 8 Projected labile concentrations for Zn, Cu, Pb and Cd to 2100 without lateral mixing and with idealised trend inputs. The following two scenarios are considered: constant pH (solid line) and pH decreasing exponentially to 0.5 units below the 2014 value (dotted line).

Projecting a decline in porewater pH to 2100 produces varying results across different metals and plots. In NIL, FYM and SS, Zn, Cu and Cd are consistently projected to be lower in 2100 as a result of lower pH. Conversely, in COM the same metals are projected to increase slightly under the declining pH trend. This is also observed for Pb in SS. In a number of further cases, notably Zn in FYM and SS, Cu in SS and Pb in FYM, the labile metal concentration under declining pH is initially projected to be higher than that under constant pH but then to drop to or below the constant pH projection in 2100. These observations can be rationalised by considering the interplay among the metal input rate, the rate of labile metal loss in leaching and the pH dependence of the ageing rate of added labile metal. In COM, the metal input rates are sufficiently high enough to project continued net metal accumulation given the soil chemistry, while in the other plots input rates are insufficient to allow continued accumulation under the changing soil chemical conditions. Increased leaching losses with declines in pH are expected to be most important for Zn and Cd. For Cu, and to an extent Pb, the distribution of organic matter between solid and solution is expected to drive partitioning and, hence, leaching losses. In all the plots except NIL, the ratio of DOC flux to SOC pool is projected to increase between 2014 and 2100. This trend will have a positive influence on metal leaching by driving an increase in the proportion of metal complexed with DOC in soil porewaters relative to the proportion adsorbed to the soil solids. In contrast to the other factors, declining pH will reduce the rate at which input metal is projected to age and so drive greater retention within the labile pool. This latter process is projected to be particularly important for Pb, but can also be seen to influence projections of the other metals, particularly in COM.

Figure 9 Projected risk characterisation ratios (RCRs) for Zn, Cu, Pb and Cd to 2100 without lateral mixing and with idealised trend inputs. The following two scenarios are considered: constant pH (solid line) and pH decreasing exponentially to 0.5 units below the 2014 value (dotted line).

Figure 9 clearly shows that the largest projected risks occur in SS and are driven by the relatively high metal inputs. In particular, it is notable that, for all metals, exceedance of the critical limit is predicted to have occurred prior to 2014 and to be maintained until 2100. The risk characterisation ratios (RCRs; the ratio of modelled to critical labile concentration) for Zn are projected to exceed 11 around the present day and then to decline under projected constant pH but to remain relatively high (∼ 9) in 2100. A similar pattern is seen for Cd, but the RCRs are smaller, approaching ∼ 1.2 around the present day and declining to ∼ 0.7 in 2100. In contrast, RCRs for Cu and Pb are predicted to continuously increase from 1949 to 2100, although the rate of increase is projected to be relatively low under the assumptions of future conditions.

The declining future pH trend results in higher projected risks for Pb in SS after 2014. The labile Pb concentration is projected to be lower after 2014 in the declining pH scenario (Fig. 8). Therefore, the increased risk must be due to the lowering of the critical limit concentration due to the decline in pH. Similar considerations can explain the slightly higher risk predicted for Zn and Cd, but the interplay among the driving factors appears more complex, since the projected labile metals are not initially lower in the declining pH scenario (Fig. 8). For Cu, the decline in pH has a relatively small influence on the critical limit concentration, and so only marginal differences in risk are projected.

Predicted risks to 2100 in NIL, FYM and COM are small and only limited to Zn; only marginal differences between the constant and declining pH scenarios are seen. Notably, exceedance for Zn is projected to occur in the future (around 2030), but RCR remains relatively small (RCR < 1.5 in 2100).

4 Conclusions

• Using metal chemistry parameters derived from independent data sets, we applied the dynamic soil metal model IDMM to simulate long-term soil metal trends in the Zurich Organic Fertilization Experiment trial, Switzerland.

• Following the calibration of inputs and accounting for lateral soil mixing, the IDMM reproduced the observed EDTA-extractable concentrations of Zn, Cu, Pb and Cd with an average overestimation by −1.66, −2.18, −4.34 and −0.05 mg kg⁻¹ for Zn, Cu, Pb and Cd, respectively.

• Using considerable amounts of data are necessary to run the model. In particular, we used modelling to fit observations and estimate future trends in the soil OM content and to estimate DOC fluxes from the topsoil.

• Estimating historic TE inputs as robustly as possible is important for the modelling. Where estimated inputs give biased predictions of soil TEs, optimisation may be useful within plausible limits if it aids in evaluation of the model.

• Supporting the modelling through mineral weathering inputs of TEs to soils is key; knowledge of metal weathering rates is generally poor, however. Optimisation of weathering rates is possible, but it must be interpreted cautiously. More research is required on the robust determination and modelling of weathering rates in the field.

• Mixing soil laterally due to ploughing appears to be a significant influence on the observed metal concentrations due to the net redistribution of TEs from plots with relatively high concentrations to those with relatively low concentrations under the specific conditions of the ZOFE experiment.

• Changing soil chemistry in the plots since the inception of the ZOFE experiment, notably the acidification and loss of soil organic matter, makes the soils more vulnerable to the ecological impacts of metals.

• Projecting metal concentrations, including under future conditions of constant climate and metal inputs, suggest that historic inputs of sewage sludge would result in present-day exceedances of threshold concentrations (critical limits) for all the TEs and that the exceedances would remain until at least 2100. Some minor exceedance of Cu and Zn critical limits would be expected by 2100 under manure and compost application.

• Applying the IDMM model to the ZOFE trial is a promising step towards understanding the key processes controlling past, current and future TE dynamics at the field scale.

• Acknowledging that, while long-term trials such as ZOFE have the advantage of being well characterised, their soil metal concentrations can be influenced by their distinctive management practice, as this example shows, is important.

• Identifying some key knowledge gaps that need to be addressed for large-scale model application – characterisation of metal weathering rates and DOC fluxes and the metal contents of organic amendments – has been achieved in this study.

• Recognising that, because of the complexity of the model data requirements and likely resulting uncertainties in predictions, large-scale application is likely to be most useful in assessing broad spatial and temporal trends in metal concentrations and risks.