A numerical simulation method and system for karst groundwater solute transport process

Abstract

This invention discloses a numerical simulation method and system for solute transport processes in karst groundwater, relating to the field of groundwater numerical simulation technology. It includes a basic data input module, a three-domain solute transport embryo generation module, a solute transport embryo trial calculation and evaluation module, a live-flow exchanged page number generation module, and a page number-driven simulation module. The basic data input module acquires the basic data set AW; the three-domain solute transport embryo generation module generates a pipe-fracture syn-transport embryo B1, a fracture-base stagnant embryo B2, and a base-returning retrograde embryo B3; the solute transport embryo trial calculation and evaluation module acquires the candidate exchanged page number P and the trial-optimized page number Pr; the live-flow exchanged page number generation module generates the live-flow exchanged page number CER; and the page number-driven simulation module controls the mass transfer between pipes, fractures, and the matrix based on the live-flow exchanged page number CER, outputting the three-domain solute transport results for the current simulation period.

Description

Technical Field

[0001] This invention relates to the field of groundwater numerical simulation technology, specifically to a numerical simulation method and system for the solute transport process in karst groundwater. Background Technology

[0002] Existing numerical simulation methods for solute transport in karst groundwater often simplify the exchange process between conduits, fissures, and the matrix to a fixed exchange coefficient or a unified exchange term, directly substituting the current exchange coefficient into the solute transport equation. While this approach can describe some conventional transport processes, under different hydraulic conditions such as high-water seasons, low-water seasons, heavy rainfall recharge periods, and receding water periods, the same exchange coefficient cannot simultaneously describe different exchange actions such as conduit-to-fissure diversion, fissure-to-matrix retention, and matrix-to-fissure release. Especially when the outlet concentration curve shows a rapid arrival before the peak, a slow decline after the peak, and a rebound in the tail section, the fixed exchange treatment easily leads to the mixed interpretation of the peak, retention, and tail sections, resulting in insufficient clarity in the current model's determination of the exchange direction in the three domains.

[0003] The aforementioned shortcomings arise primarily from the continuous variations in hydraulic driving, concentration gradient driving, and residence time among karst conduits, fissures, and the matrix across different simulation periods. Existing methods fail to preprocess these variations as intermediate results with exchange direction implications, nor do they incorporate the outlet concentration and exchange status from the previous simulation period to assess the continuity of the current exchange process. Consequently, the model may exhibit abnormal results such as peak outlet concentration arrival time shifts, misinterpretations of tail-end delays, misclassification of matrix back-release as new pollution input, and misclassification of fissure retention as ordinary dispersion. These anomalies affect the assessment of karst groundwater pollution migration paths, pollution duration, and pollution remediation timelines. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a numerical simulation method and system for the solute transport process in karst groundwater, solving the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention is implemented through the following technical solution: a numerical simulation system for the solute transport process in karst groundwater, comprising a basic data input module, a three-domain mobile embryo generation module, a mobile embryo trial calculation and evaluation module, a process-driven page number generation module, and a page number-driven simulation module; The basic data input module obtains the basic data set AW based on the head sequence of pipeline calculation nodes, the head sequence of fracture calculation nodes, the discrete model of karst channels, the water balance results, the volume of the three-domain discrete unit and the mass calculation results of the previous simulation period. The three-domain migratory embryo generation module analyzes the basic data set AW and generates tube-crack antegrade embryo B1, cleavage-based stagnant embryo B2, and basal-returning retrograde embryo B3. The free embryo trial calculation and evaluation module sorts and performs quality calculations on the tube-cracked antegrade embryo B1, the fractured base-retained embryo B2, and the base-returning reverse-release embryo B3 to obtain the candidate exchange page order P and the trial calculation preferred page order Pr. The on-the-fly page number generation module generates the on-the-fly page number CER for the current simulation period by combining the candidate page number P and the trial-and-error preferred page number Pr with the quality calculation results of the previous simulation period. The fold code-driven simulation module controls the mass transfer between pipes, fissures, and matrix based on the on-the-fly fold code CER, and outputs the three-domain solute transport results for the current simulation period.

[0006] Preferably, the basic data input module includes a three-domain source quantity extraction unit and a data normalization and shaping unit; The three-domain source quantity extraction unit extracts three-domain source quantities within the same simulation period from the existing data of the current numerical simulation system, and organizes data from different sources into the original basic quantities under the same period. Using the current simulation time period tn as the index, read the upstream node head Hc1 and downstream node head Hc2 in the pipeline calculation node head sequence, and read the fracture inlet node head Hf1 and fracture outlet node head Hf2 in the fracture calculation node head sequence. Simultaneously, the pipe section length Lc, fissure section length Lf, and pipe water-passing area Ac are read from the discrete model of the karst channel; the pipe flow rate Qc, fissure outflow rate Qfout, and fissure residence time Tf are read from the water balance results of the previous simulation period; the pipe water volume Wc, fissure water volume Wf, and matrix exchangeable water volume Wm are read from the three source quantities; and the pipe solute mass Mc, fissure solute mass Mf, and matrix solute mass Mm are read from the mass calculation results of the previous simulation period. Align the acquired data: if the data comes from the previous simulation period, then use it as the start state of the simulation period tn; When the data comes from the current discrete model, the data is used as a fixed structural quantity for the simulation time period tn. Based on the aligned data, the following parameters are calculated: pipe hydraulic gradient Jc, fracture hydraulic gradient Jf, pipe velocity Vc, pipe concentration Cc, fracture concentration Cf, matrix concentration Cm, pipe-to-fracture concentration difference Gcf, and fracture-to-matrix concentration difference Gfm; the formulas are as follows: ; ; ; ; ; ; ; ; The data normalization shaping unit processes the pipe hydraulic gradient Jc, fracture hydraulic gradient Jf, pipe velocity Vc, pipe concentration Cc, fracture concentration Cf, fracture residence time Tf, matrix concentration Cm, pipe-to-fracture concentration difference Gcf, and fracture-to-matrix concentration difference Gfm using the maximum-minimum normalization method to obtain normalized data, and then fits the data to obtain the basic data set AW.

[0007] Preferably, the three-domain mobile embryo generation module includes a three-domain directional mobile fragmentation unit and a mobile embryo compression forming unit; The three-domain directional wandering and fragmentation unit performs directional fragmentation processing on the data in the basic data group AW. Read the state of the concentration difference Gcf between the pipeline and the fracture. When the concentration difference Gcf between the pipeline and the fracture is positive, it means that the concentration in the pipeline is higher than the concentration in the fracture, and there is a concentration-driven phenomenon of solute being released from the pipeline to the fracture side during the current period. When the concentration difference Gcf between the pipeline and the fracture is not positive, the concentration push from the pipeline to the fracture side does not enter the fracture direction plate; Read the sign of the concentration difference Gfm between the fracture and the matrix. When the concentration difference Gfm between the fracture and the matrix is ​​positive, it means that the fracture concentration is higher than the matrix concentration, and there is a concentration drive of solute from the fracture to the matrix side in the current period. When the concentration difference between the fracture and the matrix, Gfm, is negative, it indicates that the matrix concentration is higher than the fracture concentration, and there is a concentration-driven phenomenon of solute returning from the matrix to the fracture side during the current period. The pipe hydraulic gradient Jc and pipe velocity Vc are used as the driving forces for pipe transport, the fracture hydraulic gradient Jf is used as the driving force for the entry of the fracture side into the matrix, and the fracture residence time Tf is used as the time support for the temporary storage or re-release process of the residence action in the fracture. Three directional slices are obtained, including pipe fracture directional slice D1, fracture-to-matrix directional slice D2, and matrix-to-return directional slice D3, as shown in the following formula: ; ; ; In the formula, This indicates taking the maximum value.

[0008] Preferably, the free embryo compression forming unit is based on the basic data set AW, and is analyzed in combination with the tube fracture direction sheet D1, the fracture base direction sheet D2, and the base return direction sheet D3; First, analyze whether the pipe fracture direction segment D1 has fast flow support in the pipe; When the pipe fracture direction segment D1 has a value, the pipe flow velocity Vc is again used as the forward traction term, so that the pipe fracture forward traction segment B1 is more inclined to reflect the lateral flow under the condition of high-speed pipe transport, rather than simply reflecting the concentration difference between the pipe and the fracture. Analyze whether the fracture-oriented sheet D2 has fracture retention support; When the fracture orientation segment D2 has a value, the fracture residence time Tf is used as the storage compression term, so that the fracture storage embryo B2 is more inclined to reflect the process of solute entering the matrix after staying in the fracture. Analyze whether the base return direction segment D3 has a dwell return support; When the base-return direction sheet D3 has a value, the crack residence time Tf is used as the return support term, so that the base-return reverse release B3 is more inclined to reflect the process of solutes already stored in the matrix returning through the crack. Finally, the formulas for the output tube-cracked antegrade embryo B1, the cleaved basal residing embryo B2, and the basal reversal-released embryo B3 are as follows: ; ; .

