A numerical simulation method for seepage and solute transport in three-dimensional fracture network based on topological pruning

Abstract

The application discloses a three-dimensional fissure network seepage and solute migration numerical simulation method based on topological pruning, which applies graph theory to analyze the topological structure of the fissure network from two scales of global and local, establishes the quantitative relationship between the fissure connectivity and the underground flow field, and accurately depicts the influence mechanism of the topological characteristics on the seepage field; based on the network topological analysis, the key fissure sub-network is extracted through the pruning method, and the high-precision numerical simulation is carried out by replacing the complete network with the key fissure sub-network, so that the calculation efficiency is significantly improved while the calculation precision is ensured. The application provides an efficient and reliable solution for the seepage and solute migration problems in the hard rock fissure medium.

Description

Technical Field

[0001] This invention relates to a numerical simulation method for flow fields in fracture networks in underground strata, specifically a numerical simulation method for seepage and solute transport in three-dimensional fracture networks based on topological pruning. Background Technology

[0002] In hard rock, fractures are the primary channels for groundwater flow and chemical migration, making them crucial for various large-scale underground engineering projects. Discrete fracture models (DFNs), a commonly used method for assessing the feasibility of large-scale projects located in fractured strata, can accurately characterize the uncertainties in seepage and mass transfer behavior caused by the heterogeneity of fracture distribution, but they are computationally expensive. On the other hand, due to the low permeability of hard rock and the heterogeneous connectivity of fracture networks, fluid flow within fractured rock strata often exhibits a distinct dominant path. Summary of the Invention

[0003] Purpose of the invention: The purpose of this invention is to provide a low-computational-cost numerical simulation method for seepage and solute transport in three-dimensional fractured networks based on topology pruning.

[0004] Technical solution: The present invention provides a numerical simulation method for seepage and solute transport in a three-dimensional fractured network based on topology pruning, comprising:

[0005] (1) Collect basic data on fracture networks and basic data on hydrogeology, and construct the target three-dimensional fracture network based on the basic data on fracture networks;

[0006] (2) Use graph theory to map the target three-dimensional fracture network into an equivalent two-dimensional undirected graph, and perform global and local network analysis;

[0007] (3) Generalize the hydrogeological conditions of the target area, construct and solve the groundwater flow model and solute transport model, and simulate the underground flow field and non-reactive transport and adsorption process of solute under different topological networks;

[0008] (4) Based on the analysis results obtained in step (2) and the simulation results obtained in step (3), determine the influence mechanism of network topology on seepage field and chemical field, and prune the target fracture network according to the key fractures to extract the key sub-networks that carry most of the material transport.

[0009] (5) Use key sub-networks to replace the complete network, solve the non-reactive solute migration model, obtain the groundwater flow field, solute breakthrough curve, solute migration field and solid solute adsorption in fractured strata, and realize efficient numerical simulation of seepage and solute transport in three-dimensional fractured networks.

[0010] Furthermore, in step (1), the basic data of the fracture network includes the fracture orientation distribution characteristics, fracture diameter distribution characteristics, and fracture density distribution characteristics.

[0011] Furthermore, the method for generating the fractured network is as follows:

[0012] The crack diameter l follows a power-law distribution, expressed as:

[0013] f(l)=βl -α

[0014] Where f(l) is the probability density function of the crack diameter distribution; β is the probability density function based on the minimum value of the crack diameter l. min and maximum value l max The obtained normalization exponent; α is the power law exponent;

[0015] The crack diameter l is represented as:

[0016]

[0017] Where R is a random number;

[0018] The fracture orientation follows the Fisher distribution:

[0019]

[0020] Where θ is the dip angle; θ0 is the average dip angle; φ is the dip direction; φ0 is the average dip direction; κ is the dispersion parameter; and sinh is the hyperbolic sine function.

[0021] Each fracture network is configured to have the same aperture.

