Rare earth elements in paddy fields from eroded granite hilly land in a southern China watershed

There are large amounts of ion-adsorption rare earth resources in the granite red soil region of southern China, and exploitation of rare earth elements (REEs) has caused serious soil erosion and soil pollution in the area. In this study, the spatial variability of soil REEs in Zhuxi watershed, Changting County, southern China, was analyzed using a geostatistics method. The analysis produced several important results: (1) The content of total rare earth elements (TREEs) in the soil samples ranged from 56.04 to 951.76 mg kg−1, with a mean value of 255.34 mg kg−1, which was higher than the background value of soil in China. The REE variables showed strong positive Ce anomalies and strong negative Eu anomalies, with mean values of 2.26 and 0.44, respectively. (2) The contents of TREEs in five subtypes of the soils were different, but they had broadly similar curves of chondrite-normalized REE patterns, with steeper patterns from La to Eu and flatter patterns from Eu to Y. (3) The spatial variability of light rare earth elements (LREEs) was mainly affected by natural factors, but the spatial variabilities of heavy rare earth elements (HREEs) and TREEs were influenced by the combination of natural factors and anthropogenic factors. Soil erosion can contribute significantly to REE migration, especially for HREEs. (4) The distribution of TREEs showed that the high content of TREEs was in the lowland of the western watershed. By comparing the distributions of TREEs in paddy fields and hilly land, we found that the area with a high content of TREEs was greater in paddy fields than in hilly land, so we deduced that REEs migrate from hilly land to the paddy field and accumulate in the soil there.


Introduction
REEs as a group of 17 elements comprising 15 lanthanides (La, Ce, Pr, Nd,Pm, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, and Lu) plus Y and Sc were defined by International Union for Pure and Applied Chemistry [1]. They have similar physical properties, chemical properties, and a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 The soil samples for REE determination were crushed with an agate mill, and the ground samples were passed through a 0.149 mm polyethylene sieve. Soil samples of 0.1g were dissolved in a Milestone Microwave Laboratory System (Multiwave 3000, Anton Paar, Austria) in a combination solution of hydrofluoric, hydrochloric, and nitric acid (HF 40%: HCl 38%: HNO 3 70% = 1:1:3). All soil samples were analyzed for 15 REEs (Sc was excluded due to its different chemical properties and Pm is a trace material in nature) by inductively coupled plasma mass spectrometry (ICP-MS, X Series 2, Thermo Scientific, USA). The Detailed analytical procedures for the REE analysis are described by our previous paper [13].

Data analysis
All statistical analyses were performed using SPSS19.0 software (SPSS Inc., Chicago, IL, USA). Normality was verified using the value of skewness and kurtosis. Natural logarithmic transformation was used to meet the assumptions of normality when the raw data did not obey a normal distribution, but the raw (untransformed) data are listed in the tables in this paper. From the analysis based on the data of the 15 REEs, we calculated six REE variables, total rare earth elements (TREEs), LREEs, HREEs, L/H, δCe, and δEu. TREEs is the sum of the 15 REEs; LREEs is the sum of La, Ce, Pr, Nd, Sm, and Eu; HREEs is the sum of Gd, Tb, Dy, Ho, Er, Tm, Yb, Lu, and Y; and L/H is the ratio of LREEs and HREEs. δCe and δEu are defined as follows: where N refers to chondrite normalization [23]. The REE contents were normalized to chondrite for the purpose of comparison [24]. The geostatistical method was used to study the spatial variability of REEs in soils of the watershed. Geostatistics is based on the theory of a regionalized variable, which is distributed in space with spatial coordinates and shows spatial autocorrelation such that samples close together in space are more alike than those farther apart [25]. The geostatistics approach consists of the following two parts: the calculation of an experimental variogram from the data and model fitting, and estimation at unsampled locations [26,27]. The semivariogram of each soil property was constructed using the following model: where γ (h) is the semi-variance for the internal distance class h, h is the lag interval, and N(h) is the total number of sample pairs for the lag interval h. Z(xi) is the measured sample value at point i, and Z(xi + h) is the measured sample value at point i + h. The semivariogram was calculated using GS+ 7.0 software (Gamma Design Software Inc., Plainwell, MI, USA) and the variability maps were created in ArcGIS 10.2 software (ESRI Inc., Redlands, CA, USA).

