Farmland Trace Metal Contamination and Management Model—Model Development and a Case Study in Central Taiwan

Abstract

In this study, the water quality of the irrigation system and concentration of trace metals in the sediments were combined to establish a farmland trace metal contamination and management model (FTM_CMM). The purpose of this model was to clarify the main sources of the trace metals that have caused the contamination of paddy soil in central Taiwan. The results of the model simulation showed that the trace metals in the paddy soil mainly came from the irrigation water and especially from the sediments in the irrigation channels. The contribution of the sediments in the irrigation channel to the individual trace metals in the paddy soil ranged from 56% to 72% as the contributions for Cr, Cu, Ni, and Zn were 72%, 68%, 56%, and 62%, respectively. The trace metal species according to their concentration in the contaminated soil ranked in the order of Zn > Cr ≈ Cu > Ni, which is about the same as in the channel sediment. During the simulation process, Cr, Cu, Ni, and Zn exceeded the control standards for farmlands in the 18th, 12th, 13th, and 17th years, respectively. This highlights that, in addition to the management of irrigation water quality, the management of trace metal contaminated sediment in adjacent canal irrigation systems is also an important part of the prevention of trace metal contamination in farmland.

1. Introduction

The harmful effects of environmental contamination were often ignored in the past because economic development was always given the first priority. However, once the problems from contamination emerge, there are often considerable social costs that must be paid. In Taiwan, agricultural soils contaminated with trace metals (TMs) are an obvious example. In past decades, due to the negligence of the competent authorities, farmland irrigation and drainage systems in Taiwan used the same channels, which caused TM-containing industrial wastewater to be discharged directly into irrigation channels, thus contaminating farmlands with TMs [1,2]. Due to the low mobility of TMs in a soil environment [3], they are not easily diffused or diluted, so they easily accumulate in surface soil. Soil contaminated with TMs in agricultural areas is likely to affect people’s health and safety once the TMs have translocated into the edible parts of crops [4,5]. TMs also affect the growth and quality of crops [6]. Over past decades, the Taiwanese government has spent billions of dollars on farmland pollution control and remediation. According to data from the Environmental Protection Agency of Taiwan (EPA Taiwan), as of the end of 2018, a total of 1138.5 hectares of farmland sites across Taiwan have been declared as contaminated. A total of approximately USD 60 million has been allocated as expenditures for soil remediation on f-armlands [7].

In order to prevent pollutants from directly entering the human body through the food chain and to safeguard people’s health, EPA Taiwan began performing soil contamination investigations in 1983 [8]. First, the farmland soil on which edible crops were grown was given priority for contamination prevention and control. Then, a systematic and nationwide survey on TM contamination in agricultural lands was launched. The survey included eight TMs and was divided into five stages that spanned from the development of a 1600-hectare grid of a large sampling area to the detailed investigation of small single-hectare sampling areas. Originally, grid sampling was the standard procedure for the soil surveillance performed by environmental protection and agricultural authorities. However, this type of study excludes the connection between the contaminants in paddy soils and those in the irrigation water. Based on the characteristics of the various sources of TMs in paddy soil, farmland contamination is closely related to the water source of the irrigation channel as TMs are transmitted through the channel. Many studies have used a kriging and mapping system to investigate the relationship between anthropogenic activities and the distribution of soil contamination [9,10,11,12]. These results have pointed out that the spatial distribution of TMs in contaminated soils is strongly related to the location of industrial plants and irrigation systems and that the most important contaminants are TMs such as chromium (Cr), copper (Cu), nickel (Ni), and zinc (Zn).

For paddy soils where the TM content exceeds the control standards announced by EPA Taiwan, statistical analysis with paired t-tests, the correlation between the TM content in the topsoil and that in the corresponding aqueduct sediments and irrigation water has been determined to be 86% and 57%, respectively [13]. Furthermore, it can be inferred that the factories in nearby areas have made a significant contribution to the occurrence of TM contamination [14]. Industrial wastewater discharged into irrigation channels and then used as irrigation water by farmers has resulted in the accumulation of TMs in farmland soils. In addition, many studies in the past have shown that farmland contamination is closely related to irrigation and drainage systems, especially when irrigation water is contaminated by the wastewater from factories [1,15].

However, there have been few past studies that have discussed the direct and cumulative effects of channel water quality and sediment TM concentrations on the TMs in paddy soil. In this study, the water quality of the irrigation system and the concentrations of TMs in the sediment were combined to establish a farmland TM contamination and management model (FTM_CMM), which may be used to understand the main sources of TMs in paddy soil and to provide direction to the development of current contamination control strategies and farmland conservation programs. As a result, this proposed screening method can effectively target a contaminated area and provide government agencies the ability to enforce farmland management, investigation, and restoration effectively.

2. Materials and Methods

2.1. Establishment of Background Information

The target paddy soils that were studied are located in central Taiwan and have been heavily contaminated with TMs, including Cr, Cu, Ni, and Zn. The irrigation water is supplied by the East–West Second Canal irrigation system, one of the major irrigation systems in Changhua County. The process of this study is shown in Figure 1 and included the following main steps:

