Multi-isotope based playa-lake groundwater pollution source apportionment system and method

Abstract

This invention provides a multi-isotope-based system and method for analyzing pollution sources in sandy lakes and groundwater, relating to the field of isotope mass spectrometry detection and tracing technology. The invention involves sub-region division, calculation of hydraulic parameters, and construction of a lake-groundwater hydraulic constraint matrix; collection of multiple water body samples to determine multi-isotope and hydrochemical indicators, constructing a fingerprint database; cross-validation of fluxes using the radon-222 mass balance method and the Darcy method, followed by sandy-adapted kinetic fractionation correction to retrieve original isotope data; identification of pollution sources and restoration of their initial isotopic composition using microbial fractionation correction; and finally, combining flux weight constraints, quantitatively solving for contribution rates using a Bayesian mixture model, and outputting the zonal analysis results after cross-validation. This invention achieves high-precision zonal quantitative analysis of pollution sources in sandy lakes and groundwater by constructing a zonal constraint matrix, coupling multi-isotope tracing, sandy-adapted fractionation correction technology, and constructing a Bayesian mixture model.

Description

Technical Field

[0001] This invention relates to the field of isotope mass spectrometry detection and tracing technology, specifically to a multi-isotope-based system and method for analyzing pollution sources in sandy lakes and groundwater. Background Technology

[0002] Sandy lake-groundwater systems are important water resources and ecological carriers in arid and semi-arid regions. Accurately identifying pollution sources is a core prerequisite for regional water environment governance and ecological protection. Among the commonly used pollution source apportionment methods, traditional hydrochemical tracing methods have low resolution. Single isotope tracing techniques are easily affected by overlapping end-member fingerprints of pollution sources, isotope fractionation due to strong evaporation in sandy areas, and fractionation effects during microbial transformation processes. Moreover, most apportionment schemes do not fully consider the spatial heterogeneity of bidirectional exchange between lakes and groundwater in the lacustrine zone and lack rigid constraints from hydrological and hydrodynamic conditions. This results in source apportionment results with strong ambiguity and insufficient quantitative accuracy, making it difficult to adapt to the complex hydrological characteristics of sandy lake-groundwater systems, which involve frequent seasonal hydraulic exchange and significant evaporation enrichment.

[0003] Current technologies for tracing water pollution mostly use single isotope tracing and traditional water chemistry methods. While these methods are quick and inexpensive, they are prone to misjudgment due to overlapping end-member fingerprints of pollution sources. Furthermore, most of these methods lack systematic isotope dynamic fractionation correction for sandy environments with strong evaporation. They also lack quantitative correction mechanisms for isotope fractionation caused by microbial transformation processes such as denitrification and sulfate reduction, resulting in distorted water isotope data that fails to accurately reflect the original characteristics of pollution sources. Consequently, the reliability of the analysis results is insufficient. Existing analysis schemes often use the entire lake as the calculation unit, without considering the spatial heterogeneity of the two-way exchange between lake and groundwater in the lakeshore zone on the composition of pollution sources. They also lack differentiated analysis of discharge zones, infiltration zones, and transition zones.

[0004] Therefore, it is necessary to provide a multi-isotope-based system and method for analyzing pollution sources in sandy lakes and groundwater to address the aforementioned problem.

[0005] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to provide a multi-isotope-based system and method for analyzing pollution sources in sandy lakes and groundwater, in order to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A multi-isotope-based source apportionment system for pollution in sandy lakes and groundwater includes the following steps: The hydrological zoning module, based on the hydrogeological data and dynamic water level monitoring data of the target area, and combined with the complementary relationship between the groundwater flow field and the lake water level, divides the area into three sub-regions, collects water level data of each sub-region, calculates the hydraulic gradient, exchange direction and interface flux, and constructs the lake-groundwater hydraulic constraint matrix. The sub-regions include the discharge zone, the infiltration zone and the transition zone. The sampling and testing module is used to collect lake water, groundwater, bottom sediment pore water and vadose zone soil solution samples from three types of sub-regions, determine the hydrogen and oxygen, nitrate nitrogen and oxygen, sulfate sulfur isotopes and radon-222 activity, simultaneously detect the concentration of anions and cations and the total dissolved solids content, and construct a multi-isotope-water chemical fingerprint database. The flux verification and inversion module establishes local atmospheric precipitation lines and regional lake evaporation lines based on hydrogen and oxygen isotope data. It calculates the groundwater inflow flux and lake water leakage in each region through the radon-222 mass balance equation, calls the constraint matrix to verify and correct the flux results, and uses the dynamic fractionation correction model to eliminate the strong evaporation isotope enrichment effect in sandy areas. It then retrieves the original isotope data of the lake water recharge source. The fractionation and reduction module is used to compare the corrected original isotope data with the characteristic values ​​of the end-member isotopes of the pollution source, identify the pollution sources of each zone by combining the water chemical ion ratio and land use information, correct the isotope fractionation caused by biotransformation through the isotope fractionation coefficient of microbial action, and restore the initial isotope composition of the pollution source. The model analysis module is used to take the initial isotope values ​​of the restored pollution sources as input, combine them with the flux weight constraints of the constraint matrix, and use a Bayesian mixture model to calculate the pollution contribution rate of each pollution source. After cross-validation, it outputs the spatial distribution and contribution ratio of the pollution sources.

[0008] Furthermore, the target area is divided into several sub-regions, and water level data for each sub-region is collected. These sub-regions include a discharge zone, an infiltration zone, and a transition zone. The method used is as follows: Hydrogeological data of the target area were retrieved, and the permeability coefficient of sediments in the lacustrine zone, the thickness of the vadose zone, and the characteristics of the regional groundwater flow field were extracted. Combined with long-term dynamic water level monitoring data of a complete hydrological year, the dynamic difference and seasonal fluctuation range between the lake water level and the groundwater level at the corresponding points in the surrounding area were calculated. If the groundwater level is higher than the lake water level throughout the year and the hydraulic gradient points towards the lake, it is designated as the discharge zone where groundwater continuously replenishes the lake. If the lake water level is higher than the groundwater level throughout the year and the hydraulic gradient points towards the aquifer, it is designated as the infiltration zone where the lake continuously infiltrates and replenishes the groundwater. If the water level difference between the dry season and the rainy season reverses, and the direction of hydraulic exchange changes with the seasons, it is designated as the transition zone where water exchange occurs in both directions between the dry and rainy seasons. In the three sub-regions of discharge zone, infiltration zone and transition zone, three-dimensional monitoring sections are set up along the hydraulic gradient direction; vertical sampling points of the lake are set up in each section, and 3-5 groundwater monitoring wells of different burial depths are set up to fully cover shallow unconfined water and medium confined water. At the same time, multi-layer soil solution samplers are set up in the vadose zone of the lake shore. Integrated online monitoring probes are deployed at each monitoring point to continuously collect water level data at each point at preset time intervals.

[0009] Furthermore, the hydraulic gradient, hydraulic exchange direction, and interface flux of each sub-region were determined, and the lake-groundwater hydraulic constraint matrix was constructed. The method used was as follows: Based on continuous water level monitoring data from each three-dimensional monitoring section, the water level values ​​of the lake and groundwater at different burial depths within the section are extracted during the same period. The ratio of the water level difference between adjacent monitoring points to the length of the seepage path is calculated along the seepage path to obtain the lateral and vertical hydraulic gradients of each sub-region. The direction of hydraulic exchange is determined based on the relationship between the lake water level and the nearshore groundwater level. The specific determination logic is as follows: when the groundwater level is higher than the lake water level throughout the year, it is determined that the groundwater is replenishing the lake; when the lake water level is higher than the groundwater level throughout the year, it is determined that the lake is infiltrating and replenishing the groundwater. For the transition zone where there is bidirectional exchange between dry and rainy seasons, the hydraulic gradient and corresponding exchange direction are calculated independently for the dry season and the rainy season. For the interface flux of each sub-region, the specific calculation method is as follows: Darcy's law is used to calculate the exchange flux per unit area of ​​the lake-groundwater interaction interface in each sub-region. A lake-groundwater hydraulic constraint matrix is ​​constructed using the discharge zone, infiltration zone, and transition zone as matrix row elements, and hydraulic gradient, hydraulic exchange direction, interface flux per unit area, total interface flux, and exchange intensity level as matrix column elements. The exchange intensity level is divided into three levels: weak, medium, and strong, based on the absolute value of the flux per unit area, and the grading threshold is set to match the hydrogeological conditions of the sandy land. This constraint matrix serves as a unified hydrological constraint benchmark for subsequent hydraulic flux verification, isotope correction, and pollution source analysis.