[0009] Preferably, the free embryo trial calculation and evaluation module includes an embryo value page sequence arrangement unit and a quality trial calculation verification unit; The embryo value page ordering unit obtains candidate exchange page order P based on the values ​​of tube-fission antegrade embryo B1, fission-based stagnant embryo B2, and basal-retrograde release embryo B3. The process is as follows: Embryo B1 (transient from tube), embryo B2 (retained at the base), and embryo B3 (reverse release from the base) are placed into the same embryo value set, and the total amount S of the three embryos is calculated using the following formula: ; When the total amount of the three embryos S is not zero, calculate the occupancy values ​​of the anterograde embryo B1, the cleavage-based residing embryo B2, and the basal-returning embryo B3 in the total amount of the three embryos S, respectively, to obtain the occupancy values ​​N1, N2, and N3 of the anterograde embryo, as follows: ; ; ; The pages are arranged in order of magnitude: pipe fracture occupancy value N1, fracture matrix occupancy value N2, and matrix reversion occupancy value N3. When the pipe crack occupancy value N1 is in the first position, the candidate swap page order P is marked as the pipe crack forward page order; When the fractured occupancy value N2 is first, the candidate exchange page order P is marked as the fractured stagnant page order. When the base-return placeholder value N3 is in the first position, the candidate swap page order P is marked as the base-return reverse page order.

[0010] If the two placeholder values ​​are the same, the current cell reads the source domain mass used by the two free embryos; for example, the tube fracture syn-transfer embryo B1 reads the tube solute mass Mc, the fracture matrix stagnant embryo B2 reads the fracture solute mass Mf, and the matrix retrograde embryo B3 reads the matrix solute mass Mm, and the free embryo with the larger source domain mass is placed first; if the sum of the three embryos S is zero, then the candidate swap page order P is the masterless free page order.

[0011] Preferably, the quality trial calculation verification unit performs a short-step quality trial calculation on each of the following: tube-cracked antegrade embryo B1, cleaved-base stagnant embryo B2, and base-returning retrograde embryo B3, based on the candidate exchange page order P, and obtains the preferred trial calculation page order Pr. The process is as follows: Read the following data: pipe fracture syn-migration embryo B1, fracture-base stagnant embryo B2, and matrix reversion embryo B3; and the solute mass Mc of the pipe, solute mass Mf of the fracture, and solute mass Mm of the matrix from the mass calculation results of the previous simulation period. Simultaneously, read the water volume Wc of the pipe, the water volume Wf of the fracture, and the exchangeable water volume Wm of the matrix. Calculate the following mass from the pipe to the fracture: Q1; from the fracture to the matrix: Q2; and from the matrix to the fracture: Q3. The formulas are as follows: ; ; ; In the formula, This indicates taking the minimum value. This represents the pipe fracture lateral migration conversion factor, with a value ranging from 0 to 1. This represents the conversion coefficient for fractured substrate storage, with a value ranging from 0 to 1. This represents the base-retroactive conversion coefficient, with a value ranging from 0 to 1; Among them, the pipe fracture lateral migration conversion coefficient Fission-induced storage conversion coefficient and base-retroactive conversion coefficient The acquisition method is as follows: first, read the tracer delivery quality, delivery time, measured sequence of outlet concentration, sequence of outlet flow rate, and sequence of historical simulation periods from the historical tracer experiment; at the same time, read the basic data group AW that has been formed in each historical simulation period, and regenerate the tube fracture syn-migration embryo B1, fracture-based stagnant embryo B2, and base-returning retrograde embryo B3 for each historical simulation period according to the processing procedure of the three-domain migratory embryo generation module.

[0012] Subsequently, based on the measured sequence of outlet concentration and the outlet flow rate, the system calculates the tracer mass discharged at the outlet for each historical simulation period, and divides the historical simulation periods into a rapid arrival phase before the peak, a gradual decline phase after the peak, and a delayed tail phase according to the shape of the outlet concentration curve. The rapid arrival phase before the peak is used to back-calculate the pipe fracture lateral migration conversion coefficient. The post-peak descent phase is used to back-calculate the conversion coefficient of fractured and stagnant reservoirs. The tail-end delay segment or tail-end rebound segment is used to back-calculate the base-reverse release conversion coefficient. .

[0013] Conversion coefficient for pipe fracture migration During the acquisition process, the system extracts the fracture acclimatization embryo B1 and the solute mass Mc of the pipeline within the rapidly arriving segment before the peak, and counts the tracer mass discharged from the outlet within the current segment. The tracer mass discharged from the outlet is compared with the fracture acclimatization base formed by the fracture acclimatization embryo B1 and the solute mass Mc of the pipeline to obtain the fracture acclimatization conversion coefficient. .

[0014] Conversion coefficient of fractured storage During the acquisition process, the system extracts the fractured embryo B2 and fracture solute mass Mf in the post-peak descent section, and statistically analyzes the delayed change in outlet mass caused by fracture retention in the current section. The delayed change in mass is compared with the amount of fractured embryo B2 and fracture solute mass Mf to obtain the fractured retention conversion coefficient. .

[0015] In the base-retroreversion conversion coefficient During the acquisition process, the system extracts the base-returning embryo B3 and matrix solute mass Mm within the tail-end delay segment or tail-end recovery segment, and statistically analyzes the change in outlet mass caused by matrix return within the current segment. The returned mass is compared with the amount of base-returning base formed by the base-returning embryo B3 and matrix solute mass Mm to obtain the base-returning conversion coefficient. .

[0016] The obtained mass values ​​Q1 from the pipe to the fracture, Q2 from the fracture to the matrix, and Q3 from the matrix to the fracture are analyzed. The residual values ​​R1, R2, and R3 of the concentration difference between the pipe and fracture are calculated, using the following formulas: Compare the values ​​of the residual values ​​R1 (difference between pipe concentration and fracture concentration), R2 (difference between fracture concentration and matrix concentration), and R3 (difference between fracture concentration and matrix concentration). When the residual value R1 of the difference between the concentration in the pipeline and the concentration in the crack is the smallest, the optimal page sequence Pr is determined by trial calculation and the pipeline crack shift page sequence is taken. When the residual value R2 of the difference between fracture concentration and matrix concentration is the minimum, the optimal page sequence Pr is calculated and the page sequence with fractured matrix stagnant storage is selected. When the residual value R3 of the difference between the fracture concentration and the matrix concentration is the minimum, the optimal page sequence Pr is calculated and the base reversion page sequence is selected.

[0017] When two residual values ​​are the same, the page order that is consistent with the candidate swap page order P is preferred; if P is an unclaimed wandering page order, the page order that makes the short-step change in the outlet concentration smaller is adopted as Pr. Preferably, the on-the-fly fold code generation module includes a page sequence homomorphism verification unit and a short-step convergence unit; First, read the page order type of the candidate swap page order P. The page order types of the candidate swap page order P include tube fracture forward migration page order, fracture base stagnant page order, base reversion reverse release page order, and unclaimed wandering page order. Then read the page order type of the trial-optimized page order Pr. The page order types of Pr include tube fracture forward migration page order, fracture base stagnant storage page order, and base reversion reverse release page order. Perform homeomorphic verification between the candidate swap page order P and the trial-optimized page order Pr: If both the candidate swap page order P and the trial-optimized page order Pr are pipe fracture forward shift page orders, then a pipe fracture diversion page seed is generated. When both the candidate swap page order P and the trial-optimized page order Pr are split-base stagnant page orders, a split-base stagnant page seed is generated. When both the candidate swap page order P and the trial-optimized page order Pr are basis-reverse page orders, a basis-reverse page seed is generated. When the candidate swap page order P is an unclaimed wandering page order, and the trial-optimized page order Pr has a definite page order, then the trial-optimized page order Pr is used to form a page seed, and the page seed is marked as a trial takeover state. If the candidate swap page order P and the trial-preferred page order Pr point to different page orders, two candidate page folding seeds are formed: the tendency page folding seed formed by the candidate swap page order P and the trial-preferred page folding seed formed by the trial-preferred page order Pr.

[0018] Preferably, the short-step convergence unit determines the number of seeds; If there is only one fold seed, the current fold seed is converted into a random exchange fold code CER. Specifically, the tube fracture shunt fold seed is converted into a tube fracture shunt fold code, the fracture basal storage fold seed is converted into a fracture basal storage fold code, and the basis reverse release fold seed is converted into a basis reverse release fold code. If there are two fold seeds, then perform short-step inheritance trial calculations on each of the two fold seeds; The short-step trial calculation reads the solute mass Mc in the pipeline, solute mass Mf in the fracture, and solute mass Mm in the matrix from the mass calculation results of the previous simulation period, as well as the water volume Wc in the pipeline, water volume Wf in the fracture, and exchangeable water volume Wm in the matrix in the three-domain discrete element volume. When the fold seed is a pipe fracture diversion fold seed, the pipe solute mass Mc and the fracture solute mass Mf are updated according to the transfer from the pipe to the fracture. When the folded seed is a cleaved folded seed, the mass of solute in the cleavage (Mf) and the mass of solute in the matrix (Mm) are updated according to the transfer from the cleavage to the matrix. When the folded seed is a base-returning folded seed, the matrix solute mass Mm and Mf and the crack solute mass Mf are updated according to the matrix returning to the crack; After completing the short-step acceptance trial calculation, the predicted short-step outlet concentration value of each fold seed is obtained. And calculate the export concentration compared to the previous simulation period. The formula for the difference in acceptance between them, Ak, is: ; In the formula, k represents any fold seed currently participating in the convergence; Read the page break code C0 from the previous simulation period and analyze its connection status with the current page break seed: If the current fold seed is the same as the fold code C0 of the previous simulation period, then the fold capacity Fk is 0; If the current fold seed continues along the sequence of tube fissure diversion, fissure base retention, and base reversion, then the fold receiving amount Fk is 1; When the current fold seed and the fold code C0 of the previous simulation period form a reverse jump, the fold carrying capacity Fk is 2; The folding convergence value Zk is synthesized from the export concentration difference Ak and the folding acceptance amount Fk, as shown in the following formula: ; For the fold convergence values ​​Zk of the two fold seeds, the fold seed with the smaller fold convergence value Zk is used to generate the shifted fold code CER for the current simulation period.