[0022] Furthermore, in step (1), the basic hydrogeological data includes structural geological maps and hydrogeological maps.

[0023] Further, in step (2), the target three-dimensional fracture network is mapped to an equivalent two-dimensional undirected graph using graph theory, including:

[0024] Treating the cracks as nodes V and the intersections between cracks as edges E, the 3D crack network is transformed into a 2D undirected graph G(V,E) as follows:

[0025] The three-dimensional fracture network F consists of N fractures f, represented as:

[0026] F={f i}, i = 1, ..., N

[0027] Define a bijective mapping Φ:

[0028] Φ:f i →v i

[0029] If there are two cracks f i with f jintersect, Then there is an edge in E connecting the corresponding vertex:

[0030]

[0031] Among them, (v i ,v j )∈E represents vertex v i and v j The edge e between ij , side e ij weight w ij Set as vertex v i and v j Length of the line of intersection between them:

[0032] w ij =Lengthf i ∩f j

[0033] Similarly, considering the flow direction, we can regard the inlet plane x0 as the source node s, and the outlet plane x L Considering the target node t, when the crack intersects with the inlet and outlet boundaries:

[0034]

[0035] Among them, e si For an edge between a source node and a non-source / sink node, e it The edges between non-source / sink nodes and the target node are thus transformed into a two-dimensional undirected graph.

[0036] Furthermore, in step (2), the global network analysis is based on three types of network topologies: the full network TN, the connected subnetwork CSN, and the dead-end-free connected subnetwork DCSN. The target fracture network is called the full network TN. When the three-dimensional fracture network is sparse, the fracture network will be divided into several parts. The fracture network connecting the source surface and the sink surface is called the connected subnetwork CSN. The CSN that does not contain the nodes corresponding to the dead-end fractures is called the dead-end-free connected subnetwork DCSN.

[0037] Furthermore, local network analysis, based on betweenness centrality BC, is defined as:

[0038]

[0039] Where BC(i) represents the betweenness centrality of node i; g jk (i) represents the number of times node i falls on the path between nodes j and k; g jk BC represents the number of times the shortest path between nodes j and k is calculated; the higher the BC value of a node, the more critical it is in the network.

[0040] Furthermore, in step (3), the governing equations of the fractured network groundwater flow model are established based on the cubic law, and the solute transport field includes advection, diffusion, and isothermal adsorption, wherein the solute adsorption equation is as follows:

[0041]

[0042] Among them, c s It is the solute adsorption capacity of the solid [mol / kg]; K d The equilibrium adsorption partition coefficient of the solute [m] 3 / mol];c s_max c is the maximum amount of solute adsorbed on the solid [mol / kg]; c is the solute concentration [mol / m]. 3 ].

[0043] Furthermore, in step (4), the pruning process of the fracture network is based on graph theory and network analysis, where the weight of the edge of the graph is assigned as the reciprocal of the length of the corresponding fracture intersection line; local network analysis is performed on the DCSN of the target fracture network, and BC is used to quantify the control capability of a single fracture on the flow field within the network; nodes with high BC values ​​are selected in the target fracture network, and Dijkstra's algorithm is used to solve the shortest seepage path network from the source node to the target node through nodes with high BC values.

[0044] Furthermore, in step (5), the non-reactive solute migration model is solved by using the finite element method in the commercial physics simulation software COMSOL Multiphysics for spatial and temporal discretization, and the Newton iteration method for solving the discretized nonlinear equations. The breakthrough curve of the solute over time is calculated at the outflow boundary surface of the target region to achieve efficient numerical simulation of seepage and solute transport in a three-dimensional fracture network.

[0045] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages:

[0046] This invention uses graph theory and network analysis to determine the influence mechanism of topological features on the seepage field and solute transport field of fracture networks. Based on betweenness centrality and Dijkstra's algorithm, fractures are pruned to obtain a high-throughput fracture subnetwork. This reduces both the complexity of the discrete fracture model network and the computational cost. When used in high-precision numerical simulations, the high-throughput fracture subnetwork significantly improves computational efficiency while maintaining computational accuracy.