REE contents and REE variables
The descriptive statistics of the REE contents and REE variables of the soils, which include mean, range, standard deviation, coefficient of variation, skewness, and kurtosis, are summarized in Table 1. All the data needed to obey a normal distribution in the subsequent semi-variance analysis, so the values of skewness and kurtosis were used to test the data for normal distribution. The value of skewness should be within the range of ±2 and the value of kurtosis should be less than 3; otherwise, it is regarded as an extreme [28]. The skewness and kurtosis results indicated that the REE contents and REE variables did not obey a normal distribution except for Ce and δEu. Therefore, we used logarithmic transformation to normalize the data. Thus, all the data successfully passed the test for normality.
The coefficient of variation values (CVs) can be used to compare the discrete degree of a property. Low CVs correspond to a spatially homogeneous distribution, and high CVs correspond to a non-homogeneous distribution [29,30]. The CVs of Pr, Nd, Sm, Eu, Gd, Tb, Dy, and Ho were greater than 100%, which indicates a high variability. The CVs of La, Ce, Er, Tm, Yb, Lu, and Y were moderate, with fluctuation only from 40.38% to 91.66%. High CVs for REE contents have also been documented by other research [30]. Previous studies have shown that CVs of heavy metals originated from natural sources are relatively low, whereas CVs of heavy metals affected by anthropogenic activities are quite high [30][31][32]. It can be concluded that exploitation activities and soil erosion influence the differentiation of REEs greatly. The CVs of the 6 REE variables were moderate, with fluctuation from 24.6% to 87.69%. The CVs of LREEs were moderate compared to those of HREEs because LREEs were largely depended on the Ce, which had the lowest CV among the 15 REEs. This indicates that HREEs were affected more strongly by anthropogenic activities than LREEs. The content of TREEs in the soil samples ranged from 56.04 to 951.76 mg kg −1 , with a mean value of 255.34 mg kg −1 . The mean value of TREEs is higher than the background value of soil in China (186.76 mg kg -1 ) [33] and higher than those in Japan (98 mg kg −1 ) [34] and Australia (105 mg kg −1 ) [35]. The content of LREEs ranged from 38.11 to 627.95 mg kg −1 , with a mean value of 178.3 mg kg −1 , and the content of HREEs ranged from 9.65 to 323.81 mg kg −1 , with a mean value of 77.04 mg kg −1 . The L/H value ranged from 0.44 to 10.73, with an average of 3.28, which shows that the content of LREEs is significantly higher than that of HREEs. This study also revealed strong positive Ce anomalies (2.26) and strong negative Eu anomalies (0.44), which is in agreement with the results of other studies in southern China, such as at Hainan Island [36] and in southern Jiangxi Province [37]. The contents of the 15 REEs tended to follow the Oddo-Harkins rule, which states that even atomic numbers are more frequent than their neighbors with odd atomic numbers. The order of the means of the 15 REE contents in the soil is Ce Similar results were reported in other studies [14,30,34].

REE patterns in different soil subtypes
The soil types in the study area include paddy soil, alluvial soil, and red soil. The area of red soil accounts for more than 74% of the total area. To distinguish the distribution of TREEs among the different soil types, the soil types were divided into subtypes. The subtypes of soil in the study area were hydromorphic paddy soil, percolated paddy soil, alluvial soil, red soil, and coarse red soil. The content of TREEs in the five soil subtypes is shown in Table 2. For comparison, the contents of TREEs were normalized to chondrite [24], as shown in Fig 2. The highest mean TREEs content of 330.36 mg kg -1 was determined in the percolated paddy soil, followed by the hydromorphic paddy soil, alluvial soil, red soil, and coarse red soil, respectively. Although there were some differences, the five subtypes of soil had similar curves of chondrite-normalized REE patterns, with steeper patterns from La to Eu and flatter patterns from Eu to Y.
The content of a TREEs pattern is linked to the soil parent material [38,39], and the parent material in this study area is the same of coarse-grained granite, so all five soil subtypes showed similar patterns. Ample evidence has shown that REEs can be mobilized during weathering processes in the tropical and subtropical monsoon climates in southern China [23,40,41], which is why the REE content presents differences between the soil subtypes in this study. Our previous study determined that REEs migrate downhill under the influence of water flow and gravity [23]. The higher TREE contents of the hydromorphic paddy soil, percolated paddy soil, and alluvial soil than the red soil and coarse red soil are probably caused by differences of altitude. Comparing the TREE contents in the red soil and the coarse red soil, we found that the red soil only 19.09 mg kg -1 less than the average of all five subtypes of soil, and the coarse red soil was much lower than the average. Considering that the coarse red soil was located in the region where soil erosion has occurred seriously in the study area [21], we can conclude that soil erosion is one of the main factors for REE migration.