• The collection of the related background information.
• The establishment of the background information.
• The model establishment and application of the GIS (Geographic Information System).
• Results and discussion.
• Review of the prevention, management, and remediation of TM-contaminated farmland.
• Conclusions and suggestions.

2.2. Model Development

The FTM_CMM uses the concept of mass balance to calculate the input and output of TMs in paddy soils. The accumulation of TMs in soil mainly comes from external sources, such as irrigation water, irrigation channel sediment, fertilizers, and atmospheric deposition (Figure 2). When irrigation water is used in a paddy field, irrigation water and sediments containing TMs are also introduced, which accumulate in the paddy soil due to soil physical or chemical adsorption and sedimentation. In addition, the use of TM-containing fertilizers and atmospheric deposition can increase the TM content of soil. In contrast, a decrease in the TM concentration in soil can occur through runoff, percolation, and plant uptake processes. The total mass of soil TMs can be written as a one-dimensional conservation equation:

M_ac,i = ΣM_in,i − ΣM_out,i

(1)

Or expressed as:

C_t,i (t) = C_0,i + C_ac,i (t)t

(2)

Among these values, M_in,i refers to the different paths for TMs of type i to be input into the soil, M_out,i refers to the different paths for TMs of type i to be released from the soil, and M_ac,i is TMs of type i that are being accumulated in the soil. C_0,i refers to the background total metal concentration of i type TMs in the soil, C_t,i (t) refers to total TM concentration of type i at the time of t, and C_ac,i (t) refers to the accumulated TM concentration rate of type i in the soil at the period time of t. Assuming that the soil properties have not changed during the simulation, the change in the mass of TMs per unit volume or weight is equivalent to the change in concentration.

2.2.1. Mass Balance Equation

Since the topsoil (plough layer) of a paddy field is regularly plowed twice a year, the TMs at a fixed depth (usually 20 cm in Taiwan) are completely mixed and the concentration is not affected by vertical depth or horizontal position. Using a homogeneous soil layer of a fixed depth (z) and bulk density (ρ) can simplify the calculation of the model. The mass balance for TMs in the homogeneous mixing soil can be expressed by the following equation [16,17,18]:

$$\mathsf{\rho }\mathrm{z}\frac{{\mathrm{dC}}_{\mathrm{t},\mathrm{i}}}{\mathrm{dt}}={\mathrm{I}}_{\mathrm{i}}-{\mathrm{U}}_{\mathrm{i}}-{\mathrm{L}}_{\mathrm{i}}-{\mathrm{R}}_{\mathrm{i}}$$

(3)

The total type i TM content of the plough layer system (topsoil), C_t,i (mg/kg), depends on the input rate at the soil surface, I_i (mg/m²/year); output rate consisting of the removal rate in the harvested products, U_i (mg/m²/year); the leaching rate in the soil lower layer, L_i (mg/m²/year); and the runoff rate at the soil surface, R_i (mg/m²/year).

2.2.2. External Inputs

The input sources of TMs (I_i) in agricultural soil may include the following: the introduction of irrigation water, I_I,i (mg/m²/year); fertilizer use, I_F,i (mg/m²/year); atmospheric deposition, I_A,i (mg/m²/year); and alluvial sediment carried from irrigation channels to paddy soil, I_S,i (mg/m²/year), which can be represented as follows:

I_i = I_I,i + I_F,i + I_A,i + I_S,i +…

(4)

2.2.3. Estimation of TMs from Irrigation Water

The sources of TMs in canal irrigation systems include the transport of irrigation water and sediment in each square meter. The annual runoff volume, Q_R (m³/m²/year), can be calculated from the field water balance expressed as

Q_R = Q_I + Q_P − ET − Q_L

(5)

where Q_I (m³/m²/year) is irrigation, Q_P (m³/m²/year) is precipitation, ET (m³/m²/year) is evapotranspiration, and Q_L (m³/m²/year) is percolation within each unit area per year. In addition, the different mass transportation paths can be expressed as the following equations, respectively:

I_I,i = Q_I × C_WI,i

(6)

R_i = Q_R × C_WR,i

(7)

where C_WI,i (mg/m³) is the TM concentration in irrigation water for type i and C_WR,i (mg/m³) is TM concentration in drainage (runoff) water for type i.

2.2.4. TM Adsorption

In soil, TMs may be combined with soil particles, the organic matter present in soil, or dissolved in the soil solution. The TMs in the water–soil–plant systems depend on their concentrations in the dissolved form. The total TM concentration can be expressed as

C_t,i = ρ K_f,i C_WI,iⁿ + θC_WI,I

(8)

where K_f,i denotes the partition coefficients (mg⁽¹⁻ⁿ⁾/kg·(m³)ⁿ) and n denotes the sorption parameter (unitless). θ denotes the moisture content of the soil (m³/m³). The former term of the above equation indicates that TMs are adsorbed in the soil, and the manner of this adsorption can be expressed by the Freundlich isothermal adsorption equation [16]. The latter term represents the phase of the TMs in the soil solution. The model assumes that under homogeneous conditions, the total concentration of trace metals in the plough layer includes solid phase and solution phase. In paddy soil, linear adsorption isotherms (n = 1) are observed for TMs [17].