[0019] Preferably, the fold code driven simulation module includes a fold code directional wandering unit and a three-domain quality write-back output unit; The page-code directional migration unit selects the unique migration embryo that enters the mass transfer calculation during the current simulation period from the tube fracture antegrade embryo B1, fracture-based stagnant embryo B2, and base-returning reverse release embryo B3, based on the process-exchange page code CER. Parsing the code pattern of the CER (Continuous Exchange Page Code): When the page number CER is the pipe fracture diversion page number, only the pipe fracture forward migration embryo B1 is allowed to participate in the calculation in the current time period. The pipe is used as the source domain and the fracture is used as the receiving domain to calculate and obtain the mass shift QER. When the time-travel exchange page number CER is the fractured stagnant page number, only the fractured stagnant embryo B2 is allowed to participate in the calculation in the current time period, and the fracture is used as the source domain and the matrix is ​​used as the receiving domain to calculate and obtain the mass migration QER. When the page number CER is the base-reverse release page number, only the base-reverse release embryo B3 is allowed to participate in the calculation in the current time period, and the matrix is ​​used as the source domain and the crack is used as the receiving domain to calculate and obtain the mass transfer quantity QER. The formula is as follows: ; ; ; The source domain is selected by the time-travel exchange fold code CER, and the selection is limited by the existing mass of the source domain; when the source domain is a pipe, the time-travel exchange fold code CER shall not exceed the pipe solute mass Mc; when the source domain is a fracture, the time-travel exchange fold code CER shall not exceed the fracture solute mass Mf; when the source domain is a matrix, the time-travel exchange fold code CER shall not exceed the matrix solute mass Mm. The three-domain mass write-back outflow unit writes back the solute mass Mc in the pipeline, the solute mass Mf in the fissure, and the solute mass Mm in the matrix based on the time-of-flight exchange (CER) page number, and outputs the three-domain solute transport results for the current simulation period. Three mutually exclusive write-back actions are performed according to the CER (Cross-Process Interchange Page Number): When the time-of-flight page number CER is the pipe fracture diversion page number, the mass shift QER is subtracted from the pipe solute mass Mc, and the mass shift QER is incorporated into the fracture solute mass Mf. The matrix solute mass Mm remains the same as at the beginning of the current simulation period. When the time-travel page number CER is the fractured matrix stagnant page number, the mass shift QER is subtracted from the fracture solute mass Mf, and the mass shift QER is incorporated into the matrix solute mass Mm. The pipeline solute mass Mc remains the mass at the beginning of the current simulation period. When the time-travel page number CER is the base-reverse release page number, the mass shift QER is subtracted from the matrix solute mass Mm, and the mass shift QER is incorporated into the fracture solute mass Mf. The pipeline solute mass Mc remains the mass at the beginning of the current simulation period. When the page number CER during the process exchange is the pipe splitting shunt page number: ; ; ; When the page number CER during the process exchange is a cleavage-dependent page number: ; ; ; When the page number CER is a base-reverse page number: ; ; ; In the formula, Indicates the mass of solute in the new pipeline. Indicates the mass of solute in the new fracture. Indicates the mass of the solute in the new matrix; Based on the solute quality of the new pipeline quality of solute in new fractures and the quality of the new matrix solute The pipe concentration Cc, fracture concentration Cf, and matrix concentration Cm at the end of the current simulation period are processed to obtain new pipe concentration nCc, new fracture concentration nCf, and new matrix concentration nCm, as shown in the following formula: ; ; ; The new pipeline concentration nCc is used as the inlet concentration for the main pipeline transport calculation, and the outlet concentration Ct for the current simulation period is obtained based on the quality update result of the discrete element at the end of the pipeline. At the same time, the new fracture concentration nCf is sent to the fracture side calculation for the next simulation period, and the new matrix concentration nCm is sent to the matrix storage or matrix back-release calculation for the next simulation period.

[0020] A numerical simulation method for solute transport processes in karst groundwater includes the following steps: Step 1: The basic data input module obtains the basic data set AW based on the pipeline calculation node head sequence, fracture calculation node head sequence, karst channel discrete model, water balance results, three-domain discrete element volume, and mass calculation results of the previous simulation period. Step 2: The three-domain migratory embryo generation module analyzes the basic data set AW to generate tube-crack syn-migratory embryo B1, cleavage-based stagnant embryo B2, and basal-retrograde release embryo B3. Step 3: The free embryo trial calculation and evaluation module sorts and performs quality calculations on the tube-cracked antegrade embryo B1, the fractured base-retained embryo B2, and the base-returning reverse-release embryo B3 to obtain the candidate exchange page order P and the trial calculation preferred page order Pr; Step 4: The on-the-fly page number generation module generates the on-the-fly page number CER for the current simulation period by combining the candidate page number P and the trial-and-error preferred page number Pr with the quality calculation results of the previous simulation period. Step 5: The fold code-driven simulation module controls the mass transfer between pipes, fissures, and matrix based on the live fold code CER, and outputs the three-domain solute transport results for the current simulation period.

[0021] This invention provides a numerical simulation method and system for solute transport processes in karst groundwater, which has the following advantages: (1) During system operation, the floating embryo trial calculation and evaluation module sorts and performs quality trials on the tube fracture antegrade embryo B1, fracture-based stagnant embryo B2, and base-returning retrograde embryo B3 to obtain the candidate exchange page order P and the trial-preferred page order Pr. Then, the follow-up exchange page code generation module combines the follow-up exchange page code of the previous simulation period and the outlet concentration of the previous simulation period to generate the follow-up exchange page code CER for the current simulation period. The current processing makes CER reflect both the floating intensity of the current simulation period and the concentration dissipation state after quality trial calculation. It also carries over the outlet concentration change of the previous simulation period, which can alleviate the problem of frequent jumps in exchange actions between adjacent simulation periods and make the concentration peak segment, stagnation segment and tail segment have a clearer connection in time.

[0022] (2) Based on the directional slice, the floating embryo compression forming unit uses the pipe flow velocity Vc as the forward traction term of the pipe fracture forward-moving embryo B1, the fracture residence time Tf as the storage compression term of the fracture base stagnant embryo B2, and the fracture residence time Tf and fracture hydraulic gradient Jf as the basis for judging the return of the base-returning reverse-release embryo B3. Thus, the pipe fracture forward-moving embryo B1, the fracture base stagnant embryo B2, and the base-returning reverse-release embryo B3 respectively reflect the three types of exchange actions in the current simulation period, providing clearer intermediate results for the generation of subsequent candidate exchange page order, trial calculation optimal page order, and follow-up exchange page number, which is helpful in explaining the phenomena such as rapid arrival before the outlet concentration peak, slow decline after the peak, and tail-end return rise.

[0023] (3) After obtaining the mass Q1 from the pipeline to the fracture, the mass Q2 from the fracture to the matrix, and the mass Q3 from the matrix to the fracture, the mass trial calculation verification unit further calculates the residual value R1 of the difference between the pipeline concentration and the fracture concentration, the residual value R2 of the difference between the fracture concentration and the matrix concentration, and the residual value R3 of the difference between the fracture concentration and the matrix concentration after matrix back-release. Based on the magnitude of the residual value, the trial calculation optimal page sequence Pr is obtained. Therefore, the exchange judgment in the current simulation period not only depends on the values ​​of the pipeline fracture antegrade embryo B1, the fracture base stagnant embryo B2, and the base back-release embryo B3, but also further examines the dissipation of the concentration difference in the three domains after each free embryo enters the mass transfer. This is helpful in distinguishing the rapid antegrade migration of the peak segment, the fracture base stagnant storage in the post-peak segment, and the matrix back-release in the tail segment.

[0024] (4) When the candidate swap page order P and the trial-preferred page order Pr point to the same page order, the corresponding page seed is directly generated; when the selected swap page order P is a masterless wandering page order and the trial-preferred page order Pr has a definite page order, the trial-preferred page order Pr forms the trial takeover state; when the candidate swap page order P and the trial-preferred page order Pr point to different page orders, both the tendency page seed and the trial-preferred page seed are retained. Through the above processing, the current system does not forcibly adopt a certain swap direction in the divergence state, but transfers the divergence state to the receiving short-step convergence unit for continued analysis, which is beneficial for handling simulation periods where the three-domain swap direction is easy to change, such as the rainstorm replenishment period and the receding period. Attached Figure Description

[0025] Figure 1 This is a schematic flowchart of a numerical simulation system for solute transport in karst groundwater according to the present invention. Figure 2 This is a schematic diagram illustrating the steps of a numerical simulation method for solute transport in karst groundwater according to the present invention. Figure 3 This is a schematic diagram of the candidate page order P determination process of the present invention; Figure 4 This is a schematic diagram of the process for determining the folding capacity Fk of the present invention. Detailed Implementation