[0010] Furthermore, a multi-isotope-water chemical fingerprint database was constructed using the following method: Using sub-region type, monitoring point number, sampling depth, and sampling time period as indexes, the isotope test results, water chemical parameters, and in-situ monitoring data of each sample are linked and stored one by one; the isotope-water chemical characteristic values ​​of typical pollution source end-members are simultaneously included to form a multi-isotope-water chemical fingerprint database covering water samples and pollution source end-members.

[0011] Furthermore, local atmospheric precipitation lines and lake evaporation lines were established, and the groundwater inflow flux and lake water leakage in each sub-region were calculated using the radon-222 mass balance equation. The method used was as follows: Based on the constructed multi-isotope-water chemical fingerprint database, atmospheric precipitation hydrogen and oxygen isotope monitoring data for three complete hydrological years in the sub-region were obtained. The local atmospheric precipitation line was obtained by fitting oxygen isotope as the abscissa and hydrogen isotope as the ordinate. The test results of hydrogen and oxygen isotopes of lake water samples in each sub-region were extracted, and the lake water evaporation lines of the discharge area, infiltration area and transition area were fitted respectively. The evaporation enrichment degree of lake water in each sub-region was quantified by the difference in slope and intercept between the evaporation line and the atmospheric precipitation line. Radon-222 mass balance equations were constructed for the three types of sub-regions, with the lake units in each sub-region as independent calculation objects, to calculate the groundwater inflow flux and lake water leakage in each sub-region.

[0012] Furthermore, the flux results were verified and corrected, and a dynamic fractionation correction model was used to eliminate the isotope enrichment effect caused by strong evaporation in sandy areas, in order to obtain the original isotope data. The method used was as follows: The constructed lake-groundwater hydraulic constraint matrix is ​​retrieved, and the interface flux calculated based on Darcy's law for each sub-region is extracted as a reference. The sub-region flux calculated by radon-222 mass balance is compared with the reference value. If the relative deviation between the two exceeds the preset threshold, the gas exchange rate and the equivalent permeability coefficient of the interface are corrected by combining the sediment particle size distribution and sandy wind speed conditions of the sub-region. The flux is then recalculated iteratively until the deviation meets the preset threshold, and the verified and corrected sub-region hydraulic exchange flux result is obtained. Based on the evaporation fractionation model, and considering the strong evaporation characteristics of sandy areas with high temperature, low relative humidity, and high wind speed, local measured meteorological parameters are introduced to correct the dynamic fractionation term. Measured values ​​of hydrogen and oxygen isotopes in lake water, mean values ​​of atmospheric precipitation isotopes, near-surface air temperature, lake surface relative humidity, and lake surface wind speed parameters are input for each sub-region. The water balance fractionation coefficient and dynamic fractionation coefficient are calculated, and a calculation model for lake water evaporation enrichment corresponding to each sub-region is constructed. The isotopic composition of the remaining water body after evaporation satisfies the Rayleigh fractionation relation. Based on the verified and corrected hydraulic exchange flux of each sub-region, the evaporation residual ratio of lake water in each sub-region is calculated by combining the lake water balance equation and substituted into the dynamic fractionation correction model to obtain the original isotopic composition of the lake water replenishment source. The dynamic fractionation correction model is based on the evaporation fractionation model.

[0013] Furthermore, the corrected original isotopic data are compared with the end-member isotopic characteristic values ​​of the pollution source, and the pollution sources in each zone are identified by combining water chemical ion ratios and land use information. The method used is as follows: Samples of five typical pollution sources in the target area, namely chemical fertilizer, human and animal excrement, domestic sewage, soil organic nitrogen, and mining wastewater, were collected and relevant indicators were tested. The characteristic value ranges of nitrate nitrogen-15, nitrate oxygen-18, and sulfate sulfur-34 isotopes of various pollution sources and the corresponding ion ratio ranges were established to form a standardized pollution source end-member characteristic set. The corrected water isotope data of each sub-region are compared with the end-member feature set. The matching degree is quantitatively determined by the isotope feature distance, and the molar concentration ratio of water feature ions is calculated to eliminate the interference of dilution and evaporation concentration. By combining the land use type and spatial distribution data of pollution sources corresponding to each sub-region, candidate sources with mismatched spatial distribution are eliminated, and the main pollution source types of each sub-region are finally determined.

[0014] Furthermore, the isotopic fractionation induced by biotransformation was corrected by the isotopic fractionation coefficient of microbial action, thus restoring the initial isotopic composition of the pollution source. The method used is as follows: By combining the measured data of dissolved oxygen and redox potential of water bodies in each sub-region, as well as the spatial decay law of nitrate and sulfate concentrations along the hydraulic gradient, it is determined whether denitrification and sulfate reduction microbial transformation processes occur in the corresponding zone; when the corresponding substrate concentration decreases along the path and the isotope value is enriched, it is determined that there is a biofractionation effect and the dominant reaction type is determined. The Rayleigh fractionation model was used to quantitatively correct the isotopic fractionation effect in the biotransformation process. Based on the correction formula and the isotopic values ​​of the water body after the reaction, the initial isotopic values ​​of the pollution source were inferred. In view of the hydrological environment characteristics of sandy land with high temperature and strong aeration of the vadose zone, in-situ culture experiments are preferred to determine the isotope enrichment coefficients of denitrification and sulfate reduction processes in the target area. When in-situ measurement is not possible, empirical values ​​from literature in sandy areas with the same climate and lithology are selected as substitutes. By substituting the measured isotope values, enrichment coefficients, and residual substrate fractions into the correction formula, the superimposed influence of biotransformation on isotope composition is eliminated, and the initial isotope characteristic dataset of the pollution sources corresponding to each sub-region is restored, which serves as the standardized endmember input parameters for subsequent quantitative source apportionment.

[0015] Furthermore, a Bayesian mixture model was used to calculate the pollution contribution rate of each pollution source. After cross-validation, the spatial distribution and contribution ratio of the pollution sources were output. The method used was as follows: The initial isotope values ​​of each pollution source endmember obtained by restoration are used as the input terms of the model endmembers, and the corrected isotope values ​​of each sub-region water body are used as the input terms of the mixed sample. The flux data of the partition interface in the lake-groundwater hydraulic constraint matrix are retrieved, and the flux ratio corresponding to each pollution migration path is converted into the prior distribution parameters of the Bayesian mixture model. Hydrological path hard constraints are applied to the range of values ​​of each pollution source contribution rate. The specific method for quantitatively solving the pollution contribution rate of each pollution source using the Bayesian mixture model is as follows: Based on the principle of isotope mass balance, a multi-isotope mixture equation is constructed, and the posterior probability distribution of the contribution rate of each pollution source is solved by Markov chain Monte Carlo sampling. The mean, standard deviation and 95% confidence interval of the contribution rate of each type of pollution source in each sub-region are output. Finally, the contribution rate results are combined with the flux data of the lake-groundwater hydraulic constraint matrix and substituted into the groundwater solute transport model to simulate the spatial distribution of pollutant concentration. The degree of fit between the simulated and measured values ​​is quantified by the coefficient of determination. The contribution ratio of each pollution source is statistically analyzed according to the discharge zone, infiltration zone, and transition zone. The final spatial distribution results of pollution sources and the contribution rate ranking list are output.

[0016] The present invention also provides a method for source apportionment of pollution in sandy lakes and groundwater based on multiple isotopes. The method is executed by the aforementioned system for source apportionment of pollution in sandy lakes and groundwater, and includes: Step 1: Based on the hydrogeological data and dynamic water level monitoring data of the target area, combined with the complementary relationship between the groundwater flow field and the lake water level, three types of sub-regions are divided, water level data of each sub-region are collected, hydraulic gradient, exchange direction and interface flux are calculated, and a lake-groundwater hydraulic constraint matrix is ​​constructed. The sub-regions include discharge zone, infiltration zone and transition zone. Step 2: Collect lake water, groundwater, bottom sediment pore water and vadose zone soil solution samples for the three types of sub-regions, determine the hydrogen and oxygen, nitrate nitrogen and oxygen, sulfate sulfur isotopes and radon-222 activity, and simultaneously detect the anion and cation concentrations and total dissolved solids content to construct a multi-isotope-water chemical fingerprint database. Step 3: Based on hydrogen and oxygen isotope data, establish local atmospheric precipitation lines and regional lake evaporation lines. Calculate the groundwater inflow flux and lake water leakage in each region using the radon-222 mass balance equation. Verify and correct the flux results using the constraint matrix. Use a dynamic fractionation correction model to eliminate the strong evaporation isotope enrichment effect in sandy areas. Invert to obtain the original isotope data of the lake water recharge source. Step 4: Compare the corrected original isotope data with the characteristic values ​​of the pollution source end-member isotopes, identify the pollution sources in each zone by combining the water chemical ion ratio and land use information, correct the isotope fractionation caused by biotransformation by the isotope fractionation coefficient of microbial action, and restore the initial isotope composition of the pollution source. Step 5: Using the restored initial isotope values ​​of the pollution sources as input, and combining the flux weight constraints of the constraint matrix, a Bayesian mixture model is used to calculate the pollution contribution rate of each pollution source. After cross-validation, the spatial distribution and contribution ratio of the pollution sources are output.

