Study on Migration Monitoring Technology of Chromium-Contaminated Site Based on Dual-Frequency Induced Polarization Method

: In order to reveal the migration process of chromium-contaminated sites, COMSOL5.6 software was used to build the initial model of the chromium pollution source and analyze the migration change characteristics over time based on the rule of groundwater movement and Convection– Dispersion equation. COMSOL provides fully coupled multiphysics modeling for most engineering ﬁelds. The results show that: the permeable layer with a high permeability coefﬁcient during the migration of chromium pollutants has a strong adsorption effect on the heavily polluted part, and it will enhance the lateral diffusion of the polluted area within a short period. Based on the migration model, the intermediate gradient and the symmetrical quadrupole sounding model are built. The variation law of apparent amplitude frequency and apparent resistivity under ﬂat and undulating terrain is analyzed based on the dual frequency IP (Induced Polarization) method. The results show that: The intermediate gradient detection is better than the symmetrical quadrupole sounding in the correspondence between the peak value of apparent amplitude frequency and the valley value of apparent resistivity. The arrangement of the bathymetric measurement point in the center of the projected edge of the pollution ﬁeld can be a better way to monitor the pollution. Monitoring with the intermediate gradient detection device and the symmetrical quadrupole sounding device creates “multiple peaks” in the curve as it passes through the valley. Arranging the power supply electrodes on the mountain frontiers on both sides of the raised peaks and synthesizes the apparent amplitude frequency and apparent resistivity curves of the pollution ﬁeld, which can effectively determine the speciﬁc orientation of the heavily polluted area of the pollution ﬁeld. This research makes theoretical additions to the migration characteristics of the Cr pollution ﬁeld. It provides technical guidance for the real-time monitoring of the pollution ﬁeld, which is of positive signiﬁcance for promoting ecological and environmental protection.


Introduction
With the rapid development of China's economy and society, the problem of soil pollution is becoming increasingly serious. In total, 82.8% of the polluted soil is contaminated with heavy metals such as chromium (Cr), cadmium (Cd), lead (Pb), copper (Cu), zinc (Zn), and nickel (Ni), and the exceedance rate of Cr is 1.1% [1]. Cr has nine valence states in soil but is mainly present in two stable valence states, trivalent and hexavalent, with hexavalent Cr being 100 times more toxic than trivalent Cr. Hexavalent chromium is highly toxic, non-degradable, soluble, mobile, concealed, and can have an adverse impact on ecological sustainability and human health [2]. When humans consume plants and animals containing hexavalent Cr, it can cause skin irritation [3], allergic dermatitis [4], DNA damage [5], and, in serious cases, cancer [6]. Soil quality is a strategic issue related to the sustainable development of the national economy and the overall progress of society.

Convection-Dispersion Modelling of Contaminated Sites
The migration process of Cr-contaminated sources is mainly influenced by convection and dispersion, and its migration process changes according to the following equation [17]: where: θ v is the volumetric water content of the soil; D ij is the tensor of dispersion coefficient; C is the concentration of hexavalent Cr in groundwater, mol/m 3 ; x i and x j are Cartesian coordinates; q i is the Darcy flow rate in the direction of x i , m/s; R k is the dissolved contaminant yield for the kth reaction out of n different reactions.

