A rare earth mine area ecological risk assessment method and system based on rare earth element migration prediction

Abstract

The application discloses a rare earth mine area ecological risk assessment method and system based on rare earth element migration prediction, and belongs to the field of environmental risk assessment. The method comprises the following steps: acquiring environmental characteristic data sets of soil profiles at different depths in a target area; performing interface coupling effect calculation on the environmental characteristic data sets to generate a multi-medium migration flux data set, wherein the interface coupling effect calculation is used to determine interface characteristic parameters evolving over time according to third characteristic data, and to determine dynamic source and sink items based on the interface characteristic parameters, and the multi-medium migration flux data set indicates the migration flux of different phase states of rare earth at each depth evolving over time; and outputting a three-dimensional risk map according to the multi-medium migration flux data set, wherein the spatial dimension of the three-dimensional risk map indicates the migration and diffusion range of the rare earth, and the chemical dimension of the three-dimensional risk map indicates the bioavailability of the rare earth. The application improves the pertinence and accuracy of risk assessment.

Description

Technical Field

[0001] This invention relates to the field of environmental risk assessment, and in particular to a method and system for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction. Background Technology

[0002] During the mining process of rare earth resources, especially in-situ leaching operations using ion adsorption, rare earth elements migrate and diffuse downstream into the soil and groundwater, posing a continuous threat to the ecological environment and human health of the mining area and its surroundings. Therefore, predicting the migration behavior of rare earth elements in the soil profile of the mining area and assessing the ecological and health risks accordingly is of significant practical importance. Existing methods for assessing rare earth migration have the following technical shortcomings:

[0003] 1) Existing methods neglect the metabolic activities of microbial communities in mining soils, or treat the interfacial electrical parameters of mineral colloids as fixed constants that do not change over time. 2) Existing methods often track rare earth solutes in soil as a single dissolved state, ignoring the two key implicit migration pathways of colloidal transport and biological complexation, and failing to describe microbially driven phase transformation reactions. 3) Existing methods typically use only the total concentration of rare earth elements as the basis for risk assessment, but a large portion of the total rare earth element concentration in soil exists in an inert, non-migratory, and non-bioabsorbable form; furthermore, light rare earth elements and heavy rare earth elements have fundamental differences in toxicology, with heavy rare earth elements being far more toxic than light rare earth elements; additionally, plant roots have an active regulatory capacity for rare earth element absorption. 4) Pure mechanistic prediction models struggle to balance physical conservation and quantitative accuracy in the field. Frequent wet and dry cycles, preferential flow caused by localized rainfall, and nonlinear succession of microbial communities leave significant systematic residuals between mechanistic predictions and actual monitoring values. Furthermore, under extreme hydrological events such as torrential rains or tailings dam leaks, colloids may be released in large quantities and suddenly, which conventional models often severely underestimate. Summary of the Invention

[0004] One of the objectives of this invention is to provide a method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction, in order to solve the problems in the existing technology that statically treat interface parameters and ignore the dynamic disturbance of microbial metabolism, track solutes only in a single dissolved state and cannot characterize multi-phase multi-channel migration, and that pure mechanism models cannot take into account both physical conservation and field quantitative accuracy.

[0005] This invention is achieved through the following technical solution: a method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction, comprising the following steps: acquiring environmental feature datasets of soil profiles at different depths, wherein the environmental feature datasets are collected in layers at different depths and include at least: first feature data for characterizing the occurrence state and flow driving conditions of water in the soil porous media, second feature data for characterizing the occurrence concentration of rare earth elements in different phases, and third feature data for characterizing the metabolic activities of microbial communities in the soil; performing interface coupling calculation on the environmental feature datasets to generate a multi-media migration flux dataset, wherein the interface coupling calculation is used to determine interface feature parameters that evolve over time based on the third feature data, and to determine dynamic source and sink terms based on the interface feature parameters, and the multi-media migration flux dataset indicates the migration flux of rare earth elements in different phases at each depth over time; outputting a three-dimensional risk map based on the multi-media migration flux dataset; the spatial dimension of the three-dimensional risk map indicates the migration and diffusion range of rare earth elements, and the chemical dimension indicates the bioavailability of rare earth elements.

[0006] Furthermore, the environmental feature dataset of soil profiles at different depths in the target area is obtained by: acquiring the environmental feature dataset through at least one of in-situ sensor arrays deployed at different depths, soil profile sampling and testing, and microbial molecular biology detection, and uploading it to the electronic device via an Internet of Things gateway.

[0007] Furthermore, the different phases of rare earth represented by the second characteristic data include at least one of free rare earth, colloidal rare earth and bio-organic complexed rare earth. The free rare earth is rare earth existing in the pore liquid in the form of hydrated cations, the colloidal rare earth is rare earth adsorbed or co-precipitated on mineral colloids, and the bio-organic complexed rare earth is rare earth coordinated with biological macromolecules.

[0008] Furthermore, the third characteristic data includes at least one of a first metabolic parameter and a second metabolic parameter, wherein the first metabolic parameter indicates the mass concentration of extracellular polymers secreted by microorganisms per unit volume of pore fluid, and the second metabolic parameter indicates the molar concentration of low molecular weight organic acids in the pore fluid.

[0009] Furthermore, the interface characteristic parameters include at least one of the surface potential and surface charge density of the mineral colloid; the interface characteristic parameters that evolve over time include: establishing a mapping relationship between the third characteristic data and the interface characteristic parameters, and converting the interface characteristic parameters from a fixed input into a state variable that evolves over time according to the mapping relationship.

[0010] Furthermore, the mapping relationship includes a first mapping relationship for the surface potential, which characterizes that: the higher the first metabolic parameter, the stronger the shielding effect on the outer electric double layer of the mineral colloid, the lower the absolute value of the surface potential, and the shielding effect tends to saturate as the first metabolic parameter increases; and the negative contribution of the second metabolic parameter to the surface potential due to deprotonation in the pore fluid, wherein the negative contribution increases as the second metabolic parameter increases and is enhanced as the pH of the environment increases relative to the degree of dissociation of the organic acid.

[0011] Furthermore, the mapping relationship also includes a second mapping relationship for the surface charge density. The second mapping relationship is used to characterize that the adsorption amount of the first metabolic parameter on the surface of the mineral colloid tends to saturate as the first metabolic parameter increases, and the contribution of the adsorption of the first metabolic parameter to the surface charge density is determined based on the adsorption amount and the effective charge valence state of the functional groups carried by the first metabolic parameter.

[0012] Further, determining the dynamic source-sink term based on the interface feature parameters includes: determining the modification interaction potential energy based on the interface feature parameters, wherein the modification interaction potential energy is used to characterize the net interaction tendency between the mineral colloid carrying rare earth elements as it approaches the inner wall of the soil particle pores, and the interaction tends to change with the separation distance; extracting the potential energy feature quantity from the modification interaction potential energy, and determining the deposition rate parameter and release rate parameter based on the preset release threshold determination rule; and using the deposition rate parameter and release rate parameter as part of the dynamic source-sink term.

[0013] Furthermore, the modified interaction potential energy is the sum of the first interaction component, the second interaction component, and the third interaction component; the first interaction component characterizes the van der Waals attraction between the mineral colloid and the soil particles; the second interaction component characterizes the double-layer repulsion formed by the like charges on their surfaces, and the second interaction component evolves with the surface potential output from the mapping relationship; the third interaction component characterizes the steric repulsion generated by the polymer layer adsorbed on the surface of the mineral colloid when it is compressed.

[0014] Furthermore, the first active component includes a delay correction term for the separation distance, which is used to reduce the determined van der Waals attraction as the separation distance increases.

[0015] Furthermore, the second active component is modified according to the system ionic strength; the greater the system ionic strength, the shorter the effective distance of the determined double-layer repulsion.

[0016] Furthermore, the third active component is activated only when the separation distance is less than twice the effective extended thickness of the adsorbed polymer layer; when activated, the third active component increases with the degree of compression of the separation distance relative to twice the thickness.

[0017] Furthermore, the potential energy characteristic quantities include: the height of the barrier formed by the repulsive effect on the modification interaction potential energy, and the first minimum depth formed on the side adjacent to the inner wall of the pore; the release threshold determination rule includes: the higher the barrier height, the smaller the determined deposition rate parameter; the greater the energy difference between the barrier height and the first minimum depth, the smaller the determined release rate parameter.

[0018] Furthermore, determining the dynamic source-sink term based on the interface feature parameters also includes: determining phase transformation components through independent and parallel first reaction channels, second reaction channels, and third reaction channels, and then coupling the determined phase transformation components, after performing mass conservation equations, as at least a part of the dynamic source-sink term to the reaction transport equation set; the reaction transport equation set is used to track the concentration fields of free rare earth, colloidal rare earth, and bio-organic complexed rare earth at each depth over time.

[0019] Furthermore, the first reaction channel is used to read the concentration of competing metal cations in the environmental feature data set, and to determine the displacement and desorption amount of the target rare earth element on the adsorption site by the competing metal cations based on the multi-component competitive adsorption relationship. The higher the concentration of the competing metal cations, the smaller the net adsorption amount of the target rare earth element determined.

[0020] Furthermore, the second reaction channel is used to determine the net complexation amount of the conversion of free rare earth to bio-organic complexed rare earth according to the first metabolic parameter and the concentration of free rare earth, based on a reversible reaction relationship. The positive complexation amount increases with the product of the first metabolic parameter and the concentration of free rare earth, and the reverse dissociation amount increases with the concentration of bio-organic complexed rare earth.

[0021] Furthermore, the third reaction channel is used to read the reduction metabolism of iron-reducing bacteria in the third feature data, determine the structural solubility of the iron-rich mineral colloid, and determine the secondary release amount of rare earth elements bound to the iron-rich mineral colloid based on the structural solubility and the rare earth mass solidified per unit mass of iron-rich mineral colloid. The secondary release amount increases with the population density of iron-reducing bacteria and is inhibited with the increase of the system redox potential.

[0022] Furthermore, the competing metal cations include at least one of aluminum ions and iron ions.

[0023] Furthermore, in the reaction transport equation set, the equation term tracking colloidal rare earths includes not only the change in the concentration of colloidal rare earths in the liquid phase, but also the change in the amount of rare earths carried by the colloids deposited on the solid phase framework. Moreover, the deposition term and the release term in this equation term are driven by the deposition rate parameter and the release rate parameter as described above, respectively.