[0017] Compared with the prior art, the beneficial effects of the present invention are: This invention significantly improves the calculation accuracy of hydraulic exchange flux in various zones of the lakeside zone in sandy areas by constructing a lake-groundwater hydraulic constraint matrix and combining the radon-222 mass balance method with Darcy's law cross-validation mechanism. At the same time, it establishes a sandy-adaptive isotope dual-order correction system, which eliminates the isotope enrichment bias caused by strong evaporation through dynamic fractionation correction and restores the initial isotopic characteristics of pollution sources through microbial fractionation correction. This effectively solves the problem of water body isotope data distortion under the special hydrological environment of sandy areas and provides real and reliable basic data support for pollution source analysis. Furthermore, this invention incorporates a Bayesian mixture model with hydrological flux weight constraints to conduct quantitative source apportionment, and utilizes hydraulic laws to narrow the model solution space, effectively reducing the ambiguity of traditional isotope source apportionment; at the same time, it conducts differentiated zoning analysis for discharge zone, infiltration zone, and transition zone, accurately matching the spatial heterogeneity of seasonal bidirectional exchange between sandy lakes and groundwater, and can output refined results on the spatial distribution and contribution ratio of pollution sources. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the system module flow of the present invention.

[0019] Figure 2 This is a schematic diagram of the overall method flow of the present invention. Detailed Implementation

[0020] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0021] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0022] Example: Please see Figure 1 A multi-isotope-based source apportionment system for pollution in sandy lakes and groundwater includes the following steps: The hydrological zoning module, based on the hydrogeological data and dynamic water level monitoring data of the target area, and combined with the complementary relationship between the groundwater flow field and the lake water level, divides the area into three sub-regions, collects water level data of each sub-region, calculates the hydraulic gradient, exchange direction and interface flux, and constructs the lake-groundwater hydraulic constraint matrix. The sub-regions include the discharge zone, the infiltration zone and the transition zone. In this embodiment, the hydrological differentiation module abandons the traditional calculation model of homogenization of the entire lake. Based on the dynamic laws of groundwater replenishment and discharge between the lake and the lake, it divides the lake into discharge zone, infiltration zone, and seasonal bidirectional exchange transition zone. This accurately matches the characteristics of strong spatial heterogeneity of hydraulic exchange in sandy lake shore zones and the reversal of exchange direction with the seasons. It effectively avoids the systematic error of the overall average calculation masking local interaction differences. This provides a spatial benchmark framework for subsequent refined pollution source analysis by zone and constructs a unified hydrological constraint benchmark for the entire process. By quantifying the hydraulic gradient, exchange direction, and interface flux of each zone, a structured hydraulic constraint matrix is ​​formed. This provides an independent hydrological verification basis for the flux calculation of the radon-222 mass balance method and also assigns flux weight constraints to the source analysis of the terminal Bayesian mixture model.

[0023] Furthermore, the target area is divided into several sub-regions, and water level data for each sub-region is collected. These sub-regions include a discharge zone, an infiltration zone, and a transition zone. The method used is as follows: Hydrogeological data of the target area were retrieved, and the permeability coefficient of sediments in the lacustrine zone, the thickness of the vadose zone, and the characteristics of the regional groundwater flow field were extracted. Combined with long-sequence dynamic water level monitoring data of a complete hydrological year, the dynamic difference and seasonal fluctuation range between the lake water level and the groundwater level at the corresponding points in the surrounding area were calculated. Collect water level monitoring data for a complete hydrological year, calculate the dynamic difference between the lake water level and the groundwater level at corresponding points in the surrounding area on a daily basis, and divide the area into sub-regions according to the following quantitative rules: If, within a hydrological year, the groundwater level is higher than the lake water level for ≥90% of the days, and the hydraulic gradient direction consistently points towards the lake, it is designated as a discharge zone; if, within a hydrological year, the lake water level is higher than the groundwater level for ≥90% of the days, and the hydraulic gradient direction consistently points towards the aquifer, it is designated as an infiltration zone; if, within a hydrological year, the water level difference between the dry season and the rainy season reverses, and the corresponding direction accounts for ≥60% of the consecutive days in both seasons, it is designated as a transition zone. Within the three sub-regions of discharge zone, infiltration zone, and transition zone, three-dimensional monitoring sections are deployed along the hydraulic gradient direction. Each section extends 50m towards the lake and 100m towards the land. Within each section, two vertical sampling points are set up in the lake area. Four groundwater monitoring wells are deployed sequentially towards the land, with well depths of 2m, 5m, 10m, and 15m, and a well spacing of 20m, completely covering shallow unconfined water and medium-level confined water. Simultaneously, a set of multi-layer soil solution samplers is deployed in the vadose zone of the lakeshore, with sampling depths of 10cm, 30cm, 60cm, and 100cm. Integrated online monitoring probes are deployed at each monitoring point, with a sampling interval of 1 hour to continuously collect water level data at each point. A manual water level check is conducted once a month to ensure that the water level monitoring error is ≤ ±1cm.

[0024] Furthermore, the hydraulic gradient, hydraulic exchange direction, and interface flux of each sub-region were determined, and the lake-groundwater hydraulic constraint matrix was constructed. The method used was as follows: Based on continuous water level monitoring data from various three-dimensional monitoring sections, the water level values ​​of the lake and groundwater at different depths within the sections are extracted. The ratio of the water level difference between adjacent monitoring points to the length of the seepage path is calculated along the seepage path to obtain the lateral and vertical hydraulic gradients of each sub-region. The formula used to calculate the hydraulic gradient is as follows: In the formula: This represents the hydraulic gradient of each sub-region. The difference in water level between adjacent monitoring points in each sub-region. The seepage path length is for each sub-region. It should be noted that the lateral hydraulic gradient is calculated using the water level difference between adjacent monitoring wells along the horizontal direction of the cross section, and the seepage path length is the horizontal straight-line distance between the two points. The vertical hydraulic gradient is calculated using the water level difference between monitoring wells at different burial depths at the same point, and the seepage path length is the vertical distance between the midpoints of the filter pipes of the two wells. The direction of hydraulic exchange is determined based on the relationship between the lake water level and the nearshore groundwater level. The specific determination logic is as follows: when the groundwater level is higher than the lake water level throughout the year, it is determined that the groundwater is recharging the lake; when the lake water level is higher than the groundwater level throughout the year, it is determined that the lake is infiltrating and recharging the groundwater. For the transition zone, the average water level during the dry season and the rainy season are extracted separately, and the hydraulic gradient and exchange direction of the two periods are calculated independently. The specific calculation method for the interface flux in each sub-region is as follows: Darcy's law is used to calculate the exchange flux per unit area of ​​the lake-groundwater interface in each sub-region, based on the following formula: in, The exchange flux per unit area of ​​each sub-region, The equivalent permeability coefficient of the interface sediments was determined through field double-ring infiltration tests. At least three test points were set up in each sub-region, and the geometric mean was taken as the partition of each sub-region. value; Total interface flux of each sub-region The calculation formula is: In the formula: Let be the area of ​​the interactive interface of the corresponding sub-region, and It is calculated by multiplying the length of the lakeshore in the sub-region by the width of the lateral interaction zone; A lake-groundwater hydraulic constraint matrix is ​​constructed using the discharge zone, infiltration zone, and transition zone as matrix row elements, and hydraulic gradient, hydraulic exchange direction, interfacial flux per unit area, total interfacial flux, and exchange intensity level as matrix column elements. The exchange intensity level is divided into three levels—weak, moderate, and strong—based on the absolute value of the flux per unit area. The grading threshold is set to match the hydrogeological conditions of the sandy land. The grading threshold is set as follows: Weak exchange (… ), China Exchange ( ), strong swap ( This constraint matrix serves as a unified hydrological constraint benchmark for subsequent hydraulic flux verification, isotope correction, and pollution source apportionment. It is the standard unit of flux in the fields of hydrogeology and groundwater dynamics, and its full name is "meter per day".