[0026] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0027] Example 1 This invention provides a numerical simulation system for solute transport processes in karst groundwater. Please refer to [link / reference]. Figures 1 to 4 It includes a basic data input module, a three-domain wandering embryo generation module, a wandering embryo trial calculation and evaluation module, a process-interchange folded number generation module, and a folded number-driven simulation module; The basic data input module obtains the basic data set AW based on the head sequence of pipeline calculation nodes, the head sequence of fracture calculation nodes, the discrete model of karst channels, the water balance results, the volume of the three-domain discrete unit and the mass calculation results of the previous simulation period. The three-domain migratory embryo generation module analyzes the basic data set AW and generates tube-crack antegrade embryo B1, cleavage-based stagnant embryo B2, and basal-returning retrograde embryo B3. The free embryo trial calculation and evaluation module sorts and performs quality calculations on the tube-cracked antegrade embryo B1, the fractured base-retained embryo B2, and the base-returning reverse-release embryo B3 to obtain the candidate exchange page order P and the trial calculation preferred page order Pr. The on-the-fly page number generation module generates the on-the-fly page number CER for the current simulation period by combining the candidate page number P and the trial-and-error preferred page number Pr with the quality calculation results of the previous simulation period. The fold code-driven simulation module controls the mass transfer between pipes, fissures, and matrix based on the on-the-fly fold code CER, and outputs the three-domain solute transport results for the current simulation period.

[0028] In this embodiment, the basic data input module first organizes the head sequences of pipeline calculation nodes, fracture calculation nodes, karst channel discrete models, water balance results, three-domain discrete unit volumes, and mass calculation results from the previous simulation period into the same simulation period, obtaining the basic data set AW. This ensures that the hydraulic, concentration, and mass states in the pipeline, fracture, and matrix are at the same data scale before entering subsequent calculations. Compared to directly calling scattered data for exchange coefficient calculation, this current processing provides a unified data entry point for the three-domain exchange judgment, reducing the bias in solute transport direction judgment caused by inconsistent data source periods.

[0029] The three-domain migration embryo generation module processes the basic data set AW into three embryos: pipe-fracture syn-migration embryo B1, fracture-matrix stagnant embryo B2, and matrix-fracture back-release embryo B3. This allows the three types of exchange processes—pipe-to-fracture diversion, fracture-to-matrix stagnant storage, and matrix-to-fracture back-release—to be entered into independent intermediate results. Compared to methods using fixed exchange coefficients or single exchange terms, the current processing can distinguish different concentration curve morphologies such as rapid pre-peak arrival, slow post-peak decline, and tail-end back-release, preventing karst groundwater solute transport simulation from mixing different three-domain exchange processes into a single transport process.

[0030] The simulation module, driven by fold-codes, controls mass transfer between pipes, fissures, and the matrix based on the live fold-code exchange (CER) system. It outputs the three-domain solute transport results for the current simulation period, ensuring that the final output includes not only the outlet concentration curve but also the exchange actions, mass transfer direction, and mass changes in the three domains during the current period. This processing helps determine whether pollutants migrate rapidly along pipes, enter fissures, and are temporarily stored in the matrix, or are released back from the matrix to the fissures. This provides clearer numerical data for analyzing the migration paths, duration, and subsequent remediation periods of karst groundwater pollution.

[0031] Example 2 Please refer to Figure 1 Specifically: the basic data input module includes a three-domain source quantity extraction unit and a data normalization and shaping unit; The three-domain source quantity extraction unit extracts three-domain source quantities within the same simulation period from the existing data of the current numerical simulation system, and organizes data from different sources into the original basic quantities under the same period. Using the current simulation time period tn as the index, read the upstream node head Hc1 and downstream node head Hc2 in the pipeline calculation node head sequence, and read the fracture inlet node head Hf1 and fracture outlet node head Hf2 in the fracture calculation node head sequence. Simultaneously, the pipe section length Lc, fissure section length Lf, and pipe water-passing area Ac are read from the discrete model of the karst channel; the pipe flow rate Qc, fissure outflow rate Qfout, and fissure residence time Tf are read from the water balance results of the previous simulation period; the pipe water volume Wc, fissure water volume Wf, and matrix exchangeable water volume Wm are read from the three source quantities; and the pipe solute mass Mc, fissure solute mass Mf, and matrix solute mass Mm are read from the mass calculation results of the previous simulation period. Align the acquired data: if the data comes from the previous simulation period, then use it as the start state of the simulation period tn; When the data comes from the current discrete model, the data is used as a fixed structural quantity for the simulation time period tn. Based on the aligned data, the following parameters are calculated: pipe hydraulic gradient Jc, fracture hydraulic gradient Jf, pipe velocity Vc, pipe concentration Cc, fracture concentration Cf, matrix concentration Cm, pipe-to-fracture concentration difference Gcf, and fracture-to-matrix concentration difference Gfm; the formulas are as follows: ; ; ; ; ; ; ; ; The data normalization shaping unit processes the pipe hydraulic gradient Jc, fracture hydraulic gradient Jf, pipe velocity Vc, pipe concentration Cc, fracture concentration Cf, fracture residence time Tf, matrix concentration Cm, pipe-to-fracture concentration difference Gcf, and fracture-to-matrix concentration difference Gfm using the maximum-minimum normalization method to obtain normalized data, and then fits the data to obtain the basic data set AW.

[0032] The three-domain free embryo generation module includes a three-domain directional free embryo fragmentation unit and a free embryo compression forming unit; The three-domain directional wandering and fragmentation unit performs directional fragmentation processing on the data in the basic data group AW. Read the state of the concentration difference Gcf between the pipeline and the fracture. When the concentration difference Gcf between the pipeline and the fracture is positive, it means that the concentration in the pipeline is higher than the concentration in the fracture, and there is a concentration-driven phenomenon of solute being released from the pipeline to the fracture side during the current period. When the concentration difference Gcf between the pipeline and the fracture is not positive, the concentration push from the pipeline to the fracture side does not enter the fracture direction plate; Read the sign of the concentration difference Gfm between the fracture and the matrix. When the concentration difference Gfm between the fracture and the matrix is ​​positive, it means that the fracture concentration is higher than the matrix concentration, and there is a concentration drive of solute from the fracture to the matrix side in the current period. When the concentration difference between the fracture and the matrix, Gfm, is negative, it indicates that the matrix concentration is higher than the fracture concentration, and there is a concentration-driven phenomenon of solute returning from the matrix to the fracture side during the current period. The pipe hydraulic gradient Jc and pipe velocity Vc are used as the driving forces for pipe transport, the fracture hydraulic gradient Jf is used as the driving force for the entry of the fracture side into the matrix, and the fracture residence time Tf is used as the time support for the temporary storage or re-release process of the residence action in the fracture. Three directional slices are obtained, including pipe fracture directional slice D1, fracture-to-matrix directional slice D2, and matrix-to-return directional slice D3, as shown in the following formula: ; ; ; In the formula, This indicates taking the maximum value.

[0033] The free embryo compression forming unit is based on the basic data set AW, and is analyzed in combination with the tube fracture direction sheet D1, the fracture base direction sheet D2, and the base return direction sheet D3. First, analyze whether the pipe fracture direction segment D1 has fast flow support in the pipe; When the pipe fracture direction segment D1 has a value, the pipe flow velocity Vc is again used as the forward traction term. Analyze whether the fracture-oriented sheet D2 has fracture retention support; When the fracture orientation segment D2 has a value, the fracture residence time Tf is used as the storage compression term; Analyze whether the base return direction segment D3 has a dwell return support; When the base return direction piece D3 has a value, the crack residence time Tf is used as the return support item; Finally, the formulas for the output tube-cracked antegrade embryo B1, the cleaved basal residing embryo B2, and the basal reversal-released embryo B3 are as follows: ; ; .

[0034] In this embodiment, the basic data input module uses a three-domain source quantity extraction unit to process the head sequence of pipeline calculation nodes, the head sequence of fracture calculation nodes, the discrete model of karst channels, the water balance results, the three-domain source quantities, and the mass calculation results of the previous simulation period into the same simulation period. This ensures that the head, structural quantities, flow rates, volume, and mass quantities have a common time period basis before entering subsequent processing. As a result, the sources of the pipeline hydraulic gradient Jc, fracture hydraulic gradient Jf, pipeline velocity Vc, three-domain concentration, and three-domain concentration difference are all clear, which can alleviate the problem of unclear exchange direction due to the direct mixing of data from different sources in existing models.

[0035] The data normalization and shaping unit normalizes the pipe hydraulic gradient Jc, fracture hydraulic gradient Jf, pipe velocity Vc, pipe concentration Cc, fracture concentration Cf, fracture residence time Tf, matrix concentration Cm, pipe-to-fracture concentration difference Gcf, and fracture-to-matrix concentration difference Gfm, forming a basic data set AW. This allows data of different dimensions to participate in three-domain exchange analysis at the same computational scale. The current processing ensures that hydraulic drive, concentration drive, and fracture residence actions are no longer entered into the model as discrete parameters, but rather as calculable quantities in the basic data set AW, which is beneficial for distinguishing different exchange states during high-water, low-water, and receding water periods.