[0024] Furthermore, after generating the multi-media migration flux dataset, the method further includes error compensation processing of the multi-media migration flux dataset. The error compensation processing includes: deriving the baseline flux at the target evaluation section, where the baseline flux is the result of the sum of the convection and dispersion components of each phase of rare earth at the target evaluation section accumulated over time; concatenating the perturbation factor data with the intermediate state of the mechanism output from the interface coupling effect calculation to form a joint feature matrix; inputting the joint feature matrix into a pre-trained nonlinear compensation model to obtain flux compensation values; and superimposing the flux compensation values ​​onto the baseline flux to correct the multi-media migration flux dataset.

[0025] Furthermore, the disturbance factor data includes at least one of the following: microenvironment pH, redox potential, and wet-dry alternation frequency.

[0026] Furthermore, the nonlinear compensation model includes a decision tree ensemble model, which is used to capture the nonlinear residual distribution between the perturbation factor data and the actual rare earth flux through the ensemble of multiple decision trees.

[0027] Furthermore, the nonlinear compensation model includes a support vector regression kernel function model, which is used to assign greater weights to the support vectors obtained through offline training when the current joint feature matrix is ​​more similar to each support vector, thereby estimating the current prediction residual and outputting the flux compensation value.

[0028] Furthermore, the output of a three-dimensional risk map includes: separating light rare earth components and heavy rare earth components from the multi-media migration flux dataset; determining a separation coefficient based on the selective adsorption degree of light and heavy rare earth components by natural mineral colloids and bio-organic composite colloids, wherein the separation coefficient indicates the degree of difference in the distribution preference of light and heavy rare earth components between the colloidal and free phases; constructing an exposure level based on the conversion flux of specific rare earth forms between activation and fixation, replacing the total amount of rare earth as the risk criterion, based on the separation coefficient; and outputting the three-dimensional risk map based on the exposure level.

[0029] Furthermore, when constructing the exposure level, a mobility penalty coefficient is assigned to the heavy rare earth components adsorbed on the bio-organic composite colloid, which is higher than that assigned to the light rare earth components adsorbed on the natural mineral colloid, in order to characterize the hidden penetration risk of the heavy rare earth components with the colloid in deep diving.

[0030] Furthermore, the light rare earth component is a rare earth element component with a relatively large ionic radius and relatively weak inner-layer complexing ability, while the heavy rare earth component is a rare earth element component with a relatively small ionic radius and relatively high charge density.

[0031] Furthermore, the three-dimensional risk map is output based on the exposure level, including: combining the exposure level, the spatial distribution of the population in the target area, and the food chain transmission relationship to convert the bioavailability of rare earth elements into an average daily intake; mapping the probability distribution of the average daily intake exceeding the safety threshold to a geographic information system to generate the three-dimensional risk map, wherein the planar coordinates of the three-dimensional risk map indicate the migration and diffusion range on the land surface, and the depth coordinates indicate the bioavailability risk at the vertical profile.

[0032] Furthermore, the determination of food chain transmission relationships includes: based on the rhizosphere regulatory parameters of ground cover crops in the target area, regulating the translocation ratio of free rare earth elements and small molecule complexed rare earth elements entering the root system, so as to convert soil liquid flux into organism uptake; wherein the rhizosphere regulatory parameters include at least one of root cation exchange capacity and rhizosphere acidification coefficient, and the rhizosphere acidification coefficient is used to characterize the extent to which root exudates cause a decrease in rhizosphere pH relative to the bulk soil.

[0033] Furthermore, when calculating the daily intake, a higher toxicity weighting is applied to the heavy rare earth components than to the light rare earth components.

[0034] Furthermore, when calculating the daily intake, a higher toxicity weighting is applied to the heavy rare earth components than to the light rare earth components.

[0035] Furthermore, the standard probability distribution is determined by the log-normal cumulative distribution function based on the assumption that the average daily intake follows a log-normal distribution.

[0036] Furthermore, the ecological risk assessment method for rare earth mining areas also includes shear correction processing:

[0037] When the rainfall infiltration intensity in the first feature data exceeds the pore water saturation critical value and triggers an extreme hydrological event, the shear correction process is performed. The shear correction process includes: determining the equivalent shear amount applied to the interface between the colloid and soil particles based on the Darcy flow velocity and volumetric water content in the first feature data; when the equivalent shear amount exceeds the critical stripping amount, determining the release amplification factor based on the degree to which the equivalent shear amount exceeds the critical stripping amount, and amplifying the release rate parameter with the release amplification factor to obtain an adaptive release rate parameter; feeding the adaptive release rate parameter back to the equation term for tracking colloidal rare earth elements, and re-determining the multi-media migration flux dataset.

[0038] Furthermore, the critical peeling amount is a decreasing function of the wet-dry alternation frequency; the higher the wet-dry alternation frequency, the lower the critical peeling amount. Moreover, the release amplification coefficient decreases with the increase of the biopolymer damping amount, which characterizes the extracellular polymer network's resistance to shear peeling.

[0039] Furthermore, the release rate parameter is amplified by the release amplification factor in a segmented manner: the release amplification factor is used only when the equivalent shear amount exceeds the critical peeling amount; otherwise, the release rate parameter remains unchanged.

[0040] Another aspect of the present invention provides an ecological risk assessment system for rare earth mining areas based on rare earth element migration prediction, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements any of the rare earth mining area ecological risk assessment methods based on rare earth element migration prediction as described above.

[0041] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0042] 1. This invention incorporates microbial metabolic characteristic parameters into an environmental characteristic dataset and establishes a mapping relationship between them and interface characteristic parameters. This transforms the interface electrical parameters, which were originally considered fixed constants, into state variables that evolve over time, thereby truly reflecting the dynamic perturbation of colloidal migration behavior by microbial activity. This overcomes the prediction inaccuracies caused by staticizing interface parameters. Simultaneously, by constructing multiple independent and parallel sub-reaction channels as dynamic source-sink terms and coupling these terms to a set of reaction transport equations for tracking free, colloidal, and bio-organic complexed rare earth elements, the invention synchronously tracks the cross-media evolution trajectory of rare earth elements in different phases, thus overcoming the systematic bias of tracking only the dissolved state.

[0043] 2. This invention constructs a bioavailability exposure level based on activation-fixed conversion flux by separating light and heavy rare earth components and applying differentiated mobility and toxicity weighting. This results in a three-dimensional risk map that combines spatial diffusion and chemobioavailability dimensions. The final assessment results can accurately reflect the portion of rare earth that can be absorbed and utilized by organisms, thereby significantly improving the targeting and accuracy of risk assessment and providing a high-precision spatial decision-making basis for targeted governance of mining areas.

[0044] 3. This invention compensates for nonlinear errors through a mechanism-data collaborative approach. It locks in fundamental laws such as mass conservation through physical mechanisms and learns field disturbance residuals using a nonlinear compensation model. This avoids both the idealization bias of pure mechanism models and the physical unconstrained drift of pure data models. Combined with adaptive shear correction logic for extreme hydrological events, it achieves physical reduction of the explosive release of colloidal rare earth elements in sudden scenarios such as rainstorms, improving robustness and quantitative accuracy in real field environments. Attached Figure Description

[0045] The accompanying drawings, which are included to provide a further understanding of embodiments of the invention and form part of this application, do not constitute a limitation thereof. In the drawings:

[0046] Figure 1 This is a flowchart of the method provided in Embodiment 1 of the present invention.

[0047] Figure 2 The diagram shows the dynamic evolution of the zeta potential of mineral colloids provided in Embodiment 1 of the present invention.

[0048] Figure 3 This is a Langmuir saturation curve provided in Embodiment 1 of the present invention.

[0049] Figure 4 The total interaction potential energy curve of MM-DLVO provided in Embodiment 1 of the present invention is shown.

[0050] Figure 5 The diagram shows the relationship between the barrier height and the deposition / release rate constant provided in Embodiment 1 of the present invention.

[0051] Figure 6 This is a schematic diagram of the exponential regulation of bioreduction release source and sink terms provided in Embodiment 1 of the present invention.

[0052] Figure 7 This is a concentration distribution map of rare earth elements along a soil depth profile provided in Embodiment 1 of the present invention.

[0053] Figure 8 This is a schematic diagram of the inhibition effect provided in Embodiment 1 of the present invention.

[0054] Figure 9 This is a comparison chart of cumulative flux prediction and actual measurement provided in Embodiment 1 of the present invention.

[0055] Figure 10 This is a comparison chart of the predicted residual distribution provided in Embodiment 1 of the present invention.

[0056] Figure 11 This is a diagram showing the evolution of interface differentiation coefficients over time and the identification of hazardous areas, provided in Embodiment 1 of the present invention.

[0057] Figure 12 This is a schematic diagram of the nonlinear enhancement effect provided in Embodiment 1 of the present invention.

[0058] Figure 13 The three-dimensional ecological and health risk map provided in Embodiment 1 of the present invention.

[0059] Figure 14 This is a schematic diagram of a colloidal rare earth burst flux pulse provided in Embodiment 1 of the present invention. Detailed Implementation

[0060] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0061] Example 1

[0062] This embodiment discloses a method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction. Figure 1 The flowchart of the method in this embodiment is shown. As can be seen from the flowchart, this embodiment includes the following steps:

[0063] Step 1: Obtain environmental feature datasets of soil profiles at different depths in the target rare earth mining area.

[0064] Among them, the target rare earth mining area refers to a specific geographical area where rare earth mining (especially ion adsorption type rare earth in-situ leaching operations) is currently underway or has been carried out, and where there is a risk of rare earth elements migrating and spreading to downstream soil and groundwater.

[0065] The environmental feature dataset refers to a set of multi-source physicochemical and biological physical quantities that influence rare earth migration in the longitudinal profile of the soil in the mining area. This dataset is collected in layers at different depths (e.g., several discrete sampling sections near the surface, subsurface, and water table), including:

[0066] The environmental feature dataset includes, but is not limited to, physical and hydrological parameters, basic data on rare earth element concentrations, and microbial metabolic feature parameters.

[0067] Physical hydrological parameters refer to physical quantities used to describe the state of water content and flow driving conditions in soil porous media, including but not limited to soil volumetric water content, Darcy velocity, saturated permeability coefficient, soil dry bulk density, critical pore water saturation value, and rainfall infiltration intensity. These parameters determine the macroscopic driving force and boundary conditions for the convection-dispersion migration of solutes in porous media.