Numerical Simulation of the Migration Characteristics of the Contaminated Site
Based on the basic law of groundwater flow and the mechanism of pollutant migration, based on the finite unit method, the global flow coupling of Darcy's law module and the porous medium dilute matter transfer module was carried out by using the COMSOL software. A 250 m × 10 m two-dimensional chromium heavy metal contamination field model was constructed (see Figure 1), with poorly permeable pulverized clay in material 1, fine sand with better permeability in material 2, and medium sand with good permeability in material 3. Only the parameters such as porosity, density, mass water content, and permeability coefficient of the material are set individually in the simulation. The resistivity of the whole model area is uniformly set to 390 Ω·m. Detailed model parameter settings are shown in Table 1. The raised part at the surface is a static storage of chromium pollution source, and the pollutant is a solution of Cr ions with relative concentration of 3 mol/m 3 , and its projected orientation on the ground is within 50 m to 100 m. The upper boundary (y = 0 m) is set to have a rainfall of 500 mm. The pressure head is set to be at x = 250 m, 5.4 m from the surface. The groundwater demarcation line is set to be at x = 0, and no fluid passes through the two sides of the demarcation line. The bottom of aquifer (y = −10 m) is set without flow. The rest of the parameter settings are shown in Table 1. In order to improve the solution accuracy and computational efficiency, the quadrilateral mesh adaptive algorithm is used to mesh the constructed model area [21]. Then, the mesh is automatically encrypted in the pollutant area and near the electrode arrangement line, and the mesh is automatically sparsed in the non-data collection area. The boundary conditions and parameters are set with reference to the numerical forward simulation in [17,22]. Appl Studies have shown that porosity  , mass water content w  , and Cr ion concentration (mg/kg)  are the main factors influencing soil resistivity [23,24], and the relationship between these and soil resistivity can be described by the Alchian equation [25]:  Studies have shown that porosity ϕ, mass water content θ w , and Cr ion concentration λ(mg/kg) are the main factors influencing soil resistivity [23,24], and the relationship between these and soil resistivity can be described by the Alchian equation [25]: where ρ 0 is the rock resistivity, Ω·m; ρ w is the pore water resistivity in the soil, Ω·m; k, m, and n are the soil property correlation coefficients. In Equation (2): where Λ is the molar conductivity of the soil solution, S·m 2 /mol; n 0 is the mass concentration of Cr ions, mol/L; θ v is the volumetric water content; λ is the concentration of Cr contaminants, mg/kg; ρ is the density of the soil medium, g/cm 3 . The empirical equation for soil resistivity derived from the literature shows that [26]: where ϕ is the porosity, θ w is the mass water content, and λ is the concentration of Cr ions in the soil (mg/kg), which are the main factors affecting the resistivity of the soil. According to Equation (3), the resistivity of the contaminated site can be calculated based on the porosity and mass water content of the contaminated site in a soil with a certain density. The migration process of the contaminated field over time was simulated using COM-SOL software. The resistivity changes were obtained for migration times of 10, 100, 200, 300, 400, and 500 days, respectively, as shown in Figure 2.
As shown in Figure 2, the migration of Cr pollutants within 10-100 days in the pulverized clay area is mainly influenced by gravity. The apparent resistivity of the pollutants is uniformly distributed in a step-like manner. The contaminated area did not undergo significant lateral migration due to the poor water permeability of the pulverized clay and the insignificant effect of convective dispersion. The ground pollutant has completely disappeared in 100 days, and the depth of pollution reaches 4.7 m, at this time, the polluted area is about 265 m 2 . As well, 100 days after the pollutant enters into the fine sand layer, it starts to spread obviously along the lateral direction, and the scope of the polluted area varies greatly in 100-200 days. The center of the pollutant migrates in the lateral direction from 50-100 to 122-127 m. The pollutant field spreads to the inner part of the aquifer, and the bottom edge has come to 9.8 m. Which indicates that the high permeability coefficient of the fine sand and the aquifer have a great influence on the migration of the pollution field. From the 200th day onwards, Cr pollutants take the pollution center (dark blue area) as the radiation source and keep the whole pollution field spreading downward and to the right at the same time spreading to the surrounding area in the form of a ladder. At this time, the pollution field has already radiated a new semi-ellipsoidal contamination area. The top of it is again approaching to the ground vertically, and its height is even more than that of the original contaminated field. On the 300th day, the heavily contaminated area had completely left the fine sand area. The whole contaminated area was between the fine sand and the medium sand area, with a tendency to continue to penetrate into the medium sand area. At this moment, the contaminated area had already reached 746 m 2 , which was 12 times as much as that on the 10th day. After 400 days of migration time, the upper right part of the pollution field has moved laterally into the middle sand area and spread laterally over a large area in the middle sand area. Meanwhile, the rest of the field continues to migrate downward in the aquifer. Then, 500 days later, the middle sand area and its vicinity had encompassed nearly half of the polluted area, which indicates that the permeable layer with a high permeability coefficient has a very good adsorption effect on the pollutants. At this juncture, the pollution field has already migrated out of the model range and spreads to a deeper and wider area. Therefore, it is very necessary to control the migration of pollution sources as early as possible. In Figure 2, the overall contaminated area gradually increases with time, but the dark blue heavily contaminated area is gradually decreasing. Which is due to the continuous diffusion of Cr-contaminated solutes to the surrounding area, resulting in a slight decrease in the concentration of the center of the contamination. The above shows that if there are factors such as water-rich soil or permeable layer under the ground that can accelerate the diffusion of pollution, it will lead to a rapid increase in the area of heavy metal contamination, and even produce multiple heavy contamination main areas.
where 0  is the rock resistivity, Ω·m; w  is the pore water resistivity in the soil, Ω·m; k, m, and n are the soil property correlation coefficients. In Equation (2): , where Λ is the molar conductivity of the soil solution, S·m 2 /mol; 0 n is the mass concentration of Cr ions, mol/L; v  is the volumetric water content;  is the concentration of Cr contaminants, mg/kg;  is the density of the soil medium, g/cm 3 . The empirical equation for soil resistivity derived from the literature shows that [26]: where  is the porosity, w  is the mass water content, and  is the concentration of Cr ions in the soil (mg/kg), which are the main factors affecting the resistivity of the soil. According to Equation (3), the resistivity of the contaminated site can be calculated based on the porosity and mass water content of the contaminated site in a soil with a certain density.
The migration process of the contaminated field over time was simulated using COMSOL software. The resistivity changes were obtained for migration times of 10, 100, 200, 300, 400, and 500 days, respectively, as shown in Figure 2. As shown in Figure 2, the migration of Cr pollutants within 10-100 days in the pulverized clay area is mainly influenced by gravity. The apparent resistivity of the pollutants is uniformly distributed in a step-like manner. The contaminated area did not undergo significant lateral migration due to the poor water permeability of the pulverized clay and the insignificant effect of convective dispersion. The ground pollutant has completely dis-

Comparative Analysis of Numerical Forward Simulation and Soil Trench Modeling Experiments
Considering the practical factors, such as the inability to control the infiltration of chromium solution in the soil and the irregularity of the chromium-contaminated area, it is very difficult to build a physical experimental model of Cr pollutant transport process in soil tanks. The paper adopts the COMSOL software for numerical forward simulation of the chromium contamination field migration process of electrometallurgical monitoring. In this section, the reliability and accuracy of the numerical orthogonal simulation are verified by static soil trench physical modeling experiments.