Spatial variability of REEs
The REE contents in the soil may have been affected by intrinsic or extrinsic sources that cannot be discriminated by descriptive statistics. Thus, the spatial correlation structure of the REEs was explored by semivariogram. This method not only considered the randomness of the data, but also the spatial structure characteristics of the data. The best variogram model and parameters of REE contents and REE variables are listed in Table 3. They were fitted with spherical, exponential, linear, and Gaussian semivariograms with determination coefficient values (R 2 ) ranging from 0.449 to 0.958. The R 2 of all fitting results is greater than 0.3, indicates that the fitting results is better [42,43].
The ratio of nugget effect (C 0 ) over sill (C 0 + C) is the nugget variance, which expresses the percent of total semi-variance [30]. It was used to judge the spatial dependency of the REE parameters and provide a quantitative basis for interpolating unsampled locations. If the ratio was equal to or lower than 25%, the variables were considered strongly dependent; if it was between 25% and 75%, the variables were considered moderately dependent; and if it was greater than 75%, the variables were considered weakly dependent [44,45]. In this study, the semivariograms indicated strong spatial dependence for REEs such as La, Pr, Nu, Sm, Eu, and LREEs but indicated moderate spatial dependence for Ce, Tb, Dy, Ho, Er, Dy, Yb, Lu, Y, HREEs, TREEs, L/H, and δEu. Only Gd and δCe showed weak spatial dependence. The strong spatial dependence of REEs in the soil may be the result of natural factors, such as strong pedogenesis [39]. Thus, we can see that the 6 LREEs except for Ce, present strong spatial dependence. With the variables of LREEs presenting strong spatial dependence, we concluded that the LREEs in the soil were affected mainly by natural factors. Moderate spatial dependence indicates that anthropogenic factors changed the soil texture spatial correlation through activities such as farming, management, practices, industrial production, and other human activities  [46,47]. The 9 HREEs, except for Gd, presented moderate spatial dependence, which indicates that HREEs were more affected by anthropogenic factors than LREEs. The TREEs were influenced by a combination of natural and anthropogenic factors [48]. In the case of Gd, as the transition element between LREEs to HREEs, the difference between sampling sites may have weakened the spatial correlation. δCe also presented weak spatial dependence caused by a complicated oxidation-deoxidation environment.
The range values showed large variability among the REEs, which can be a useful principle for mapping [49,50]. The results indicate that the spatial correlations (range) of the REE contents and REE variables vary widely from 510 m (Eu) to 5586 m (HREEs). The different ranges of the spatial dependence among the REEs may be a result of the parent material, erosiondeposition factors, or topography [48]. A large range indicates that the observed values are affected over greater distance by other values of the parameter compared to the parameters with smaller ranges [51]. Thus, the range of 5586 m for HREEs implies that HREEs values influenced neighboring values over greater distances than did values of other soil variables.