2.2.5. Bottom Sediments

When farmland is in the process of diversion irrigation, the irrigation water may washout bottom sediment that contains TMs, I_S,i (mg/m²/year), in the channel into the farmland. The TMs in the sediment will also enter the paddy soil, and b (unitless) is the specific thickness ratio or the weight ratio of the bottom sediment washout in the paddy soil per year. C_S,i (mg/kg) is the content of TMs in the bottom sediment. It is assumed that the mass and concentration of the sediment will not change during the process of transporting to the paddy field, and the boundary conditions of the plough layer will not be affected. This can be expressed as

I_S,i = b × C_S,i

(9)

2.2.6. Atmospheric Deposition and Fertilizers

Wet and dry atmospheric deposition may also contain TMs. Its flux can be expressed as I_A,i (mg/m²/year). It is related to the amount of atmospheric deposition, D_A (kg/m²/year) and the content of TMs in the deposition (or dust fall), C_A,i (mg/kg). D_A refers to the average weight of the deposition in each square meter area of paddy soil per year. In order to simplify the model, atmospheric deposition is represented by dust fall (the particle size is about 10–100 μm). When dust fall enters the topsoil, it is regarded as completely mixed, and its properties are the same as that of the topsoil and do not consider the proportion of the dissolved constituents in the falling dust.

I_A,i = D_A × C_A,i

(10)

Fertilizers are often applied during cultivation, and TMs in fertilizers may also accumulate in the soil. Considering that the cumulative flux of fertilizer use, I_F,i (mg/m²/year), is related to annual fertilization, W_F (kg/m²/year), and different TM contents, C_F,i (mg/kg), then this can be expressed as

I_F,i = ΣW_F × C_F,i

(11)

where ΣW_F means the total average weight of different fertilization applied in each square meter area of soil in each year. Based on individual fertilizer usage and TM content, the total amount of TM accumulation in the soil can be obtained. Because different chemical fertilizers are mainly used in the study area, it is assumed that all the fertilizers are 100% soluble.

2.2.7. Leaching

Leaching, L_i (mg/m²/year), is one TM output pathway from soil. It always relates to the rate of percolation, v (m/year); the thickness of the soil layer, z (m); and the TM concentration in the soil, C_t,i (mg/kg). When the TMs have a low concentration, the concentration gradient in the different soil layers and depths can be ignored. The fixed depth of the soil may be assumed to be uniform throughout the plough layer. The leaching process can be described simply by the convective–dispersive transport equation [17]:

L_i = θ v C_t,i/K_d,i

(12)

where θ denotes the moisture content of soil (m³/m³) and v denotes the rate of percolation (m/year). Linear isotherm adsorption (n = 1) is used in this model, and k_d,i are partition coefficients.

2.2.8. Plant Uptake

Plant uptake, U_i (mg/m²/year), is the amount of TMs removed from the soil by plants that are harvested. Plant uptake depends on the total soil TM concentration and can be expressed as follows [19]:

$${\mathrm{U}}_{\mathrm{i}}={\mathrm{k}}_{\mathrm{c},\mathrm{i}}{\mathrm{C}}_{\mathrm{t},\mathrm{i}}^{\mathrm{n}}$$

(13)

where k_c,i is a coefficient relating the uptake with the yield of the plant (Y_i, kg/m²/year) and the bioconcentration factor (BCF_i), which is the relationship between the TM concentration in the harvested plant (C_p,i) and the total soil TM concentration (C_t,i). C_p,i represents the TM concentration in the harvested plant parts (mg/kg in dry matter based) [20]. In order to simplify the calculation, it is assumed that the BCF value of each TM will not change with a change in concentration in the soil, that is, the ratio of C_p/C_t is assumed to be constant during the simulation.

$${\mathrm{k}}_{\mathrm{c},\mathrm{i}}={\mathrm{Y}}_{\mathrm{i}}\times {\mathrm{BCF}}_{\mathrm{i}}$$

(14)

BCF_i = C_p,i/C_t,i

(15)

$${\mathrm{U}}_{\mathrm{i}}={\mathrm{Y}}_{\mathrm{i}}\times \frac{{\mathrm{C}}_{\mathrm{p},\mathrm{i}}}{{\mathrm{C}}_{\mathrm{t},\mathrm{i}}}{\mathrm{C}}_{\mathrm{t},\mathrm{i}}^{\mathrm{n}}$$

(16)

If it is a linear relationship, then n = 1

U_i = Y_i × C_p,i = Y_i × BCF_i × C_t,i

(17)

2.2.9. Governing Equations

Finally, if the thickness of the topsoil is fixed in the homogeneous mixing soil, the external input equation can be obtained as follows:

I_i = (Q_I × C_WI,i + D_A × C_A,i + ΣW_F × C_F,i) − Q_R × C_WR,i

(18)

Combining the above equations, the governing equation can be expressed as

ρz dC_t,i/dt = (Q_I × C_WI,i + D_A × C_A,i + ΣW_F × C_F,i − Q_R × C_WR,i) − Y_i × BCF_i × C_t,i – θ × v/K_d,i × C_t,i

(19)

or simplified as