Experimental Program Design
The test area (6 m × 2 m) was selected and the average resistivity of the surface soil in the tested area was 537 Ω·m. The specification of the chromium-contaminated body was 50 cm × 40 cm × 20 cm as shown in Figure 3. The bottom layer of 15-20 cm is dry clay with a resistivity of about 759 Ω·m, which serves to ensure that the chromium solution stored in the fine sand will not be lost by downward infiltration. The middle layer of 5-15 cm is the chromium-contaminated area, which is made of fine river sand with high water absorption and permeability. By grinding fine river sand to a particle size of 0.25-0.35 mm (see Figure 4) and compacting it after passing through a soil sieve, the final layer was made. After the fine sand was filled and compacted, 2 mol/m 3 of Na 2 CrO 4 solution (yellow liquid in the beaker in Figure 4) was uniformly placed on the surface and left for 1 h to make the chromium solution penetrate and mix with the fine sand. The soil at the top of the pollutant was sandy soil from the experimental area, and this type of sandy soil was placed at the top, 0-5 cm. After the preparation of the chromium-contaminated body was completed, it was determined that the average resistivity of the chromium-contaminated body at the depth of 5~15 cm was about 93 Ω·m.

Comparative Analysis of Numerical Forward Simulation and Soil Trench Modeling Experiments
Considering the practical factors, such as the inability to control the infiltration of chromium solution in the soil and the irregularity of the chromium-contaminated area, it is very difficult to build a physical experimental model of Cr pollutant transport process in soil tanks. The paper adopts the COMSOL software for numerical forward simulation of the chromium contamination field migration process of electrometallurgical monitoring. In this section, the reliability and accuracy of the numerical orthogonal simulation are verified by static soil trench physical modeling experiments.

Experimental Program Design
The test area (6 m × 2 m) was selected and the average resistivity of the surface soil in the tested area was 537 Ω·m. The specification of the chromium-contaminated body was 50 cm × 40 cm × 20 cm as shown in Figure 3. The bo om layer of 15-20 cm is dry clay with a resistivity of about 759 Ω·m, which serves to ensure that the chromium solution stored in the fine sand will not be lost by downward infiltration. The middle layer of 5-15 cm is the chromium-contaminated area, which is made of fine river sand with high water absorption and permeability. By grinding fine river sand to a particle size of 0.25-0.35 mm (see Figure 4) and compacting it after passing through a soil sieve, the final layer was made. After the fine sand was filled and compacted, 2 mol/m 3 of 4 2 CrO Na solution (yellow liquid in the beaker in Figure 4) was uniformly placed on the surface and left for 1 h to make the chromium solution penetrate and mix with the fine sand. The soil at the top of the pollutant was sandy soil from the experimental area, and this type of sandy soil was placed at the top, 0-5 cm. After the preparation of the chromium-contaminated body was completed, it was determined that the average resistivity of the chromium-contaminated body at the depth of 5~15 cm was about 93 Ω·m.
The physical experimental model of the constructed soil trench is shown in Figure 5. The chromium-contaminated body was placed in the center of the area, and the burial depth was from 5 cm to 25 cm. The middle gradient is used to lay the electrodes, A and B are the power supply electrodes, and the pole distance is 6 m. M and N are the measuring electrodes, and the pole distance is 10 cm. YDZ75(A) parallel DC electrometer is used for the detection, and the transmi er module transmits the 100-mA current as the excitation, and the receiver module is used for the data acquisition and processing.    The physical experimental model of the constructed soil trench is shown in Figure 5. The chromium-contaminated body was placed in the center of the area, and the burial depth was from 5 cm to 25 cm. The middle gradient is used to lay the electrodes, A and B are the power supply electrodes, and the pole distance is 6 m. M and N are the measuring electrodes, and the pole distance is 10 cm. YDZ75(A) parallel DC electrometer is used for the detection, and the transmitter module transmits the 100-mA current as the excitation, and the receiver module is used for the data acquisition and processing.

Numerical forward Simulation of the Dual-Frequency Excitation Method
When using the dual frequency IP method for contamination site detection, the IP parameter apparent amplitude frequency Fs is calculated from the potential difference between the two measurement points and is expressed as [15]: denote the difference between the high and low frequency potentials measured between the measurement electrodes MN when the high and low frequency currents are supplied, respectively, mV.
The other IP parameter, apparent resistivity ρs, is calculated by the following equation [15]: where K is the electrode device factor; Δ VMN is the potential difference between the two measurement points of the measuring electrode when powered by high frequency current, V; I is the high frequency supply current, A; AM, AN, BM, and BN are the four electrode spacings, m.

Numerical forward Simulation of the Dual-Frequency Excitation Method
When using the dual frequency IP method for contamination site detection, the IP parameter apparent amplitude frequency F s is calculated from the potential difference between the two measurement points and is expressed as [15]: where ∆V H and ∆V L denote the difference between the high and low frequency potentials measured between the measurement electrodes MN when the high and low frequency currents are supplied, respectively, mV. The other IP parameter, apparent resistivity ρ s , is calculated by the following equation [15]: where K is the electrode device factor; ∆V MN is the potential difference between the two measurement points of the measuring electrode when powered by high frequency current, V; I is the high frequency supply current, A; AM, AN, BM, and BN are the four electrode spacings, m.