[0047] This invention provides an efficient and reliable solution to the problems of seepage and solute transport in fractured hard rock media. Attached Figure Description

[0048] Figure 1This is a flowchart of the numerical simulation method for seepage and solute transport in a three-dimensional fracture network based on topology pruning provided in this embodiment of the invention.

[0049] Figure 2 (a) is the target fracture network TN in this embodiment of the invention. Figure 2 (b) is the graphical representation of TN in the embodiment of the present invention, where the coloring represents the average flow velocity on each fracture;

[0050] Figure 3 (a) is the CSN of the target fracture network in an embodiment of the present invention. Figure 3 (b) is the graphical representation of CSN in the embodiment of the present invention, where the coloring represents the average flow velocity on each fracture.

[0051] Figure 4 (a) is the DCSN of the target fracture network in an embodiment of the present invention. Figure 4 (b) is the graphical representation of DCSN in the embodiment of the present invention, where the coloring represents the average flow velocity on each fracture.

[0052] Figure 5 It is the BTC of Se(IV) in the solute transport model of CSN and DCSN in the embodiments of the present invention;

[0053] Figure 6 This is the adsorption of Se(IV) in CSN in the embodiments of the present invention (t=20a) and its graphical characterization;

[0054] Figure 7 This is a graphical representation of the DCSN in this embodiment of the invention, where the coloring represents the BC value of each crack;

[0055] Figure 8 It is the BTC of Se(IV) in the solute transport model of the subnetworks obtained by CSN, DCSN and different pruning schemes in the embodiments of the present invention. Detailed Implementation

[0056] The invention will now be further described with reference to the accompanying drawings.

[0057] like Figure 1 As shown, this embodiment of the invention provides a numerical simulation method for seepage and solute transport in a three-dimensional fracture network based on topology pruning, comprising the following steps:

[0058] (1) Collect basic data on fracture networks and basic data on hydrogeology, and construct the target three-dimensional fracture network based on the basic data on fracture networks;

[0059] Basic hydrogeological data includes structural geological maps and hydrogeological maps.

[0060] The basic data for fracture networks include fracture orientation distribution characteristics, fracture diameter distribution characteristics, and fracture density distribution characteristics. The method for generating fracture networks is as follows:

[0061] The crack diameter l follows a power-law distribution, expressed as:

[0062] f(l)=βl -α

[0063] Where f(l) is the probability density function of the crack diameter distribution; β is the probability density function based on the minimum value of the crack diameter l. min and maximum value l max The obtained normalization exponent; α is the power law exponent;

[0064] The crack diameter l is represented as:

[0065]

[0066] Where R is a random number;

[0067] The fracture orientation follows the Fisher distribution:

[0068]

[0069] Where θ is the dip angle; θ0 is the average dip angle; φ is the dip direction; φ0 is the average dip direction; κ is the dispersion parameter; and sinh is the hyperbolic sine function.

[0070] Each fracture network is configured to have the same aperture.

[0071] (2) Use graph theory to map the target three-dimensional fracture network into an equivalent two-dimensional undirected graph, and perform global and local network analysis;

[0072] Using graph theory, the target 3D fracture network is mapped to an equivalent 2D undirected graph, including:

[0073] Treating the cracks as nodes V and the intersections between cracks as edges E, the 3D crack network is transformed into a 2D undirected graph G(V,E) as follows:

[0074] The three-dimensional fracture network F consists of N fractures f, represented as:

[0075] F={f i}, i = 1, ..., N

[0076] Define a bijective mapping Φ:

[0077] Φ:f i →v i

[0078] If there are two cracks f i with fj intersect, Then there is an edge in E connecting the corresponding vertex:

[0079]