[0068] The basic data on rare earth element concentrations refer to the initial occurrence concentration information of various rare earth elements in the target profile, including but not limited to the concentration of free rare earth elements existing in the pore liquid as hydrated cations, the concentration of colloidal rare earth elements adsorbed or co-precipitated on mineral colloids, and the concentration of bio-organic complexed rare earth elements coordinated with biomacromolecules. They can be further distinguished according to light rare earth elements (LREEs, such as lanthanum and cerium) and heavy rare earth elements (HREEs, such as yttrium and ytterbium).

[0069] Microbial metabolic characteristic parameters refer to physical quantities used to characterize the intensity of metabolic activity and the level of metabolites of microbial communities (e.g., *Thiobacillus ferrooxidans*, iron-reducing bacteria, etc.) in mining soil. In this embodiment, microbial metabolic characteristic parameters include at least the secretion of extracellular polymeric substances (EPS) and the concentration of low molecular weight organic acids, and may further include the population density of iron-reducing bacteria, the ratio of polysaccharide / protein components of specific chain lengths in EPS, and the reduction metabolic rate constant of iron-reducing bacteria, etc.

[0070] Among them, the extracellular polymeric substance (EPS) secretion amount refers to the mass concentration of extracellular polymeric substances secreted by microorganisms in a unit volume of pore liquid. The charged functional groups such as carboxyl and phosphate groups carried by EPS can not only change the interfacial electrical properties of mineral colloids, but also directly coordinate with rare earth ions to form stable bio-organic colloids.

[0071] The concentration of low molecular weight organic acids refers to the molar concentration of small molecule organic acids such as oxalic acid and citric acid in pore fluid. Their ability to undergo deprotonation reactions in soil pore fluid and contribute additional negative charges to the system is regulated by the relative relationship between their dissociation constant and the ambient pH.

[0072] Understandably, traditional rare earth migration assessment methods often completely ignore the aforementioned microbial metabolic parameters, or treat only the interfacial electrical parameters of mineral colloids as fixed constants that do not change over time. This is the fundamental reason for their inaccurate predictions in bioactive mining areas. The reason why this embodiment incorporates microbial metabolic parameters into the environmental feature dataset is to use "bio-disturbance," a key driving factor that is persistent in real mining areas but has long been overlooked, as the boundary input for subsequent dynamic calculations.

[0073] In this embodiment, the environmental feature dataset can be obtained by a combination of methods, such as in-situ sensor arrays deployed at different depths in the mining area, soil profile sampling and testing, and microbial molecular biology detection, and then uploaded to the electronic device where this system is located via an Internet of Things gateway.

[0074] Step 2: Perform quaternary coupling interface interaction calculations on the environmental feature dataset to generate a multi-media rare earth migration flux dataset.

[0075] Among them, the quaternary coupling interface interaction calculation refers to the process of simultaneously solving the deposition and release behavior of rare earth element-carrying colloids in soil pores, as well as the cross-media evolution of three phases of rare earth elements in porous media, including free state, colloidal state, and bio-organic complexed state, by comprehensively considering the coupling interaction between the four components: mineral colloids, pore fluid, soil particles, and microbial metabolites.

[0076] The multi-media rare earth migration flux dataset refers to the set of mass fluxes obtained through the above calculations, describing the evolution of different phases of rare earth elements at various depth sections of the target profile over time. It serves as the direct data basis for subsequent ecological and health risk assessment. In this embodiment, the calculation of the quaternary coupling interface effect specifically includes two interconnected core steps:

[0077] First, based on microbial metabolic characteristic parameters, the total microscopic interaction potential energy of the interface between microbially modified mineral colloids and soil particles is calculated, and this is used to quantify the critical deposition rate and release rate constant of rare earth element-carrying colloids in soil pores. Second, the critical deposition rate and release rate constant are used as dynamic source and sink terms and input into a preset reactive solute transport master equation to obtain a multi-media rare earth migration flux dataset. The following will provide a detailed explanation of the two core steps.

[0078] I. Calculation of the total microscopic interaction potential energy at the interface between mineral colloids and soil particles under microbial modification.

[0079] In this embodiment, in order to ensure that the potential energy calculation can accurately reflect the dynamic perturbation of colloidal migration behavior by microbial activity, thereby overcoming the prediction inaccuracies caused by the traditional DLVO model treating interface parameters as static constants, the calculation specifically includes the following sub-steps:

[0080] Sub-step 2.1: Establish a deterministic mapping matrix between the concentration of microbial metabolites and the characteristic parameters of the colloidal interface.

[0081] The mapping matrix is ​​used to establish a deterministic mapping relationship between EPS secretion, low molecular weight organic acid concentration, and the change in zeta potential and surface functional group density (surface charge density) of mineral colloids, thus making the zeta potential, which was originally considered a fixed constant, more stable. With surface charge density Instead of fixed inputs, they are transformed into state variables that evolve over time.

[0082] The physical design of this mapping relationship is as follows: After EPS macromolecules are adsorbed onto the colloidal surface, their own neutral or weakly charged functional groups form a shielding layer around the particles, compressing the effective thickness of the electric double layer and thus reducing the absolute value of the detectable surface potential. This shielding effect is significantly concentration-dependent; the higher the EPS concentration, the stronger the shielding. Simultaneously, low-molecular-weight organic acids are deprotonated in the pore liquid, contributing additional negative charges to the system. The amount of this contribution depends on the concentration of the organic acid and the relative pH of the environment to the dissociation constant of the organic acid. The level of pH, when the environmental pH is higher than During this period, organic acids undergo significant deprotonation, with negative charges contributing the most.

[0083] For example, in this embodiment, the Zeta potential of the modified mineral colloid surface can be calculated using the following formula:

[0084]

[0085] in, for Zeta potential (mV) on the surface of the microbially modified mineral colloid. The initial Zeta potential (mV) of natural mineral colloids when undisturbed by microorganisms. for Mass concentration of EPS in pore fluid at any given time (mg / L); The molar concentration (mol / L) of low molecular weight organic acids. The value is the EPS charge shielding coefficient (L / mg). The larger the value, the stronger the compression effect of the double layer per unit EPS concentration. Contribution factor to deprotonation of organic acids (mV·L / mol); denoted as the dissociation constant of the corresponding organic acid. Figure 2 The dynamic evolution curve of the Zeta potential of the microbially modified mineral colloid in this embodiment is shown.

[0086] It is understandable that the first term in the above formula adopts an exponential decay form, reflecting the nonlinear saturation characteristics of EPS coating, that is, when the EPS concentration is extremely low, the potential is approximately equal to the initial value. When the EPS concentration is sufficiently high, the exponential term approaches zero, indicating that the shielding effect tends to saturate and the potential no longer decreases with concentration. The second term takes a linear form, approximately describing the linear cumulative contribution of organic acids to surface charge within a specific pH range. The superposition of these two terms unifies the two opposing but simultaneous physical mechanisms of shielding de-positioning and deprotonation into a single evolutionary equation.

[0087] After determining the dynamic evolution of the Zeta potential, it is necessary to further track the surface charge density of the colloidal particles. The physical basis for this change lies in the fact that EPS adsorption onto the colloidal surface not only alters the measurable potential but also directly injects additional charge into the surface through its own charged functional groups such as carboxyl and phosphate groups. This injection process exhibits typical isothermal adsorption characteristics. The number of adsorption sites on the colloidal surface is limited, and as the EPS concentration increases, the adsorption gradually approaches saturation. Therefore, it can be described using the Langmuir isotherm. Multiplying this by the effective charge valence state of the EPS functional groups yields the net contribution of EPS adsorption to the surface charge density.

[0088] For example, in this embodiment, the surface charge density of the modified colloid can be calculated using the following formula:

[0089]

[0090] in, The surface charge density of the modified colloid (C / m²) is the surface charge density of the colloid. The initial charge density is (C / m²). This represents the maximum adsorption capacity of EPS on the colloidal surface (mol / m²), indicating the upper limit of the surface adsorption site capacity. The value represents the Langmuir adsorption equilibrium constant (L / mg), which reflects the affinity between EPS and the colloidal surface. A higher value indicates that a higher adsorption coverage can be achieved at a lower concentration. The effective charge valence state of the EPS functional groups determines the surface charge increment per unit molar adsorption amount. Figure 3 The Langmuir saturation curve of surface charge density caused by EPS adsorption in this embodiment is shown.

[0091] It is understandable that the Langmuir fractional term in the above formula is... When it approaches zero, it approaches zero. Much larger The time approaches one, which precisely reflects the physical picture of the adsorption sites being gradually filled until saturation; through the two mapping equations established in sub-step 2.1, the static interface feature parameters in the traditional DLVO model are endowed with the ability to evolve in real time with microbial activities, providing dynamic and realistic boundary conditions that reflect biological disturbances for subsequent potential energy calculations.

[0092] Sub-step 2.2: Update the double-layer repulsion energy calculation formula in the DLVO theoretical model using the change in Zeta potential output by the mapping matrix, and combine it with the dynamic ion intensity to correct the van der Waals attraction energy calculation formula, thereby introducing steric repulsion energy and constructing the microbial modified DLVO (MM-DLVO) framework.

[0093] Among them, the total microscopic interaction potential is used to answer whether the net force between a mineral colloidal particle coated with EPS and the inner wall of the soil medium pores is attraction or repulsion; the former causes the colloid to be deposited on the pore wall and captured, while the latter causes the colloid to continue to migrate through the matrix with the pore fluid.

[0094] Understandably, traditional DLVO theory only considers van der Waals attraction and double-layer repulsion. However, when the colloidal surface is coated with biopolymers such as EPS, the polymer chains extend outward on the particle surface, forming an elastic polymer brush structure. When the distance between the colloid and the pore walls decreases to less than twice the extension length of the polymer chains, the polymer layers undergo spatial overlap and compression, generating strong steric repulsion. Traditional DLVO theory, precisely because it neglects this factor, severely overestimates the sedimentation tendency of colloids in systems coated with EPS. Therefore, this embodiment introduces a third term into the total potential energy: steric repulsion energy. .

[0095] In this embodiment, the colloid and soil particulate medium are connected by thermal energy. Total interaction potential energy in units It is determined by three parts: the van der Waals attraction energy, the electric double-layer repulsion energy, and the steric hindrance repulsion energy.