Introduction of the Cole-Cole Model in Dual Frequency IP Method
The Cole-Cole model has been widely used in geological exploration [27]. Liu Haorui et al. have verified the feasibility of the Cole-Cole model for characterizing the complex resistivity of soils through experiments on Cr-contaminated soils [28]. The Cole-Cole model is one of the models commonly used by researchers to describe the excitation effect. In this study, the Cole-Cole complex resistivity model is introduced in the numerical simulation, which fully considers the effects of contaminants concentration, soil particle homogeneity, porosity, and denseness of the contaminants on the excitation effect of the contaminated site. Cole-Cole can be expressed in the frequency domain as [14]: where ω is the angular frequency, rad/s; ρ 0 is the zero-frequency resistivity, Ω·m, which corresponds to the DC resistivity in the time domain; m is the polarization rate; τ is the time constant; c (0.1 ≤ c ≤ 0.6) is the frequency correlation coefficient. Where ρ 0 characterizes the electrical conductivity; m characterizes the strength of the excitation effect; τ is a parameter characterizing the position of the frequency band where the IP spectrum changes; and c characterizes the degree of dithering of the IP spectrum. The electrical conductivity and relative dielectric constant were set in the material details column of the contaminated site. The conductivity value is equal to the reciprocal of the Cole-Cole model's zero frequency resistivity, the simulation's settings are shown in Figure 6. i is an imaginary unit; rho stands for ρ 0 ; tao stands for τ; pi stands for π; 2 × pi × freq stands for ω. The relative permittivity ε r is related to the polarization rate m in the Cole-Cole model as [16]: ulation, which fully considers the effects of contaminants concentration, soil particle homogeneity, porosity, and denseness of the contaminants on the excitation effect of the contaminated site. Cole-Cole can be expressed in the frequency domain as [14]: where  is the angular frequency, rad/s; ρ0 is the zero-frequency resistivity, Ω·m, which corresponds to the DC resistivity in the time domain; m is the polarization rate; τ is the time constant; c (0.1 ≤ c ≤ 0.6) is the frequency correlation coefficient. Where ρ0 characterizes the electrical conductivity; m characterizes the strength of the excitation effect; τ is a parameter characterizing the position of the frequency band where the IP spectrum changes; and c characterizes the degree of dithering of the IP spectrum. The electrical conductivity and relative dielectric constant were set in the material details column of the contaminated site. The conductivity value is equal to the reciprocal of the Cole-Cole model's zero frequency resistivity, the simulation's se ings are shown in Figure 6. i is an imaginary unit; rho stands for ρ0; tao stands for τ; pi stands for π; 2 × pi × freq stands for  . The relative permi ivity εr is related to the polarization rate m in the Cole-Cole model as [16]: The simulation schematic in COMSOL is shown in Figure 7. The background resistivity value is set to 537 Ω·m. The chromium contaminant 0-5 cm layer is set to 537 Ω·m. The middle layer 5-15 cm is set to 93 Ω·m. The bo om layer 15-20 cm is set to 759 Ω·m. The model specification was set to 1:1 for the experiment.

Acquisition of Numerical Simulation Results
In the simulation, an ultra-low frequency sinusoidal current was used as the excitation. The frequency ratio S = fH/fL = 13 was selected for the dual frequency current, with a high frequency of 4 Hz and a low frequency of 0.308 Hz. The current amplitude is 100 mA.
The measured potential will have real and imaginary parts since the Cole-Cole complex resistivity model is introduced in the numerical forward simulation. In this case, the imaginary part should be removed, and the real part of the potential signal should be utilized to calculate the apparent amplitude frequency. In Figure 7, there are 60 measurement points from left to right, which corresponds to 60 groups of high and low frequency potential values. Taking measurement points 1-20 as an example, the apparent amplitude frequency Fs is calculated according to 20 groups of high and low potential differences, and Table 2 is obtained as follows. Where represents the potential difference between points M and N when high (low) frequency current is supplied. Significant digits

Acquisition of Numerical Simulation Results
In the simulation, an ultra-low frequency sinusoidal current was used as the excitation. The frequency ratio S = f H /f L = 13 was selected for the dual frequency current, with a high frequency of 4 Hz and a low frequency of 0.308 Hz. The current amplitude is 100 mA.
The measured potential will have real and imaginary parts since the Cole-Cole complex resistivity model is introduced in the numerical forward simulation. In this case, the imaginary part should be removed, and the real part of the potential signal should be utilized to calculate the apparent amplitude frequency. In Figure 7, there are 60 measurement points from left to right, which corresponds to 60 groups of high and low frequency potential values. Taking measurement points 1-20 as an example, the apparent amplitude frequency F s is calculated according to 20 groups of high and low potential differences, and Table 2 is obtained as follows. Where ∆V H (∆V L ) represents the potential difference between points M and N when high (low) frequency current is supplied. Significant digits of ∆V H and ∆V L were set to five in COMSOL and take the numerical precision to four decimal places. Similarly, when the device coefficient K value at 1-20 measurement points is calculated, the apparent resistivity can be calculated according to the ∆V MN obtained from numerical simulation, and the apparent resistivity under 1-20 measurement points is made as Table 3:

Analysis of Experimental Results
In the numerical forward simulation, the potential value is calculated based on the Cole-Cole model introduced by the dual-frequency excitation method, combined with Equations (4) and (5) to calculate the apparent amplitude frequency and apparent resistivity, respectively. The calculated values together with the experimentally measured values make the apparent resistivity and apparent amplitude frequency characteristic curves as shown in Figure 8. It can be found that the shape of the apparent amplitude frequency F s curve in the figure is similar to the inverted apparent resistivity ρ s curve, and there is a good correspondence at the same measurement point. In order to further analyze the accuracy of COMSOL simulation, an error analysis is required (see Table 4). Because the pollution field in the experiment and simulation are symmetric structures, it is sufficient to take 1-30 measurement points for the error analysis. Where error = |(simulation value − experimental value)/experimental value|. As can be seen from Figure 8 and Table 4, at the edge of the monitoring area (measurement points 1-20), since the soil was untreated here in the experiment, differences in soil particle properties would have perturbed the apparent resistivity values. The Fs and ρs errors are both less than 3%, which indicates that COMSOL also has a fairly high-precision detection effect for the uncontaminated area. The Fs and ρs between the experimental and simulated values reached more than 15% in the measurement points 21-27. This is because the contaminants in the simulation are only fixed in the set orientation and In order to further analyze the accuracy of COMSOL simulation, an error analysis is required (see Table 4). Because the pollution field in the experiment and simulation are symmetric structures, it is sufficient to take 1-30 measurement points for the error analysis. Where error = |(simulation value − experimental value)/experimental value|. As can be seen from Figure 8 and Table 4, at the edge of the monitoring area (measurement points 1-20), since the soil was untreated here in the experiment, differences in soil particle properties would have perturbed the apparent resistivity values. The F s and ρ s errors are both less than 3%, which indicates that COMSOL also has a fairly high-precision detection effect for the uncontaminated area. The F s and ρ s between the experimental and simulated values reached more than 15% in the measurement points 21-27. This is because the contaminants in the simulation are only fixed in the set orientation and do not infiltrate into the surrounding area. In contrast, in the experiment, the chromium solution inside the fine sand is not completely dried and volatilized, and it will still infiltrate the surrounding soil. The apparent resistivity ρ s decreases gradually when the MN measurement points are gradually approaching the chromium-contaminated body. The errors of F s and ρ s between the experimental and simulated values are both less than 5% at the projection position of the chromium-contaminated body (measurement points 28-30). In conclusion, monitoring the chromium pollution field using the Cole-Cole model introduced based on the dual frequency IP method is feasible.

Numerical Simulation of Dual Frequency Induced Polarization Monitoring of Contaminated Sites in Different Terrains
The most commonly used profiling device in the dual frequency IP method is the intermediate gradient detection device, which has flexible pole placement and reduces the number of moves of the supply electrodes [29]. The data obtained from the sideline measurements can be fully interpreted for the target anomaly. The symmetrical quadrupole sounding device requires a small supply current [29]. However, it can obtain a large potential difference, is more advantageous in complex terrain conditions, and is often used as the main device for dual frequency IP method sounding. Given both advantages, this paper will use an intermediate gradient detection device and a symmetric quadrupole sounding device to perform numerical forward simulations of contaminant site migration detection. In the simulation, an ultra-low frequency sinusoidal current was used as the excitation and the frequency ratio S = f H /f L = 13 was selected for the dual frequency current, with a high frequency of 4 Hz and a low frequency of 0.308 Hz and a current amplitude of 0.25 A.
The numerical forward simulation is based on the migration model of the contaminated field from 10 to 400 days in 1.2. The distribution of the contaminated field at 500 days is similar to that at 400 days, but the field has already migrated out of the model at 500 days, so it is not taken into consideration. Since the Cr contamination field is mostly a complex geological condition with undulating terrain, two geological conditions, flat and undulating terrain, need to be considered. Numerical forward simulations are carried out using the intermediate gradient detection device and a symmetric quadrupole sounding device, respectively.