[0025] The sampling and testing module is used to collect lake water, groundwater, bottom sediment pore water, and vadose zone soil solution samples from three types of sub-regions, and to determine the activities of hydrogen and oxygen, nitrate nitrogen and oxygen, sulfate sulfur isotopes, and radon-222. Simultaneously, it detects the concentration of anions and cations and the total dissolved solids content, and constructs a multi-isotope-water chemical fingerprint database.

[0026] In this embodiment, four types of samples were collected from the three-dimensional monitoring points in the three-dimensional monitoring areas of the discharge area, infiltration area and transition area: lake water, groundwater, bottom sediment pore water and vadose zone soil solution. The sampling frequency was one systematic sampling each during the dry season and rainy season each year, continuously covering no less than one complete hydrological year. A set of parallel samples was set up for every 10 samples, and blank samples were collected simultaneously in the field and transported for quality control.

[0027] All sample tests were performed in accordance with national or industry standard methods, as detailed below: Hydrogen and oxygen isotopes ( The wavelength scanning cavity ring-down spectroscopy method was used for testing, with the Vienna standard average seawater as the reference standard. The required testing accuracy was as follows: , ; Nitrate nitrogen-15, nitrate oxygen-18 isotopes ( After being converted to nitrous oxide using denitrifying bacteria, the nitrous oxide was analyzed by isotope mass spectrometry. The required accuracy for the analysis was as follows: , ; Sulfate sulfur-34 isotopes: Elemental analysis-isotope mass spectrometry was used for testing, with the Vienna Diablo Valley meteorite as a reference standard. The testing accuracy was [not specified]. ; Radon-222 activity: determined on-site using a portable radon analyzer via the air extraction method. The unit of measurement is [unit missing]. Detection limit ; The main anions and cations were tested using ion chromatography, and the cations were simultaneously calibrated using inductively coupled plasma atomic emission spectrometry. The detection limit was... Total dissolved solids were determined by gravimetric method and simultaneously verified by conductivity value.

[0028] Furthermore, a multi-isotope-water chemical fingerprint database was constructed using the following method: Using sub-region type, monitoring point number, sampling depth, and sampling time period as indexes, the isotope test results, water chemical parameters, and in-situ monitoring data of each sample are linked and stored one by one; the isotope-water chemical characteristic values ​​of typical pollution source end-members are also included, forming a multi-isotope-water chemical fingerprint database covering water samples and pollution source end-members. All data must undergo quality control verification before being entered, and abnormal data with a relative deviation of more than 10% for parallel samples are removed to ensure the reliability of the data in the database.

[0029] The flux verification and inversion module establishes local atmospheric precipitation lines and regional lake evaporation lines based on hydrogen and oxygen isotope data. It calculates the groundwater inflow flux and lake water leakage in each region using the radon-222 mass balance equation, verifies and corrects the flux results using the constraint matrix, and uses a dynamic fractionation correction model to eliminate the strong evaporation isotope enrichment effect in sandy areas. The module then retrieves the original isotope data of the lake water recharge source.

[0030] In this embodiment, the interface flux calculated by Darcy's law in the hydraulic constraint matrix is ​​used as an independent hydrological benchmark. It is compared and verified with the isotope tracer flux obtained by the radon-222 mass balance method. The calculation deviation is narrowed by iteratively correcting the gas exchange coefficient and the interface equivalent permeability coefficient parameters. At the same time, for the hydrological characteristics of sandy land with high temperature and strong evaporation, a dynamic fractionation correction model adapted to sandy land is constructed based on Rayleigh fractionation theory. The remaining proportion of lake water evaporation is calculated based on the verified hydraulic flux. The evaporation enrichment effect is deduced in reverse to restore the original isotope composition of the lake water recharge source, forming a closed-loop logic of "flux verification supporting isotope correction and isotope results feeding back into hydrological understanding".

[0031] Furthermore, local atmospheric precipitation lines and lake evaporation lines were established, and the groundwater inflow flux and lake water leakage in each sub-region were calculated using the radon-222 mass balance equation. The method used was as follows: Based on the constructed multi-isotope-water chemical fingerprint database, atmospheric precipitation hydrogen and oxygen isotope monitoring data for three complete hydrological years in the sub-region were obtained. The local atmospheric precipitation line was obtained by fitting oxygen isotope as the abscissa and hydrogen isotope as the ordinate. The test results of hydrogen and oxygen isotopes of lake water samples in each sub-region were extracted, and the lake water evaporation lines of the discharge area, infiltration area and transition area were fitted respectively. The evaporation enrichment degree of lake water in each sub-region was quantified by the difference in slope and intercept between the evaporation line and the atmospheric precipitation line. Radon-222 mass balance equations were constructed for the three types of sub-regions, with each sub-region's lake unit as an independent calculation object. The groundwater inflow flux and lake water leakage in each sub-region were calculated using the following formulas: In the formula, This represents the interface exchange flux between sub-regions per unit time. Positive values ​​indicate groundwater recharge to the lake, while negative values ​​indicate lake infiltration recharge to groundwater. To correspond to the radon-222 activity in the lake water in each sub-region, The radioactive decay constant of radon-222 The radon-222 gas exchange rate at the water-gas interface. The activity of radon-222 in the near-surface atmosphere. The volume of lake water corresponding to each sub-region. The background activity of radon-222 in groundwater in the corresponding sub-regions; The calculation result for the discharge zone is positive, and the output is the groundwater inflow flux into the lake; the calculation result for the infiltration zone is negative, and the absolute value is used to output the lake water leakage; the transition zone is calculated independently based on the monitoring data of the dry season and the rainy season, corresponding to the two types of fluxes exchanged in both directions.

[0032] Furthermore, the flux results were verified and corrected, and a dynamic fractionation correction model was used to eliminate the isotope enrichment effect caused by strong evaporation in sandy areas, in order to obtain the original isotope data. The method used was as follows: Retrieve the interfacial flux per unit area calculated based on Darcy's law for each subregion in the lake-groundwater hydraulic constraint matrix. As a reference standard, the total volumetric flux calculated using the radon-222 mass balance equation is... Divide by the lake-groundwater interface area of ​​the corresponding sub-region to obtain the radon-222 flux per unit area. ; The flux of the sub-region calculated by radon-222 mass balance is compared with the benchmark value. If the relative deviation between the two exceeds the preset threshold, the gas exchange rate and the equivalent permeability coefficient of the interface are corrected by combining the sediment particle size distribution of the sub-region and the sandy wind speed conditions. The flux is recalculated iteratively until the deviation meets the preset threshold, and the verified and corrected hydraulic exchange flux of the sub-region is obtained. The preset allowable relative deviation threshold between the two is ±20%, and the relative deviation is calculated using the following formula: In the formula: This represents the relative deviation between the flux of the sub-region and the corresponding baseline value. like Then, the two types of core parameters are corrected according to the following rules, and the flux is iteratively recalculated until the deviation meets the threshold. The method used is as follows: To address the deviation in gas exchange parameters caused by high wind speeds in sandy areas, a Schmidt number-normalized wind speed model is used for correction. The corrected formula is: In the formula: The gas exchange coefficient corresponding to a Schmidt number of 600 is given by the daily average wind speed at a height of 10m above the lake surface. calculate, , The Schmidt number of -222 in water is calculated from the measured water temperature. The wind speed index is used, taking -2 / 3 for wind speeds < 3 m / s and -1 / 2 for wind speeds ≥ 3 m / s. The corrected value is then... Substitute into the radon-222 mass balance formula and recalculate. and ; To address the Darcy flux deviation caused by uncertainties in sediment parameters, and based on measured data of sediment grain size distribution in sub-regions, the equivalent permeability coefficient was corrected using the Hazen empirical formula. The formula used is as follows: In the formula: The effective grain size of the sediment. The Hazen empirical coefficient ranges from 80 to 120 for sandy sediments in sandy areas, with higher values ​​indicating higher grain uniformity. The corrected value will be... Substitute Darcy's law into the baseline flux to recalculate. After completing a single round of parameter correction, the relative deviation is recalculated. If it still exceeds the threshold, the correction iteration is repeated. The iteration terminates at the specified time, and the verified and corrected sub-region hydraulic exchange flux results are output.

[0033] Based on the evaporative fractionation model, and considering the strong evaporation characteristics of sandy areas with high temperature, low relative humidity, and high wind speed, local measured meteorological parameters are introduced to correct the dynamic fractionation term. The formula used is as follows: In the formula: The relative humidity of the lake surface. The dynamic fractionation coefficient; The final total isotope enrichment coefficient is obtained as follows: In the formula: The vapor-liquid-water equilibrium fractionation coefficient representing hydrogen and oxygen isotopes; Input the measured values ​​of hydrogen and oxygen isotopes in the lake water, the mean isotope values ​​of atmospheric precipitation, the near-surface air temperature, the relative humidity of the lake surface, and the wind speed parameters of the lake surface for each sub-region. Calculate the water balance fractionation coefficient and the dynamic fractionation coefficient, and construct a calculation model for lake water evaporation enrichment for each sub-region. The isotopic composition of the water remaining after evaporation satisfies the Rayleigh fractionation relation. The water balance equation for the sub-region is: In the formula: The daily precipitation over the lake surface in the sub-region. The lateral outflow of lake water in the sub-region The vertical seepage of lake water in the sub-region. This refers to the daily evaporation intensity of the lake surface. The evaporation surplus ratio is the proportion of water remaining after evaporation to the total replenishment water volume, calculated using the following formula: ; Based on the verified and corrected hydraulic exchange flux of each sub-region, the evaporation residual ratio of lake water in each sub-region is calculated by combining the lake water balance equation and substituted into the dynamic fractionation correction model to obtain the original isotopic composition of the lake water replenishment source. The dynamic fractionation correction model is based on the evaporation fractionation model. The isotopic composition of the water remaining after evaporation follows the Rayleigh fractionation relation. By transforming the formula, the original isotopic values ​​of the replenishment source can be obtained by reverse calculation. In the formula: These are the measured isotope values ​​corresponding to the lake water in the sub-region. The original isotopic values ​​of the lake water recharge source obtained by inversion, The remaining proportion of evaporation in each subregion is calculated by substituting the hydrogen and oxygen isotope parameters, and finally the original hydrogen and oxygen isotope data of the supply source after eliminating the strong evaporation enrichment effect in each subregion are obtained.

[0034] It should be noted that the reason for correcting the gas exchange rate and the equivalent permeability coefficient parameters when the preset threshold is exceeded is as follows: In the radon-222 mass balance calculation, the gas exchange rate at the water-gas interface directly determines the accuracy of the radon-222 atmospheric emission flux calculation. Sandy lakes are generally open and unobstructed, with high wind speeds and drastic seasonal fluctuations. General empirical values ​​cannot match the local wind field characteristics, which can easily lead to distortion in the emission flux calculation and ultimately reduce the reliability of the radon method flux. Correcting the gas exchange rate by combining measured wind speed with the Schmidt number model can eliminate the parameter deviation caused by local strong winds. The core parameter for calculating flux using Darcy's law is the sediment permeability coefficient. However, the sandy sediments in the lacustrine zone are highly heterogeneous with large spatial differences in particle size distribution. The permeability coefficient obtained from the initial survey has limited representativeness, and applying it across the entire range will produce systematic errors. By using measured sediment particle size distribution data and combining the Hazen empirical formula to correct the equivalent permeability coefficient, the Darcy method benchmark flux can be made to better fit the actual hydrogeological conditions of the zone.

[0035] It should be noted that the reason for eliminating the evaporation enrichment effect and restoring the original isotopic composition of the lake water source is as follows: During the natural evaporation of water, lighter isotopes such as hydrogen and oxygen are more easily released into the atmosphere with water vapor, while heavier isotopes remain more in the remaining water. This leads to a systematic enrichment shift in the hydrogen and oxygen isotope values ​​of the lake water compared to the source water. The higher the evaporation intensity, the more significant the enrichment deviation. Sandy areas are characterized by high temperature, low relative humidity, and high wind speed, resulting in a much higher evaporation intensity than humid areas, leading to a very prominent isotope enrichment effect. Directly using the measured isotope values ​​of the lake water for analysis would bring significant problems. System bias: Subsequent pollution source identification requires comparing the water body isotopes with the original isotopic fingerprints of the pollution source endmembers and the replenishment endmembers. The feature values ​​of each type of endmember are the original values ​​that are not affected by evaporation. If the lake water data is superimposed with the evaporation enrichment effect, it will directly lead to fingerprint matching errors and misjudgment of the pollution source type and the proportion of the replenishment source. Furthermore, the quantitative calculation of the terminal Bayesian mixture model is based on the isotopic mass balance principle, which requires that the input water body isotope values ​​represent the true characteristics after multi-source mixing, rather than the distorted values ​​after superimposed evaporation and fractionation. Otherwise, it will directly lead to the calculation results of the pollution source contribution rate deviating from the true situation.

[0036] The fractionation and reduction module is used to compare the corrected original isotope data with the characteristic values ​​of the end-member isotopes of the pollution source, identify the pollution sources of each zone by combining the water chemical ion ratio and land use information, correct the isotope fractionation caused by biotransformation through the isotope fractionation coefficient of microbial action, and restore the initial isotope composition of the pollution source.

[0037] In this embodiment, the module constructs a three-level progressive pollution identification system of "multi-isotope quantitative matching - hydrochemical ratio corroboration - spatial distribution verification". It quantifies the end-member matching degree with multi-isotope characteristic distance, eliminates concentration interference caused by dilution and evaporation concentration by combining the characteristic ion ratio of chloride ion normalization, and eliminates candidate sources with spatial distribution mismatch by using land use information. It effectively solves the problems of overlapping end-member fingerprints and strong subjectivity in traditional single isotope source identification, significantly improves the accuracy and reliability of qualitative determination of pollution sources in each zone, and clarifies the scope of objects for subsequent targeted correction and quantitative analysis.

[0038] Furthermore, the corrected original isotopic data are compared with the end-member isotopic characteristic values ​​of the pollution source, and the pollution sources in each zone are identified by combining water chemical ion ratios and land use information. The method used is as follows: Samples of five typical pollution sources in the target area, namely chemical fertilizer, human and animal excrement, domestic sewage, soil organic nitrogen, and mining wastewater, were collected and relevant indicators were tested. The characteristic value ranges of nitrate nitrogen-15, nitrate oxygen-18, and sulfate sulfur-34 isotopes of various pollution sources and the corresponding ion ratio ranges were established to form a standardized pollution source end-member characteristic set. The corrected water isotope data for each sub-region are compared with the endmember feature set. The Euclidean distance is used to calculate the comprehensive isotope feature distance, and the fingerprint matching degree between the sample and various pollution sources is quantitatively determined. The formula used is as follows: In the formula: The distance is defined as the isotopic composite characteristic distance between water samples in each sub-region and the corresponding pollution source end-member. , , These are the isotopic values ​​of nitrate nitrogen-15, nitrate oxygen-18, and sulfate sulfur-34 in the water sample to be identified. , , These are the isotopic characteristic values ​​of nitrate nitrogen-15, nitrate oxygen-18, and sulfate sulfur-34, respectively, for the corresponding pollution end-members; The molar concentration ratio of characteristic ions in water bodies is calculated to eliminate the interference of dilution and evaporation concentration. The specific method is as follows: the molar ratio of nitrate / chloride ions and the molar ratio of sulfate / chloride ions are selected as auxiliary discrimination indicators. Chloride ions are conservative ions and are not affected by biological transformation and evaporation fractionation. They can be normalized to eliminate the interference of concentration fluctuations. The calculated ratio is compared with the corresponding ratio range in the endmember feature set. Candidate pollution sources with completely mismatched ratio ranges are eliminated to further narrow down the source type range. By combining the land use type and spatial distribution data of pollution sources corresponding to each sub-region, candidate sources with mismatched spatial distribution are eliminated, and the main pollution source types of each sub-region are finally determined.

[0039] Furthermore, the isotopic fractionation induced by biotransformation was corrected by the isotopic fractionation coefficient of microbial action, thus restoring the initial isotopic composition of the pollution source. The method used is as follows: Based on the measured data of dissolved oxygen (DO) and oxidation-reduction potential (ORP) of each sub-region water body, as well as the spatial decay pattern of nitrate and sulfate concentrations along the hydraulic gradient, the dominant microbial transformation process was determined, and the quantitative determination rules are as follows; When DO < 2 mg / L, ORP < 200 mV, and the nitrate concentration decreases by ≥ 30% along the seepage path, and nitrate nitrogen-15 / oxygen-18 isotopes are enriched simultaneously, a significant denitrification process is determined to exist. When the water body has DO < 1 mg / L, ORP < -100 mV, and the sulfate concentration decreases by ≥ 20% along the seepage path, and the sulfate sulfur-34 isotope is significantly enriched, a significant sulfate reduction process is determined to exist. For sub-regions identified as exhibiting significant biofractionation, the Rayleigh fractionation model was used to quantitatively correct the isotopic fractionation effect of the biotransformation process. The initial isotopic values ​​at the time of pollution source input were inferred from the measured isotopic values ​​of the water body after the reaction. The correction formula is as follows: in, These are the initial isotopic values ​​when pollution sources in each sub-region input water bodies. The values ​​are the measured isotope values ​​of the water bodies in each sub-region after evaporation correction. The isotopic enrichment coefficients corresponding to the microbial processes are: for denitrification, the enrichment coefficients are for nitrate nitrogen-15 and nitrate oxygen-18, respectively; and for sulfate reduction, the enrichment coefficient is for sulfate sulfur-34, respectively. The remaining substrate mole fraction represents the proportion of substrate remaining after biotransformation relative to the initial substrate, and the remaining substrate fraction... Calculations were performed using the conservative ion normalization method; In response to the hydrological characteristics of sandy areas, such as high temperature and strong aeration in the vadose zone, in-situ culture experiments were preferred to determine the isotope enrichment coefficients of the target areas. The in-situ, light-protected culture bottle method was used, where in-situ water samples were collected from the corresponding sub-regions, a quantitative amount of substrate was added, and the samples were cultured in the dark. Substrate concentration and isotope values ​​were monitored periodically, and the enrichment coefficients were calculated using a Rayleigh model. When in-situ measurement was not possible, empirical values ​​from literature were selected for sandy areas with similar climates and lithologies. Typical values ​​for sandy areas were: nitrogen enrichment coefficient during denitrification (-20‰ to -5‰), oxygen enrichment coefficient (-10‰ to -3‰), and sulfur enrichment coefficient during sulfate reduction (-25‰ to -5‰). By substituting the measured isotope values, enrichment coefficients, and residual substrate fractions into the correction formula, the superimposed influence of biotransformation on isotope composition is eliminated, and the initial isotope characteristic dataset of the pollution sources corresponding to each sub-region is restored, which serves as the standardized endmember input parameters for subsequent quantitative source apportionment.

[0040] The model analysis module is used to take the initial isotope values ​​of the restored pollution sources as input, combine them with the flux weight constraints of the constraint matrix, and use a Bayesian mixture model to calculate the pollution contribution rate of each pollution source. After cross-validation, it outputs the spatial distribution and contribution ratio of the pollution sources.

[0041] In this embodiment, the module transforms the flux weights of the lake-groundwater hydraulic constraint matrix into prior distribution parameters of a Bayesian mixture model. This overcomes the limitations of traditional isotope source apportionment, which relies solely on pure mathematical fitting of fingerprint data. By applying rigid constraints to the range of pollution source contribution rates based on hydrodynamic laws, the model's solution space is effectively narrowed, significantly reducing the ambiguity of source apportionment in scenarios with similar isotope endmember fingerprints and greatly improving the accuracy of quantitative calculations. Simultaneously, by solving the posterior probability distribution of contribution rates through Markov chain Monte Carlo sampling, the module outputs multi-dimensional statistical results including mean, standard deviation, and 95% confidence interval, which is more suitable for the uncertainty characteristics of pollution migration processes in sandy lake-groundwater systems.

[0042] Furthermore, a Bayesian mixture model was used to calculate the pollution contribution rate of each pollution source. After cross-validation, the spatial distribution and contribution ratio of the pollution sources were output. The method used was as follows: Using the initial isotope values ​​of each pollution source endmember obtained from the restoration as the model endmember input, and the corrected isotope values ​​of each sub-region's water body as the mixed sample input, the zonal interface flux data in the lake-groundwater hydraulic constraint matrix are retrieved. The flux proportions corresponding to each pollution migration path are converted into prior distribution parameters of the Bayesian mixture model. Hydrological path hard constraints are applied to the range of values ​​for the contribution rate of each pollution source. The formula for calculating the prior distribution parameters is as follows: In the formula: For the first The prior distribution parameters corresponding to a type of pollution source; the larger the parameter value, the higher the prior contribution weight of that pollution source. For the first The interface fluxes corresponding to the main migration paths of pollution sources are taken from the lake-groundwater hydraulic constraint matrix. For example, groundwater-carried pollution sources correspond to groundwater inflow flux into the lake, and lake surface infiltration pollution sources correspond to lake water leakage flux. The total number of candidate pollution sources. This is the prior intensity scaling factor, with a value range of 1-10. The default value for sandy scenarios is 5. It is used to balance the rigidity of hydrological constraints with the degrees of freedom of isotopic data, and to avoid excessive constraints from obscuring the true isotopic characteristics. The specific method for quantitatively solving the pollution contribution rate of each pollution source using a Bayesian mixture model is as follows: A multi-isotope mixture equation is constructed based on the isotope mass balance principle. The posterior probability distribution of the contribution rate of each pollution source is solved through Markov chain Monte Carlo sampling. The mean, standard deviation, and 95% confidence interval of the contribution rate of each type of pollution source in each sub-region are output. The core isotope mass balance equation in this process is: in, For each sub-region, the first The measured values ​​of the isotopic indicators in the water body correspond to three indicators: nitrate nitrogen-15, nitrate oxygen-18, and sulfate sulfur-34. For the first The contribution rate of pollution sources of this type For the first The first type of pollution source end-unit corresponding to the The initial values ​​of the isotope indexes are taken from the calibration results of the fractionation and reduction module. For the first The random error term of the model for each index follows a normal distribution with a mean of 0, and includes end-member isotope testing error and spatial heterogeneity error. The Bayesian mixture model is solved using Markov chain Monte Carlo sampling. Specifically, three independent Markov chains are used for parallel sampling, with a total of 100,000 iterations. The first 50,000 iterations are discarded as a burn-out period, and the remaining samples are sampled every 10 steps to reduce autocorrelation. The convergence diagnostic statistic from the Markov chain Monte Carlo sampling method is employed. To determine the convergence of the model, when all parameters are convergent... The model is considered convergent when all values ​​are less than 1.01; the final output is the posterior mean, standard deviation and 95% confidence interval of the contribution rate of each type of pollution source in each sub-region. In view of the bidirectional exchange between dry and rainy seasons in the transition zone, the prior parameters for interface flux calculation for the corresponding time period are substituted into the model and solved independently for each time period to match the differences in pollution migration paths in different seasons. Finally, the contribution rate results, combined with the flux data from the lake-groundwater hydraulic constraint matrix, were substituted into the groundwater solute transport model to simulate the spatial distribution of pollutant concentrations. Based on the contribution rate of each pollution source and the corresponding path flux, the pollutant inflow load of various pollution sources in each sub-region was calculated using the following formula: In the formula: For the first Pollutant load entering the lake from pollution sources of this type For the first Baseline concentrations of pollutants from pollution sources of this type. For the first Interfacial fluxes corresponding to the main migration paths of pollution sources; Substituting the load results into the groundwater solute transport model in the lakeside zone, the pollutant concentration values ​​at each monitoring point were simulated. The coefficient of determination was used to quantify the degree of fit between the simulated and measured values. The formula for calculating the coefficient of determination is as follows: in, As the coefficient of determination, This represents the measured concentration of the pollutant. To simulate concentrations in the model, A pass threshold is set for the average measured concentration. If the threshold is met, the quantitative fitting effect is considered acceptable. Based on the discharge zone, infiltration zone, and transition zone, the contribution ratio of each pollution source is statistically analyzed, and the final spatial distribution results and contribution rate ranking list of pollution sources are output. The pollution source emission list of the target area is retrieved, and the annual emission load of each type of pollution source in each sub-region is calculated and ranked from high to low. The contribution rate ranking output by the model is compared with the emission load ranking. When the ranking consistency (i.e., the coefficient of determination) is ≥70%, the result is determined to conform to the objective law of regional pollution distribution. If both verifications meet the requirements, the source resolution result is deemed valid; otherwise, the prior distribution coefficient is adjusted backtracking and the solution is iterated again. The results are statistically analyzed separately for three sub-regions: discharge zone, infiltration zone, and transition zone. The transition zone is output separately for dry season and rainy season. The contribution percentage, standard deviation, and 95% confidence interval of each type of pollution source in each sub-region are quantitatively output. Combined with the spatial zoning boundary of the lakeside zone, a spatial distribution map of the contribution percentage of pollution sources is generated. A list of pollution sources in each zone is generated in descending order of contribution rate, and the dominant pollution type is marked. Finally, a complete analysis of pollution sources is output.

[0043] Please see Figure 2 The present invention also provides a method for source apportionment of pollution in sandy lakes and groundwater based on multiple isotopes. The method is executed by the aforementioned system for source apportionment of pollution in sandy lakes and groundwater, and includes: Step 1: Based on the hydrogeological data and dynamic water level monitoring data of the target area, combined with the complementary relationship between the groundwater flow field and the lake water level, three types of sub-regions are divided, water level data of each sub-region are collected, hydraulic gradient, exchange direction and interface flux are calculated, and a lake-groundwater hydraulic constraint matrix is ​​constructed. The sub-regions include discharge zone, infiltration zone and transition zone. Step 2: Collect lake water, groundwater, bottom sediment pore water and vadose zone soil solution samples for the three types of sub-regions, determine the hydrogen and oxygen, nitrate nitrogen and oxygen, sulfate sulfur isotopes and radon-222 activity, and simultaneously detect the anion and cation concentrations and total dissolved solids content to construct a multi-isotope-water chemical fingerprint database. Step 3: Based on hydrogen and oxygen isotope data, establish local atmospheric precipitation lines and regional lake evaporation lines. Calculate the groundwater inflow flux and lake water leakage in each region using the radon-222 mass balance equation. Verify and correct the flux results using the constraint matrix. Use a dynamic fractionation correction model to eliminate the strong evaporation isotope enrichment effect in sandy areas. Invert to obtain the original isotope data of the lake water recharge source. Step 4: Compare the corrected original isotope data with the characteristic values ​​of the pollution source end-member isotopes, identify the pollution sources in each zone by combining the water chemical ion ratio and land use information, correct the isotope fractionation caused by biotransformation by the isotope fractionation coefficient of microbial action, and restore the initial isotope composition of the pollution source. Step 5: Using the restored initial isotope values ​​of the pollution sources as input, and combining the flux weight constraints of the constraint matrix, a Bayesian mixture model is used to calculate the pollution contribution rate of each pollution source. After cross-validation, the spatial distribution and contribution ratio of the pollution sources are output.

[0044] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0045] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0046] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0047] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A multi-isotope-based source apportionment system for pollution in sandy lakes and groundwater, characterized in that, The specific steps include: The hydrological zoning module is used to divide the target area into three sub-regions based on hydrogeological data and dynamic water level monitoring data, combined with the complementary relationship between groundwater flow field and lake water level. It collects water level data of each sub-region, calculates hydraulic gradient, exchange direction and interface flux, and constructs lake-groundwater hydraulic constraint matrix. The sub-regions include discharge zone, infiltration zone and transition zone. The sampling and testing module is used to collect lake water, groundwater, bottom sediment pore water and vadose zone soil solution samples from three types of sub-regions, determine the hydrogen and oxygen, nitrate nitrogen and oxygen, sulfate sulfur isotopes and radon-222 activity, simultaneously detect the concentration of anions and cations and the total dissolved solids content, and construct a multi-isotope-water chemical fingerprint database. The flux verification and inversion module is used to establish local atmospheric precipitation lines and regional lake evaporation lines based on hydrogen and oxygen isotope data. It calculates the groundwater inflow flux and lake water leakage in each region through the radon-222 mass balance equation, calls the constraint matrix to verify and correct the flux results, and uses the dynamic fractionation correction model to eliminate the strong evaporation isotope enrichment effect in sandy areas. It then retrieves the original isotope data of the lake water recharge source. The fractionation and reduction module is used to compare the corrected original isotope data with the characteristic values ​​of the end-member isotopes of the pollution source, identify the pollution sources of each zone by combining the water chemical ion ratio and land use information, correct the isotope fractionation caused by biotransformation through the isotope fractionation coefficient of microbial action, and restore the initial isotope composition of the pollution source. The model analysis module is used to take the initial isotope values ​​of the restored pollution sources as input, combine them with the flux weight constraints of the constraint matrix, and use a Bayesian mixture model to calculate the pollution contribution rate of each pollution source. After cross-validation, it outputs the spatial distribution and contribution ratio of the pollution sources.