[0096]

[0097] in, This represents the separation distance (nm) between the colloidal surface and the soil medium surface. The physical meaning and calculation logic of the three partial energies are explained below.

[0098] For Van der Waals attraction There exists a van der Waals dispersion force between any two objects due to instantaneous dipole fluctuations, the magnitude of which can be given by the Hamaker method. It should be noted that when the interparticle spacing is large (typically exceeding approximately 5 nm), the finite propagation speed of the electromagnetic field leads to a weaker actual attractive force than classically predicted; this is the so-called "retardation effect." This embodiment introduces a correction factor proposed by Gregory (1981). The delay effect is approximated in a single formula, where the characteristic attenuation wavelength is... Typically, 100nm is used.

[0099] For example, in this embodiment, the van der Waals attraction energy can be calculated using the following formula:

[0100]

[0101] in, is the equivalent hydrodynamic radius (nm) of the mineral colloid. The effective Hamaker constant (J) for the three-component system of colloid (1)-water (2)-soil medium (3) is determined by the sign of the van der Waals force, which is usually positive in common mineral-water-mineral systems, i.e., it is attractive. The characteristic attenuation wavelength (nm) is given. It is understandable that the negative sign before the formula indicates... It is always negative (attraction energy) and increases with separation distance. It decreases and then rapidly increases.

[0102] For double-layer repulsion energy When two surfaces carrying the same charge approach each other, their surrounding diffuse double layers overlap, causing an increase in osmotic pressure and a decrease in the entropy of the counterion configuration, thus producing a repulsive effect. Under moderate surface potential conditions, this repulsive potential can be described by a logarithmic approximation.

[0103] For example, in this embodiment, the double-layer repulsion energy can be calculated using the following formula:

[0104]

[0105] in, is the dielectric constant of the medium (F / m); , These are the surface potentials (mV) of the inner walls of the pores in mineral colloids and porous soil media, respectively. The former is provided in real time by the dynamic mapping equation of S2.1, which is a concrete implementation of the design intention of "updating the double-layer repulsion energy by using the change in Zeta potential output by the mapping matrix". The Debye sieve parameters (m⁻¹) are related to the system ionic strength. Proportional It is understandable that the higher the concentration of electrolytes such as ammonium salts in the eluent, the greater the ionic strength, the thinner the electric double layer is compressed, and the shorter the range of the repulsive potential energy. This is the physical essence of why colloids are more likely to deposit under high ionic strength conditions, which is also the inherent meaning of combining dynamic ionic strength correction in this step.

[0106] For spatial steric repulsion energy This item only applies when the distance between the two surfaces is less than twice the elongation length of the EPS polymer chain (i.e., The steric hindrance term is activated only when the distance between the polymer layers on both sides is greater than this critical value. When the distance between the layers is greater than this critical value, the polymer layers on both sides have not yet come into contact, and the steric hindrance term is zero. Once the compression range is entered, the elastic repulsion energy is described by a parabolic approximation that is proportional to the square of the compression amount. This is a simplified form of the Alexander-deGennes model under the weak compression limit.

[0107] For example, in this embodiment, the steric repulsion energy can be calculated using the following formula:

[0108]

[0109] in, The effective extension thickness (nm) of the EPS polymer layer adsorbed on the colloidal surface directly determines the "range" of the steric hindrance effect. The number of polymer chains adsorbed per unit area (m⁻²). The coverage of the polymer on the colloidal surface (dimensionless, value from 0 to 1). and The product of these factors together determines the "stiffness" of steric repulsion. The product of the Boltzmann constant and the thermodynamic temperature, as the natural unit of energy, reflects the essential contribution of thermal fluctuations to the degrees of freedom of polymer chain configuration. It is understandable that the values ​​within the parentheses... Item in When the repulsive force approaches zero, it approaches one (the repulsive force is strongest). Approaching When the time approaches zero (the steric hindrance effect just disappears), the entire transition process is continuous and smooth.

[0110] By By incorporating the total potential energy, the MM-DLVO framework can accurately predict that even during periods of vigorous microbial growth and high EPS secretion, colloidal particles can still maintain high stability and avoid sedimentation by relying on the elastic protection of the "polymer brush," thereby achieving long-distance migration. This physical essence, which cannot be reflected in traditional DLVO models, can be quantitatively characterized in this embodiment.

[0111] Sub-step 2.3: The modified double-layer repulsion energy is superimposed with the van der Waals attraction energy, and the steric hindrance repulsion energy is also superimposed to generate the total interaction potential energy curve of colloidal-soil particles that reflects the microbial modification effect.

[0112] The total interaction potential energy curve refers to the curve with the separation distance as the reference. The x-axis represents the total potential energy. The curve is plotted with the vertical axis as the ordinate. This curve is obtained by superimposing the three components from sub-step 2.2 at each separation distance, and fully depicts the energy topography experienced by the colloid as it approaches the pore wall.

[0113] It is understandable that, since the double-layer repulsion energy and steric hindrance repulsion energy in sub-step 2.2 both evolve in real time with the metabolic characteristic parameters of microorganisms, the total interaction potential energy curve is not a static curve, but a dynamic curve that changes synchronously with the intensity of microbial activity in the mining area. Figure 4 The graph shows the total interaction potential energy curve of the MM-DLVO with the superposition of the three component energies in this embodiment.

[0114] Sub-step 2.4: Extract the barrier height and the first minimum depth from the total interaction potential energy curve, and calculate the critical deposition rate and release rate constant according to the preset colloid kinetic release threshold determination rule.

[0115] Among them, the barrier height (denoted as) The first minimum depth (denoted as ) refers to the highest peak value of the repulsive energy on the total interaction potential energy curve. The colloid can only approach the pore walls and deposit after crossing this barrier. (Usually a negative value) refers to the energy well bottom formed by the potential energy curve on the side immediately adjacent to the pore wall, where the deposited colloid is trapped.

[0116] Critical deposition rate and release rate constant refer to two kinetic rate constants that can be directly applied to the macroscopic transport equations by converting the aforementioned continuous potential energy curve at the microscale. These constants represent the deposition rate constants of colloids from the liquid phase to the solid phase. (s⁻¹), and the release rate constant of the deposited colloid from the solid phase back to the liquid phase. (s⁻¹). The physical basis for this transformation is Kramers transition state theory—analogizing the process of colloids crossing a potential barrier as a thermally activated transition of particles on a potential energy surface.

[0117] For the deposition rate constant The process of colloids diffusing from a bulk solution to a solid surface and overcoming a potential barrier can be described as a flux problem driven by Brownian diffusion and modulated by the potential energy curve. Its reciprocal (i.e., the resistance encountered during deposition) is proportional to the effective transit time of the colloid under the influence of the potential energy curve, and can be obtained by integrating along the distance direction with respect to the Boltzmann weights.

[0118] For example, in this embodiment, the colloidal deposition rate constant can be calculated using the following formula:

[0119]

[0120]

[0121] in, The Brownian diffusion coefficient (m² / s) of the colloid can be obtained from the Stokes-Einstein equation. The given value reflects the inherent migration capacity of colloidal thermal motion; the integral term in the denominator is the potential energy-weighted spatial average drag. It can be understood that when the potential barrier... The higher the value of the integrand, the larger its exponential peak at the potential barrier, and the larger its integral value. The smaller the value, the better—this aligns perfectly with the physical intuition that "the higher the potential barrier, the more difficult it is for the colloid to deposit." It should be noted that this integral is typically solved numerically in practical calculations, and its value is directly updated in real-time by S2.3. The curve provides that the potential energy that evolves with microbial activity can be reflected in the deposition rate in real time.

[0122] For the release rate constant For a colloid deposited in the first minimum potential well to escape back to the liquid phase, it must overcome the energy difference between the bottom of the potential well and the top of the potential barrier. According to the Kramers / Arrhenius thermal activation theory, the escape rate is proportional to the characteristic vibrational frequency and is suppressed by an exponential Boltzmann factor.

[0123] For example, in this embodiment, the colloidal release rate constant can be calculated using the following formula:

[0124]

[0125] in, The characteristic vibrational frequency (s⁻¹) of the colloid in the potential well represents the frequency at which the colloid attempts to escape due to thermal motion at the bottom of the potential well. The height of the first potential barrier on the potential energy curve (in terms of...) (in units) The depth of the first minimum value (in) (The unit is usually negative). It's understandable that the values ​​within the parentheses... If the value is negative, its absolute value is the effective activation energy barrier; the higher the barrier, the greater the activation energy barrier. The smaller the value, the more difficult it is for the deposited colloids to spontaneously desorb, and the system exhibits irreversible deposition; conversely, when EPS modification raises the steric hindrance repulsion barrier and simultaneously shallows the first minimum value, The colloids significantly increase in size, transforming into a reversible adsorption mode. This is precisely the kinetic mechanism by which microbial activity promotes colloid migration, and also the physical meaning of the colloid kinetic release threshold determination rule in sub-step 2.4. Figure 5 The diagram illustrates the relationship between the barrier height and the deposition / release rate constant in this embodiment.

[0126] At this point, the data channel from microscopic colloidal physics (potential energy curves) to macroscopic transport equations (kinetic rate constants) has been fully established. and This parameter will be used as a parameter that evolves dynamically over time and will be fed into the following master equation governing reactive solute transport in real time.

[0127] 2. The critical deposition rate and release rate constant are used as dynamic source and sink terms input into the master equation governing reactive solute transport.

[0128] In this embodiment, to realistically depict the multi-form, multi-pathway transport and dynamic phase transformation of rare earth elements in the soil profile of the mining area from a mechanistic perspective, the main governing equation for reactive solute transport is extended to a three-phase multi-channel coupled partial differential equation system. Dynamic source and sink terms are constructed through calculations of three independent and parallel sub-reaction channels. The following will first explain the calculation logic of the three sub-reaction channels, and then explain the construction of the coupled transport equation system:

[0129] (a) Calculation of three independent and parallel sub-reaction channels.

[0130] Understandably, traditional convection-dispersion equations treat soil solutes as a single dissolved state, completely ignoring the two key implicit pathways of colloidal transport and biocomplexation, and failing to describe microbially driven phase transformation reactions. This leads to systematic prediction biases in the highly complex biogeochemical environment of rare earth mining areas. Therefore, this embodiment sets up the following three independent and parallel sub-reaction pathways:

[0131] 1) Competitive adsorption channel: This channel reads aluminum ions from the input water chemistry data in real time. With iron ions The concentration was dynamically calculated based on a surface complexation model to determine the displacement and desorption rate of the target rare earth element at the adsorption sites on the colloidal surface. The physical background is that in industrial scenarios involving in-situ leaching of rare earth minerals, the injection of the leaching solution leads to the accumulation of a large amount of [unclear - possibly a specific substance or process] in the pore fluid. and These two trivalent cations compete fiercely with rare earth elements for adsorption sites on mineral surfaces, significantly reducing the effective adsorption capacity of rare earth elements. Since the simple single-component Langmuir equation cannot describe this substitution relationship, this embodiment employs a multi-component competitive Langmuir equation. This equation incorporates the competition relationship between the three cations for adsorption sites into the same denominator, and quantitatively describes the degree of competition by weighting the values ​​using their respective affinity constants.