Numerical Simulation of Intermediate Gradient Detection Device Detection
Taking the migration model of the pollution field on the 100th day as an example, the intermediate gradient detection model of the pollution field in different terrains was constructed as shown in Figure 9a,b. Where the height of each valley bump in Figure 9b was 2 m, and the width of each valley bump was 50 m. The migration detection model of the pollution field for other days was the same as that in Figure 9. In order to get more accurate excitation anomaly data, the AB electrodes were set at the two ends of the wiring, and the measurement electrode MN pole spacing was set to 5 m. The spacing of each point number in Figure 9a,b was 5 m (the horizontal distance of each measurement point in the case of undulating terrain was 5 m). The total number of measurement points was 50. The variation rules of apparent amplitude frequency and apparent resistivity at different topographies of the polluted field from 10 to 400 days were obtained, as shown in Figures 10 and 11, respectively. the measurement electrode MN pole spacing was set to 5 m. The spacing of each point number in Figure 9a,b was 5 m (the horizontal distance of each measurement point in the case of undulating terrain was 5 m). The total number of measurement points was 50. The variation rules of apparent amplitude frequency and apparent resistivity at different topographies of the polluted field from 10 to 400 days were obtained, as shown in Figure 10 and Figure 11, respectively. As shown in Figure 10a, the apparent amplitude frequency on the 10th day of contaminated field migration jumps from 1.5% at point No.8 to 13.6% at point No.11 when the terrain is leveled. This is because there is a large concentration of Cr contaminants retained on the ground on the 10th day, which makes the surface contaminated area and the seepage into the subsurface portion produce excitation anomalies at the same time, As shown in Figure 10a, the apparent amplitude frequency on the 10th day of contaminated field migration jumps from 1.5% at point No.8 to 13.6% at point No.11 when the terrain is leveled. This is because there is a large concentration of Cr contaminants retained on the ground on the 10th day, which makes the surface contaminated area and As shown in Figure 10a, the apparent amplitude frequency on the 10th day of contaminated field migration jumps from 1.5% at point No.8 to 13.6% at point No.11 when the terrain is leveled. This is because there is a large concentration of Cr contaminants retained on the ground on the 10th day, which makes the surface contaminated area and the seepage into the subsurface portion produce excitation anomalies at the same time, resulting in the jumps of apparent amplitude frequency at the boundaries of the source of the contamination. The depths of the center of the heavily polluted area were 2.2 m, 3.9 m, 5.2 m, and 5.5 m at the 100th, 200th, 300th, and 400th days, respectively. The corresponding peak values of the apparent amplitude frequency were 10.27%, 8.50%, 7.58%, and 6.30%, which indicated that the peak value of the apparent amplitude frequency decreased gradually with time when the pollution field spread downward. This is because as the depth of the pollution field migrates, the value of the secondary potential that can be detected gradually decreases, and the excitation effect anomaly has become less obvious. The projected orientation of the center of the heavily contaminated area on the ground at days 100, 200, 300, and 400 was 75 m, 130 m, 153 m, and 178 m, respectively. The corresponding point numbers of the peaks of the apparent amplitude frequency were 15, 25, 30, and 35, respectively, which indicated that the point numbers of the peaks of the apparent amplitude frequency gradually increased when the contaminated field migrated to the right. As the time increases, the angle of the apparent amplitude frequency curve also increases, indicating that the contaminated area of the detected soil area is increasing. On the 200th day, due to the good permeability of the fine sand, under the action of convective dispersion, the contaminated area spreads significantly in the transverse direction. The part that moves out of the fine sand area is subjected to the action of the aquifer and the pressure head in the longitudinal direction, which leads to the appearance of double peaks in the apparent amplitude frequency at the point number of 15-25 at this time. The peak on the left side represents the part of the heavy contamination that is retained in the fine sand area (the low-resistivity part corresponds to the dark blue color in Figure 2). The peak on the right side represents the part of the heavily contaminated area (the black-blue semi-ellipse part in Figure 2c).
As shown in Figure 10b, when the terrain is undulating, the overall trend of the apparent amplitude frequency curve change rule is similar to when the terrain is flat. However, the smoothness of the curve is worse than that when the terrain is flat. Thus, the judgment of the main body of the contamination will be affected by the undulation of the ground. It can also be seen from the figure that on the 10th and 100th days, the pollution field was closer to the ground. The peaks on the top of the polluted area were just above it, and at this time, the pollution field showed a trend of longitudinal migration with no obvious horizontal diffusion. At this time, the apparent amplitude frequency curve still showed a "peak" for a heavily polluted area. Which indicated that the effect of the raised mountain peaks on the shape of the apparent amplitude frequency curve was not obvious. It only weakened the value of the apparent amplitude frequency monitoring. On the 200th, 300th, and 400th days, the valley locations (point numbers 20, 30, and 40) caused significant fluctuations in the apparent amplitude frequency curves, which resulted in a "multi-peak" effect in the apparent amplitude frequency curves. Avoiding the valley locations during the monitoring is a better way to analyze the number and location of heavily polluted areas from the curve patterns.
From Figure 11a, it can be seen that the valleys of the apparent resistivity curves at the 100th, 200th, 300th, and 400th days correspond to the measurement point numbers 15, 25, 30, and 36, respectively. In contrast to Figure 10a, the valleys of the apparent resistivity curves in Figure 11a correspond to the peaks of the apparent amplitude frequency curves, almost the same as the orientations of the heavily polluted areas. In Figure 11b, the apparent resistivity curves are affected by the undulation of the ground, and the phenomenon of "multiple valleys" occurs in a heavily polluted area. This shows that it is not possible to accurately determine the orientation of the heavily polluted area only based on the apparent resistivity map. The numerical values of the measurement points where valleys occur increase with the migration of the pollution source to the right for both topographies. All the valleys of the apparent resistivity curves are near 30-60 Ω·m, which indicates that the downward migration of the pollution source does not have a significant effect on the valleys of the apparent resistivity curves.