$$\mathsf{\rho }\mathrm{z}\frac{{\mathrm{dC}}_{\mathrm{t},\mathrm{i}}}{\mathrm{dt}}={\mathrm{I}}_{\mathrm{i}}-{\mathrm{Y}}_{\mathrm{i}}\times {\mathrm{BCF}}_{\mathrm{i}}\times {\mathrm{C}}_{\mathrm{t},\mathrm{i}}-\mathsf{\theta }\times \mathsf{\nu }\times \frac{{\mathrm{C}}_{\mathrm{t},\mathrm{i}}}{{\mathrm{K}}_{\mathrm{d},\mathrm{i}}}$$

(20)

The analytical solution of the above equation is as follows:

C_t,i = (Q_I × C_WI,i + D_A ×C_A,i + ΣW_F × C_F,i − Q_R × C_WR,i)/ρz × (1 − exp(−Y_i × BCF_i − θ v/K_d,i)/ρz × t))/((−Y_i × BCF_i − θ v/K_d,i)/ρz) + exp(−Y_i × BCF_i − θ v/K_d,i)/ρz × t) × C_0,i

(21)

where C_0,i is for the initial TM concentration in soil for each TM. In order to simplify the description of the equation, the following formulas are expressed as

I_i = (Q_I × C_WI,i + D_A × C_A,i + ΣW_F × C_F,i − Q_R × C_WR,i)/ρz

(22)

k_u,i = Y_i × BCF_i/ρz

(23)

k_l,i = θ × v/K_d,i/ρz

(24)

where k_u,i is equal to the plant uptake rate constant (year) and k_l,i is equal to the leaching rate constant (per year). The governing equation can be simplified as

dC_t,i/dt = I_i − k_u,i × C_t,i − k_l,i × C_t,i

(25)

The analytical solution of the equation can simplify the description as follows:

$${\mathrm{C}}_{\mathrm{t},\mathrm{i}}=\frac{{\mathrm{I}}_{\mathrm{i}}\left(1-\exp \left(-{\mathrm{k}}_{\mathrm{u},\mathrm{i}}\mathrm{t}-{\mathrm{k}}_{\mathrm{l},\mathrm{i}}\mathrm{t}\right)\right)}{{\mathrm{k}}_{\mathrm{u},\mathrm{i}}+{\mathrm{k}}_{\mathrm{l},\mathrm{i}}}+\exp \left(-{\mathrm{k}}_{\mathrm{u},\mathrm{i}}\mathrm{t}-{\mathrm{k}}_{\mathrm{l},\mathrm{i}}\mathrm{t}\right){\mathrm{C}}_{0,\mathrm{i}}$$

(26)

2.2.10. Combined Bottom Sediments

The sediment in an irrigation channel may enter a paddy field due to water diversion, flood erosion, or other factors. When a paddy field is in the process of diversion irrigation, the irrigation water may washout bottom sediment containing TMs, I_S,i (mg/m²/year), in the channel into the paddy soil, and b (unitless) is the specific thickness or the weight ratio of the bottom sediment washout in the soil per year. C_S,i (mg/kg) is the content of TMs in the bottom sediment. The flux in TMs via sediment transportation can be expressed as I_S,i (mg/m²/year):

I_S,i = b × C_S,i

(27)

If the transportation of sediment is included in the mass balance equation and the topsoil is of a fixed depth, then the governing equation is as follows:

C_T,i = C_t,i × (1 − b) + b × I_S,i = C_t,i × (1 − b) + b × C_S,i

(28)

$${\mathrm{C}}_{\mathrm{T},\mathrm{i}}=\left(\frac{{\mathrm{I}}_{\mathrm{i}}\left(1-\exp \left(-{\mathrm{k}}_{\mathrm{u},\mathrm{i}}\mathrm{t}-{\mathrm{k}}_{\mathrm{l},\mathrm{i}}\mathrm{t}\right)\right)}{{\mathrm{k}}_{\mathrm{u},\mathrm{i}}+{\mathrm{k}}_{\mathrm{l},\mathrm{i}}}+\exp \left(-{\mathrm{k}}_{\mathrm{u},\mathrm{i}}\mathrm{t}-{\mathrm{k}}_{\mathrm{l},\mathrm{i}}\mathrm{t}\right){\mathrm{C}}_{0,\mathrm{i}}\right)\times \left(1-\mathrm{b}\right)+\mathrm{b}\times {\mathrm{C}}_{\mathrm{S},\mathrm{i}}$$

(29)

where b is the weight ratio of the sediment alluvial to the total weight of the soil in a unit square meter per year. C_T,i is the target year soil total trace metal concentration after combining sediment.

2.2.11. Case Study