[0080] Among them, (v i ,v j )∈E represents vertex v i and v j The edge e between ij , side e ij weight w ij Set as vertex v i and v j Length of the line of intersection between them:

[0081] w ij =Lengthf i ∩f j

[0082] Similarly, considering the flow direction, we can regard the inlet plane x0 as the source node s, and the outlet plane x L Considering the target node t, when the crack intersects with the inlet and outlet boundaries:

[0083]

[0084] Among them, e si For an edge between a source node and a non-source / sink node, e it The edges between non-source / sink nodes and the target node are thus transformed into a two-dimensional undirected graph.

[0085] Global network analysis is based on three types of network topologies: the complete network (TN), the connected subnetwork (CSN), and the dead-end-free connected subnetwork (DCSN). The target fracture network is called the complete network (TN). When the 3D fracture network is sparse, it is divided into several parts. The fracture network connecting the source and sink surfaces is called the connected subnetwork (CSN); the CSN that does not contain nodes corresponding to dead-end fractures is called the dead-end-free connected subnetwork (DCSN).

[0086] Local network analysis, based on betweenness centrality (BC), is defined as follows:

[0087]

[0088] Where BC(i) represents the betweenness centrality of node i; g jk (i) represents the number of times node i falls on the path between nodes j and k; g jkBC is the number of times the shortest path between nodes j and k is found; BC measures the influence of a node in the network, and the higher the BC value of a node, the more critical it is in the network.

[0089] (3) Generalize the hydrogeological conditions of the target area, construct and solve the groundwater flow model and solute transport model, and simulate the underground flow field and non-reactive transport and adsorption process of solute under different topological networks;

[0090] The governing equations of the fractured network groundwater flow model are established based on the cubic law. The solute transport field includes advection, diffusion, and isothermal adsorption, with the solute adsorption equation as follows:

[0091]

[0092] Among them, c s It is the solute adsorption capacity of the solid [mol / kg]; K d The equilibrium adsorption partition coefficient of the solute [m] 3 / mol];c s_max c is the maximum amount of solute adsorbed on the solid [mol / kg]; c is the solute concentration [mol / m]. 3 ].

[0093] (4) Based on the analysis results obtained in step (2) and the simulation results obtained in step (3), determine the influence mechanism of network topology on seepage field and chemical field, and prune the target fracture network according to the key fractures to extract the key sub-networks that carry most of the material transport.

[0094] The pruning process of the fracture network is based on graph theory and network analysis, where the weights of the graph edges are assigned the reciprocal of the length of the corresponding fracture intersection line. Local network analysis is performed on the DCSN of the target fracture network, and BC quantification is used to assess the control capability of a single fracture on the flow field. Nodes with high BC values ​​within the target fracture network are selected, and Dijkstra's algorithm is used to find the shortest seepage path network from the source node to the target node through nodes with high BC values. The pruned network exhibits the largest pressure gradient, indicating the likelihood of dominant flow.

[0095] (5) Use key sub-networks to replace the complete network, solve the non-reactive solute migration model, obtain the groundwater flow field, solute breakthrough curve, solute migration field and solid solute adsorption in fractured strata, and realize efficient numerical simulation of seepage and solute transport in three-dimensional fractured networks.

[0096] To solve the non-reactive solute migration model, the finite element method in the commercial physics simulation software COMSOL Multiphysics was used for spatial and temporal discretization, and the Newton-Raphson iteration method was used to solve the discretized nonlinear equations. The solute breakthrough curve over time was calculated at the outflow boundary surface of the target region to achieve efficient numerical simulation of seepage and solute transport in a three-dimensional fractured network.

[0097] This invention provides a numerical simulation method for seepage and solute transport in a three-dimensional fracture network based on topology pruning. It uses network analysis and Dijkstra's algorithm to extract a high-pass quantum network of three-dimensional fractures, which greatly reduces the complexity of discrete fracture models and alleviates the computational burden on underground flow and chemical fields. This provides an efficient and reliable solution for seepage and solute transport problems in hard rock fractured media.