Numerical Simulation of Symmetrical Quadrupole Sounding Device Detection
The symmetrical quadrupole sounding detection model of the pollution field with different topography is constructed as shown in Figure 12a,b. The migration detection model of the pollution field for other days is the same as that in Figure 12. In order to have enough area to arrange the A electrodes to the left, the model range was increased on the left side of the pollution field, which was increased from 250 m to 300 m. In order to make the model still reflect the excitation level at a deep enough position, the bottom of the model was set to be grounded. There are 601 points on the two types of terrain, each with a horizontal spacing of 0.5 m. The horizontal spacing of the measuring electrodes MN is 1 m, and O is the midpoint of MN. When a certain measuring point is used for the sounding work on the two types of terrain, the A and B electrodes are gradually increased towards the two ends. In the paper, under both methods, 75 sets of AB positions were taken and then 75 sets of potential values were measured separately. Then, the horizontal spacing of AB is gradually increased from 2 m to 76 m. In the simulation detection, three orientations of measurement points were arranged under the conditions of flat terrain and undulating terrain: The apparent amplitude frequency in Figure 13 and the apparent resistivity graph in Figure 14 both adopt logarithmic coordinates with 1/2 of the horizontal net distance of AB electrodes as the horizontal axis, which can better reflect the anomalies of the polluted field buried in the middle and deep. When AB/2 is taken as the horizontal axis, the measured depth tends to be the same as the actual depth within the 0-100 m sounding section [30,31]. The apparent amplitude frequency in Figure 13 and the apparent resistivity graph in Figure 14 both adopt logarithmic coordinates with 1/2 of the horizontal net distance of AB electrodes as the horizontal axis, which can better reflect the anomalies of the polluted field buried in the middle and deep. When AB/2 is taken as the horizontal axis, the measured depth tends to be the same as the actual depth within the 0-100 m sounding section [30,31].
The apparent amplitude frequency curves obtained from the sounding when leveling the terrain are shown in Figure 13a-c. From which it can be seen that, on the 10th day of the migration of the contaminated field, the peaks of the apparent amplitude frequency of the extended line measured at the measurement points 1-3 are 7.82%, 4.88%, and 1.97%, respectively. This indicates that a larger excitation anomaly can be detected by arranging the measurement points right above the center of the contaminated field. The more by the edge of the arrangement of the measurement points, the more the excitation anomaly is not obvious. On the 10th day of the migration, the depth of the contaminated field was very shallow, and the curves under all the measurement points showed a two-layer pattern. The apparent amplitude frequency F s decreased with the increase of the pole distance AB, according to which the depth of the contaminated field could be estimated. When the contaminated field migrates downward to leave the surface (100-400 days), the apparent amplitude frequency curve obtained at the measurement point 1 still shows a two-layer pattern. The F s increases with the increase of the pole distance AB and the depth of the contaminated field at the top of the field can be estimated according to the curve information. As shown in Figure 13b,c, the apparent amplitude frequency curves detected at measurement points 2 and 3 show a three-layer pattern from 100 to 400 days, regardless of the distance of the contaminated field from the ground. Taking 100 days as an example, at this time, the amplitude of the apparent amplitude frequency excitation anomaly at measurement point 2 is about 90% of that at measurement point 1, and that at measurement point 3 is about 34% of that at measurement point 1. However, the apparent amplitude frequency amplitude of the two measurements is much lower than that at point 1. The horizontal axis coordinates corresponding to the raised portion of the apparent amplitude frequency curves at measurement points 2 and 3 can be used to better estimate the depths of the contaminated field at the top and bottom of the field. The advantage of point 2 over point 3 is that the anomaly contrast of the apparent amplitude-frequency curve obtained from point 2 is greater when the depth of migration of the portion with the largest contamination concentration increases. Anomaly contrast refers to the visibility of the abnormality relative to the background value, defined as the ratio of the abnormal maximum value to the normal field value.
The apparent amplitude frequency curves obtained by a symmetrical quadrupole sounding device in undulating terrain are shown in Figure 13d-f. The overall trend of the apparent amplitude frequency F s curves is still similar to that of flat terrain. However, there are slight fluctuations due to the influence of terrain undulation. The apparent frequency F s excitation anomalies detected at point 1 in Figure 13d decrease with the increase of the number of days. However, the curves are all in a two-layer pattern, and it is not possible to judge the burial depth at the bottom of the contaminated field. The curves from point 2 in Figure 13e can more accurately estimate the depths at the top and bottom of the contaminated site, and its anomaly contrast is larger than that of point 3 in Figure 13d.
Through the above analysis, measurement point 2 is the best detection among the three measurement points in both terrains. The apparent resistivity change curve corresponding to measurement point 2 (see Figure 14) is now made for further analysis. As shown in Figure 14, the topographic relief will cause fluctuations in the apparent resistivity value under a symmetrical quadrupole sounding device. However, it will not affect the valley of the apparent resistivity. This means that the apparent resistivity at the center of (e) (f) The apparent amplitude frequency curves obtained from the sounding when leveling the terrain are shown in Figure 13a-c. From which it can be seen that, on the 10th day of the migration of the contaminated field, the peaks of the apparent amplitude frequency of the extended line measured at the measurement points 1-3 are 7.82%, 4.88%, and 1.97%, respectively. This indicates that a larger excitation anomaly can be detected by arranging the measurement points right above the center of the contaminated field. The more by the edge of the arrangement of the measurement points, the more the excitation anomaly is not obvious. On the 10th day of the migration, the depth of the contaminated field was