In order to identify the distribution of TMs among the farmlands irrigated by the East–West Second and Third irrigation canals with high contamination potential within Changhua County in central Taiwan (Figure 3), EPA Taiwan carried out a study entitled “The Investigation Project for National Agricultural Land with High Potential of Heavy Metal Contamination (IV)” from 2016 to 2017. After collecting and analyzing the relevant documents on the environmental factors (such as pH, water quality), related historical survey records, irrigation system information, surrounding geological characteristics, and spatial–temporal distribution data, it was concluded that the farmland’s TM contamination was closely related to the irrigation system. In general, the irrigation system has been proved a major path for contaminant migration and a cause of agricultural soil contamination [21]. Obviously, poor management of the irrigation system leads to an indiscriminate mix of irrigation water and discharged wastewater from the manufacturing industries, which is the main cause of farmland TMs contamination. In order to understand and confirm the real cause of pollution, this study used mathematical methods to verify whether the results of past investigations are correct. The relevant parameters and data analysis of the model in case study area were as follows: The data of the geographic information system was established from the National Land Surveying and Mapping Center, Ministry of the Interior, R.O.C. The target year of the simulation was set at 40 years, mainly because Taiwan’s economy began to take off rapidly in the 1970s and 1980s. There are a large number of metal plating plants in this study area, and TM contaminants are discharged into irrigation channels along with the discharges. The study area is the main rice production area in the western plain of Taiwan. In the past few decades, chemical fertilizers were the main fertilizer used. In the simulation process, the cumulative amount of TMs was calculated based on the annual average fertilizer usage and converted into the model calculation based on the TM content in the fertilizer. The content of TMs in fertilizers refers to the related research results [22]. The air deposition (dust fall) annual average value was taken from the study area from 2006 to 2015 [23,24], and according to its TM content [25], it was converted into a cumulative milligram per square meter per year for the input data of this model. The consumptive characteristic of irrigation water for paddy fields was evaluated from Chen et al. [26]. The mean value (including the upstream, middle, and downstream channel) of irrigation water quality (C_WI,i) and sediment metal concentration (C_S,i) were taken from the more complete survey results in 2006 [27] to represent the mean value of the past decades. The water budget of each component paddy field was taken from the study results of Fan et al. [28] during 2004–2009, including irrigation water (Q_I), precipitation (Q_P), evapotranspiration (ET), runoff (Q_R), and vertical percolation (Q_L) [28,29]. Considering that the paddy field irrigation rate is fast during irrigation period, in order to simplify the calculation process, this study regards the water quality of irrigation water and runoff water as the same. As there is no relevant research report on paddy field partition coefficients (k_d,i) in Taiwan, the relevant parameters of this study refer to the recommended values of the U.S. EPA [30]. There are two harvests a year in the study area, and the annual yield of rice grain (Y_g) is about 1.1 kg/m²/year [31]. The accumulation ability of rice for different trace metals was significantly different [32]. The mean BCF (bioconcentration factor) values of Cu, Ni, and Zn in rice were taken from the study results of Kunhikrishnan et al. [33]. Since there were no BCF relevant data on Cr, this study referred to the study results of Chen in 2000 [2] and used the half value of Ni as the calculated value for the input data of this model. All relevant data and parameters are summarized in

Table 1. Default values of model parameters and definition for simulation in paddy soil.

Model Parameters	Definition	Assigned Values	Units	Sources
t	Simulation period	0–40	year	Case study
C0,i	The background mass concentration of i type TM in the soil	Cr: 97.2	mg/kg	Case study [10]
		Cu: 85.9
		Ni: 108.4
		Zn: 215.5
Ct,i	The total soil TM concentration in year t (without combined sediments)		mg/kg	Case study
CT,i	The total soil TM concentration in year t (combined sediments)		mg/kg	Case study
z	Soil depth	0.2	m	Case study
ρ	Soil bulk density	1300	kg/m3	Case study
ΣWF	The total average weight of different fertilization applied in each year	0.2	kg/m2/year	[22]
CF,i	Concentration of TMs in fertilizer (mean)	Cr: 14.6	mg/kg	[22]
		Cu: 3.4
		Ni: 10.2
		Zn: 56.1
DA	Air deposition (dust fall)	0.033	kg/m2/year	[23,24]
CA,i	Concentration of TMs in air deposition (mean)	Cr: 141	mg/kg	[25]
		Cu: 908
		Ni: 262
		Zn: 5360
QI	Irrigation water	3.384	m3/m2/year	[26]
CWI,i	Irrigation water quality (mean)	Cr: 38	mg/m3	[27]
		Cu: 165
		Ni: 172
		Zn: 676
QR	Runoff water	1.928	m3/m2/year	[28]
QP	Precipitation	1.851	m3/m2/year	[28]
ET	Evapotranspiration	2.931	m3/m2/year	[28]
CWR,i	Runoff water quality	Cr: 38	mg/m3	Case study
		Cu: 165
		Ni: 172
		Zn: 676
QL	Percolation water	0.376	m3/m2/year	[29]
v	Percolation rate	0.376	m/year	[29]
θ	Moisture content of soil	0.27	m3/m3	Case study
CS,i	Sediment TM concentration (mean)	Cr: 756.8	mg/kg	[27]
		Cu: 757.0
		Ni: 594.8
		Zn: 1802.8
Kd,i	Partition coefficients	Cr: 7.9	kg/m3	[30]
		Cu: 0.5
		Ni: 1.3
		Zn: 1.3
Yi	Uptake with the yield of the plant		mg/m2/year	Case study
Yg	Yield of rice grain	1.1	kg/m2/year	[31]
BCFi	Bioconcentration factor: TM concentration ratio in the harvested plant part (based on dry matter) to soil.	Cr: 0.017	unitless	[32] [2] [33]
		Cu: 0.399
		Ni: 0.033
		Zn: 0.473