[0098] Here is a specific example.

[0099] This example demonstrates a numerical simulation method for seepage and solute transport in a three-dimensional fractured network based on topology pruning. It comprehensively considers the topological properties of the fractured network and mainly includes two parts: 1) global and local topological analysis of the fractured network based on graph theory; 2) numerical model of seepage and solute transport in a three-dimensional fractured network based on topology pruning.

[0100] 1) Global and local topology analysis of fracture networks based on graph theory methods

[0101] This invention selects the non-reactive migration and adsorption process of Se(IV) in a high-level radioactive waste disposal project as an example. This example is a common engineering application of DFNs and has certain practical significance for actual engineering. The specific discrete fracture network parameters are shown in Table 1. In this example, the fracture region is a cube with a side length of 10m, and a fracture network containing 151 discrete fractures is constructed.

[0102] Table 1 Parameters of Discrete Fragment Network

[0103]

[0104] Where φ0 is the average dip direction; θ0 is the average dip angle; and l0 is the average fracture diameter.

[0105] Figures 2 to 4Three network topologies, TN, CSN, and DCSN, and their graphical representations are presented. Except for fractures intersecting the source and sink planes, node positions are the projections of the fracture plate centers onto the z=0 plane. Fractures intersecting the inlet and outlet are considered source and target nodes, respectively. Source and sink nodes are represented by equilateral and inverted triangles, respectively, positioned at x=0m and x=10m based on the y-axis coordinates of their corresponding fracture plate centers. To more intuitively represent the seepage within the network, nodes are colored with the average flow velocity on each fracture in the high-fidelity model.

[0106] The outbound traffic of the TN, CSN, and DCSN network topologies is 1.92 × 10⁻⁶. -2 m 3 / d、1.94×10 -2 m 3 / d and 1.88×10 -2 m 3 / d, containing 151, 108, and 77 fracture plates respectively. That is, although the DCSN only contains 51% of the target fracture network, it carries the vast majority of flow transmission. The edges between nodes represent the intersection lines of fractures, and the width of the edge comes from the length of the intersection line. Comparing the average flow velocities of nodes in the three networks reveals that the fracture plates with dominant seepage are basically the same. In the TN, the local average flow of fracture plates not connected to the source-sink surface is approximately 0. In other words, only fracture subnetworks connected to the source-sink surface can effectively induce seepage; the contribution of other fracture plates and fracture subnetworks to seepage can be ignored. In the TN and CSN, the average flow within most dead-end fractures is approximately zero, but some fractures exhibit internal seepage due to pressure from adjacent fractures. In summary, there are clearly dominant channels within the fracture network, mainly those connected to the source-sink surface. Dead-end fractures contribute negligibly to the overall seepage, but local flow exists within them.

[0107] like Figure 5 As shown, the penetration curves (BTC) of CSN and DCSN peak at 11.5a and 11.4a, respectively, with peak values ​​of 7.07 × 10⁻⁶. -5 mol / (m 3 ·d) and 7.13×10 -5 mol / (m 3 ·d). Both BTC exhibit a distinct tail, but as dead-end gaps are incorporated into the network, CSN's BTC is biased towards a later time, with a smaller peak amplitude and a longer tail. Figure 6 The graph shows the adsorption of Se(IV) in the CSN at t=20a, with the colored areas representing the amount of solute adsorbed on the fissures. The results indicate that even in fissure-free networks, there are dominant channels for solute transport.

[0108] 2) Numerical model of seepage and solute transport in a three-dimensional fracture network based on topology pruning