Determination of the Location of the Heavily Contaminated Area
In summary, it can be seen that: the integrated apparent amplitude frequency and apparent resistivity curve of the pollution field can be more effective in determining the specific location of the heavily polluted area of the pollution field. In a symmetrical quadrupole sounding device, the use of the edge of the pollution field to arrange the measurement point can be a better way to make a judgment on the top and bottom of the pollution field depth. This section takes the migration model of the pollution field on the 100th day as an example; synthesizing the monitoring curves of the apparent amplitude frequency and apparent resistivity of the pollution field and makes a comprehensive profile and sounding map (see Figure 15). (Sounding maps for both terrains were monitored using measurement point 2).
On the 100th day of pollutant migration, the exact position of the contaminated area in the ground projection is at 50-100 m, and the depth under flat terrain is at 0.1-4.7 m. With the depth of the top of the more contaminated area at 1.5 m, at which time the depth of the heavily contaminated area ranges from 1.5-4.7 m. The depth of the contaminated area under undulating terrain is at 2.1-6.7 m (with the height of the peak at 2 m). The depth of the top of the heavily polluted area is 3.5 m, at which time the depth of the heavily polluted area ranges from 3.5 to 6.7 m. From Figure 15, under the flat terrain, the horizontal projection orientation of the Cr heavy pollution area is located at 54.2-91.8 m, and the error is 4.2-8.2 m compared with the model, and the coincidence with the exact orientation is 75.2%. The sounding orientation is located at 1.9-4.3 m, and the error is 0.4 m, and the coincidence with the exact orientation is 75.0%. In the undulating terrain, the horizontal projection bearing is located at 50. rupole sounding device, the use of the edge of the pollution field to arrange the measurement point can be a be er way to make a judgment on the top and bo om of the pollution field depth. This section takes the migration model of the pollution field on the 100th day as an example; synthesizing the monitoring curves of the apparent amplitude frequency and apparent resistivity of the pollution field and makes a comprehensive profile and sounding map (see Figure 15). (Sounding maps for both terrains were monitored using measurement point 2). To summarize, the methods to improve the monitoring effect of the heavily polluted area are as follows: a.
Arrange the measuring points at the edge of the polluted field during the bathymetric sounding. b.
Combine the apparent amplitude frequency and apparent resistivity curve to make a comprehensive determination. c.
Arrange power supply electrodes on the mountain frontier on both sides.
The analysis on point c is as follows: the electric field lines of the power supply electrodes imported into the earth will be affected by the mountain area to produce a certain degree of convergence. This can better induce the polluted area to produce the polarization phenomenon and make the monitoring of the secondary potential, which can have better detection effect than that in the flat ground.

Discussion
The significance of this work is as follows: chromium pollutants are extremely harmful to nature, and analyzing their migration characteristics and monitoring their migration status in real time is of positive significance for carrying out soil management. The article has pointed out the superiority of using software to numerically simulate the migration of contaminated fields in the Section 5.
In the study of this paper, the author simulated the migration of chromium contamination field with time using COMSOL software. Based on the model of the contamination field obtained after the migration, electrical monitoring with an intermediate gradient detection device and the symmetrical quadrupole sounding device was carried out. On this basis, different topographic structures were introduced to analyze the effects of the migration of the pollution field on the apparent resistivity and apparent amplitude frequency curves.
First of all, the author used numerical forward simulation to simulate the migration of the contamination field because it is very difficult to implement experiments for validation. Difficulties include: The impossibility in practice of controlling the penetration of chromium solutions into the soil. The often irregular nature of chromium-contaminated areas. The fact that the spread of chromium-contaminated areas is subject to many uncertainties, and so on. Fortunately, although it is not possible to simulate the dynamic migration of chromium contamination field, the author has managed to implement electro-methodic monitoring of static chromium-contaminated areas, which has a positive significance in demonstrating the numerical simulation accuracy of our electro-methodic monitoring.
In future research, we advocate for the introduction of more influencing factors into the simulation of chromium contamination field migration. We will adapt more exploratory devices for electro-methodological monitoring, and the pursuit of a systematic analysis on a three-dimensional level. This is of positive significance for realizing the regional monitoring of chromium pollution.
In conclusions, despite some limitations, our study provides a large number of numerical simulations and experimental evidence, which illustrates that the intermediate gradient detection device and the symmetrical quadrupole sounding device are excellent for detecting chromium pollution field migration.

Conclusions
(1) Factors that can accelerate the spread of pollution, such as the presence of a permeable layer with a high permeability coefficient or a water-rich soil layer under the ground, which can lead to a rapid increase in the size of the area contaminated with heavy metals or even to the creation of multiple main areas of heavy pollution. Specifically, the permeable layer with a high permeability coefficient has a strong adsorption effect on the heavily polluted part. It will enhance the lateral diffusion of the polluted area in a short period. (2) The peaks of the apparent amplitude frequency curves corresponded well with the valleys of the apparent resistivity when the Intermediate Gradient detection device was used for detection on flat terrain. It was determined that the combination of the two curves gave a better picture of the orientation of the heavily contaminated area. When detecting undulating terrain, the valleys will lead to the phenomenon of "multiple peaks" of the apparent amplitude frequency curve in a heavily polluted area. (3) In the symmetrical quadrupole sounding device, the placement of the bathymetric survey points at the center of the edges of the contaminated field projections during soundings provides a better estimation of the depths at the top and bottom of the heavily contaminated areas of the field based on the apparent amplitude frequency curves. Consistent with the mid-staircase profile, an increase in the migration depth of the contaminated field causes a significant decrease in the peak apparent amplitude frequency. However, it has no significant effect on the apparent resistivity valley. (4) Influenced by the raised mountain peaks, the arrangement of power supply electrodes on the mountain frontiers on both sides of them allows for more efficient and precise monitoring and judgment of heavily polluted areas. Institutional Review Board Statement: Not applicable. This study did not involve humans or animals, so this statement does not apply to this study.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data in this article is publicly shareable. All relevant data are within the paper. All the data in this paper were obtained by numerical simulation through COMSOL, and then ORIGIN was used for chart drawing. (1) The setting of simulation parameters of pollution field has been described in Section 2. Figures 1 and 2 are obtained by numerical simulation through COMSOL. (2) The simulation parameter Settings in Sections 3 and 4 have been described in the paper, where the values of m, τ, and c are as follows: m = 0.03, τ = 1 s, c = 0.25. After parameter setting, the data obtained can be calculated according to the formula in the paper, and the chart in the paper can be obtained.