2. The multi-isotope-based source apportionment system for pollution sources in sandy lakes and groundwater according to claim 1, characterized in that, The target area was divided into several sub-regions, and water level data were collected for each sub-region. These sub-regions included a discharge zone, an infiltration zone, and a transition zone. The method used was as follows: Hydrogeological data of the target area were retrieved, and the permeability coefficient of sediments in the lacustrine zone, the thickness of the vadose zone, and the characteristics of the regional groundwater flow field were extracted. Combined with long-term dynamic water level monitoring data of a complete hydrological year, the dynamic difference and seasonal fluctuation range between the lake water level and the groundwater level at the corresponding points in the surrounding area were calculated. If the groundwater level is higher than the lake water level throughout the year and the hydraulic gradient points towards the lake, it is designated as the discharge zone where groundwater continuously replenishes the lake. If the lake water level is higher than the groundwater level throughout the year and the hydraulic gradient points towards the aquifer, it is designated as the infiltration zone where the lake continuously infiltrates and replenishes the groundwater. If the water level difference between the dry season and the rainy season reverses, and the direction of hydraulic exchange changes with the seasons, it is designated as the transition zone where water exchange occurs in both directions between the dry and rainy seasons. In the three sub-regions of discharge zone, infiltration zone and transition zone, three-dimensional monitoring sections are set up along the hydraulic gradient direction; vertical sampling points of the lake are set up in each section, and 3-5 groundwater monitoring wells of different burial depths are set up to fully cover shallow unconfined water and medium confined water. At the same time, multi-layer soil solution samplers are set up in the vadose zone of the lake shore. Integrated online monitoring probes are deployed at each monitoring point to continuously collect water level data at each point at preset time intervals.

