Three-dimensional fracture network seepage and solute transport numerical simulation method based on topological pruning

Abstract

The invention discloses a three-dimensional fracture network seepage and solute transport numerical simulation method based on topological pruning, which comprises the following steps of: analyzing a fracture network topological structure from global and local scales by applying a graph theory, establishing a quantitative relation between fracture connectivity and an underground flow field, and accurately describing an influence mechanism of topological characteristics on a seepage field; on the basis of network topology analysis, the key fracture sub-network is extracted through a pruning method, the key fracture sub-network is used for replacing a complete network for high-precision numerical simulation, and the calculation efficiency is remarkably improved while the calculation precision is guaranteed. The invention provides an efficient and reliable solution for the problems of seepage and solute transport in a hard rock fractured medium.

Description

Technical Field

[0001] The present invention relates to a numerical simulation method for a fracture network flow field in an underground stratum, and in particular to a numerical simulation method for seepage and solute transport in a three-dimensional fracture network based on topological pruning. Background Art

[0002] In hard rock, fractures are the primary pathways for groundwater flow and chemical migration, making them crucial for large-scale underground projects. Discrete fracture models (DFNs), a common method for assessing the feasibility of large-scale projects located in fractured formations, can accurately characterize the uncertainty in seepage and mass transfer caused by the heterogeneity of fracture distribution, but they are computationally expensive. Furthermore, due to the low permeability of hard rock and the heterogeneous connectivity of fracture networks, fluid flow within fractured rock formations often exhibits distinct, dominant paths. Summary of the Invention

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

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

[0005] (1) Collect basic data on fracture networks and hydrogeological data, 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 solute non-reactive transport and adsorption process 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 mechanism by which the network topology affects the seepage field and the chemical field, and prune the target fracture network according to the key fractures to extract the key sub-network that carries most of the material migration;

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

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

[0011] Furthermore, the fracture network is generated as follows:

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

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

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

[0015] The crack diameter l is expressed as:

[0016]

[0017] Where R is a random number;

[0018] The fracture occurrence follows the Fisher distribution:

[0019]

[0020] Where θ is the inclination; θ0 is the average inclination; φ is the dip; φ0 is the average dip; κ is the dispersion parameter; sinh is the hyperbolic sine function;

[0021] Sets the fractures in each fracture network to have the same opening.

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

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

[0024] Considering the fractures as nodes V and the intersections between fractures as edges E, the three-dimensional fracture network is converted into a two-dimensional undirected graph G(V,E) as follows:

[0025] The three-dimensional fracture network F consists of N fractures f, which can be expressed as:

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

[0027] Define a bijective mapping Φ:

[0028] Φ:f i →v i

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

[0030]

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

[0032] w ij =Lengthf i ∩f j

[0033] Similarly, considering the flow direction, the inlet plane x0 is regarded as the source node s, and the outlet plane x L Considered as the target node t, when the crack intersects the entry and exit boundaries:

[0034]

[0035] Among them, e si is the edge between the source node and the non-source sink node, e it is the edge between the non-source node and the target node; at this point, the three-dimensional fracture network is converted into a two-dimensional undirected graph.

[0036] Furthermore, in step (2), the global network analysis is based on three topological networks, namely the full network TN, the connected subnetwork CSN, and the dead-end-free connected subnetwork DCSN, where 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, and 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 is based on betweenness centrality (BC), which is defined as:

[0038]

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

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

[0041]

[0042] Among them, c s is the solute adsorption capacity of the solid [mol / kg]; K d is the equilibrium adsorption distribution coefficient of the solute [m 3 / mol]; c s_max is the maximum adsorption capacity of the solute 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 methods and network analysis, where the weights of the edges of the graph are assigned the inverse of the lengths of the corresponding fracture intersections; a local network analysis is performed on the DCSN of the target fracture network, and the control ability of a single fracture in the network over the flow field is quantified using BC; nodes with high BC values ​​in the target fracture network are selected, and the Dijkstra algorithm is used to solve the shortest seepage path network from the source node to the target node through the nodes with high BC values.

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

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