[0109] Based on the global and local analysis results of the fractured network, pruning is based on the DCSN of the target network, and the importance of network nodes is quantified by BC. The results are as follows: Figure 7 As shown, the BC values ​​of the fractures exhibit strong heterogeneity, with only a small number of fractures being relatively important. By using Dijkstra's algorithm to obtain the subnetwork of the shortest path set through important fractures, maximum network pruning can be achieved while ensuring hydraulic connectivity. In summary, this example designs five different pruning schemes. The pruned networks are subnetworks of the shortest path sets of the fracture sets with the top 10%, 30%, 50%, 70%, and 90% BC values ​​in the DCSN, respectively. High-fidelity simulations of seepage and solute transport were performed on the pruned networks.

[0110] The seepage simulation results are shown in Table 2. As the number of high-BC fractures in the pruning scheme increases, the number of fractures in the pruning network also increases. The flow rate shows a logarithmic growth trend relative to the number of fractures, indicating that mass transport in the fracture network is mainly concentrated in the dominant channels. When the pruning network incorporates more high-BC fractures, it exhibits better consistency with the original network, but the efficiency decreases accordingly.

[0111] BTC simulation of solute transport, such as Figure 8 As shown, all networks exhibit a clear tail in their BTC peak values. The 10% and 30% pruned networks show a significant reduction in BTC peak value, demonstrating that these two pruning schemes did not retain enough critical gaps, leading to larger errors in the flow and chemical fields. From 50% to 90% pruned networks, the BTC peak value increased from 6.82 × 10⁻⁶. -5 mol / (m 3 ·d) decreased to 6.69×10 -5 mol / (m 3 ·d), and the time of appearance is delayed by 0.1a. At the same time, the tail of BTC is also significantly elongated, proving that non-critical gaps are included in the pruned network. Therefore, the 50% simplified solution is the most cost-effective simplified solution because it obtains a higher proportion of dominant channels within the subnetwork.

[0112] Table 2. High-fidelity simulation results of seepage in DCSN, TN, and different pruned networks.

[0113]

[0114] In summary, network analysis and Dijkstra's algorithm can effectively prune the entire network and extract key subnetworks. Using the pruned network to establish a high-fidelity seepage and solute transport model, the key seepage and solute transport processes of the original network are restored while significantly reducing unnecessary computation. Therefore, this pruning method provides an efficient and reliable solution to seepage and solute transport problems in hard rock fractured media.

Claims

1. A numerical simulation method for seepage and solute transport in a three-dimensional fractured network based on topology pruning, characterized in that, include: (1) Collect basic data on fracture networks and basic data on hydrogeology, and construct the target three-dimensional fracture network based on the basic data on fracture networks; (2) Use graph theory to map the target three-dimensional fracture network into an equivalent two-dimensional undirected graph, and perform global and local network analysis; (3) Generalize the hydrogeological conditions of the target area, construct and solve the groundwater flow model and solute transport model, and simulate the underground flow field and non-reactive transport and adsorption process of solute under different topological networks; (4) Based on the analysis results obtained in step (2) and the simulation results obtained in step (3), determine the influence mechanism of network topology on seepage field and chemical field, and prune the target fracture network according to the key fractures to extract the key sub-networks that carry most of the material transport. (5) Use key sub-networks to replace the complete network, solve the non-reactive solute migration model, obtain the groundwater flow field, solute breakthrough curve, solute migration field and solid solute adsorption in fractured strata, and realize efficient numerical simulation of seepage and solute transport in three-dimensional fractured networks.

2. The numerical simulation method according to claim 1, characterized in that, In step (1), the basic data of the fracture network includes the fracture orientation distribution characteristics, fracture diameter distribution characteristics, and fracture density distribution characteristics.

3. The numerical simulation method according to claim 2, characterized in that, The method for generating a fractured network is as follows: The crack diameter l follows a power-law distribution, expressed as: f(l)＝βl -α Where f(l) is the probability density function of the crack diameter distribution; β is the probability density function based on the minimum value of the crack diameter l. min and maximum value l max The obtained normalization exponent; α is the power law exponent; The crack diameter l is represented as: Where R is a random number; The fracture orientation follows the Fisher distribution: Where θ is the dip angle; θ0 is the average dip angle; φ is the dip direction; φ0 is the average dip direction; κ is the dispersion parameter; and sinh is the hyperbolic sine function. Each fracture network is configured to have the same aperture.