3. The multi-isotope-based sandy lake-groundwater pollution source apportionment system according to claim 2, characterized in that, The hydraulic gradient, hydraulic exchange direction, and interface flux of each sub-region were determined, and the lake-groundwater hydraulic constraint matrix was constructed using the following method: Based on continuous water level monitoring data from each three-dimensional monitoring section, the water level values ​​of the lake and groundwater at different burial depths within the section are extracted during the same period. The ratio of the water level difference between adjacent monitoring points to the length of the seepage path is calculated along the seepage path to obtain the lateral and vertical hydraulic gradients of each sub-region. The direction of hydraulic exchange is determined based on the relationship between the lake water level and the nearshore groundwater level. The specific determination logic is as follows: when the groundwater level is higher than the lake water level throughout the year, it is determined that the groundwater is replenishing the lake; when the lake water level is higher than the groundwater level throughout the year, it is determined that the lake is infiltrating and replenishing the groundwater. For the transition zone where there is bidirectional exchange between dry and rainy seasons, the hydraulic gradient and corresponding exchange direction are calculated independently for the dry season and the rainy season. For the interface flux of each sub-region, the specific calculation method is as follows: Darcy's law is used to calculate the exchange flux per unit area of ​​the lake-groundwater interaction interface in each sub-region. A lake-groundwater hydraulic constraint matrix is ​​constructed using the discharge zone, infiltration zone, and transition zone as matrix row elements, and hydraulic gradient, hydraulic exchange direction, interface flux per unit area, total interface flux, and exchange intensity level as matrix column elements. The exchange intensity level is divided into three levels: weak, medium, and strong, based on the absolute value of flux per unit area. The classification threshold is set to match the hydrogeological conditions of the sandy land. This constraint matrix serves as a unified hydrological constraint benchmark for subsequent hydraulic flux verification, isotope correction, and pollution source analysis.