[0046] This paper uses graph theory and network analysis to determine how topological characteristics influence the seepage and solute transport fields of a fracture network. Based on betweenness centrality and the Dijkstra algorithm, fractures are pruned to create a high-throughput fracture subnetwork. This approach reduces both the complexity of the discrete fracture model network and computational cost. When used in high-precision numerical simulations, this high-throughput fracture subnetwork significantly improves computational efficiency while maintaining accuracy.

[0047] The present invention provides an efficient and reliable solution to the problems of seepage and solute migration in hard rock fracture media. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

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

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

[0053] Figure 6 The adsorption of Se(IV) on CSN in the embodiment of the present invention (t=20a) and its graphical representation;

[0054] Figure 7 is a graphical representation of the DCSN in an embodiment of the present 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 embodiment of the present invention. DETAILED DESCRIPTION

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

[0057] like Figure 1 As shown, an embodiment of the present invention provides a method for numerical simulation of 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 hydrogeological data, and construct the target three-dimensional fracture network based on the basic data on fracture networks;

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

[0060] The basic data of the fracture network include the fracture occurrence distribution characteristics, fracture diameter distribution characteristics and fracture density distribution characteristics. The fracture network is generated as follows:

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

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

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

[0064] The crack diameter l is expressed as:

[0065]

[0066] Where R is a random number;

[0067] The fracture occurrence follows the Fisher distribution:

[0068]

[0069] Where θ is the inclination; θ0 is the average inclination; φ is the dip; φ0 is the average dip; κ is the dispersion parameter; sinh is the hyperbolic sine function;

[0070] Sets the fractures in each fracture network to have the same opening.

[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] Graph theory is used to map the target 3D fracture network into an equivalent 2D undirected graph, including:

[0073] Considering the fractures as nodes V and the intersections between fractures as edges E, the three-dimensional fracture network is converted into a two-dimensional undirected graph G(V,E) as follows:

[0074] The three-dimensional fracture network F consists of N fractures f, which can be expressed as:

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

[0076] Define a bijective mapping Φ:

[0077] Φ:f i →v i

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

[0079]

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

[0081] w ij =Lengthf i ∩f j

[0082] Similarly, considering the flow direction, the inlet plane x0 is regarded as the source node s, and the outlet plane x L Considered as the target node t, when the crack intersects the entry and exit boundaries:

[0083]

[0084] Among them, e si is the edge between the source node and the non-source sink node, e it is the edge between the non-source node and the target node; at this point, the three-dimensional fracture network is converted into a two-dimensional undirected graph.

[0085] Global network analysis is based on three topological networks: 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, 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 is based on Betweenness centrality (BC), which is defined as:

[0087]

[0088] Where BC(i) is the betweenness centrality of node i; g jk (i) is the number of times node i falls on the path between nodes j and k; g jkis the number of shortest paths between nodes j and k; BC measures the influence of a node in the network. 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 solute non-reactive transport and adsorption process under different topological networks;

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

[0091]

[0092] Among them, c s is the solute adsorption capacity of the solid [mol / kg]; K d is the equilibrium adsorption distribution coefficient of the solute [m 3 / mol]; c s_max is the maximum adsorption capacity of the solute 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 mechanism by which the network topology affects the seepage field and the chemical field, and prune the target fracture network according to the key fractures to extract the key sub-network that carries most of the material migration;

[0094] The fracture network pruning process is based on graph theory and network analysis, where graph edges are weighted by the inverse of the length of the corresponding fracture intersection. A local network analysis is performed on the DCSN of the target fracture network, using the BC value to quantify the control ability of individual fractures within the network over the flow field. Nodes with high BC values ​​within the target fracture network are selected, and the Dijkstra algorithm is used to find the shortest flow path from the source node to the target node through nodes with high BC values. The pruned network has the largest pressure gradient, indicating the highest probability of preferential flow.

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

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

[0097] An embodiment of the present invention provides a numerical simulation method for seepage and solute transport in a three-dimensional fracture network based on topological pruning. It uses network analysis and the Dijkstra algorithm to extract a three-dimensional fracture high-pass quantum network, greatly reducing the complexity of the discrete fracture model and alleviating the computational load of the underground flow and chemical fields. This provides an efficient and reliable solution to the problems of seepage and solute transport in hard rock fracture media.

[0098] A specific example is given below.