b	Sediment ratio, the ratio of sediment alluvial weight to the total weight of soil in the unit square meter per year	0.005–0.02	unitless	Case study
		(0.0141)
Ii	The input sources of TMs		mg/m2/year	Case study
II,i	The input sources of TMs by irrigation water		mg/m2/year	Case study
IF,i	The input sources of TMs by fertilizers		mg/m2/year	Case study
IA,i	The input sources of TMs by air deposition		mg/m2/year	Case study
IS,i	The input sources of TMs by sedimentation		mg/m2/year	Case study
Ui	Plant uptake (rice) is the amount of TMs removed from the soil by plants		mg/m2/year	Case study
Li	TM output pathway from soil by leaching		mg/m2/year	Case study
Ri	TM output pathway from soil by runoff		mg/m2/year	Case study

. Those parameters were converted into the annual mean value or constant input to the model for simulation, while ignoring the influence of other external environmental factors (such as temperature, pH, soil type, etc.). Soil sampling and analysis were completed by the EPA project [21]. First, 1300 sets of soil samples were screened by portable XRF (X-ray fluorescence spectrometer) in the studied irrigation area using the grid method, and then samples with XRF screening results exceeding the regulatory standards were selected, and the aqua regia digestion method was used for the full analysis of TMs. Among them, there were 119 sampling sites in area (a) and 174 sampling sites in area (b). The detection limits of the detector (710 ICP-OE) for Cr, Cu, Ni, and Zn were 1.55, 1.38, 1.24, and 1.52 mg/kg, respectively. The main contaminants in the agricultural lands were Cr, Cu, Ni, and Zn with an average aqua regia soluble concentration between 296 and 1056 mg/kg (Figure 4). The average concentrations from the 119 sampling sites for Cr, Cu, Ni, and Zn in study area (a) were 410, 464, 348, and 1024 mg/kg, respectively, and from the 174 sampling sites in study area (b) they were 390, 296, 430, and 1056 mg/kg, respectively. Since 1991, the remediation of farmlands contaminated with TMs in Changhua County has been performed continuously through many remediation projects funded by central and local governments. The above two groups of paddy field soil contamination data belong to two different projects, and the remediation work began in 2019.

2.2.12. Model Parameters

The required parameters include soil specific gravity, irrigation water quality, irrigation water volume, field water consumption, sediment TM concentration, fertilizer consumption and characteristics, air deposition and TM content, crop yield and TM uptake rate, and other factors. The related parameters were assigned as shown in Table 1.

2.2.13. Model Calibration and Verification

Surveying data from 2016 to 2017 [21] were used for model calibration and verification. The data included irrigation channel water quality, sediment TM concentrations, paddy soil TM concentrations, and other data. The model calibration and verification used mean absolute percentage error (MAPE) as a measure of quality for this model [34]:

$$\mathrm{MAPE}=\frac{1}{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left|\frac{{\mathrm{A}}_{\mathrm{i}}-{\mathrm{F}}_{\mathrm{i}}}{{\mathrm{A}}_{\mathrm{i}}}\right|$$

(30)

where A_i is the actual value and F_i is the forecast value.

2.2.14. Model Limits

This model does not consider the horizontal transport and adsorption of TMs in a paddy field, and it is limited to assessing the accumulation of TMs in the soil at a limited depth (in the upper arable soil horizon, plough layer, or topsoil, 0–20 cm) per unit area. In order to simplify the calculation of the model, some parameters were taken from long-term survey or observational data, such as atmospheric deposition, fertilizer use, harvesting, irrigation water quality, evaporation, infiltration, and other data. The model ignored the synergistic effects in the process of leaching or uptake of metals and also ignored the impact factors of changes in oxidation–reduction potential, microbial activity, soil/solution pH, etc., on the efficiency of soil absorption/desorption of trace metals.

2.3. Software

In this study, Microsoft Excel 2019 and ArcGIS Desktop 10 were used as data collation and plotting software, and Microsoft Excel 2019 was used as data input and numerical calculation software.

3. Results

3.1. Model Calibration and Verification

In Taiwan, a series of surveys on TM-contaminated soils and background concentrations in agricultural land began being conducted in 1982 [2]. The background concentrations of TMs in the paddy soils of the study area of the current study are shown in

Table 2. The background concentrations of TMs (mg/kg) in the agricultural soils of this study and the results of the model calibration and verification.

Item	Cr	Cu	Ni	Zn	Remarks
C0,i	97.2	85.9	108.4	215.5	Background concentration of each TM in the topsoil [10]
CT,i	391.7	392.6	352.9	1000.4	Prediction for the target year of each TM concentration
Cm,i	409.9 (119)	463.5 (119)	348.1 (119)	1023.8 (119)	Measured TM mean concentrations in the contaminated soil of two different areas [21]. (Numbers in brackets are expressed as sample size)
	390.3 (174)	430.3 (174)	295.8 (174)	1056.2 (174)
Calibration	4%	15%	1%	2%	Model evaluation criterion, MAPE:
Verification	0%	9%	19%	5%	<10% high prediction
					10~20% good prediction
					20~50% reasonable prediction
					>50% inaccurate prediction

and are based on the survey results over years [10]. The background average concentrations of Cr, Cu, Ni, and Zn were 97.2 mg/kg, 85.9 mg/g, 108.4 mg/kg, and 215.5 mg/kg, respectively, in the study area. In order to evaluate the credibility and reliability of the model simulation results, we use the model evaluation criterion MAPE for testing the model accuracy of this study. The results of the calibration and verification (Table 2) show that Cu and Ni are in the good and high prediction range, and the other TMs meet the high prediction range. In addition, according to the results of calibration and verification, the optimal ratio for sediment alluvial weight is about 1.4–1.5% (Figure 5). The average concentration of TMs in the channel sediments in the study area for Cr, Cu, Ni, and Zn was 756.8, 757.0, 594.8, and 1802.8 mg/kg, respectively [27].

3.2. Model Simulation Analysis

In the results of this study, as shown in Figure 6, Figure 7 and Figure 8, it was found that the TMs in the paddy soil have mainly come from the sediment in the irrigation channel. For the paddy soil in the study area, the sediment in the channel acts as a TM pool and is one of the main sources for many of the TMs in the contaminated paddy soil. The TM species according to their concentration in the contaminated soil rank in the order of Zn > Cr ≈ Cu > Ni, which is about the same as in the channel sediment.

With the simulation process, the first TM that exceeded the control standard for farmland (200 mg/kg) was Cu, which occurred at about the 12th year, and in the 40th year, it was almost close to the control standard for normal lands (400 mg/kg) (

Table 3. The control and monitoring standards for the concentrations of TMs in the soil of normal lands and farmlands as announced by the Environmental Protection Agency (EPA) of Taiwan.

Item		As	Cd	Cr	Cu	Hg	Ni	Pb	Zn
		mg/kg
Normal land	Control standard a	60	20	250	400	20	200	2000	2000
	Monitoring standard b	30	10	175	220	10	130	1000	1000
Farmland	Control standard a	60	5	250	200	5	200	500	600
	Monitoring standard b	30	2.5	175	120	2	130	300	260

^a Soil pollution control standard refers to the soil pollution control limits set to prevent the deterioration of soil pollution. ^b Soil pollution monitoring standard refers to the concentration of pollutants that are required to be monitored for soil pollution based on soil pollution prevention purposes.

). The next was Ni, which surpassed the control standard for normal lands at 200 mg/kg in the 13th year. The accumulation rate of Cr was almost the same as that of Cu. Owing to the higher standard value of regulation, it exceeded the control standard for normal lands with 250 mg/kg in the 18th year. Zn exceeded the control standard for farmlands (600 mg/kg) at about the 17th year, and the simulated value was about 1000 mg/kg in the 40th year. Since Zn has a higher standard value for normal land (2000 mg/kg), its concentration by the end of the simulation period was only about half of the control standard.

For the monitoring standard of farmlands, the concentrations of Zn, Ni, Cu, and Cr exceeded the monitoring standard for farmland in the second, third, fourth, and ninth year, respectively. Except for Cr, the concentrations of the other three TMs exceeded the monitoring standards for farmlands within five years. This stimulation result reveals that the competent authorities should regularly implement soil sampling and monitoring projects to prevent the contamination of paddy soil.

3.3. Irrigation Water without Sediment Transportation

According to the results of previous surveys on irrigation channels, the long-term average concentrations of TMs in irrigation water can meet the standards for irrigation water quality (

Table 4. Irrigation water quality during the simulation period in this study area.

Item	Cd	Cr	Cu	Hg	Ni	Pb	Zn
	mg/L
Mean	0.003	0.038	0.165	0.000	0.172	0.053	0.676
SD	0.002	0.032	0.215	0.000	0.099	0.081	1.259
Maximum	0.005	0.102	0.653	0.000	0.324	0.147	3.780
Minimum	0.000	0.004	0.016	0.000	0.048	0.000	0.094
Standard a	0.01	0.1	0.2	0.002	0.2	0.1	2

^a Standards for irrigation water quality as announced by the Council of Agriculture, Taiwan.

). The results of the model simulation also show that with long-term farming (Figure 7), concentrations of TMs are still below monitoring standards, and there is no problem with the sustainable use of the farmland. However, once the TMs in the irrigation water exceed irrigation water standards for decades, the paddy soil may still potentially become contaminated and exceed control standards, such as with Zn. A related study has also been found in India [35].

3.4. Irrigation Channel Sediment

When TMs are released into a water environment, the dilution effect of the water and its chemical composition can change the total concentration of TMs and the relative content of their different chemical species. They can also interact with the suspended solids and sediment in the water. As a result, the distribution of TMs in the water phase and solid phase changes and reactions occur, including diffusion, precipitation, adsorption, sedimentation, and more, and TMs accumulate in the irrigation channel sediment.

The TMs in a water environment can be roughly divided into solid matter and dissolved matter. The solid matter in water, including suspended solids and sediment, plays an important role in the transmission and distribution of TMs in a water environment. Sediment may enter a paddy field and contaminate paddy soils through diversion irrigation, flooding, or farmers dredging the sediment. In Figure 8, when comparing the situation with and without sediment, the difference in the final concentration of TM can be found, and the amount of change in concentration is the degree of influence of the sediment. The simulation results of this study have shown that 56–72% of the TMs in the contaminated paddy soil were a result of the sedimentation process in the sediment of the irrigation channel. The individual concentrations of TMs in the contaminated paddy soil resulting from the sediment in the irrigation channel for Cr, Cu, Ni, and Zn were 72%, 68%, 56%, and 62%, respectively, as shown in Figure 8. This also shows that, consistent with the survey results, most of the TMs accumulate at the inlet of the irrigated paddy field [13].

3.5. Plant Uptake

To protect the air quality in Taiwan, farmers are prohibited from burning rice straw and weeds in the field. Most of the rice straw is chopped by harvesters and then directly decomposed into fertilizers by plowing and burying in the paddy field. Therefore, the plant uptake removal of TMs in the paddy field mainly comes from the rice grains at harvest. In general, the TMs that are more likely to accumulate in rice are relatively more likely to be removed [33,36]. Related studies have pointed out that Zn and Cu are trace essential elements required for plant growth and are more easily absorbed by rice. Ni and Cr are less likely to be absorbed by plants and are more stable in soil [37,38]. The simulation results show that through rice harvesting, the removal rate (percent per year) of TMs in the paddy soil follows Zn > Cu > Ni > Cr, with rates of 4%, 3%, 0.3%, and 0.1%, respectively, as shown in Figure 8.

3.6. Model Sensitivity Analysis

To understand the influence of each main parameter on the analysis results of this model, this study conducted a sensitivity analysis on several main parameters and sources. When each parameter was increased by 50%, the significance of the impact on the estimated value of the target year was evaluated. Based on the default data for the study area, the sensitive factors of this model were sediment concentration (C_S), bioconcentration factor (BCF), irrigation water quality (C_WI), atmospheric deposition concentration (C_A), percolation velocity (ν), and fertilization concentration (C_F). From the analysis results (Figure 9), it can be seen that the influence of sediment concentration (C_S) is the highest, and the range of influence is from 36% to 42%. The variety of different TMs mainly comes from the influence of the concentration ratios of the TMs in the sediment and irrigation water and the higher the concentration ratio, the greater the influence. The effect of BCF is between −0.1% and −2.0%, which ranges for Cr, Ni, Cu, and Zn in this order. For Zn and Cu, which are essential elements, the cumulative concentrations in plants are higher than those of Cr and Ni. The removal rates of Zn and Cu for the contaminated soil are also higher. Another influencing factor is the concentration of TMs in the contaminated soil: the higher the concentration of the TMs in the soil, the higher the amount absorbed by the plants. Relatively, Zn and Cu are also TMs with higher concentrations in the contaminated paddy soil. The order of influence from the irrigation water is Zn, Ni, Cu, and Cr, which is slightly different from the order of the sediment. The main reason is that the concentration of Ni in the irrigation channel is higher than that of Cu and the concentration of Ni more easily becomes higher in the paddy soil because of its relatively higher concentration in the environment. Ni also forms oxide precipitates less easily than Cu, so its mobility is higher than that of Cu. In this study, since all of the TM levels in the irrigation water were lower than the irrigation water quality standards, the impact of irrigation water is not so significant for the contamination of paddy soil, and its impact range is between 1.0% and 6.7%.

3.7. Review on the Prevention, Management, and Remediation of TM-Contaminated Farmlands

In Taiwan, the current prevention and control strategy for TM-contaminated farmlands focuses on remediation. Through the use of data ascertained by well-designed investigations of contaminated farmland, government authorities are able to announce and conduct a variety of plans for the control and evaluation of farmland remediation. However, potential sources for the contaminants of farmland have been neglected. Therefore, even after the remediation of TM-contaminated farmland is complete, there still exists the possibility for recontamination, thus losing the original remediation benefits.

The current TM contamination problem in farmland is serious, and the competent authorities have adopted a method for remediating contaminated farmland in stages. However, most of the contaminated sites have been rectified and deregulated. At the same time, the strict inspection and control of the wastewater discharge from surrounding high-polluting factories are being performed to avoid further contamination of farmland with TMs. However, most of the contaminated sediments in the irrigation channels have not been properly treated. The current control strategy only focuses on monitoring the quality of discharge and irrigation water. To be able to avoid the recontamination of farmland with TMs, the management of the surrounding irrigation channels, especially of sediment contaminated with TMs, should take priority over the remediation of the farmland.

4. Conclusions

In this study, we combined the water quality of the irrigation system and the concentrations of TMs in the sediment to establish a farmland TM contamination and management model (FTM_CMM). The model is simple and easy to use and is suitable for the long-term evaluation of paddy field management strategies. The simulation shows that the TM content in the topsoil and the corresponding TM content in the irrigation channel sediment are 56% to 72%. This highlights that, in addition to the management of irrigation water quality, the management of TM-contaminated sediment in adjacent canal irrigation systems is also an important part of the prevention of TM contamination in farmland.