[0036] The three-domain directional migration segmentation unit, based on the sign states of the pipe-to-fracture concentration difference Gcf and the fracture-to-matrix concentration difference Gfm, segments three directions: pipe fracture release, fracture-to-matrix transfer, and matrix return. It then combines the pipe hydraulic gradient Jc, pipe velocity Vc, fracture hydraulic gradient Jf, and fracture residence time Tf to obtain pipe fracture direction segment D1, fracture-to-matrix direction segment D2, and matrix return direction segment D3. This current processing prevents the three-domain exchange process from being represented by a fixed exchange coefficient in a general way. Instead, it first breaks it down into intermediate segments with directional meaning, allowing rapid pipe diversion, fracture retention, and matrix return to be separately included in subsequent simulation analyses.

[0037] Example 3 Please refer to Figure 3 Specifically: the free embryo trial calculation and evaluation module includes an embryo value page sequence arrangement unit and a quality trial calculation and verification unit; The embryo value page ordering unit obtains candidate exchange page order P based on the values ​​of tube-fission antegrade embryo B1, fission-based stagnant embryo B2, and basal-retrograde release embryo B3. The process is as follows: Embryo B1 (transient from tube), embryo B2 (retained at the base), and embryo B3 (reverse release from the base) are placed into the same embryo value set, and the total amount S of the three embryos is calculated using the following formula: ; When the total amount of the three embryos S is not zero, calculate the occupancy values ​​of the anterograde embryo B1, the cleavage-based residing embryo B2, and the basal-returning embryo B3 in the total amount of the three embryos S, respectively, to obtain the occupancy values ​​N1, N2, and N3 of the anterograde embryo, as follows: ; ; ; The pages are arranged in order of magnitude: pipe fracture occupancy value N1, fracture matrix occupancy value N2, and matrix reversion occupancy value N3. When the pipe crack occupancy value N1 is in the first position, the candidate swap page order P is marked as the pipe crack forward page order; When the fractured occupancy value N2 is first, the candidate exchange page order P is marked as the fractured stagnant page order. When the base-return placeholder value N3 is in the first position, the candidate swap page order P is marked as the base-return reverse page order.

[0038] The quality trial calculation verification unit performs a short-step quality trial calculation on each of the following embryos: tube-cracked antegrade embryo B1, cleaved base-retained embryo B2, and base-returning retrograde embryo B3, based on the candidate exchange page order P, and obtains the optimal page order Pr. The process is as follows: Read the following data: pipe fracture syn-migration embryo B1, fracture-base stagnant embryo B2, and matrix reversion embryo B3; and the solute mass Mc of the pipe, solute mass Mf of the fracture, and solute mass Mm of the matrix from the mass calculation results of the previous simulation period. Simultaneously, read the water volume Wc of the pipe, the water volume Wf of the fracture, and the exchangeable water volume Wm of the matrix. Calculate the following mass from the pipe to the fracture: Q1; from the fracture to the matrix: Q2; and from the matrix to the fracture: Q3. The formulas are as follows: ; ; ; In the formula, This indicates taking the minimum value. Indicates the pipe fracture lateral migration conversion factor. Indicates the conversion coefficient of fractured substrate storage. Indicates the base-retroactive conversion coefficient; The obtained mass values ​​Q1 from the pipe to the fracture, Q2 from the fracture to the matrix, and Q3 from the matrix to the fracture are analyzed. The residual values ​​R1, R2, and R3 of the concentration difference between the pipe and fracture are calculated, using the following formulas: Compare the values ​​of the residual values ​​R1 (difference between pipe concentration and fracture concentration), R2 (difference between fracture concentration and matrix concentration), and R3 (difference between fracture concentration and matrix concentration). When the residual value R1 of the difference between the concentration in the pipeline and the concentration in the crack is the smallest, the optimal page sequence Pr is determined by trial calculation and the pipeline crack shift page sequence is taken. When the residual value R2 of the difference between fracture concentration and matrix concentration is the minimum, the optimal page sequence Pr is calculated and the page sequence with fractured matrix stagnant storage is selected. When the residual value R3 of the difference between the fracture concentration and the matrix concentration is the minimum, the optimal page sequence Pr is calculated and the base reversion page sequence is selected.

[0039] In this embodiment, the embryo value page sequence arrangement unit places the pipe-to-fracture syn-mold embryo B1, the fracture-to-matrix stagnant embryo B2, and the matrix-to-fracture retro-release embryo B3 into the same embryo value set. It then obtains the candidate exchange page sequence P using the sum of the three embryo values ​​S and the pipe-to-fracture occupancy values ​​N1, fracture-to-matrix occupancy values ​​N2, and matrix-to-fracture occupancy values ​​N3, allowing the three types of exchange actions to enter the comparison process as placeholders. Compared to directly using fixed exchange coefficients in the transport equations, this current processing can first distinguish between the three states: pipe-to-fracture syn-mold, fracture-to-matrix stagnant, and matrix-to-fracture retro-release, making the dominant exchange direction clearer during the current simulation period.

[0040] The candidate swap page order P is not directly used as the final swap result, but rather as the entry point for the quality trial calculation verification unit. This allows the pipe-fracture syn-migration embryo B1, the fracture-base stagnant embryo B2, and the matrix-returning retrograde embryo B3 to each participate in a short-step quality trial calculation. The pipe-to-fracture mass Q1, the fracture-to-matrix mass Q2, and the matrix-to-fracture mass Q3 represent the pipe-to-fracture mass, the fracture-to-matrix mass, and the matrix-to-fracture mass, respectively. A source domain mass restriction method is used to ensure that any calculated mass does not exceed the existing solute mass in the current source domain. This current processing prevents the mass transfer from exceeding the current domain's memory capacity when the syn-migration embryo value is large, ensuring that the three-domain mass transfer process remains consistent with the actual solute content in the pipe, fracture, and matrix.

[0041] In this embodiment, through the dual output of candidate exchange page order P and trial-calculated preferred page order Pr, the current module forms two judgment criteria for the subsequent on-the-fly exchange page number (CER) generation process: candidate exchange page order P reflects the occupancy status of the three types of mobile embryos in the current simulation period, and trial-calculated preferred page order Pr reflects the residual concentration status of the three types of mobile embryos after mass calculation. Compared with the existing methods that directly fix the exchange direction or select the exchange item based on only a single parameter, the current processing can explain the phenomena such as rapid arrival before the peak of the outlet concentration, slow decline after the peak, and tail-end rebound within the simulation period in segments, reducing the risk of matrix back-release being misjudged as new pollution input and fissure stagnation being misjudged as ordinary diffusion.

[0042] Example 4 Please refer to Figure 4 Specifically: the on-the-fly page number generation module includes a page sequence homomorphism verification unit and a short-step convergence unit; First, read the page order type of the candidate swap page order P. The page order types of the candidate swap page order P include tube fracture forward migration page order, fracture base stagnant page order, base reversion reverse release page order, and unclaimed wandering page order. Then read the page order type of the trial-optimized page order Pr. The page order types of Pr include tube fracture forward migration page order, fracture base stagnant storage page order, and base reversion reverse release page order. Perform homeomorphic verification between the candidate swap page order P and the trial-optimized page order Pr: If both the candidate swap page order P and the trial-optimized page order Pr are pipe fracture forward shift page orders, then a pipe fracture diversion page seed is generated. When both the candidate swap page order P and the trial-optimized page order Pr are split-base stagnant page orders, a split-base stagnant page seed is generated. When both the candidate swap page order P and the trial-optimized page order Pr are basis-reverse page orders, a basis-reverse page seed is generated. When the candidate swap page order P is an unclaimed wandering page order, and the trial-optimized page order Pr has a definite page order, then the trial-optimized page order Pr is used to form a page seed, and the page seed is marked as a trial takeover state. If the candidate swap page order P and the trial-preferred page order Pr point to different page orders, two candidate page folding seeds are formed: the tendency page folding seed formed by the candidate swap page order P and the trial-preferred page folding seed formed by the trial-preferred page order Pr.

[0043] The short-step convergence unit determines the number of seeds; If there is only one fold seed, the current fold seed is converted into a random exchange fold code CER. Specifically, the tube fracture shunt fold seed is converted into a tube fracture shunt fold code, the fracture basal storage fold seed is converted into a fracture basal storage fold code, and the basis reverse release fold seed is converted into a basis reverse release fold code. If there are two fold seeds, then perform short-step inheritance trial calculations on each of the two fold seeds; The short-step trial calculation reads the solute mass Mc in the pipeline, solute mass Mf in the fracture, and solute mass Mm in the matrix from the mass calculation results of the previous simulation period, as well as the water volume Wc in the pipeline, water volume Wf in the fracture, and exchangeable water volume Wm in the matrix in the three-domain discrete element volume. When the fold seed is a pipe fracture diversion fold seed, the pipe solute mass Mc and the fracture solute mass Mf are updated according to the transfer from the pipe to the fracture. When the folded seed is a cleaved folded seed, the mass of solute in the cleavage (Mf) and the mass of solute in the matrix (Mm) are updated according to the transfer from the cleavage to the matrix. When the folded seed is a base-returning folded seed, the matrix solute mass Mm and Mf and the crack solute mass Mf are updated according to the matrix returning to the crack; After completing the short-step acceptance trial calculation, the predicted short-step outlet concentration value of each fold seed is obtained. And calculate the export concentration compared to the previous simulation period. The formula for the difference in acceptance between them, Ak, is: ; In the formula, k represents any fold seed currently participating in the convergence; Read the page break code C0 from the previous simulation period and analyze its connection status with the current page break seed: If the current fold seed is the same as the fold code C0 of the previous simulation period, then the fold capacity Fk is 0; If the current fold seed continues along the sequence of tube fissure diversion, fissure base retention, and base reversion, then the fold receiving amount Fk is 1; When the current fold seed and the fold code C0 of the previous simulation period form a reverse jump, the fold carrying capacity Fk is 2; The folding convergence value Zk is synthesized from the export concentration difference Ak and the folding acceptance amount Fk, as shown in the following formula: ; For the fold convergence values ​​Zk of the two fold seeds, the fold seed with the smaller fold convergence value Zk is used to generate the shifted fold code CER for the current simulation period.