[0099] This example develops a numerical simulation method for seepage and solute transport in three-dimensional fracture networks based on topological pruning, taking into account the topological properties of the fracture network. It mainly includes two parts: 1) global and local topological analysis of the fracture network based on graph theory methods; 2) a numerical model for seepage and solute transport in three-dimensional fracture networks based on topological pruning.

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

[0101] This paper uses 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 projects. The specific discrete fracture network parameters are shown in Table 1. The fracture region in this example is a cube with a side length of 10m, and a fracture network consisting of 151 discrete fractures was constructed.

[0102] Table 1 Discrete fracture network parameters

[0103]

[0104] Among them, φ0 is the average inclination; θ0 is the average dip; l0 is the average crack diameter.

[0105] Figures 2 to 4Three topological networks, TN, CSN, and DCSN, and their graphical representations are shown. Except for fractures intersecting the source-sink plane, node locations are the projections of the fracture slice center 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, and are placed at x = 0 m and x = 10 m, based on the y-axis coordinates of their corresponding fracture slice centers. To more intuitively represent the percolation within the network, nodes are colored by the average flow velocity across each fracture in the high-fidelity model.

[0106] The egress traffic of the three topology networks of TN, CSN and DCSN 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 slices, respectively. This means that although the DCSN only contains 51% of the target fracture network, it carries the vast majority of flow. Edges between nodes represent fracture intersections, and the edge width is derived from the length of the intersection. Comparing the average flow velocities of nodes across the three networks reveals that the fracture slices with dominant seepage are essentially identical. In the TN, the local average flow rate of fracture slices that are not connected to the source / sink interface is approximately zero. In other words, only fracture subnetworks connected to the source / sink interface can effectively induce seepage, and the contribution of other fracture slices and subnetworks to seepage can be ignored. The average flow rate within most dead-end fractures in both the TN and CSN is approximately zero, but some fractures experience internal seepage due to pressure from adjacent fractures. In summary, there are distinct dominant pathways within the fracture network, primarily those connecting to the source / sink interface. Dead-end fractures contribute minimally to the overall seepage, but localized flow within them does occur.

[0107] like Figure 5 As shown in the figure, the penetration curves (BTC) of CSN and DCSN reach their peak values ​​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 BTCs show a distinct tail, but as dead-end rifts are incorporated into the network, CSN’s BTC shifts to a later time, with a smaller peak amplitude and a longer tail. Figure 6 The figure shows the Se(IV) adsorption on the CSN at t = 20a and its graphical representation, where the color represents the amount of solute adsorbed on the cracks. The results show that there are also advantageous pathways for solute transport in the dead-end crack network.

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

[0109] Based on the global and local analysis results of the crack network, pruning is based on the DCSN of the target network and the importance of network nodes is quantified by BC. Figure 7 As shown, the BC values ​​of fractures are highly heterogeneous, with only a small number of fractures being significant. By using the Dijkstra algorithm to obtain a subnetwork consisting of the shortest paths through important fractures, we can achieve maximum network pruning while ensuring hydraulic connectivity. In summary, this example designed five different pruning schemes. The pruned networks consist of subnetworks consisting of the shortest paths through the top 10%, 30%, 50%, 70%, and 90% of the fracture sets with the top BC values ​​in the DCSN, respectively. High-fidelity simulations of seepage and solute transport were then performed on these 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 pruned network also increases. The flow rate exhibits a logarithmic growth trend relative to the number of fractures, indicating that mass transfer in the fracture network is primarily concentrated in the dominant channels. When the pruned network incorporates more high-BC fractures, it exhibits better consistency with the original network, but its efficiency decreases.

[0111] BTC for solute transport simulation Figure 8 As shown in Figure 2, the BTC of all networks exhibits a significant tail. The BTC peaks of the 10% and 30% pruned networks are significantly reduced, proving that these two pruning schemes do not retain enough key cracks, resulting in large errors in the flow and chemical fields. From 50% to 90% pruned networks, the BTC peak decreases from 6.82×10 -5 mol / (m 3 d) decreased to 6.69×10 -5 mol / (m 3 d), and its appearance is delayed by 0.1a. Simultaneously, the tail of BTC is significantly lengthened, demonstrating that non-critical cracks are incorporated into the pruned network. Therefore, the 50% reduction scheme is the most cost-effective, as it yields a higher proportion of advantageous channels within the obtained subnetwork.

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