4. The multi-isotope-based source apportionment system for sandy lake-groundwater pollution according to claim 3, characterized in that, The method used to construct the multi-isotope-water chemical fingerprint database is as follows: Using sub-region type, monitoring point number, sampling depth, and sampling time period as indexes, the isotope test results, water chemical parameters, and in-situ monitoring data of each sample are linked and stored one by one; the isotope-water chemical characteristic values ​​of typical pollution source end-members are simultaneously included to form a multi-isotope-water chemical fingerprint database covering water samples and pollution source end-members.

5. The multi-isotope-based source apportionment system for sandy lake-groundwater pollution according to claim 4, characterized in that, Local atmospheric precipitation lines and lake evaporation lines were established, and the groundwater inflow flux and lake water leakage in each sub-region were calculated using the radon-222 mass balance equation. The method used was as follows: Based on the constructed multi-isotope-water chemical fingerprint database, atmospheric precipitation hydrogen and oxygen isotope monitoring data for three complete hydrological years in the sub-region were obtained. The local atmospheric precipitation line was obtained by fitting oxygen isotope as the abscissa and hydrogen isotope as the ordinate. The test results of hydrogen and oxygen isotopes of lake water samples in each sub-region were extracted, and the lake water evaporation lines of the discharge area, infiltration area and transition area were fitted respectively. The evaporation enrichment degree of lake water in each sub-region was quantified by the difference in slope and intercept between the evaporation line and the atmospheric precipitation line. Radon-222 mass balance equations were constructed for the three types of sub-regions, with the lake units in each sub-region as independent calculation objects, to calculate the groundwater inflow flux and lake water leakage in each sub-region.