[0044] In this embodiment, the page sequence homeomorphic verification unit reads the candidate exchange page sequence P and the trial-preferred page sequence Pr separately and performs homeomorphic verification, so that the exchange judgment in the current simulation period is simultaneously constrained by the three-domain migrating embryo occupancy status and the mass trial calculation result. Compared with the method of directly determining the exchange direction based solely on the larger value among the tube fracture anterograde embryo B1, the fractured basal stagnant embryo B2, and the basal retrograde release embryo B3, the current processing can distinguish between two situations: "strong migration tendency" and "closer to concentration change after entering mass transfer", so that the determination of tube fracture splitting, fractured basal stagnant, and basal retrograde release no longer depends on a single value.

[0045] When two folding seeds exist, the short-step convergence unit performs short-step convergence trials separately, and obtains the predicted short-step outlet concentration value for each folding seed based on the mass calculation results of the previous simulation period. Then compare it with the export concentration of the previous simulation period. The difference in acceptance is calculated. The current processing ensures that the CER of the switching page is generated not only based on the wandering embryo and mass trial results of the current time period, but also constrained by the change in outlet concentration in the previous time period. This can alleviate the problem of abrupt changes in the outlet concentration curve caused by repeated switching of the switching code in adjacent simulation time periods.

[0046] The short-step convergence unit further reads the previous simulation period's on-the-fly exchange fold code C0, and uses the fold acceptance amount Fk to distinguish the continuation, sequential advancement, or reverse jump state between the current fold seed and C0. Then, it synthesizes the fold convergence value Zk by combining the outlet concentration acceptance difference Ak and the fold acceptance amount Fk, and selects the fold seed with the smaller Zk value to generate the on-the-fly exchange fold code CER. Thus, the on-the-fly exchange fold code CER can simultaneously express the current exchange direction, the concentration dissipation state after mass calculation, and the exchange acceptance state between previous and subsequent simulation periods. This is beneficial for distinguishing between rapid pre-peak pipe fissure splitting, post-peak fissure base retention, and tail-end base backflow, reducing the possibility of matrix backflow being misjudged as new pollution input and fissure retention being misjudged as ordinary diffusion.

[0047] Example 5 Please refer to Figure 1 Specifically: the fold code driven simulation module includes a fold code directional wandering unit and a three-domain quality write-back output unit; The page-code directional migration unit selects the unique migration embryo that enters the mass transfer calculation during the current simulation period from the tube fracture antegrade embryo B1, fracture-based stagnant embryo B2, and base-returning reverse release embryo B3, based on the process-exchange page code CER. Parsing the code pattern of the CER (Continuous Exchange Page Code): When the page number CER is the pipe fracture diversion page number, only the pipe fracture forward migration embryo B1 is allowed to participate in the calculation in the current time period. The pipe is used as the source domain and the fracture is used as the receiving domain to calculate and obtain the mass shift QER. When the time-travel exchange page number CER is the fractured stagnant page number, only the fractured stagnant embryo B2 is allowed to participate in the calculation in the current time period, and the fracture is used as the source domain and the matrix is ​​used as the receiving domain to calculate and obtain the mass migration QER. When the page number CER is the base-reverse release page number, only the base-reverse release embryo B3 is allowed to participate in the calculation in the current time period, and the matrix is ​​used as the source domain and the crack is used as the receiving domain to calculate and obtain the mass transfer quantity QER. The formula is as follows: ; ; ; The three-domain mass write-back outflow unit writes back the solute mass Mc in the pipeline, the solute mass Mf in the fissure, and the solute mass Mm in the matrix based on the time-of-flight exchange (CER) page number, and outputs the three-domain solute transport results for the current simulation period. Three mutually exclusive write-back actions are performed according to the CER (Cross-Process Interchange Page Number): When the time-of-flight page number CER is the pipe fracture diversion page number, the mass shift QER is subtracted from the pipe solute mass Mc, and the mass shift QER is incorporated into the fracture solute mass Mf. The matrix solute mass Mm remains the same as at the beginning of the current simulation period. When the time-travel page number CER is the fractured matrix stagnant page number, the mass shift QER is subtracted from the fracture solute mass Mf, and the mass shift QER is incorporated into the matrix solute mass Mm. The pipeline solute mass Mc remains the mass at the beginning of the current simulation period. When the time-travel page number CER is the base-reverse release page number, the mass shift QER is subtracted from the matrix solute mass Mm, and the mass shift QER is incorporated into the fracture solute mass Mf. The pipeline solute mass Mc remains the mass at the beginning of the current simulation period. When the page number CER during the process exchange is the pipe splitting shunt page number: ; ; ; When the page number CER during the process exchange is a cleavage-dependent page number: ; ; ; When the page number CER is a base-reverse page number: ; ; ; In the formula, Indicates the mass of solute in the new pipeline. Indicates the mass of solute in the new fracture. Indicates the mass of the solute in the new matrix; Based on the solute quality of the new pipeline quality of solute in new fractures and the quality of the new matrix solute The pipe concentration Cc, fracture concentration Cf, and matrix concentration Cm at the end of the current simulation period are processed to obtain new pipe concentration nCc, new fracture concentration nCf, and new matrix concentration nCm, as shown in the following formula: ; ; ; The new pipeline concentration nCc is used as the inlet concentration for the main pipeline transport calculation, and the outlet concentration Ct for the current simulation period is obtained based on the quality update result of the discrete element at the end of the pipeline. At the same time, the new fracture concentration nCf is sent to the fracture side calculation for the next simulation period, and the new matrix concentration nCm is sent to the matrix storage or matrix back-release calculation for the next simulation period.

[0048] In this embodiment, the fold code directional migration unit uses the on-the-fly fold code CER as the exchange control entry for the current simulation period, and selects one migration embryo from the following three types: pipe fracture syn-migration embryo B1, fracture-based storage embryo B2, and matrix retrograde release embryo B3, to participate in the mass transfer calculation. This current processing ensures that only one type of exchange action—pipe fracture diversion, fracture-based storage, or matrix retrograde release—is executed within the same simulation period, preventing the three types of exchange actions from canceling each other out after entering the calculation simultaneously, thus making the mass transfer direction between the pipe, fracture, and matrix clearer.

[0049] The mass transfer quantity QER is jointly obtained based on the shifted page number CER, the selected free embryo, and the solute mass of the source domain. A min function is used to limit the mass transfer quantity QER of the shifted page number to not exceed the existing solute mass of the source domain. Therefore, when the pipe is the source domain, the mass transfer quantity does not exceed the pipe solute mass Mc; when the fracture is the source domain, the mass transfer quantity does not exceed the fracture solute mass Mf; and when the matrix is ​​the source domain, the mass transfer quantity does not exceed the matrix solute mass Mm, ensuring that the mass transfer process in the three domains conforms to the solute inventory state during the current simulation period.

[0050] The three-domain mass write-back outflow unit performs mutually exclusive write-back operations based on the CER (Continuous Exchange Page Code), generating new pipe solute mass nMc, new fracture solute mass nMf, and new matrix solute mass nMm, respectively. The current processing ensures that pipe deduction and fracture incorporation during pipe fracture splitting, fracture deduction and matrix incorporation during fracture matrix stagnation, and matrix deduction and fracture incorporation during matrix reversion all have clear computational paths. This allows the previously generated page code results to be transformed into three-domain mass results that can continue to participate in transport calculations within the current simulation period.

[0051] In this embodiment, the system re-obtains the new pipeline solute mass nMc, the new fracture solute mass nMf, and the new matrix solute mass nMm based on the new pipeline solute mass nCc, the new fracture solute mass nMf, and the new matrix solute mass nMm. nCc is then used to calculate the outlet concentration Ct for the current simulation period. The new fracture concentration nCf and the new matrix concentration nCm are respectively fed into the fracture-side calculation and matrix storage or matrix back-release calculation for the next simulation period. Therefore, the outlet concentration Ct is no longer solely calculated from a fixed exchange coefficient, but is driven by the current period's pipe fracture diversion, fracture-based matrix storage, or matrix back-release actions. This helps explain phenomena such as rapid arrival before the outlet concentration peak, slow decline after the peak, and a tail-end rise, and reduces the possibility of matrix back-release being misjudged as new pollution input and fracture storage being misjudged as ordinary diffusion.