4. The numerical simulation method according to claim 1, characterized in that, In step (1), the basic hydrogeological data includes structural geological maps and hydrogeological maps.

5. The numerical simulation method according to claim 1, characterized in that, In step (2), graph theory is used to map the target three-dimensional fracture network into an equivalent two-dimensional undirected graph, including: Treating the cracks as nodes V and the intersections between cracks as edges E, the 3D crack network is transformed into a 2D undirected graph G(V,E) as follows: The three-dimensional fracture network F consists of N fractures f, represented as: F={f i },i=1,……,N Define a bijective mapping Φ: F:f i →v i If there are two cracks f i with f j intersect, Then there is an edge in E connecting the corresponding vertex: Among them, (v i ,v j )∈E represents vertex v i and v j The edge e between ij , side e ij weight w ij Set as vertex v i and v j Length of the line of intersection between them: w ij =Lengthf i ∩f j Similarly, considering the flow direction, we can regard the inlet plane x0 as the source node s, and the outlet plane x L Considering the target node t, when the crack intersects with the inlet and outlet boundaries: Among them, e si For an edge between a source node and a non-source / sink node, e it The edges between non-source / sink nodes and the target node are thus transformed into a two-dimensional undirected graph.

6. The numerical simulation method according to claim 5, characterized in that, In step (2), the global network analysis is based on three types of network topologies: the full network TN, the connected subnetwork CSN, and the dead-end-free connected subnetwork DCSN. The target fracture network is called the full network TN. When the three-dimensional fracture network is sparse, the fracture network will be divided into several parts. The fracture network connecting the source surface and the sink surface is called the connected subnetwork CSN. The CSN that does not contain the nodes corresponding to the dead-end fractures is called the dead-end-free connected subnetwork DCSN.

7. The numerical simulation method according to claim 6, characterized in that, Local network analysis, based on betweenness centrality BC, is defined as follows: Where BC(i) represents the betweenness centrality of node i; g jk (i) represents the number of times node i falls on the path between nodes j and k; g jk BC represents the number of times the shortest path between nodes j and k is calculated; the higher the BC value of a node, the more critical it is in the network.

8. The numerical simulation method according to claim 1, characterized in that, In step (3), the governing equations of the fractured network groundwater flow model are established based on the cubic law. The solute transport field includes advection, diffusion, and isothermal adsorption, among which the solute adsorption equation is as follows: Among them, c s It is the solute adsorption capacity of the solid [mol / kg]; K d The equilibrium adsorption partition coefficient of the solute [m] 3 / mol];c s_max c is the maximum amount of solute adsorbed on the solid [mol / kg]; c is the solute concentration [mol / m]. 3 ].

9. The numerical simulation method according to claim 7, characterized in that, In step (4), the pruning process of the fracture network is based on graph theory and network analysis, where the weight of the edge of the graph is assigned as the reciprocal of the length of the corresponding fracture intersection line; local network analysis is performed on the DCSN of the target fracture network, and BC is used to quantify the control capability of a single fracture on the flow field within the network; nodes with high BC values ​​are selected in the target fracture network, and Dijkstra's algorithm is used to solve the shortest seepage path network from the source node to the target node through nodes with high BC values.

10. The numerical simulation method according to claim 1, characterized in that, In step (5), the non-reactive solute migration model is solved by using the finite element method in the commercial physics simulation software COMSOL Multiphysics for spatial and temporal discretization, and the Newton iteration method for solving the discretized nonlinear equations. The breakthrough curve of solute over time is calculated at the outflow boundary surface of the target region to achieve efficient numerical simulation of seepage and solute transport in a three-dimensional fracture network.

Images