[0132] For example, in this embodiment, competing adsorption source and sink terms It can be calculated using the following formula:

[0133]

[0134] in, Soil dry bulk density (g / cm³); The maximum adsorption capacity (mg / g) of rare earth elements on porous media surfaces is limited by soil mineral type and specific surface area. , , The Langmuir affinity constants (L / mg) for rare earth elements, aluminum ions, and iron ions are respectively, and their relative magnitudes determine their competitive positions. The concentration of free rare earth ions (mg / L); , This represents the real-time concentration (mol / L) of free competing metal cations in the environment. It's understandable that this term takes the total derivative with respect to time because the adsorption amount dynamically changes with the environmental concentration; only its transient change represents the true source-sink intensity in the convection-dispersion equation. In the denominator... and The larger the value, the smaller the net adsorption capacity of rare earth elements, meaning... Decrease (strength of adsorption sink), free state The proportion retained in the liquid phase actually increases, which is the quantitative mechanism by which competing ions promote rare earth migration.

[0135] 2) Bio-complexation channel: This channel calculates the mass percentage of stable bio-organic colloids formed by direct coordination between EPS and rare earth elements, based on EPS secretion levels and the polysaccharide / protein ratio of specific chain lengths. The physical background is that the carboxyl and phosphate functional groups abundant in EPS have a strong coordination affinity for trivalent rare earth ions, capable of converting free rare earth elements into soluble EPS-rare earth complexes. This process can be described by second-order reversible reaction kinetics: the forward complexation rate is proportional to the product of EPS concentration and free rare earth concentration, while the reverse dissociation rate is proportional to the complex concentration; the difference between the two is the net complexation flux.

[0136] Exemplarily, in this embodiment, biological complex source-sink term It can be calculated using the following formula:

[0137]

[0138] in, The soil volumetric water content (dimensionless) indicates that the reaction occurs in the liquid phase pore liquid and the reaction rate needs to be converted to the mass flux per unit volume of soil. The positive rate constant for biological complexation (L / (mg·s)) reflects the coordination reaction activity between EPS functional groups and rare earth ions. The rate constant for complex dissociation (s⁻¹) reflects the thermodynamic stability of EPS-rare earth complexes. The smaller the value, the more stable the complex and the more difficult it is to dissociate. This refers to the concentration of rare earth elements in biologically organically complexed forms (mg / L). It is understandable that when... Greater than hour, A positive value indicates a net transfer of rare earth elements from a free state to a biologically organically complexed state.

[0139] 3) Biodegradation and Reduction Release Channel: This channel reads the reduction metabolic rate constant of iron-reducing bacteria, calculates the structural dissolution rate of iron-rich mineral colloids, and then derives the secondary release rate of rare earth elements bound to these iron-rich mineral colloids. The physical background is that in an anaerobic microenvironment (low redox potential)... (Lower or negative values) Iron-reducing bacteria can use ferric iron in iron-rich mineral colloids as a terminal electron acceptor for anaerobic respiration, reducing it to soluble ferrous iron, thereby causing the iron oxides that form the colloidal framework to dissolve. Simultaneously, rare earth elements originally fixed on the colloid are released back into the liquid phase, forming a strong internal source of free rare earth elements. The rate of this process is controlled by three factors: the density of iron-reducing bacteria, the amount of available iron-rich colloidal substrate, and the redox potential.

[0140] Exemplarily, in this embodiment, the bioreduction release source and sink... It can be calculated using the following formula:

[0141]

[0142]

[0143] in, is the specific activity reaction rate constant of iron-reducing bacteria (g / (cells·s)); for The density of iron-reducing bacteria in the soil at any given time (cells / g) directly reflects the first-order driving effect of biomass on the reaction rate. The mass (mg / g) of reducible iron-rich mineral colloids in the soil framework represents the upper limit of substrate availability. Real-time redox potential (mV); is the redox sensitivity coefficient (mV⁻¹); This represents the percentage of solidified rare earth elements (mg / g) in a unit mass of iron-rich colloid, used as a stoichiometric factor to convert the dissolved amount of colloid into the released amount of rare earth elements. It is understandable that when... When the value is positive (under oxidation conditions), the exponential term When the value is much smaller than one, the activity of iron-reducing bacteria is strongly inhibited; when When the value is reduced to negative (strong reducing conditions), the exponential term increases significantly, and the activity of iron-reducing bacteria is fully released, which is consistent with the biological fact that iron-reducing bacteria can only carry out dissimilar iron reduction under anaerobic conditions. Figure 6 The redox potential for bioreduction release source and sink term is shown in this embodiment. A diagram illustrating index regulation.

[0144] It should be noted that the above three sub-reaction channels are independent of each other and calculated in parallel. They are based on three mechanistic pathways: competitive adsorption, biological complexation, and reduction-release, and constitute the three source components of the dynamic source-sink term.

[0145] (ii) Generation of dynamic source and sink terms and construction of coupled transport control equations.

[0146] After obtaining the calculation results of the three sub-reaction channels, this embodiment further performs mass conservation equations to construct a time-varying set of third-order partial differential terms, which are then coupled to the main control equation for reactive solute transport.

[0147] Specifically, this set of coupled transport governing equations tracks the spatiotemporal evolution of three phases of rare earth elements: the concentration of free rare earth ions... (mg / L) represents rare earth elements existing in the pore liquid as hydrated cations that can directly participate in adsorption-desorption reactions; concentration of colloidal rare earth complexes. (mg / L) represents the rare earth elements adsorbed or co-precipitated on suspended mineral colloids and migrating with the colloids; the concentration of rare earth elements in bio-organic complexed form. (mg / L) represents a stable complex formed with biomacromolecules such as EPS through coordination bonds.

[0148] Channel 1 is used to track the spatiotemporal evolution of free rare earth ions: its governing equations follow a standard convection-dispersion framework and include three source-sink terms, namely the adsorption flux to the soil solid phase. (Negative sign, representing the sink), the reaction flux of EPS complexing to generate organic complexes. (Negative sign, representing the pool) and the rare earth flux released back into the liquid phase when iron-reducing bacteria decompose iron-rich colloids. (Use the positive sign as the source):

[0149]

[0150] ,

[0151] Channel 2 is used to track the spatiotemporal evolution of suspended colloidal rare earth elements: the left side of this equation includes not only the time derivative of the concentration of the colloidal state in the liquid phase, but also the time derivative of the amount of rare earth elements carried by the colloids deposited on the solid framework, to reflect the "two-phase" characteristics of the colloids; the colloidal deposition and release terms on the right side are the rate constants obtained from the dynamic calculation in S2.4 mentioned above. , Directly driven:

[0152]

[0153] ,

[0154] Channel 3 is used to track the spatiotemporal evolution of bio-organic complexed rare earth elements: their formation depends entirely on the chemical chelation of free rare earth elements by EPS, source item The corresponding sink in Channel 1 is a mirror image of the channel, thus ensuring the overall quality conservation of the system.

[0155] ,

[0156] in, Soil volumetric water content (dimensionless). Soil dry bulk density (g / cm³); Darcy velocity (m / s) represents the macroscopic convection driving force of the pore fluid; , , The values ​​represent the hydrodynamic dispersion coefficients (m² / s) for the three phases, which combine the contributions of molecular diffusion and mechanical dispersion. The amount of colloids deposited per unit mass on the surface of the soil skeleton (mg / g). The rare earth content (mg / mg) adsorbed on the deposited colloid. These represent the adsorption source and sink terms of colloidal rare earth elements on the solid phase.

[0157] Understandably, the three equations described above constitute a set of one-dimensional finite element reactive solute transport master equations. The displacement-desorption rate, bio-organic colloid mass ratio, and secondary release rate calculated from the three sub-reaction channels, after being combined with the equations based on mass conservation, serve as a time-varying set of third-order partial differential source-sink terms, which are coupled to the hydrodynamic convection-dispersion control node of this set of equations. Since the concentration change in any channel will affect other channels through the source-sink terms, the system can simultaneously track the cross-media evolution trajectory of free rare earth ions and colloidal rare earths with different carrier types.

[0158] By solving the above set of equations simultaneously, the concentration field of each phase of rare earth at each depth section of the target profile as a function of time can be obtained, and then a multi-media rare earth migration flux dataset can be synthesized.

[0159] III. Mechanism - Nonlinear error compensation under the data collaborative driving mechanism.

[0160] In this embodiment, after obtaining the multi-medium rare earth migration flux dataset, this step may also include a step of nonlinear error compensation using a mechanism-data collaborative driving mechanism.

[0161] It is understandable that even if the aforementioned physical mechanism equations are structurally sophisticated, the multiple random disturbances objectively present in the field mining environment—such as frequent alternations between wet and dry conditions leading to sudden changes in water content, activation of preferential flow channels due to localized rainfall, and nonlinear succession of microbial communities at different times—will still leave a significant systematic residual between the mechanistic predictions and actual monitoring values. Directly ignoring these residuals will result in a continuous deviation in the output cumulative flux prediction; conversely, incorporating all uncertainties into the mechanistic equations will lead to severe overdeterminism of the model parameters, resulting in unacceptable computational costs. Therefore, this embodiment adopts a mechanism-data collaborative driving strategy: the physical mechanism equations are used to lock in fundamental physical laws such as mass conservation and reaction stoichiometry to ensure the accuracy of the qualitative trend and magnitude range of the predictions; simultaneously, a deterministic nonlinear kernel regression model is introduced to specifically learn and compensate for the residuals of the mechanistic model under nonlinear disturbances in the field.