Spatial distribution of REEs
To express the REE distribution clearly, six REE variables were chosen for spatial interpolation by Kriging. The spatial distributions of these REE variables are shown in Fig 3. It is useful to identify areas with high contents of REEs and assess the possible variation of REEs in the watershed.
The distribution map of LREEs presents an irregular trend in different directions and some discontinuous high content areas located in the central part of the watershed. On the whole, areas with a high content of LREEs are located in the west and northwest. The map of HREEs distribution shows an increasing trend from east to west and northwest, except for the small plot of low content in the central part of the watershed. The distribution map of TREEs presents a trend similar to the distribution trend of HREEs. The topography of the watershed is tilted from east to west, and a previous study found that REEs can migrate downhill to lower sites by water flow and gravity [23], so the high REE content area is located in the lowland of the west. The central part of the watershed is the area where soil erosion was serious [21]. The distribution map of L/H shows a decreasing trend from east to west, except for the central part of the watershed. The soil erosion caused migration of REEs and preferential release of HREEs resulting in a greater depletion of HREE than LREE [52,53]. Thus, we further demonstrated that soil erosion caused migration of HREEs more than LREEs.
The depletion or enrichment of Ce and Eu usually occurs in nature due to their oxidation state and mobility under different oxidation-reduction conditions [54]. The distribution of δCe shows a decreasing trend from the central part of the area to the east and to the west. A strongly positive Ce anomaly area (δCe >2.5) is located in the north-central part of the watershed, and it may be influenced by weathering processes [4] and soil erosion. Ce 3+ oxidized to Ce 4+ under the exogenous environment was stored in the surface soil stably, and soil erosion caused losses of other REEs producing positive Ce anomaly development in this area [54]. The distribution of δEu shows a decreasing trend from southeast to northwest. The Eu anomaly is affected mainly by parent material, and the development process of the soil is also a process of Eu loss increase [30]. Thus, because the parent material originates form coarse-grained biotitic granite, the variation of negative Eu anomaly was small, ranging from 0.20 to 0.60. The strongly negative Eu anomaly area is located in the northwestern part of the watershed where the degree of soil maturation is high. The Ce positive anomaly and Eu negative anomaly indicate that differentiation occurred between Ce, Eu, and other REEs in the weathering process of the parent rock [30,55]. Rare earth elements in paddy fields from eroded granite hilly land

Distribution of TREEs between paddy fields and hilly land
From the analysis above, we concluded that REEs migrated due to soil erosion. To verify whether REEs accumulated in lowland, where paddy fields are largely located, we compared the TREE distribution between paddy fields and hilly land. The areas of different classes of TREEs are shown in Table 4.
The Table 4 shows that the area of TREE content lower than 180 mg kg -1 (background value of soil in China is 187.60 mg kg -1 ) [33] accounted for more than 41.14% of the distribution in the hilly land but for only 6.05% in the paddy fields. Additionally, in the class of 180-230 mg kg -1 (background value of soil in Fujian Province is 223.47 mg kg -1 ) [13], the area in the hilly land accounted for 20.16%, which means that about 61.30% of the area of TREE content in the hilly land is lower than the Fujian background value. However, the area of TREE content in paddy fields higher than the Fujian background value accounted for more than 66.83% of the total area. The percentage of high content classes of TREEs was greater in paddy fields than in hilly land. The parent material originated from the coarse-grained biotitic granite in the watershed, so the distribution of REEs in soil derived from the parent material should be similar in different regions. The difference of REEs in soil between hilly land and paddy fields further indicates that REEs migrate in the watershed. Combining these results with the previous findings, we suggest that soil erosion can significantly contribute to the migration of REEs in the soil. The REEs migrated from hilly land to the paddy fields in the watershed and then accumulated there.

Summary and conclusions
The content of REEs in the soil of the granite area in southern China was higher than the average value in China and higher than the average value in Fujian Province. The trend of REEs in soil followed the decreasing order of Ce > Y > La > Nd > Pr >Gd > Dy > Sm > Yb > Er > Ho > Tb > Eu > Lu > Tm and obeyed to the Oddo-Harkins rule, presenting strong positive Ce anomalies and strong negative Eu anomalies, which is similar to the results of other studies in southern China. Although the contents of TREEs in the five subtypes soil were different, they had broadly similar curves of chondrite-normalized REE patterns, indicating that the REE pattern was affected by parent material.
The nugget variance indicated that the spatial variability of LREEs was mainly affected by natural factors, while the spatial variabilities of HREEs and TREEs were influenced by the combination of natural factors and anthropogenic factors. Soil erosion can significantly contribute to REE migration, especially for HREEs. The distribution of TREEs showed that the high-TREE content area was located in the lowland of the western watershed. Additionally, the percentages of high content classes of TREEs were greater in paddy fields than in hilly land, so we conclude that REEs can migrate from hilly land to the paddy fields in the watershed.