[0052] Example 6 A numerical simulation method for solute transport processes in karst groundwater; please refer to [reference needed]. Figure 2 Specifically, it includes the following steps: Step 1: The basic data input module obtains the basic data set AW based on the pipeline calculation node head sequence, fracture calculation node head sequence, karst channel discrete model, water balance results, three-domain discrete element volume, and mass calculation results of the previous simulation period. Step 2: The three-domain migratory embryo generation module analyzes the basic data set AW to generate tube-crack syn-migratory embryo B1, cleavage-based stagnant embryo B2, and basal-retrograde release embryo B3. Step 3: The free embryo trial calculation and evaluation module sorts and performs quality calculations on the tube-cracked antegrade embryo B1, the fractured base-retained embryo B2, and the base-returning reverse-release embryo B3 to obtain the candidate exchange page order P and the trial calculation preferred page order Pr; Step 4: The on-the-fly page number generation module generates the on-the-fly page number CER for the current simulation period by combining the candidate page number P and the trial-and-error preferred page number Pr with the quality calculation results of the previous simulation period. Step 5: The fold code-driven simulation module controls the mass transfer between pipes, fissures, and matrix based on the live fold code CER, and outputs the three-domain solute transport results for the current simulation period.

[0053] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and variations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended technical solutions and their equivalents.

Claims

1. A numerical simulation system for solute transport processes in karst groundwater, characterized in that: It includes a basic data input module, a three-domain wandering embryo generation module, a wandering embryo trial calculation and evaluation module, a process-interchange folding code generation module, and a folding code-driven simulation module; The basic data input module obtains the basic data set AW based on the head sequence of pipeline calculation nodes, the head sequence of fracture calculation nodes, the discrete model of karst channels, the water balance results, the volume of the three-domain discrete unit and the mass calculation results of the previous simulation period. The three-domain migratory embryo generation module analyzes the basic data set AW and generates tube-crack antegrade embryo B1, cleavage-based stagnant embryo B2, and basal-returning retrograde embryo B3. The free embryo trial calculation and evaluation module sorts and performs quality calculations on the tube-cracked antegrade embryo B1, the fractured base-retained embryo B2, and the base-returning reverse-release embryo B3 to obtain the candidate exchange page order P and the trial calculation preferred page order Pr. The on-the-fly page number generation module generates the on-the-fly page number CER for the current simulation period by combining the candidate page number P and the trial-and-error preferred page number Pr with the quality calculation results of the previous simulation period. The fold code-driven simulation module controls the mass transfer between pipes, fissures, and matrix based on the on-the-fly fold code CER, and outputs the three-domain solute transport results for the current simulation period.

2. The numerical simulation system for solute transport processes in karst groundwater according to claim 1, characterized in that: The basic data input module includes a three-domain source quantity extraction unit and a data normalization and shaping unit; The three-domain source quantity extraction unit extracts three-domain source quantities within the same simulation period from the existing data of the current numerical simulation system, and organizes data from different sources into the original basic quantities under the same period. Using the current simulation time period tn as the index, read the upstream node head Hc1 and downstream node head Hc2 in the pipeline calculation node head sequence, and read the fracture inlet node head Hf1 and fracture outlet node head Hf2 in the fracture calculation node head sequence. Simultaneously, the pipe section length Lc, fissure section length Lf, and pipe water-passing area Ac are read from the discrete model of the karst channel; the pipe flow rate Qc, fissure outflow rate Qfout, and fissure residence time Tf are read from the water balance results of the previous simulation period; the pipe water volume Wc, fissure water volume Wf, and matrix exchangeable water volume Wm are read from the three source quantities; and the pipe solute mass Mc, fissure solute mass Mf, and matrix solute mass Mm are read from the mass calculation results of the previous simulation period. Align the acquired data: if the data comes from the previous simulation period, then use it as the start state of the simulation period tn; When the data comes from the current discrete model, the data is used as a fixed structural quantity for the simulation time period tn. Based on the aligned data, the following parameters are calculated: pipe hydraulic gradient Jc, fracture hydraulic gradient Jf, pipe velocity Vc, pipe concentration Cc, fracture concentration Cf, matrix concentration Cm, pipe-to-fracture concentration difference Gcf, and fracture-to-matrix concentration difference Gfm; the formulas are as follows: ; ; ; ; ； ； ； ； The data normalization shaping unit processes the pipe hydraulic gradient Jc, fracture hydraulic gradient Jf, pipe velocity Vc, pipe concentration Cc, fracture concentration Cf, fracture residence time Tf, matrix concentration Cm, pipe-to-fracture concentration difference Gcf, and fracture-to-matrix concentration difference Gfm using the maximum-minimum normalization method to obtain normalized data, and then fits the data to obtain the basic data set AW.

3. The numerical simulation system for solute transport processes in karst groundwater according to claim 2, characterized in that: The three-domain free embryo generation module includes a three-domain directional free embryo fragmentation unit and a free embryo compression forming unit; The three-domain directional wandering and fragmentation unit performs directional fragmentation processing on the data in the basic data group AW. Read the state of the concentration difference Gcf between the pipeline and the fracture. When the concentration difference Gcf between the pipeline and the fracture is positive, it means that the concentration in the pipeline is higher than the concentration in the fracture, and there is a concentration-driven phenomenon of solute being released from the pipeline to the fracture side during the current period. When the concentration difference Gcf between the pipeline and the fracture is not positive, the concentration push from the pipeline to the fracture side does not enter the fracture direction plate; Read the sign of the concentration difference Gfm between the fracture and the matrix. When the concentration difference Gfm between the fracture and the matrix is ​​positive, it means that the fracture concentration is higher than the matrix concentration, and there is a concentration drive of solute from the fracture to the matrix side in the current period. When the concentration difference between the fracture and the matrix, Gfm, is negative, it indicates that the matrix concentration is higher than the fracture concentration, and there is a concentration-driven phenomenon of solute returning from the matrix to the fracture side during the current period. The pipe hydraulic gradient Jc and pipe velocity Vc are used as the driving forces for pipe transport, the fracture hydraulic gradient Jf is used as the driving force for the entry of the fracture side into the matrix, and the fracture residence time Tf is used as the time support for the temporary storage or re-release process of the residence action in the fracture. Three directional slices are obtained, including pipe fracture directional slice D1, fracture-to-matrix directional slice D2, and matrix-to-return directional slice D3, as shown in the following formula: ； ； ； In the formula, This indicates taking the maximum value.

4. The numerical simulation system for solute transport processes in karst groundwater according to claim 3, characterized in that: The free embryo compression forming unit is based on the basic data set AW, and is analyzed in combination with the tube fracture direction sheet D1, the fracture base direction sheet D2, and the base return direction sheet D3. First, analyze whether the pipe fracture direction segment D1 has fast flow support in the pipe; When the pipe fracture direction segment D1 has a value, the pipe flow velocity Vc is again used as the forward traction term. Analyze whether the fracture-oriented sheet D2 has fracture retention support; When the fracture orientation segment D2 has a value, the fracture residence time Tf is used as the storage compression term; Analyze whether the base return direction segment D3 has a dwell return support; When the base return direction piece D3 has a value, the crack residence time Tf is used as the return support item; Finally, the formulas for the output tube-cracked antegrade embryo B1, the cleaved basal residing embryo B2, and the basal reversal-released embryo B3 are as follows: ； ； 。 5. The numerical simulation system for solute transport processes in karst groundwater according to claim 4, characterized in that: The free embryo trial calculation and evaluation module includes an embryo value page sequence arrangement unit and a quality trial calculation verification unit; The embryo value page ordering unit obtains candidate exchange page order P based on the values ​​of tube-fission antegrade embryo B1, fission-based stagnant embryo B2, and basal-retrograde release embryo B3. The process is as follows: Embryo B1 (transient from tube), embryo B2 (retained at the base), and embryo B3 (reverse release from the base) are placed into the same embryo value set, and the total amount S of the three embryos is calculated using the following formula: ; When the total amount of the three embryos S is not zero, calculate the occupancy values ​​of the anterograde embryo B1, the cleavage-based residing embryo B2, and the basal-returning embryo B3 in the total amount of the three embryos S, respectively, to obtain the occupancy values ​​N1, N2, and N3 of the anterograde embryo, as follows: ； ； ； The pages are arranged in order of magnitude: pipe fracture occupancy value N1, fracture matrix occupancy value N2, and matrix reversion occupancy value N3. When the pipe crack occupancy value N1 is in the first position, the candidate swap page order P is marked as the pipe crack forward page order; When the fractured occupancy value N2 is first, the candidate exchange page order P is marked as the fractured stagnant page order. When the base-return placeholder value N3 is in the first position, the candidate swap page order P is marked as the base-return reverse page order.

6. The numerical simulation system for solute transport processes in karst groundwater according to claim 5, characterized in that: The quality trial calculation verification unit performs a short-step quality trial calculation on each of the following embryos: tube-cracked antegrade embryo B1, cleaved base-retained embryo B2, and base-returning retrograde embryo B3, based on the candidate exchange page order P, and obtains the optimal page order Pr. The process is as follows: Read the following data: pipe fracture syn-migration embryo B1, fracture-base stagnant embryo B2, and matrix reversion embryo B3; and the solute mass Mc of the pipe, solute mass Mf of the fracture, and solute mass Mm of the matrix from the mass calculation results of the previous simulation period. Simultaneously, read the water volume Wc of the pipe, the water volume Wf of the fracture, and the exchangeable water volume Wm of the matrix. Calculate the following mass from the pipe to the fracture: Q1; from the fracture to the matrix: Q2; and from the matrix to the fracture: Q3. The formulas are as follows: ； ； ； In the formula, This indicates taking the minimum value. Indicates the pipe fracture lateral migration conversion factor. Indicates the conversion coefficient of fractured substrate storage. Indicates the base-retroactive conversion coefficient; The obtained mass values ​​Q1 from the pipe to the fracture, Q2 from the fracture to the matrix, and Q3 from the matrix to the fracture are analyzed. The residual values ​​R1, R2, and R3 of the concentration difference between the pipe and fracture are calculated, using the following formulas: Compare the values ​​of the residual values ​​R1 (difference between pipe concentration and fracture concentration), R2 (difference between fracture concentration and matrix concentration), and R3 (difference between fracture concentration and matrix concentration). When the residual value R1 of the difference between the concentration in the pipeline and the concentration in the crack is the smallest, the optimal page sequence Pr is determined by trial calculation and the pipeline crack shift page sequence is taken. When the residual value R2 of the difference between fracture concentration and matrix concentration is the minimum, the optimal page sequence Pr is calculated and the page sequence with fractured matrix stagnant storage is selected. When the residual value R3 of the difference between the fracture concentration and the matrix concentration is the minimum, the optimal page sequence Pr is calculated and the base reversion page sequence is selected.