[0162] It should be noted that the nonlinear kernel regression model can be implemented by decision tree ensemble algorithms such as random forest, or by the support vector regression (SVR) kernel function in this embodiment. The two are completely consistent in their design idea of ​​capturing the nonlinear residual distribution law between multiple environmental disturbance factors and the actual rare earth flux through machine learning. This embodiment uses the SVR kernel function as an example for explanation.

[0163] Nonlinear error compensation can specifically include the following logic:

[0164] First, the mechanistic intermediate state feature vector output at the end of the integration step of the master equation governing reactive solute transport is extracted, and the target evaluation section (depth is...) is derived. The baseline theoretical cumulative flux at point () is given. At any given time, the rare earth mass flux passing through this section is determined by the sum of the convection and dispersion terms of each of the three channels. Integrating this sum over time yields the baseline theoretical cumulative flux.

[0165] For example, in this embodiment, the baseline theoretical cumulative flux at the target section can be calculated using the following formula:

[0166]

[0167] ,

[0168] In this context, the two terms within each set of parentheses represent the convective flux (velocity-driven mass transport) and diffusion flux (concentration gradient-driven diffusion transport) of the corresponding phase, respectively; subscripts... This indicates that only the value at the target depth section is taken; the external integration over time accumulates the instantaneous flux density into the total mass flux. It is understandable that... This represents the theoretically required mass of rare earth elements to traverse the cross section under idealized and precisely parameterized conditions, serving as a benchmark for subsequent residual compensation.

[0169] Secondly, multiple environmental disturbance factors measured in the experiment, including microenvironmental pH and redox potential, will be considered. The frequency of wet-dry alternation is concatenated with the eigenvectors of intermediate states of the mechanism to form a joint feature matrix. The design idea is that the prediction bias of the mechanism model is not random noise, but rather there is a complex nonlinear dependency between the current state of the system (pore liquid concentration, water content, pH, Eh) and historical environmental events (frequency of wet-dry alternation). If this relationship can be found, it is possible to predict in advance whether the mechanism model will be overestimated or underestimated under the current system state, and how much the bias will be, and then make corrections accordingly.

[0170] For example, in this embodiment, the joint feature matrix (high-dimensional feature vector) used for residual prediction can be defined as:

[0171]

[0172] ,

[0173] in, , The free and colloidal rare earth concentrations output by the mechanism model at the target section are called the "mechanism intermediate state feature vector". This represents the measured moisture content. , As a real-time geochemical state indicator; For the front The historical wet-dry cycle frequency within each assessment period is used to capture the cumulative effect of historical events on the current residuals. This is the transpose symbol for a matrix.

[0174] Next, the joint feature matrix is ​​input into a pre-trained nonlinear kernel compensation model. A kernel function ensemble algorithm is used to capture the nonlinear residual distribution between multiple environmental disturbance factors and the actual rare earth flux, outputting the corresponding flux compensation residual value. This embodiment employs a deterministic SVR model with a Gaussian radial basis function kernel (RBF-Kernel). Based on historical 'feature vector-measured residual' sample pairs accumulated in the offline training set, the aforementioned nonlinear mapping relationship is learned, thereby enabling online estimation of the predicted residual at the current moment.

[0175] For example, in this embodiment, the residual compensation amount predicted by the SVR kernel function can be calculated using the following formula:

[0176]

[0177] in, and For each of the Lagrange multipliers (all non-negative), the sign of the difference between them determines the sign of the first multiplier. Whether the residual contributed by each support vector is positive or negative compensation is determined by offline training on historical observation residual data; For the first Support vectors are the feature vector sample points selected in the training set to represent the kernel function expansion. The width parameter (m⁻²) of the RBF kernel controls the radiation radius of each support vector's influence. The larger the value, the sharper the kernel function, the more refined the model's fit to local states, but the weaker its generalization ability. The bias constant is determined through training; The total number of support vectors is determined by the amount of training data and the regularization parameter. It is understandable that the RBF kernel function in the current state... With support vectors The larger the similarity (the smaller the Euclidean distance), the larger the value. This reflects the interpolation logic of residual similarity under similar historical states, which is the mathematical basis for kernel methods to learn nonlinear manifolds.

[0178] Finally, the flux compensation residual value is superimposed on the theoretical baseline flux calculated by the master equation of reactive solute transport to correct the final multi-media rare earth migration flux dataset under extreme mining conditions.

[0179] For example, in this embodiment, the final output after dual-drive correction can be calculated using the following formula:

[0180]

[0181] Understandable, The physical mechanism model provides a qualitatively accurate prediction benchmark, ensuring that the prediction results do not violate basic physical laws such as the conservation of mass. The residual patterns of nonlinear perturbations in the field are learned online by a nonlinear kernel model, ensuring the quantitative accuracy of the prediction results in real field environments. The two are interdependent and each performs its own function, avoiding both the idealization bias of pure mechanistic models and the unconstrained drift in physics of pure data models. Figure 7 The concentration distribution of rare earth elements in three different phases along the soil depth profile in this embodiment is shown. Figure 8 This diagram illustrates the inhibitory effect of multi-component competitive adsorption on the effective adsorption amount of rare earth elements in this embodiment. Figure 9 This embodiment shows a comparison of cumulative flux prediction and actual measurement using both mechanism-data dual-driven and pure mechanism approaches. Figure 10 The diagram shows a comparison of the prediction residual distribution before and after the mechanism-data dual-drive compensation in this embodiment.

[0182] Step 3: Based on the multi-media rare earth migration flux dataset, analyze the multi-scale morphological distribution characteristics of rare earth elements during the migration process, and output a three-dimensional ecological and health risk map that includes the spatial dimension of migration and diffusion radius and the chemical dimension of biological effectiveness.

[0183] In this embodiment, to overcome the distortion of risk assessment caused by traditional methods that only use the total concentration of rare earth elements as the basis for risk assessment, ignore the essential differences in toxicology between light and heavy rare earth elements (heavy rare earth elements, such as yttrium and ytterbium, have much stronger neurotoxicity than light rare earth elements such as lanthanum and cerium), and also ignore the active regulatory ability of plant roots to absorb rare earth elements, this step may include the following sub-steps:

[0184] Sub-step 3.1: Separate the light rare earth (LREE) components and heavy rare earth (HREE) components from the multi-media rare earth migration flux dataset.

[0185] Among them, the light rare earth (LREEs) component refers to the rare earth element components represented by lanthanum and cerium, which have relatively large ionic radii and relatively weak inner-layer complexing ability; the heavy rare earth (HREEs) component refers to the rare earth element components represented by yttrium and ytterbium, which have gradually decreased ionic radii due to the lanthanide contraction effect and have high charge density.

[0186] It is understandable that separating light and heavy rare earth elements is a prerequisite for subsequently characterizing their differentiation during migration and applying differentiated toxicity weights.

[0187] Sub-step 3.2: Based on the thermodynamic constants of selective adsorption of light and heavy rare earth elements by natural mineral colloids and bio-organic composite colloids, calculate the separation coefficients of light and heavy rare earth elements at different media interfaces.

[0188] Among them, the interface differentiation coefficient (denoted as ) It is used to quantitatively characterize the degree of difference in the partition preference of light rare earth elements and heavy rare earth elements between colloidal and free phases.

[0189] The physical basis of this differentiation phenomenon lies in the fact that the ionic radii of light rare earth elements (LREEs) and heavy rare earth elements gradually decrease due to the 4f electron orbital contraction effect (lanthanide contraction), leading to a differentiation in their coordination chemistry. Heavy rare earth elements, due to their higher charge density, have a stronger affinity for oxygen-containing ligands in EPS and for the surfaces of colloidal minerals, and tend to accumulate in the colloidal phase; while light rare earth elements are relatively more retained in the free liquid phase. In this embodiment, the distribution ratio of light rare earth elements between the colloidal and free phases, divided by the distribution ratio of heavy rare earth elements between the same two phases, is defined as the interfacial differentiation coefficient.

[0190] For example, in this embodiment, the interface differentiation coefficient can be calculated using the following formula:

[0191]

[0192] in, and These represent the sum of the concentrations (mg / L) of each light rare earth element in the colloidal phase (including the sum of mineral colloids and bio-organic colloids) and the free dissolved phase, respectively; the corresponding denominators are... and This corresponds to the amounts of heavy rare earth elements. It is understandable that when... When the value is greater than one, it indicates that light rare earth elements are more enriched in the colloidal phase than heavy rare earth elements, meaning that colloids have a selective carrying advantage for light rare earth elements; when When the value is less than one, it indicates that heavy rare earth elements tend to bind to colloids, and colloidal migration will preferentially carry heavy rare earth elements. It will serve as the core weighting factor in subsequent calculations of plant availability concentration. Figure 11 The evolution of the interface differentiation coefficient over time and the danger zone identification diagram in this embodiment are shown.

[0193] The logic setting for the differentiation coefficient in this embodiment is as follows: Since light rare earth elements have a larger ionic radius and weaker inner-layer complexing ability than heavy rare earth elements, when the system iterates through the bioavailability exposure level in step S3, it assigns a higher mobility weighted penalty coefficient to heavy rare earth elements adsorbed on organic EPS media than to light rare earth elements adsorbed on inorganic kaolinite colloids. This accurately characterizes the hidden penetration risk of heavy rare earth elements in deep diving. That is, although heavy rare earth elements tend to bind to colloids, they may migrate long distances with pore fluid through colloids protected by polymer brushes, thus forming a hidden risk in deep diving that is not easily detected by conventional total indicators.

[0194] Sub-step 3.3: Based on the differentiation coefficient, the total rare earth content assessment index in the soil is abandoned, and instead a bioavailability exposure level based on the activation-fixation conversion flux of specific rare earth forms is constructed.

[0195] Among them, the bioavailability exposure level refers to the index system that abandons the traditional approach of using the total amount of rare earth elements in the soil as the risk criterion, and instead uses the conversion flux of specific forms of rare earth elements between activation (release to a free state or small molecule complex state that can be absorbed by organisms) and fixation (adsorption, deposition, stable complexation) as the basis to classify the exposure risk of different spatial locations.