6. The multi-isotope-based source apportionment system for sandy lake-groundwater pollution according to claim 5, characterized in that, The flux results were verified and corrected, and a dynamic fractionation correction model was used to eliminate the isotope enrichment effect caused by strong evaporation in sandy areas, in order to obtain the original isotope data. The method used was as follows: The constructed lake-groundwater hydraulic constraint matrix is ​​retrieved, and the interface flux calculated based on Darcy's law for each sub-region is extracted as a reference. The sub-region flux calculated by radon-222 mass balance is compared with the reference value. If the relative deviation between the two exceeds the preset threshold, the gas exchange rate and the equivalent permeability coefficient of the interface are corrected by combining the sediment particle size distribution and sandy wind speed conditions of the sub-region. The flux is then recalculated iteratively until the deviation meets the preset threshold, and the verified and corrected sub-region hydraulic exchange flux result is obtained. Based on the evaporation fractionation model, and considering the strong evaporation characteristics of sandy areas with high temperature, low relative humidity, and high wind speed, local measured meteorological parameters are introduced to correct the dynamic fractionation term. Measured values ​​of hydrogen and oxygen isotopes in lake water, mean values ​​of atmospheric precipitation isotopes, near-surface air temperature, lake surface relative humidity, and lake surface wind speed parameters are input for each sub-region. The water balance fractionation coefficient and dynamic fractionation coefficient are calculated, and a calculation model for lake water evaporation enrichment corresponding to each sub-region is constructed. The isotopic composition of the remaining water body after evaporation satisfies the Rayleigh fractionation relation. Based on the verified and corrected hydraulic exchange flux of each sub-region, the evaporation residual ratio of lake water in each sub-region is calculated by combining the lake water balance equation and substituted into the dynamic fractionation correction model to obtain the original isotopic composition of the lake water replenishment source. The dynamic fractionation correction model is based on the evaporation fractionation model.

7. The multi-isotope-based sandy lake-groundwater pollution source apportionment system according to claim 6, characterized in that, The corrected original isotopic data were compared with the characteristic values ​​of pollution source end-member isotopes, and combined with water chemical ion ratios and land use information to identify pollution sources in each zone. The method used was as follows: Samples of five typical pollution sources in the target area, namely chemical fertilizer, human and animal excrement, domestic sewage, soil organic nitrogen, and mining wastewater, were collected and relevant indicators were tested. The characteristic value ranges of nitrate nitrogen-15, nitrate oxygen-18, and sulfate sulfur-34 isotopes of various pollution sources and the corresponding ion ratio ranges were established to form a standardized pollution source end-member characteristic set. The corrected water isotope data of each sub-region are compared with the end-member feature set. The matching degree is quantitatively determined by the isotope feature distance, and the molar concentration ratio of water feature ions is calculated to eliminate the interference of dilution and evaporation concentration. By combining the land use type and spatial distribution data of pollution sources corresponding to each sub-region, candidate sources with mismatched spatial distribution are eliminated, and the main pollution source types of each sub-region are finally determined.

8. The multi-isotope-based source apportionment system for sandy lake-groundwater pollution according to claim 7, characterized in that, The method used to correct the isotopic fractionation induced by biotransformation by adjusting the isotopic fractionation coefficient of microbial action, thereby restoring the initial isotopic composition of the pollution source, is as follows: By combining the measured data of dissolved oxygen and redox potential of water bodies in each sub-region, as well as the spatial decay law of nitrate and sulfate concentrations along the hydraulic gradient, it is determined whether denitrification and sulfate reduction microbial transformation processes occur in the corresponding zone; when the corresponding substrate concentration decreases along the path and the isotope value is enriched, it is determined that there is a biofractionation effect and the dominant reaction type is determined. The Rayleigh fractionation model was used to quantitatively correct the isotopic fractionation effect in the biotransformation process. Based on the correction formula and the isotopic values ​​of the water body after the reaction, the initial isotopic values ​​of the pollution source were inferred. In view of the hydrological environment characteristics of sandy land with high temperature and strong aeration of the vadose zone, in-situ culture experiments are preferred to determine the isotope enrichment coefficients of denitrification and sulfate reduction processes in the target area. When in-situ measurement is not possible, empirical values ​​from literature in sandy areas with the same climate and lithology are selected as substitutes. By substituting the measured isotope values, enrichment coefficients, and residual substrate fractions into the correction formula, the superimposed influence of biotransformation on isotope composition is eliminated, and the initial isotope characteristic dataset of the pollution sources corresponding to each sub-region is restored, which serves as the standardized endmember input parameters for subsequent quantitative source apportionment.

9. A multi-isotope-based source apportionment system for pollution sources in sandy lakes and groundwater according to claim 8, characterized in that, A Bayesian mixture model was used to calculate the pollution contribution rate of each pollution source. After cross-validation, the spatial distribution and contribution ratio of the pollution sources were output. The method used was as follows: The initial isotope values ​​of each pollution source endmember obtained by restoration are used as the input terms of the model endmembers, and the corrected isotope values ​​of each sub-region water body are used as the input terms of the mixed sample. The flux data of the partition interface in the lake-groundwater hydraulic constraint matrix are retrieved, and the flux ratio corresponding to each pollution migration path is converted into the prior distribution parameters of the Bayesian mixture model. Hydrological path hard constraints are applied to the range of values ​​of each pollution source contribution rate. The specific method for quantitatively solving the pollution contribution rate of each pollution source using the Bayesian mixture model is as follows: Based on the principle of isotope mass balance, a multi-isotope mixture equation is constructed, and the posterior probability distribution of the contribution rate of each pollution source is solved by Markov chain Monte Carlo sampling. The mean, standard deviation and 95% confidence interval of the contribution rate of each type of pollution source in each sub-region are output. Finally, the contribution rate results are combined with the flux data of the lake-groundwater hydraulic constraint matrix and substituted into the groundwater solute transport model to simulate the spatial distribution of pollutant concentration. The degree of fit between the simulated and measured values ​​is quantified by the coefficient of determination. The contribution ratio of each pollution source is statistically analyzed according to the discharge zone, infiltration zone, and transition zone. The final spatial distribution results of pollution sources and the contribution rate ranking list are output.

10. A method for source apportionment of pollution in sandy lakes and groundwater based on multiple isotopes, characterized in that, The method is performed by the multi-isotope-based sandy lake-groundwater pollution source apportionment system according to any one of claims 1-9, including: Step 1: Based on the hydrogeological data and dynamic water level monitoring data of the target area, combined with the complementary relationship between the groundwater flow field and the lake water level, three types of sub-regions are divided, water level data of each sub-region are collected, hydraulic gradient, exchange direction and interface flux are calculated, and a lake-groundwater hydraulic constraint matrix is ​​constructed. The sub-regions include discharge zone, infiltration zone and transition zone. Step 2: Collect lake water, groundwater, bottom sediment pore water and vadose zone soil solution samples for the three types of sub-regions, determine the hydrogen and oxygen, nitrate nitrogen and oxygen, sulfate sulfur isotopes and radon-222 activity, and simultaneously detect the anion and cation concentrations and total dissolved solids content to construct a multi-isotope-water chemical fingerprint database. Step 3: Based on hydrogen and oxygen isotope data, establish local atmospheric precipitation lines and regional lake evaporation lines. Calculate the groundwater inflow flux and lake water leakage in each region using the radon-222 mass balance equation. Verify and correct the flux results using the constraint matrix. Use a dynamic fractionation correction model to eliminate the strong evaporation isotope enrichment effect in sandy areas. Invert to obtain the original isotope data of the lake water recharge source. Step 4: Compare the corrected original isotope data with the characteristic values ​​of the pollution source end-member isotopes, identify the pollution sources in each zone by combining the water chemical ion ratio and land use information, correct the isotope fractionation caused by biotransformation by the isotope fractionation coefficient of microbial action, and restore the initial isotope composition of the pollution source. Step 5: Using the restored initial isotope values ​​of the pollution sources as input, and combining the flux weight constraints of the constraint matrix, a Bayesian mixture model is used to calculate the pollution contribution rate of each pollution source. After cross-validation, the spatial distribution and contribution ratio of the pollution sources are output.

Images