7. The numerical simulation system for solute transport processes in karst groundwater according to claim 6, characterized in that: The on-the-fly folding code generation module includes a page sequence homomorphism verification unit and a short-step convergence unit; First, read the page order type of the candidate swap page order P. The page order types of the candidate swap page order P include tube fracture forward migration page order, fracture base stagnant page order, base reversion reverse release page order, and unclaimed wandering page order. Then read the page order type of the trial-optimized page order Pr. The page order types of Pr include tube fracture forward migration page order, fracture base stagnant storage page order, and base reversion reverse release page order. Perform homeomorphic verification between the candidate swap page order P and the trial-optimized page order Pr: If both the candidate swap page order P and the trial-optimized page order Pr are pipe fracture forward shift page orders, then a pipe fracture diversion page seed is generated. When both the candidate swap page order P and the trial-optimized page order Pr are split-base stagnant page orders, a split-base stagnant page seed is generated. When both the candidate swap page order P and the trial-optimized page order Pr are basis-reverse page orders, a basis-reverse page seed is generated. When the candidate swap page order P is an unclaimed wandering page order, and the trial-optimized page order Pr has a definite page order, then the trial-optimized page order Pr is used to form a page seed, and the page seed is marked as a trial takeover state. If the candidate swap page order P and the trial-preferred page order Pr point to different page orders, two candidate page folding seeds are formed: the tendency page folding seed formed by the candidate swap page order P and the trial-preferred page folding seed formed by the trial-preferred page order Pr.

8. The numerical simulation system for solute transport processes in karst groundwater according to claim 7, characterized in that: The short-step convergence unit determines the number of seeds: If there is only one fold seed, the current fold seed is converted into a random exchange fold code CER. Specifically, the tube fracture shunt fold seed is converted into a tube fracture shunt fold code, the fracture basal storage fold seed is converted into a fracture basal storage fold code, and the basis reverse release fold seed is converted into a basis reverse release fold code. If there are two fold seeds, then perform short-step inheritance trial calculations on each of the two fold seeds; The short-step trial calculation reads the solute mass Mc in the pipeline, solute mass Mf in the fracture, and solute mass Mm in the matrix from the mass calculation results of the previous simulation period, as well as the water volume Wc in the pipeline, water volume Wf in the fracture, and exchangeable water volume Wm in the matrix in the three-domain discrete element volume. When the fold seed is a pipe fracture diversion fold seed, the pipe solute mass Mc and the fracture solute mass Mf are updated according to the transfer from the pipe to the fracture. When the folded seed is a cleaved folded seed, the mass of solute in the cleavage (Mf) and the mass of solute in the matrix (Mm) are updated according to the transfer from the cleavage to the matrix. When the folded seed is a base-returning folded seed, the matrix solute mass Mm and Mf and the crack solute mass Mf are updated according to the matrix returning to the crack; After completing the short-step acceptance trial calculation, the predicted short-step outlet concentration value of each fold seed is obtained. And calculate the export concentration compared to the previous simulation period. The formula for the difference in acceptance between them, Ak, is: ; In the formula, k represents any fold seed currently participating in the convergence; Read the page break code C0 from the previous simulation period and analyze its connection status with the current page break seed: If the current fold seed is the same as the fold code C0 of the previous simulation period, then the fold capacity Fk is 0; If the current fold seed continues along the sequence of tube fissure diversion, fissure base retention, and base reversion, then the fold receiving amount Fk is 1; When the current fold seed and the fold code C0 of the previous simulation period form a reverse jump, the fold carrying capacity Fk is 2; The folding convergence value Zk is synthesized from the export concentration difference Ak and the folding acceptance amount Fk, as shown in the following formula: ； For the fold convergence values ​​Zk of the two fold seeds, the fold seed with the smaller fold convergence value Zk is used to generate the shifted fold code CER for the current simulation period.

9. A numerical simulation system for solute transport processes in karst groundwater according to claim 8, characterized in that: The fold code driven simulation module includes a fold code directional wandering unit and a three-domain quality write-back output unit; The page-code directional migration unit selects the unique migration embryo that enters the mass transfer calculation during the current simulation period from the tube fracture antegrade embryo B1, fracture-based stagnant embryo B2, and base-returning reverse release embryo B3, based on the process-exchange page code CER. Parsing the code pattern of the CER (Continuous Exchange Page Code): When the page number CER is the pipe fracture diversion page number, only the pipe fracture forward migration embryo B1 is allowed to participate in the calculation in the current time period. The pipe is used as the source domain and the fracture is used as the receiving domain to calculate and obtain the mass shift QER. When the time-travel exchange page number CER is the fractured stagnant page number, only the fractured stagnant embryo B2 is allowed to participate in the calculation in the current time period, and the fracture is used as the source domain and the matrix is ​​used as the receiving domain to calculate and obtain the mass migration QER. When the page number CER is the base-reverse release page number, only the base-reverse release embryo B3 is allowed to participate in the calculation in the current time period, and the matrix is ​​used as the source domain and the crack is used as the receiving domain to calculate and obtain the mass transfer quantity QER. The formula is as follows: ； ； ； The three-domain mass write-back outflow unit writes back the solute mass Mc in the pipeline, the solute mass Mf in the fissure, and the solute mass Mm in the matrix based on the time-of-flight exchange (CER) page number, and outputs the three-domain solute transport results for the current simulation period. Three mutually exclusive write-back actions are performed according to the CER (Cross-Process Interchange Page Number): When the time-of-flight page number CER is the pipe fracture diversion page number, the mass shift QER is subtracted from the pipe solute mass Mc, and the mass shift QER is incorporated into the fracture solute mass Mf. The matrix solute mass Mm remains the same as at the beginning of the current simulation period. When the time-travel page number CER is the fractured matrix stagnant page number, the mass shift QER is subtracted from the fracture solute mass Mf, and the mass shift QER is incorporated into the matrix solute mass Mm. The pipeline solute mass Mc remains the mass at the beginning of the current simulation period. When the time-travel page number CER is the base-reverse release page number, the mass shift QER is subtracted from the matrix solute mass Mm, and the mass shift QER is incorporated into the fracture solute mass Mf. The pipeline solute mass Mc remains the mass at the beginning of the current simulation period. When the page number CER during the process exchange is the pipe splitting shunt page number: ； ； ； When the page number CER during the process exchange is a cleavage-dependent page number: ； ； ； When the page number CER is a base-reverse page number: ； ； ； In the formula, Indicates the mass of solute in the new pipeline. Indicates the mass of solute in the new fracture. Indicates the mass of the solute in the new matrix; Based on the solute quality of the new pipeline quality of solute in new fractures and the quality of the new matrix solute The pipe concentration Cc, fracture concentration Cf, and matrix concentration Cm at the end of the current simulation period are processed to obtain new pipe concentration nCc, new fracture concentration nCf, and new matrix concentration nCm, as shown in the following formula: ； ； ； The new pipeline concentration nCc is used as the inlet concentration for the main pipeline transport calculation, and the outlet concentration Ct for the current simulation period is obtained based on the quality update result of the discrete element at the end of the pipeline. At the same time, the new fracture concentration nCf is sent to the fracture side calculation for the next simulation period, and the new matrix concentration nCm is sent to the matrix storage or matrix back-release calculation for the next simulation period.

10. A numerical simulation method for solute transport processes in karst groundwater, applied to the numerical simulation system for solute transport processes in karst groundwater as described in any one of claims 1 to 9, characterized in that: Includes the following steps: Step 1: The basic data input module obtains the basic data set AW based on the pipeline calculation node head sequence, fracture calculation node head sequence, karst channel discrete model, water balance results, three-domain discrete element volume, and mass calculation results of the previous simulation period. Step 2: The three-domain migratory embryo generation module analyzes the basic data set AW to generate tube-crack syn-migratory embryo B1, cleavage-based stagnant embryo B2, and basal-retrograde release embryo B3. Step 3: The free embryo trial calculation and evaluation module sorts and performs quality calculations on the tube-cracked antegrade embryo B1, the fractured base-retained embryo B2, and the base-returning reverse-release embryo B3 to obtain the candidate exchange page order P and the trial calculation preferred page order Pr; Step 4: The on-the-fly page number generation module generates the on-the-fly page number CER for the current simulation period by combining the candidate page number P and the trial-and-error preferred page number Pr with the quality calculation results of the previous simulation period. Step 5: The fold code-driven simulation module controls the mass transfer between pipes, fissures, and matrix based on the live fold code CER, and outputs the three-domain solute transport results for the current simulation period.

Images