[0113]

[0114] Overall, network analysis and the Dijkstra algorithm effectively pruned the entire network and extracted key subnetworks. Using the pruned network, a high-fidelity seepage and solute transport model was established, eliminating unnecessary computational effort while restoring the key seepage and solute transport processes in the original network. Therefore, this pruning method provides an efficient and reliable solution to the seepage and solute transport problems in hard rock fractured media.

Claims

1. A numerical simulation method for seepage and solute transport in three-dimensional fracture networks based on topological pruning, characterized in that: include: (1) Collect basic data on fracture networks and hydrogeological data, 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 solute non-reactive transport and adsorption process under different topological networks; (4) Based on the analysis results obtained in step (2) and the simulation results obtained in step (3), determine the mechanism by which the network topology affects the seepage field and the chemical field, and prune the target fracture network according to the key fractures to extract the key sub-network that carries most of the material migration; (5) The key subnetwork is used to replace the complete network to solve the non-reactive solute migration model, obtain the groundwater flow field, solute breakthrough curve, solute migration field and solid solute adsorption in the fractured formation, and realize the efficient numerical simulation of seepage and solute transport in the three-dimensional fracture network.

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

3. The numerical simulation method according to claim 2, characterized in that: The fracture network is generated as follows: The crack diameter l follows a power law distribution, which can be expressed as: f(l)＝βl -α Where f(l) is the probability density function of the crack diameter distribution; β is the minimum value l of the crack diameter. min and the maximum value l max The normalized exponent obtained; α is the power law exponent; The crack diameter l is expressed as: Where R is a random number; The fracture occurrence follows the Fisher distribution: Where θ is the inclination; θ0 is the average inclination; φ is the dip; φ0 is the average dip; κ is the dispersion parameter; sinh is the hyperbolic sine function; Sets the fractures in each fracture network to have the same opening.

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

5. The numerical simulation method according to claim 1, characterized in that: In step (2), the target three-dimensional fracture network is mapped into an equivalent two-dimensional undirected graph using graph theory, including: Considering the fractures as nodes V and the intersections between fractures as edges E, the three-dimensional fracture network is converted into a two-dimensional undirected graph G(V,E) as follows: The three-dimensional fracture network F consists of N fractures f, which can be expressed as: F={f i },i=1,……,N Define a bijective mapping Φ: F:f i →v i If the two cracks f i With f j intersect, Then there is an edge in E connecting the corresponding vertices: Among them, (v i ,v j )∈E represents vertex v i and v j The edge between ij , side e ij The weight w ij Set to vertex v i and v j The length of the intersection between: w ij =Lengthf i ∩f j Similarly, considering the flow direction, the inlet plane x0 is regarded as the source node s, and the outlet plane x L Considered as the target node t, when the crack intersects the entry and exit boundaries: Among them, e si is the edge between the source node and the non-source sink node, e it is the edge between the non-source node and the target node; at this point, the three-dimensional fracture network is converted 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 topological networks: 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 is based on betweenness centrality (BC), which is defined as: Where BC(i) is the betweenness centrality of node i; g jk (i) is the number of times node i falls on the path between nodes j and k; g jk is the number of the shortest path between nodes j and k; 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 equation of the fracture network groundwater flow model is established based on the cubic law. The solute transport field includes advection, diffusion, and isothermal adsorption, where the solute adsorption equation is as follows: Among them, c s is the solute adsorption capacity of the solid [mol / kg]; K d is the equilibrium adsorption distribution coefficient of the solute [m 3 / mol]; c s_max is the maximum adsorption capacity of the solute 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 the inverse of the length of the corresponding fracture intersection line; a local network analysis is performed on the DCSN of the target fracture network, and the control ability of a single fracture in the network over the flow field is quantified using BC; nodes with high BC values ​​in the target fracture network are selected, and the Dijkstra algorithm is used to solve the shortest seepage path network from the source node to the target node through the 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 solute breakthrough curve over time is calculated at the outflow boundary surface of the target area to achieve efficient numerical simulation of seepage and solute transport in three-dimensional fracture networks.

Images