[0196] It is understandable that a large portion of the total rare earth elements in the soil exists in an inert, non-migratory, and non-bioavailable form. Including all of these in the risk assessment would seriously overestimate the actual harm. This embodiment uses activation-fixation conversion flux as the core criterion and combines it with the mobility-weighted penalty for heavy rare earth elements in sub-step 3.2 to ensure that the exposure level truly reflects the portion of rare earth elements that can be absorbed and utilized by organisms, thereby significantly improving the relevance and accuracy of the risk assessment.

[0197] Sub-step 3.4: Combining the exposure level with the population spatial distribution and food chain enrichment and transmission attenuation equation of the target area, the bioavailability flux is converted into the average daily intake of rare earth elements in a specific valence state in humans, and the probability distribution of the average daily intake exceeding the safety threshold is mapped to the geographic information system to generate a three-dimensional ecological and health risk map.

[0198] The food chain enrichment and transport attenuation equation describes the attenuation law of rare earth elements from the soil liquid phase through crop roots, from absorption and enrichment to final transfer to the human body. In this embodiment, the execution conditions of the equation include: loading the rhizosphere microenvironment acidification coefficient of a specific crop based on the differences in root cation exchange capacity (CEC) characteristics of the ground cover crop varieties in the target mining area; and using the acidification coefficient to regulate the transport ratio of free rare earth ions and small-molecule organic complexed rare earth elements entering the root cortical cells, thereby achieving a deterministic conversion from soil liquid phase flux to organismal intake.

[0199] It is understandable that the concentration of rare earth elements (REEs) in the soil liquid phase is not equivalent to the effective dose actually absorbed by plants. Crop roots influence REE absorption through two active regulatory mechanisms: first, the root cation exchange capacity determines the maximum adsorption capacity of REE cations on the root surface; second, roots secrete organic acids to actively acidify the rhizosphere microenvironment (lowering the pH relative to the bulk soil pH by a certain margin), making the rhizosphere pH lower than the bulk pH, thereby increasing the solubility and availability of REE cations. Simultaneously, based on the aforementioned toxicological differences, when summarizing the effective doses of light and heavy REEs, a higher toxicity-weighted penalty must be applied to heavy REEs.

[0200] For example, in this embodiment, the bioavailable equivalent concentration that can be absorbed by crops into the food chain can be calculated using the following formula:

[0201]

[0202]

[0203]

[0204] in, The value represents the bioavailable equivalent concentration (mg / kg) that can be absorbed by crops and enter the food chain; the first term in square brackets is the rhizosphere regulation coefficient, where... It provides the basic absorption capacity (cmol / kg). The exponential term describes the enhancing effect of rhizosphere acidification on bioavailability. The lower the effective pH, the larger the exponential term and the higher the bioavailability. This relationship reflects the nonlinear characteristics of the effect of pH on ion activity through an exponential function. The root acidity sensitivity coefficient reflects the intensity of a specific crop species' response to rhizosphere acidification. The pH of the soil itself, This represents the extent of rhizosphere acidification. and These are the toxicological weighting coefficients for light and heavy rare earth elements, respectively, reflecting the assessment logic that heavy rare earth elements are more toxic. An example can be set as follows: This involves multiplying the bioeffective dose of heavy rare earth elements by 4.5 times and converting it into an equivalent hazard dose of light rare earth elements (it should be noted that this is an exemplary setting; in actual applications, calibration must be performed based on the specific toxicological data of the heavy rare earth elements). It is understandable that... The introduction of this ensures that the differentiation effect of colloidal carrying bias towards heavy rare earth elements in sub-step 3.2 is correctly amplified and included in the risk assessment. Figure 12 A schematic diagram illustrating the nonlinear enhancement effect of rhizosphere acidification on bioavailability in this embodiment is shown.

[0205] After obtaining the bioavailable equivalent concentration, it is further converted into the actual human intake dose. By using the biotransfer coefficient of rare earth elements from the soil to the edible parts of crops, as well as the average daily dietary intake and body weight of the population, the average daily intake (ADI) of contaminated agricultural products consumed by surrounding residents through diet can be calculated.

[0206] For example, in this embodiment, the average daily intake of a human body can be calculated using the following formula:

[0207]

[0208] in, The biotransfer coefficient (dimensionless) of rare earth elements from soil to edible parts of crops reflects the differences in the enrichment capacity of different crop species for rare earth elements. The target population's daily crop intake (kg / d); Average human body weight (kg); The unit is mg / (kg·d), which is the daily intake per unit body weight.

[0209] In obtaining After considering the time series, Influenced by various random factors such as spatial heterogeneity of soil, batch differences in agricultural products, and individual differences in population intake, its statistical distribution often approximates a log-normal distribution. This embodiment will... Input a log-normal cumulative distribution function and assess whether it exceeds the national security critical threshold. (Exemplary example) The probability of exceeding the standard is calculated, and a three-dimensional topological risk map is formed in space. , The coordinates represent the planar migration and diffusion radius on the Earth's surface. Axial depth represents the bioavailability risk level at the vertical profile below ground.

[0210] For example, in this embodiment, the probability of exceeding the three-dimensional limit in space can be calculated by the following formula:

[0211]

[0212] in, It is the standard normal cumulative distribution function; and They are respectively The mean and standard deviation can be estimated from historical monitoring data or Monte Carlo sampling. It is understandable that when a certain spatial point... The mean is exactly equal to hour, The value equals 0.5, meaning the probability of exceeding the limit is 50%; when When significantly higher than the threshold, When the value approaches one, it indicates that the position is almost certainly exceeding the limit; when When significantly below the threshold, A value approaching zero indicates that the location is safe. The generated 3D map simultaneously presents the spatial footprint (planar distribution) of pollutant diffusion and the biotoxicity depth of specific valence states (longitudinal profile), changing the past practice of crudely determining risk based solely on one-dimensional physical boundaries or two-dimensional planar total concentration, and providing a high-precision spatial decision-making basis for identifying priority areas for targeted governance in mining areas. Figure 13 A three-dimensional ecological and health risk map is shown in this embodiment.

[0213] This step also includes adaptive correction logic for extreme hydrological event scenarios:

[0214] The aforementioned steps, under normal hydrodynamic conditions, can fully describe the migration and risk evolution of rare earth elements. However, in extreme hydrological events such as torrential rains or tailings dam leaks, the pore fluid velocity can increase by several orders of magnitude in a short period. The resulting high-speed fluid shear force is sufficient to forcibly strip away colloidal particles already deposited and anchored on the pore walls, triggering a sudden, large-scale release of colloids, manifested as a steep peak in the cumulative flux. The colloidal release rate constant calibrated under normal conditions... Under these conditions, the release intensity will be severely underestimated, leading to inaccurate risk prediction. Soil that has experienced multiple cycles of drying and wetting weakens the physical bond between colloidal particles and the matrix due to repeated shrinkage-expansion cycles (microcracks develop on the particle surface), and its critical peeling stress also decreases, making it more sensitive to rainstorm events. Therefore, this embodiment further includes adaptive correction logic for extreme hydrological event scenarios:

[0215] When the monitoring system detects that rainfall infiltration intensity in the hydrological feature dataset exceeds the pore water saturation threshold, leading to extreme wet-dry cycles, a strong shear hydrodynamic response branch is triggered. Within this branch, the pore water velocity gradient variable is forcibly injected into the calculation of the total interaction potential energy curve under microbial modification. This recalibrates the EPS hydration film rupture threshold caused by extreme hydrological scouring, thereby correcting the sudden flux peak caused by the explosive release of colloidal rare earth elements during heavy rainfall scenarios. This branch is driven by real-time fluid shear force, modulated by historical wet-dry cycles, and resisted by EPS biomass cohesion, reconstructing the colloidal release rate under extreme events in real time. Its specific logic is as follows:

[0216] First, the macroscopic Darcy velocity in the porous medium is converted into an estimate of the equivalent wall shear force applied to the colloidal-matrix interface. At the Darcy scale, the average velocity within the pores can be approximated by dividing the Darcy velocity by the volumetric water content. Combined with hydrodynamic viscosity, the equivalent shear strength of the pore walls can then be estimated.

[0217] For example, in this embodiment, the pore-scale equivalent wall shear force can be estimated by the following formula:

[0218]

[0219] in, Equivalent wall shear force (Pa); The dynamic viscosity of water is (Pa·s). Real-time Darcy velocity (m / s); Real-time volumetric moisture content (dimensionless). Let be the saturated permeability coefficient of the medium (m / s). It should be noted that this expression is an engineering approximation formula used to qualitatively capture the amplification effect of extreme flow velocity events on shear force levels at the Darcy scale, rather than a rigorously derived exact solution at the pore scale from the Navier-Stokes equations. It has reasonable applicability in engineering early warning applications.

[0220] Secondly, based on the calculated shear force, it is further determined whether this shear force exceeds the critical threshold that the colloid can resist in its current state, and the release amplification factor corresponding to the degree of exceedance is quantified. Among these factors, the critical peel stress is... It is the historical frequency of alternation between dry and wet seasons. The decreasing function, the more times it undergoes alternating wet and dry conditions, the more... The lower the temperature, the weaker the physical anchoring ability of the soil to the colloids; at the same time, the EPS network provides additional elastic cohesion to the colloids due to its viscoelastic characteristics. It acts as a buffer layer to resist shearing and stripping. The larger the value, the smaller the amplification of the release corresponding to the same shear force exceedance.

[0221] For example, in this embodiment, the shear force damping correction factor can be calculated using the following formula:

[0222]

[0223] in, The shear force-induced release rate amplification factor (dimensionless); molecule This is the excess amount (Pa) of shear force relative to the critical threshold. Corrections are only meaningful when this value is positive (i.e., after the shear force exceeds the critical value); denominator The elastic damping cohesion (Pa) provided to the EPS network serves to normalize the excess and modulate the amplification factor. The richer the EPS, the greater the excess under the same excess. The smaller the value, the more the amplification effect is suppressed by the cohesive force of EPS. This is understandable. It should be calibrated in the laboratory or in the field through dry-wet cycle colloidal peeling experiments to establish its functional relationship with the decreasing frequency of dry-wet alternation.

[0224] Furthermore, the aforementioned correction coefficients are embedded into the normal release rate constant in the form of a piecewise function to achieve an adaptive response to extreme hydrodynamic events. The correction is activated only when the shear force exceeds the critical threshold, thus avoiding the introduction of unnecessary disturbances at normal flow rates.

[0225] For example, in this embodiment, the adaptive colloid release rate under rainstorm dynamics correction can be calculated by the following formula:

[0226]

[0227] in, The base release rate constant (s⁻¹) is calculated under normal hydrodynamic conditions (i.e., S2.4). A value greater than 1 indicates the amplification factor of the release rate by extreme shear events, and its value increases with... The exponential growth of the increase reflects the nonlinear strengthening effect of high-speed fluid shear stripping.

[0228] at last, After calculation at each time step, the colloidal release term will be immediately fed back in real time to replace the colloidal release term in the aforementioned channel two transport equation, triggering an adjustment to the colloidal concentration field. The amount of solid deposition was recalculated, and then the system of equations with mass conservation was updated in channel one. The amount of driving substrate ultimately enables the entire three-channel equation system to respond synchronously to extreme hydrodynamic events, recalculate the flux mutation at the target cross section, and thus obtain a complete physical reconstruction of the evolution trajectory of the sudden increase in pollution flux under the scenario of a prolonged drought followed by a sudden rainstorm within the model framework. Figure 14 A schematic diagram of the explosive flux pulse of colloidal rare earth elements under an extreme rainstorm event is shown in this embodiment.

[0229] Example 2

[0230] This embodiment discloses an ecological risk assessment system for rare earth mining areas based on rare earth element migration prediction.

[0231] Specifically, this system can be integrated into electronic devices, such as terminals and servers. Terminals can be mobile phones, tablets, industrial field data acquisition gateways, laptops, or personal computers; servers can be single servers, server clusters consisting of multiple servers, or cloud computing platforms carrying IoT monitoring data for the mining area.

[0232] In this embodiment, the rare earth mining area ecological risk assessment system based on rare earth element migration prediction can also be integrated into multiple electronic devices. For example, the system can be integrated into multiple servers, with each server undertaking sub-tasks such as interface potential energy calculation, migration equation solving, error compensation, and risk map rendering, to collaboratively realize the rare earth mining area ecological risk assessment method based on rare earth element migration prediction in Embodiment 1 of this application.

[0233] In this embodiment, the server can also be implemented in the form of a terminal.

[0234] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction, characterized in that, The ecological risk assessment methods for rare earth mining areas include: Obtain environmental feature datasets from soil profiles at different depths. The environmental feature dataset is collected in layers at different depths and includes at least: First characteristic data used to characterize the occurrence state and flow driving conditions of water in soil porous media. Secondary characteristic data used to characterize the occurrence concentration of rare earth elements in different phases. Third characteristic data used to characterize the metabolic activities of microbial communities in soil; Interface coupling effects are calculated on the environmental feature dataset to generate a multi-media migration flux dataset, wherein, The interface coupling calculation is used to determine the interface feature parameters that evolve over time based on the third feature data. And based on the interface feature parameters, dynamic source and sink terms are determined. The multi-media migration flux dataset indicates the migration flux of different phase rare earth elements at each depth over time; Based on the multi-media migration flux dataset, output a three-dimensional risk map; The spatial dimension of the three-dimensional risk map indicates the migration and diffusion range of rare earth elements, while the chemical dimension indicates their bioavailability.

2. The method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction according to claim 1, characterized in that, The third feature data includes at least: At least one of the first metabolic parameter and the second metabolic parameter. The first metabolic parameter indicates the mass concentration of extracellular polymers secreted by microorganisms per unit volume of pore fluid. The second metabolic parameter indicates the molar concentration of low molecular weight organic acids in the pore fluid; The interface characteristic parameters include at least one of the surface potential and surface charge density of the mineral colloid; The interface feature parameters that evolve over time include: Establish a mapping relationship between the third feature data and the interface feature parameters, and based on the mapping relationship, convert the interface feature parameters from fixed inputs into state variables that evolve over time.

3. The method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction according to claim 2, characterized in that, The mapping relationship includes a first mapping relationship for the surface potential. The first mapping relationship is used to characterize that: the higher the first metabolic parameter, the stronger the shielding effect on the outer electric double layer of the mineral colloid, the lower the absolute value of the surface potential, and the shielding effect tends to saturate as the first metabolic parameter increases; And the negative contribution of the second metabolic parameter to the surface potential due to deprotonation in the pore fluid, the negative contribution increasing with the increase of the second metabolic parameter and increasing with the increase of the environmental pH relative to the degree of dissociation of the organic acid; The mapping relationship also includes a second mapping relationship for the surface charge density. The second mapping relationship is used to characterize that the adsorption amount of the first metabolic parameter on the surface of the mineral colloid tends to saturate as the first metabolic parameter increases, and the contribution of the adsorption of the first metabolic parameter to the surface charge density is determined based on the adsorption amount and the effective charge valence state of the functional groups carried by the first metabolic parameter.

4. The method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction according to claim 1, characterized in that, The determination of dynamic source and sink items based on the interface feature parameters includes: The modification interaction potential energy is determined based on the interface characteristic parameters. The modification interaction potential energy is used to characterize the net interaction tendency between the mineral colloid carrying rare earth elements as it approaches the inner wall of the soil particle pores, and how this interaction changes with the separation distance. Extract the potential energy characteristic quantity from the modified interaction potential energy, and determine the deposition rate parameter and release rate parameter from the potential energy characteristic quantity according to the preset release threshold determination rule; The deposition rate parameter and release rate parameter are included as part of the dynamic source-sink term; The modified interaction potential energy is the sum of the first action component, the second action component, and the third action component; The first active component characterizes the van der Waals attraction between the mineral colloid and the soil particles; The second action component characterizes the repulsive effect of the double layer formed by the same charge on the surfaces of the two objects, and the second action component evolves with the surface potential output by the mapping relationship; The third action component characterizes the steric repulsion effect generated when the polymer layer adsorbed on the surface of the mineral colloid is compressed.

5. The method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction as described in claim 1 or 4, characterized in that, The step of determining dynamic source and sink items based on the interface feature parameters further includes: Phase transition components were determined using the first, second, and third reaction channels, which are independent and parallel to each other. The determined phase transformation components are then combined with the mass conservation equations and used as at least a part of the dynamic source-sink terms, coupled to the reaction transport equations. The reaction transport equations are used to track the concentration fields of free rare earth elements, colloidal rare earth elements, and bio-organic complexed rare earth elements at each depth over time.

6. The method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction according to claim 5, characterized in that, The first reaction channel is used to read the concentration of competing metal cations in the environmental feature dataset, and to determine the displacement and desorption amount of the target rare earth element on the adsorption site by the competing metal cations based on the multi-component competitive adsorption relationship. The higher the concentration of the competing metal cations, the smaller the net adsorption amount of the target rare earth element determined. The second reaction channel is used to determine the net complexation amount of the conversion of free rare earth to bio-organic complexed rare earth according to the first metabolic parameter and the concentration of free rare earth, based on a reversible reaction relationship. The positive complexation amount increases with the product of the first metabolic parameter and the concentration of free rare earth, and the reverse dissociation amount increases with the concentration of bio-organic complexed rare earth. The third reaction channel is used to read the reduction metabolism of iron-reducing bacteria in the third feature data, determine the structural solubility of the iron-rich mineral colloid, and, based on the structural solubility and the mass of rare earth elements solidified per unit mass of iron-rich mineral colloid, determine the secondary release amount of rare earth elements bound to the iron-rich mineral colloid. The amount of secondary release increases with the population density of iron-reducing bacteria and is inhibited by the increase of the system's redox potential.

7. The method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction according to claim 1, characterized in that, After generating the multi-media migration flux dataset, the process also includes error compensation processing of the multi-media migration flux dataset. The error compensation process includes: The baseline flux at the target evaluation section is derived. The baseline flux is the result of the sum of the convective and diffuse components of each phase of rare earth at the target evaluation section accumulated over time. The intermediate state of the mechanism calculated by the perturbation factor data and the interface coupling effect is concatenated into a joint feature matrix; The joint feature matrix is ​​input into a pre-trained nonlinear compensation model to obtain the flux compensation value. The flux compensation value is superimposed on the baseline flux to correct the multi-media migration flux dataset.

8. The method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction according to claim 1, characterized in that, The output three-dimensional risk map includes: Light rare earth components and heavy rare earth components were separated from the multi-media migration flux dataset; Based on the selective adsorption degree of light rare earth components and heavy rare earth components by natural mineral colloids and bio-organic composite colloids, the separation coefficient is determined. The separation coefficient indicates the degree of difference in the distribution preference of light rare earth components and heavy rare earth components between the colloidal phase and the free phase. Based on the differentiation coefficient, an exposure level is constructed based on the conversion flux of specific rare earth elements between activation and fixation, to replace the total amount of rare earth elements as the risk criterion. Output the three-dimensional risk map based on the exposure level; When constructing the exposure level, a mobility penalty coefficient is assigned to the heavy rare earth components adsorbed on the bio-organic composite colloid, which is higher than that assigned to the light rare earth components adsorbed on the natural mineral colloid, in order to characterize the hidden penetration risk of the heavy rare earth components with the colloid in deep diving.

9. The method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction according to claim 1, characterized in that, The ecological risk assessment method for rare earth mining areas also includes shear correction processing: When the rainfall infiltration intensity in the first feature data exceeds the pore water saturation critical value and triggers an extreme hydrological event, the shear correction process is performed. The shearing correction process includes: Based on the Darcy flow rate and volumetric water content in the first feature data, the equivalent shear amount applied to the interface between the colloid and soil particles is determined. When the equivalent shear amount exceeds the critical peel amount, a release amplification factor is determined based on the degree to which the equivalent shear amount exceeds the critical peel amount, and the release rate parameter is amplified by the release amplification factor to obtain an adaptive release rate parameter. The adaptive release rate parameter is fed back into the equation for tracking colloidal rare earth elements to redetermine the multi-media migration flux dataset. The critical peeling amount is a decreasing function of the wet-dry alternation frequency; the higher the wet-dry alternation frequency, the lower the critical peeling amount. Furthermore, the release amplification factor decreases with the increase of the biopolymer damping amount, which characterizes the extracellular polymer network's resistance to shear peeling.

10. An ecological risk assessment system for rare earth mining areas based on rare earth element migration prediction, characterized in that, The rare earth mining area ecological risk assessment system includes: processor; The memory stores a computer program, which, when executed by a processor, implements the method for ecological risk assessment of rare earth mining areas based on rare earth element migration prediction as described in any one of claims 1 to 9.

Images