Method and system for simulating transport of reactive solute in groundwater in three-dimensional fracture network

Abstract

Disclosed in the present invention is a method for simulating the transport of a reactive solute in groundwater in a three-dimensional fracture network, which method is used for simulating a reactive transport pattern of a solute in groundwater in a fracture medium. The method can take into consideration the detailed depiction of the transport process of a reactive solute in a fracture and the hydraulic connection between the fracture and surrounding rock, thereby affecting the spatial distribution of a water flow field and a chemical field in a fracture network; and a water flow field and chemical field coupling model in a three-dimensional fracture network is constructed to simulate the hydraulic connection between the fracture and the surrounding rock, so as to study the transport process of a reactive solute in groundwater in the fracture network. The method provided in the present invention comprises: taking into consideration the transport of a reactive solute in groundwater in a fracture and the hydraulic connection between the fracture and surrounding rock; constructing a model for the process of coupling a water flow field and a chemical field; and carrying out numerical simulation on the transport of the reactive solute in groundwater in a three-dimensional fracture network. The simulation technique provided in the present invention can provide a scientific and rational theoretical basis for the study of space-time distribution of reactive transport of a target solute in groundwater in a fracture network.

Description

A method and system for simulating reactive solute transport in groundwater in a three-dimensional fracture network Technical Field

[0001] The present invention relates to a three-dimensional fracture network groundwater reactive solute migration simulation method and system, belonging to the field of earth science and engineering. Background Art

[0002] In areas with well-developed fracture networks, the transport of reactive solutes in aquifers is highly complex. Research on the transport and transformation of solutes in bedrock fracture water is closely related to the assessment and prediction of fracture water contamination, as well as the risk assessment of geological disposal of nuclear waste and greenhouse gases. Quantitatively describing the transport of reactive solutes in fractured media has long been a hot topic and a challenge in international water environment and hydrogeology research. Unlike the transport mechanism of inert solutes, reactive solute transport is a complex process that requires simultaneous consideration of water flow, solute transport, and chemical reactions. This complexity is particularly pronounced in highly heterogeneous fractured media. Reactive solute transport in fractured media exhibits "anomalous" transport that deviates from the traditional reactive convection-diffusion equation. This is manifested by premature arrival of products, tailing, and lower-than-predicted average concentrations. This leads to practical problems such as prolonged pollutant remediation cycles and lower-than-expected remediation efficiency. Recent studies both domestically and internationally have demonstrated that the fractured rock matrix serves as an important hydrogeochemical buffer zone and significantly retards solute transport in the surrounding fractured media. Considering the geochemical reactions and mass exchange between these two regions is key to accurately characterizing the overall buffering capacity of fractured groundwater in a region and the related hydrogeochemical evolution processes.

[0003] However, current research primarily studies reactive transport processes in fractured media at small scales through laboratory experiments and models, with limited quantitative characterization of reactive solute transport processes in fractured groundwater at large field scales. The proposed method can consider the influence of fractures and fractured rock mass on groundwater solute transport, effectively characterizing the spatiotemporal variations of reactive solute transport in groundwater within a three-dimensional fracture network. This provides key regional information for groundwater pollution monitoring, thereby better guiding field work on the prevention and control of pollutant migration and diffusion in fractured areas, as well as their subsequent remediation and treatment.

[0004] Summary of the Invention

[0005] Purpose of the invention: In view of the deficiencies of the above-mentioned prior art, the purpose of the present invention is to provide a method and system for simulating the migration of reactive solutes in groundwater in a three-dimensional fracture network. This method can take into account the migration of reactive solutes in fractured groundwater and the hydraulic connection between fractures and surrounding rock masses, construct a model of coupled water flow field-chemical field processes, and carry out numerical simulation of the migration of reactive solutes in groundwater in a three-dimensional fracture network.

[0006] Technical solution: The present invention proposes a method for simulating reactive solute transport in groundwater in a three-dimensional fracture network, which includes the following steps:

[0007] Step SS1, collecting geological structure data and local hydrogeological parameters of the fractured plots in the target study area;

[0008] Step SS2: Based on the geological structure data collected in step SS1, the data is modeled and parameterized, the lithologic parameters of the formation are set, a discrete fracture network is generated, and the plane fracture grid is upscaled to a three-dimensional cubic grid, mapped into an equivalent continuous porous medium, and divided into a three-dimensional fracture body grid to establish a geological structure model of the study area;

[0009] Step SS3: Based on the geological structure model constructed in step SS2, the initial hydraulic head conditions, boundary conditions, and permeability coefficients of fractures and fracture bodies of the model are determined according to the collected hydrogeological parameters, and the spatiotemporal distribution of the groundwater level in the study area is solved to obtain the initial groundwater flow field.

[0010] Step SS4, based on the initial groundwater flow field obtained in step SS3, setting the pollutant physicochemical parameters and source-sink phase conditions to determine the spatiotemporal distribution of solute concentration in the fractured groundwater;

[0011] In step SS5, based on the solute transport model of step SS4, relevant chemical reaction parameters are set, and the basic equations of conservation of mass and energy, the seepage field equations, and the solute transport control equations are used to couple the water flow field and the chemical field to solve the results of reactive solute transport in the fractured medium driven by hydrodynamic and hydrogeochemical reactions, thereby completing the numerical simulation of reactive solute transport in groundwater in a three-dimensional fracture network.

[0012] Furthermore, in step SS1, the data include the following for the target study area: (1) hydrogeological conditions, stratum lithology and thickness, burial depth, thickness and spatial distribution of aquifers and aquicludes, recharge and drainage conditions of hydrogeological units, groundwater ion content, groundwater flow rate, and hydrogeochemical characteristics; (2) fracture medium conditions, the overall number and spatial distribution of fractures, the occurrence, opening, development depth, roughness of fractures, physicochemical properties of fracture fillings and fracture geological bodies, and spatial distribution of fracture network permeability coefficient.

[0013] Furthermore, in step SS2, the three-dimensional fracture network model is constructed as follows:

[0014] The characteristic parameters of the three-dimensional fracture network geological structure model are determined by using the discrete fracture network parameter control equations and the model upscaling grid subdivision control equations;

[0015] (1) Discrete fracture network parameter control equations

[0016] a. The direction of the normal vector of the crack surface is characterized by the von Mises-Fisher distribution function;

[0017] in, is the crack average normal vector matrix; κ is the density coefficient; T represents the matrix transpose; I is the identity matrix;

[0018] b. The crack radius R adopts a truncated power law distribution, defined by the following probability density function;

[0019] Among them, R u and R l are the upper and lower bounds of the crack radius, respectively, and α is a constant coefficient;

[0020] c. The crack opening φ is determined by the exponential function of the crack radius:

[0021] φ=F·R γ (3)

[0022] Among them, F and γ are parameter values ​​related to water permeability;

[0023] d. Fracture permeability k is calculated based on fracture opening φ;

[0024] (2) Model upscaling and meshing control equations

[0025] a. The porosity n of the fractured medium is calculated based on the opening of the fractures in each grid;

[0026] Where l is the length of the mesh edge, N is the total number of cracks passing through the mesh, and φ i represents the opening of the i-th crack in the target grid, which is calculated by formula (3);

[0027] b. Fracture medium permeability k f Assign values ​​based on the permeability of each fracture in the reference grid;

[0028] Among them, k i The permeability of the i-th crack in the target grid is calculated by formula (4).

[0029] Furthermore, in step SS4, the pollutant parameters are set to include: groundwater background ion components and concentrations, pollution source ion components and concentrations, and source-sink boundary conditions; in step SS5, the chemical reaction parameters include: groundwater ion component reaction equations and reaction equilibrium constants, mineral components and their volume fractions, reaction rate constants, and mineral phase reaction specific surface areas.

[0030] Furthermore, in step SS5, the groundwater reactive solute transport model is constructed as follows:

[0031] The seepage field equation, solute transport control equation, mass and energy conservation equation are used to calculate the spatiotemporal evolution of target pollutants in the model driven by hydrodynamic and hydrogeochemical effects.

[0032] a. The seepage field equation is:

[0033] Where n is porosity; s is saturation; ρ is fluid density; w is infiltration recharge; k is absolute permeability; k r is the relative permeability; μ is the fluid viscosity coefficient; P is the pressure; g is the acceleration of gravity; z is the position head;

[0034] b. The solute transport control equation is derived based on the law of conservation of mass and is expressed as:

[0035] Among them, C j is the concentration of j particles; C i ′ is the concentration of the secondary species i particle corresponding to the j particle; υ ji is the stoichiometric number of all liquid phase reactions involving particles j and i; N i is the total number of particles of minor species i; q is the groundwater velocity; D is the molecular diffusion tensor; Q j is the source and sink term of particle j; I m is the mineral reaction rate; jm ′ is the stoichiometric number of particle j related to the mineral phase reaction; M is the total number of particles j involved in the mineral phase reaction; particle j represents the particle with the largest relative concentration in the liquid phase, and the corresponding minor species i particle represents the particle j ionized in the liquid phase or combined with other group particles, for example, particle j is Ca 2+ When the corresponding secondary species i particles include CaCO3(aq), CaHCO3 + 、Ca(OH) + 、CaCl + , CaCl2(aq), CaSO4(aq), etc.; j particle is HCO3 - When the corresponding secondary species i particles include CO3 2- , CO2(aq), etc.

[0036] c. Mineral phase reaction rate I m Calculated by the law of conservation of mass;

[0037] in, is the mineral volume fraction, is the molar volume of the mineral;

[0038] d. According to the law of conservation of mass and energy of equilibrium reaction, the concentration of j particles of the same type is C j The corresponding secondary species i particle concentration C i Quantitative relationship of ′:

[0039] Among them, K i is the reaction equilibrium constant of minor species i, γ j and γ i are the activity coefficients of particle j and its corresponding secondary species particle i, respectively.

[0040] The present invention also proposes a three-dimensional fracture network groundwater reactive solute transport simulation system, which includes the following modules:

[0041] Model parameter collection module, used to collect geological structure data and local hydrogeological parameters of fractured plots in the target study area;

[0042] The 3D fracture network construction module is used to perform modeling and parameterization processing on the collected geological structure data, set the lithologic parameters of the formation, generate a discrete fracture network, and upscale the plane fracture grid into a 3D cubic grid, mapping it into an equivalent continuous porous medium, and then dividing it into a 3D fracture body grid to establish a geological structure model of the study area.

[0043] The groundwater flow module is used to determine the initial head conditions, boundary conditions, and permeability coefficients of fractures and fracture bodies based on the established geological structure model and the collected hydrogeological parameters, and to solve the spatiotemporal distribution of groundwater levels in the study area to obtain the initial groundwater flow field.

[0044] The solute transport module is used to set the physical and chemical parameters of pollutants and the source-sink phase conditions based on the initial groundwater flow field obtained, and to determine the spatiotemporal distribution of solute concentration in fractured groundwater;

[0045] The fracture reaction transport module is used to set relevant chemical reaction parameters based on the step-by-step solute transport model. According to the basic equations of conservation of mass and energy, the seepage field equations, and the solute transport control equations, the flow field-chemical field coupling is used to solve the results of reactive solute transport in fractured media driven by hydrodynamic and hydrogeochemical reactions, and complete the numerical simulation of reactive solute transport in groundwater in a three-dimensional fracture network.

[0046] Furthermore, the data include the following for the target study area: (1) hydrogeological conditions, stratum lithology and thickness, burial depth, thickness and spatial distribution of aquifers and aquicludes, recharge and drainage conditions of hydrogeological units, groundwater ion content, groundwater flow rate, and hydrogeochemical characteristics; (2) fracture medium conditions, the overall number and spatial distribution of fractures, fracture occurrence, opening, development depth, roughness, physical and chemical properties of fracture fillings and fracture geological bodies, and spatial distribution of fracture network permeability coefficient.

[0047] Furthermore, the three-dimensional fracture network model is constructed as follows:

[0048] The characteristic parameters of the three-dimensional fracture network geological structure model are determined by using the discrete fracture network parameter control equations and the model upscaling grid subdivision control equations;

[0049] (1) Discrete fracture network parameter control equations

[0050] a. The direction of the normal vector of the crack surface is characterized by the von Mises-Fisher distribution function;

[0051] in, is the crack average normal vector matrix; κ is the density coefficient; T represents the matrix transpose; I is the identity matrix;

[0052] b. The crack radius R adopts a truncated power law distribution, defined by the following probability density function;

[0053] Among them, R u and R l are the upper and lower bounds of the crack radius, respectively, and α is a constant coefficient;

[0054] c. The crack opening φ is determined by the exponential function of the crack radius:

[0055] φ=F·R γ (3)

[0056] Among them, F and γ are parameter values ​​related to water permeability;

[0057] d. Fracture permeability k is calculated based on fracture opening φ;

[0058] (2) Model upscaling and meshing control equations

[0059] a. The porosity n of the fractured medium is calculated based on the opening of the fractures in each grid;

[0060] Where l is the length of the mesh edge, N is the total number of cracks passing through the mesh, and φ irepresents the opening of the i-th crack in the target grid, which is calculated by formula (3);

[0061] b. Fracture medium permeability k f Assign values ​​based on the permeability of each fracture in the reference grid;

[0062] Among them, k i The permeability of the i-th crack in the target grid is calculated by formula (4).

[0063] Furthermore, the pollutant parameters are set to include: groundwater background ion components and concentrations, pollution source ion components and concentrations, and source-sink boundary conditions; in step SS5, the chemical reaction parameters include: groundwater ion component reaction equations and reaction equilibrium constants, mineral components and their volume fractions, reaction rate constants, and mineral phase reaction specific surface areas.

[0064] Furthermore, the groundwater reactive solute transport model is constructed as follows:

[0065] The seepage field equation, solute transport control equation, mass and energy conservation equation are used to calculate the spatiotemporal evolution of target pollutants in the model driven by hydrodynamic and hydrogeochemical effects.

[0066] a. The seepage field equation is:

[0067] Where n is porosity; s is saturation; ρ is fluid density; w is infiltration recharge; k is absolute permeability; k r is the relative permeability; μ is the fluid viscosity coefficient; P is the pressure; g is the acceleration of gravity; z is the position head;

[0068] b. The solute transport control equation is derived based on the law of conservation of mass and is expressed as:

[0069] Among them, C j is the concentration of j particles; C i ′ is the concentration of the secondary species i particle corresponding to the j particle; υ ji is the stoichiometric number of all liquid phase reactions involving particles j and i; N i is the total number of particles of minor species i; q is the groundwater velocity; D is the molecular diffusion tensor; Q j is the source and sink term of particle j; I m is the mineral reaction rate; jm ′ is the stoichiometric number of particle j related to the mineral phase reaction; M is the total number of particles j involved in the mineral phase reaction; particle j represents the particle with the largest relative concentration in the liquid phase, and the corresponding minor species i particle represents the particle j ionized in the liquid phase or combined with other group particles, for example, particle j is Ca2+ When the corresponding secondary species i particles include CaCO3(aq), CaHCO3 + 、Ca(OH) + 、CaCl + , CaCl2(aq), CaSO4(aq), etc.; j particle is HCO3 - When the corresponding secondary species i particles include CO3 2- , CO2(aq), etc.

[0070] c. Mineral phase reaction rate I m Calculated by the law of conservation of mass;

[0071] in, is the mineral volume fraction, is the molar volume of the mineral;

[0072] d. According to the law of conservation of mass and energy of equilibrium reaction, the concentration of j particles of the same type is C j The corresponding secondary species i particle concentration C i Quantitative relationship of ′:

[0073] Among them, K i is the reaction equilibrium constant of minor species i, γ j and γ i are the activity coefficients of particle j and its corresponding secondary species particle i, respectively.

[0074] Beneficial effects: Compared with the prior art, the technical solution of the present invention has the following beneficial effects:

[0075] The three-dimensional fracture network groundwater reactive solute transport model in the present invention can take into account the influence of fractures and fractured rock masses on groundwater solute transport, quantitatively simulate fracture water solute transport and related hydrogeochemical processes in large and medium-sized sites, effectively characterize the temporal and spatial variation patterns of groundwater reactive solute transport in three-dimensional fracture networks, provide key regional information for groundwater pollution monitoring, and better guide on-site work on the prevention and control of pollutant migration and diffusion in fracture-developed areas and their subsequent remediation and treatment. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] Figure 1 is a flow chart of the numerical simulation method for reactive solute transport in groundwater in a three-dimensional fracture network;

[0077] Figure 2 is a conceptual model diagram of reactive solute transport in groundwater in a three-dimensional fracture network;

[0078] Figure 3 is a diagram of a three-dimensional discrete fracture network model;

[0079] Figure 4 is a diagram of a three-dimensional fracture network equivalent continuous porous medium geological structure model;

[0080] Figure 5 shows the spatial distribution of the overall permeability of the three-dimensional fracture network;

[0081] Figure 6 shows the spatial distribution of the water pressure field in the three-dimensional fracture network as a whole (a) and in the fractures (b);

[0082] Figure 7 shows the spatiotemporal distribution of pH in a three-dimensional fracture network;

[0083] Figure 8 shows the spatial distribution of the concentrations of the main groundwater components in the three-dimensional fracture network in the 10th year of simulation. DETAILED DESCRIPTION

[0084] The present invention will be further described below in conjunction with the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solutions of the present invention and are not intended to limit the scope of protection of the present invention.

[0085] As shown in FIG1 , the present invention proposes a method for simulating reactive solute transport in groundwater in a three-dimensional fracture network, which includes the following steps:

[0086] Step SS1, collecting geological structure data and local hydrogeological parameters of the fractured plots in the target study area;

[0087] Step SS2: Based on the geological structure data collected in step SS1, the data is modeled and parameterized, the lithologic parameters of the formation are set, a discrete fracture network is generated, and the plane fracture grid is upscaled to a three-dimensional cubic grid, mapped into an equivalent continuous porous medium, and divided into a three-dimensional fracture body grid to establish a geological structure model of the study area;

[0088] Step SS3: Based on the geological structure model constructed in step SS2, the initial hydraulic head conditions, boundary conditions, and permeability coefficients of fractures and fracture bodies of the model are determined according to the collected hydrogeological parameters, and the spatiotemporal distribution of the groundwater level in the study area is solved to obtain the initial groundwater flow field.

[0089] Step SS4, based on the initial groundwater flow field obtained in step SS3, setting the pollutant physicochemical parameters and source-sink phase conditions to determine the spatiotemporal distribution of solute concentration in the fractured groundwater;

[0090] In step SS5, based on the solute transport model of step SS4, relevant chemical reaction parameters are set, and the basic equations of conservation of mass and energy, the seepage field equations, and the solute transport control equations are used to couple the water flow field and the chemical field to solve the results of reactive solute transport in the fractured medium driven by hydrodynamic and hydrogeochemical reactions, thereby completing the numerical simulation of reactive solute transport in groundwater in a three-dimensional fracture network.

[0091] Furthermore, in step SS1, the data include the following for the target study area: (1) hydrogeological conditions, stratum lithology and thickness, burial depth, thickness and spatial distribution of aquifers and aquicludes, recharge and drainage conditions of hydrogeological units, groundwater ion content, groundwater flow rate, and hydrogeochemical characteristics; (2) fracture medium conditions, the overall number and spatial distribution of fractures, the occurrence, opening, development depth, roughness of fractures, physicochemical properties of fracture fillings and fracture geological bodies, and spatial distribution of fracture network permeability coefficient.

[0092] Furthermore, in step SS2, the three-dimensional fracture network model is constructed as follows:

[0093] The characteristic parameters of the three-dimensional fracture network geological structure model are determined by using the discrete fracture network parameter control equations and the model upscaling grid subdivision control equations;

[0094] (1) Discrete fracture network parameter control equations

[0095] a. The direction of the normal vector of the crack surface is characterized by the von Mises-Fisher distribution function;

[0096] in, is the crack average normal vector matrix; κ is the density coefficient; T represents the matrix transpose; I is the identity matrix;

[0097] b. The crack radius R adopts a truncated power law distribution, defined by the following probability density function;

[0098] Among them, R u and R l are the upper and lower bounds of the crack radius, respectively, and α is a constant coefficient;

[0099] c. The crack opening φ is determined by the exponential function of the crack radius:

[0100] φ=F·R γ (3)

[0101] Among them, F and γ are parameter values ​​related to water permeability;

[0102] d. Fracture permeability k is calculated based on fracture opening φ;

[0103] (2) Model upscaling and meshing control equations

[0104] a. The porosity n of the fractured medium is calculated based on the opening of the fractures in each grid;

[0105] Where l is the length of the mesh edge, N is the total number of cracks passing through the mesh, and φi represents the opening of the i-th crack in the target grid, which is calculated by formula (3);

[0106] b. Fracture medium permeability k f Assign values ​​based on the permeability of each fracture in the reference grid;

[0107] Among them, k i The permeability of the i-th crack in the target grid is calculated by formula (4).

[0108] Furthermore, in step SS4, the pollutant parameters are set to include: groundwater background ion components and concentrations, pollution source ion components and concentrations, and source-sink boundary conditions; in step SS5, the chemical reaction parameters include: groundwater ion component reaction equations and reaction equilibrium constants, mineral components and their volume fractions, reaction rate constants, and mineral phase reaction specific surface areas.

[0109] Furthermore, in step SS5, the groundwater reactive solute transport model is constructed as follows:

[0110] The seepage field equation, solute transport control equation, mass and energy conservation equation are used to calculate the spatiotemporal evolution of target pollutants in the model driven by hydrodynamic and hydrogeochemical effects.

[0111] a. The seepage field equation is:

[0112] Where n is porosity; s is saturation; ρ is fluid density; w is infiltration recharge; k is absolute permeability; k r is the relative permeability; μ is the fluid viscosity coefficient; P is the pressure; g is the acceleration of gravity; z is the position head;

[0113] b. The solute transport control equation is derived based on the law of conservation of mass and is expressed as:

[0114] Among them, C j is the concentration of j particles; C i ′ is the concentration of the secondary species i particle corresponding to the j particle; υ ji is the stoichiometric number of all liquid phase reactions involving particles j and i; N i is the total number of particles of minor species i; q is the groundwater velocity; D is the molecular diffusion tensor; Q j is the source and sink term of particle j; I m is the mineral reaction rate; jm′ is the stoichiometric number of particle j related to the mineral phase reaction; M is the total number of particles j involved in the mineral phase reaction; particle j represents the particle with the largest relative concentration in the liquid phase, and the corresponding minor species i particle represents the particle j ionized in the liquid phase or combined with other group particles, for example, particle j is Ca 2+ When the corresponding secondary species i particles include CaCO3(aq), CaHCO3 + 、Ca(OH) + 、CaCl + , CaCl2(aq), CaSO4(aq), etc.; j particle is HCO3 - When the corresponding secondary species i particles include CO3 2- , CO2(aq), etc.

[0115] c. Mineral phase reaction rate I m Calculated by the law of conservation of mass;

[0116] in, is the mineral volume fraction, is the molar volume of the mineral;

[0117] d. According to the law of conservation of mass and energy of equilibrium reaction, the concentration of j particles of the same type is C j The corresponding secondary species i particle concentration C i Quantitative relationship of ′:

[0118] Among them, K i is the reaction equilibrium constant of minor species i, γ j and γ i are the activity coefficients of particle j and its corresponding secondary species particle i, respectively.

[0119] The present invention also proposes a three-dimensional fracture network groundwater reactive solute transport simulation system, which includes the following modules:

[0120] Model parameter collection module, used to collect geological structure data and local hydrogeological parameters of fractured plots in the target study area;

[0121] The 3D fracture network construction module is used to perform modeling and parameterization processing on the collected geological structure data, set the lithologic parameters of the formation, generate a discrete fracture network, and upscale the plane fracture grid into a 3D cubic grid, mapping it into an equivalent continuous porous medium, and then dividing it into a 3D fracture body grid to establish a geological structure model of the study area.

[0122] The groundwater flow module is used to determine the initial head conditions, boundary conditions, and permeability coefficients of fractures and fracture bodies based on the established geological structure model and the collected hydrogeological parameters, and to solve the spatiotemporal distribution of groundwater levels in the study area to obtain the initial groundwater flow field.

[0123] The solute transport module is used to set the physical and chemical parameters of pollutants and the source-sink phase conditions based on the initial groundwater flow field obtained, and to determine the spatiotemporal distribution of solute concentration in fractured groundwater;

[0124] The fracture reaction transport module is used to set relevant chemical reaction parameters based on the step-by-step solute transport model. According to the basic equations of conservation of mass and energy, the seepage field equations, and the solute transport control equations, the flow field-chemical field coupling is used to solve the results of reactive solute transport in fractured media driven by hydrodynamic and hydrogeochemical reactions, and complete the numerical simulation of reactive solute transport in groundwater in a three-dimensional fracture network.

[0125] Furthermore, the data include the following for the target study area: (1) hydrogeological conditions, stratum lithology and thickness, burial depth, thickness and spatial distribution of aquifers and aquicludes, recharge and drainage conditions of hydrogeological units, groundwater ion content, groundwater flow rate, and hydrogeochemical characteristics; (2) fracture medium conditions, the overall number and spatial distribution of fractures, fracture occurrence, opening, development depth, roughness, physical and chemical properties of fracture fillings and fracture geological bodies, and spatial distribution of fracture network permeability coefficient.

[0126] Furthermore, the three-dimensional fracture network model is constructed as follows:

[0127] The characteristic parameters of the three-dimensional fracture network geological structure model are determined by using the discrete fracture network parameter control equations and the model upscaling grid subdivision control equations;

[0128] (1) Discrete fracture network parameter control equations

[0129] a. The direction of the normal vector of the crack surface is characterized by the von Mises-Fisher distribution function;

[0130] in, is the crack average normal vector matrix; κ is the density coefficient; T represents the matrix transpose; I is the identity matrix;

[0131] b. The crack radius R adopts a truncated power law distribution, defined by the following probability density function;

[0132] Among them, R u and R l are the upper and lower bounds of the crack radius, respectively, and α is a constant coefficient;

[0133] c. The crack opening φ is determined by the exponential function of the crack radius:

[0134] φ=F·R γ (3)

[0135] Among them, F and γ are parameter values ​​related to water permeability;

[0136] d. Fracture permeability k is calculated based on fracture opening φ;

[0137] (2) Model upscaling and meshing control equations

[0138] a. The porosity n of the fractured medium is calculated based on the opening of the fractures in each grid;

[0139] Where l is the length of the mesh edge, N is the total number of cracks passing through the mesh, and φ i represents the opening of the i-th crack in the target grid, which is calculated by formula (3);

[0140] b. Fracture medium permeability k f Assign values ​​based on the permeability of each fracture in the reference grid;

[0141] Among them, k i The permeability of the i-th crack in the target grid is calculated by formula (4).

[0142] Furthermore, the pollutant parameters are set to include: groundwater background ion components and concentrations, pollution source ion components and concentrations, and source-sink boundary conditions; in step SS5, the chemical reaction parameters include: groundwater ion component reaction equations and reaction equilibrium constants, mineral components and their volume fractions, reaction rate constants, and mineral phase reaction specific surface areas.

[0143] Furthermore, the groundwater reactive solute transport model is constructed as follows:

[0144] The seepage field equation, solute transport control equation, mass and energy conservation equation are used to calculate the spatiotemporal evolution of target pollutants in the model driven by hydrodynamic and hydrogeochemical effects.

[0145] a. The seepage field equation is:

[0146] Where n is porosity; s is saturation; ρ is fluid density; w is infiltration recharge; k is absolute permeability; k r is the relative permeability; μ is the fluid viscosity coefficient; P is the pressure; g is the acceleration of gravity; z is the position head;

[0147] b. The solute transport control equation is derived based on the law of conservation of mass and is expressed as:

[0148] Among them, C j is the concentration of j particles; C i ′ is the concentration of the secondary species i particle corresponding to the j particle; υ ji is the stoichiometric number of all liquid phase reactions involving particles j and i; N i is the total number of particles of minor species i; q is the groundwater velocity; D is the molecular diffusion tensor; Q j is the source and sink term of particle j; I m is the mineral reaction rate; jm ′ is the stoichiometric number of particle j related to the mineral phase reaction; M is the total number of particles j involved in the mineral phase reaction; particle j represents the particle with the largest relative concentration in the liquid phase, and the corresponding minor species i particle represents the particle j ionized in the liquid phase or combined with other group particles, for example, particle j is Ca 2+ When the corresponding secondary species i particles include CaCO3(aq), CaHCO3 + 、Ca(OH) + 、CaCl + , CaCl2(aq), CaSO4(aq), etc.; j particle is HCO3 - When the corresponding secondary species i particles include CO3 2- , CO2(aq), etc.

[0149] c. Mineral phase reaction rate I m Calculated by the law of conservation of mass;

[0150] in, is the mineral volume fraction, is the molar volume of the mineral;

[0151] d. According to the law of conservation of mass and energy of equilibrium reaction, the concentration of j particles of the same type is C j The corresponding secondary species i particle concentration C i Quantitative relationship of ′:

[0152] Among them, K i is the reaction equilibrium constant of minor species i, γ j and γ i are the activity coefficients of particle j and its corresponding secondary species particle i, respectively.

[0153] Example 1: The embodiment of the present invention takes the established benchmark model as an example.

[0154] 1) Conceptual Model

[0155] The simulation area of ​​this example is a rectangle with a length of 150m and a width of 80m. The thickness of the aquifer is 30m. A three-dimensional discrete fracture network is generated within the simulation area. The fracture statistical characteristics are set to truncated power law distribution. The fracture length is limited to 10-40m. The fracture density parameter P is 0. 32 The value is set to 1.1 (i.e. the sum of all crack areas within a unit volume is 1.1m 2 ), the main fracture group strikes 90° and dips 90°. The discrete fracture network model generates a total of 253 connected effective fractures, which are mapped into an equivalent continuous porous medium model and divided into 75 rows × 40 columns × 15 layers, with a total of 45,000 effective units. The water flow in the simulation area is from west to east. The east and west boundaries are set as constant head boundaries, with a head difference of 0.03m on both sides. The north and south boundaries are set as zero flow boundaries, and the aquifer bottom plate is a water-proof boundary. Rainfall infiltration and water evaporation are not considered in the simulation area. The model rock matrix is ​​a homogeneous isotropic water-bearing medium, and its permeability, porosity, and tortuosity are set to 1.0×10 -20 , 5.0×10 -3 , 0.2, and the longitudinal, transverse and vertical diffusions are 1m, 0.001m and 0.001m respectively; the parameters of the fracture part are set as an equivalent continuous porous medium of heterogeneous anisotropy according to the physical properties such as shape, size, opening and roughness of the discrete fracture network.

[0156] The simulated area was initially balanced by slightly alkaline water (pH = 8) and limestone matrix, and then acidic water (pH = 5) flowed in from the west at a fixed concentration. The main mineral of the entire three-dimensional fracture model rock mass is calcite, and its initial volume fraction is 1.0×10 -5 The initial concentrations of the main ions in the acidic water and model groundwater are shown in Table 1, and the main reaction parameters are shown in Table 2. The calcite dissolution reaction rate constant in the fractured rock mass is 1.0×10 -6 mol / m 2 s, specific surface area is 1.0m 2 / m 3 The total stress period of the simulation is set to 10 years.

[0157] Table 1 Main ion parameters of the model in the embodiment of the invention (unit: mol / L)

[0158] Table 2 Main reaction parameters of the model in the embodiment of the invention

[0159] 2) Determine the control equation

[0160] (1) Discrete fracture network parameter control equations

[0161] a. The direction of the normal vector of the crack surface is characterized by the von Mises-Fisher distribution function:

[0162] in, is the crack average normal vector matrix; κ is the density coefficient; T represents the matrix transpose; I is the identity matrix;

[0163] b. The crack radius R adopts a truncated power law distribution, defined by the following probability density function:

[0164] Among them, R u and R l are the upper and lower bounds of the crack radius respectively; α is a constant coefficient;

[0165] c. The crack opening φ is determined by the exponential function of the crack radius:

[0166] φ=F·R γ (3)

[0167] Among them, F and γ are empirical parameter values ​​related to water permeability;

[0168] d. Fracture permeability k is calculated based on fracture opening φ:

[0169] (2) Model upscaling and meshing control equations:

[0170] a. The porosity n of the fractured medium is calculated based on the opening of the fractures in each grid:

[0171] Where l is the length of the mesh edge; N is the total number of cracks passing through the mesh; φ i represents the opening of the i-th crack in the target grid (calculated by formula (3));

[0172] b. Fracture medium permeability k f Assign values ​​based on the permeability of each fracture in the reference grid:

[0173] Where N is the total number of cracks passing through the mesh; k i represents the permeability of the i-th fracture in the target grid (calculated by formula (4)).

[0174] (3) Governing equations of groundwater reactive solute transport model:

[0175] a. The seepage field equation is:

[0176] Where n is porosity; s is saturation; ρ is fluid density; w is infiltration recharge; k is absolute permeability; k r is the relative permeability; μ is the fluid viscosity coefficient; P is the pressure; g is the acceleration of gravity; z is the position head;

[0177] b. The solute transport control equation is derived based on the law of conservation of mass and is expressed as:

[0178] Among them, C j is the concentration of j particles; C i ′ is the concentration of the secondary species i particle corresponding to the j particle; υ ji is the stoichiometric number of all liquid phase reactions involving particles j and i; N i is the total number of particles of minor species i; q is the groundwater velocity; D is the molecular diffusion tensor; Q j is the source and sink term of particle j; I m is the mineral reaction rate; jm ′ is the stoichiometric number of particle j related to the mineral phase reaction; M is the total number of particles j involved in the mineral phase reaction; particle j represents the particle with the largest relative concentration in the liquid phase, and the corresponding minor species i particle represents the particle j ionized in the liquid phase or combined with other group particles, for example, particle j is Ca 2+ When the corresponding secondary species i particles include CaCO3(aq), CaHCO3 + 、Ca(OH) + 、CaCl + , CaCl2(aq), CaSO4(aq), etc.; j particle is HCO3 - When the corresponding secondary species i particles include CO3 2- , CO2(aq), etc.

[0179] c. Mineral phase reaction rate I m Calculated by the law of conservation of mass;

[0180] in, is the mineral volume fraction, is the molar volume of the mineral;

[0181] d. According to the law of conservation of mass and energy of equilibrium reaction, the concentration of j particles of the same type is C j The corresponding secondary species i particle concentration C i Quantitative relationship of ′:

[0182] Among them, K i is the reaction equilibrium constant of minor species i, γ j and γ iare the activity coefficients of particle j and its corresponding secondary species particle i, respectively.

[0183] 3) Analysis of numerical simulation results of a three-dimensional fracture network groundwater reactive solute transport model

[0184] Figure 6 shows the distribution of groundwater pressure within a fracture network. It can be seen that the groundwater flow field at shallow depths is affected by the distribution of fractures, resulting in significant nonuniformity in the water pressure field in the upstream region. However, this effect decreases with increasing depth, indicating that the spatial distribution of water gradients in medium- to deep-seated fractures is similar to that of groundwater in saturated porous media.

[0185] Figure 7 shows the spatiotemporal distribution of pH values ​​within the fracture network during the first eight years of a simulated scenario of sustained leakage of acidic water at a constant concentration. Comparing the pH distribution of the fractured groundwater in the second year of the simulation with that in the fourth, sixth, and eighth years reveals that the migration and diffusion efficiency of acidic water is greatest when the acidic water is introduced into the fractured groundwater for a short time. However, after the second year, as the distance traveled by the acidic water increases, its migration rate gradually decreases. This is likely due to the fact that as the acidic water migrates through the fractured medium, its contact area with the limestone in the fractured groundwater increases, increasing the buffering effect of the fractured rock mass on the acidic water and effectively hindering its migration. Furthermore, the figure shows a significant preferential flow of acidic water in areas with dense longitudinal intersecting fractures. This is due to the clustering of water-conducting fractures along the descending hydraulic gradient, which increases the permeability of these areas. These patterns are consistent with real-world phenomena, demonstrating the robustness of the model for characterizing reactive solute transport in three-dimensional fractured groundwater.

[0186] Figure 8 shows the spatial distribution of the concentrations of the main groundwater components in the three-dimensional fracture network in the 10th year of simulation. As can be seen from the figure, the maximum distance of the acidic water impact is 92m, mainly concentrated in the fracture medium. In addition, as the acidic water continues to migrate downstream, calcite is also continuously buffered and dissolved, resulting in a large range of Ca behind the acid front. 2+ plume; the buffering effect also involves the original high concentration of HCO3 - , HCO3 in the acidic region - While being reacted, it is also reacted with low concentration HCO3 in acidic water. - displacement, resulting in the overall HCO3 - The free ion concentration decreases. By comparing Figure 7, we can see that the solute migration of free ions in the fractured groundwater in the simulation area is mainly located in the fractured medium area with greater water conductivity; while the chemical composition of the rock matrix remains basically unchanged, with only slight changes at the fractured medium interface. However, this has a significant impact on the hydrogeochemical reactions in the entire fractured groundwater, which is also an important link that cannot be ignored in guiding the remediation and treatment of fractured groundwater.

[0187] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the technical principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A method for simulating reactive solute transport in groundwater in a three-dimensional fracture network, characterized in that: The method comprises the following steps: Step SS1, collecting geological structure data and local hydrogeological parameters of the fractured plots in the target study area; Step SS2, based on the geological structure data collected in step SS1, the data is modeled and parameterized, the lithological parameters of the formation are set, a discrete fracture network is generated, and the plane fracture grid is scaled to a three-dimensional cubic grid and mapped into an equivalent continuous porous medium, and divided into a three-dimensional fracture body grid to establish a geological structure model of the study area; Step SS3, based on the geological structure model constructed in step SS2, the initial water head condition, boundary condition, fracture and fracture body permeability coefficient of the model are determined according to the collected hydrogeological parameters, the spatiotemporal distribution of the groundwater level in the study area is solved, and the initial groundwater flow field is obtained; Step SS4, based on the initial groundwater flow field obtained in step SS3, setting the pollutant physicochemical parameters and source-sink phase conditions, and determining the spatiotemporal distribution of solute concentration in fractured groundwater; Step SS5, based on the solute transport model of step SS4, set the relevant chemical reaction parameters, and solve the results of reactive solute transport in the fracture medium driven by hydrodynamic and hydrogeochemical reactions according to the basic equations of conservation of mass and energy, the seepage field equations, and the solute transport control equations to couple the water flow field and the chemical field, and complete the numerical simulation of reactive solute transport in groundwater in a three-dimensional fracture network.

2. A method for simulating reactive solute transport in groundwater in a three-dimensional fracture network according to claim 1, characterized in that: In step SS1, the data include the following: (1) hydrogeological conditions, stratigraphic lithology and thickness, burial depth, thickness and spatial distribution of aquifers and aquicludes, recharge and drainage conditions of hydrogeological units, groundwater ion content, groundwater flow rate, and hydrogeochemical characteristics of the target study area; (2) fracture medium conditions, the overall number and spatial distribution of fractures, the occurrence, opening, development depth, roughness, physical and chemical properties of fracture fillings and fracture geological bodies, and spatial distribution of fracture network permeability coefficients.

3. A method for simulating reactive solute transport in groundwater in a three-dimensional fracture network according to claim 1, characterized in that: In step SS2, the three-dimensional fracture network model is constructed as follows: The characteristic parameters of the three-dimensional fracture network geological structure model are determined by using the discrete fracture network parameter control equation and the model upscaling grid generation control equation. (1) Discrete fracture network parameter control equation a. The direction of the normal vector of the crack surface is characterized by the von Mises-Fisher distribution function; in, is the crack average normal vector matrix; κ is the density coefficient; T represents the matrix transpose; I is the unit matrix; b. The crack radius R adopts a truncated power law distribution, defined by the following probability density function; Among them, R u and R l are the upper and lower bounds of the crack radius, respectively, and α is a constant coefficient; c. The crack opening φ is determined by the exponential function of the crack radius; φ=F·R γ (3) Among them, F and γ are parameter values ​​related to water permeability; d. Fracture permeability k is calculated based on fracture opening φ; (2) Model upscaling and meshing control equations a. The porosity n of the fractured medium is calculated based on the opening of the fractures in each grid; Where l is the length of the mesh edge, N is the total number of cracks passing through the mesh, and φ i represents the opening of the i-th crack in the target grid, which is calculated by formula (3); b. Fracture medium permeability k f The permeability of each fracture in the reference grid is assigned; Among them, k i The permeability of the i-th fracture in the target grid is calculated by formula (4).

4. A method for simulating reactive solute transport in groundwater in a three-dimensional fracture network according to claim 1, characterized in that: In step SS4, the pollutant parameters are set including: groundwater background ion components and concentrations, pollution source ion components and concentrations, and source-sink boundary conditions; in step SS5, the chemical reaction parameters include: groundwater ion component reaction equation and its reaction equilibrium constant, mineral components and their volume fractions, reaction rate constants, and mineral phase reaction specific surface area.

5. A method for simulating reactive solute transport in groundwater in a three-dimensional fracture network according to claim 1, characterized in that: In step SS5, the groundwater reactive solute transport model is constructed as follows: The seepage field equation, solute transport control equation, mass and energy conservation equation are used to calculate the temporal and spatial evolution of the target pollutants in the model driven by hydrodynamic and hydrogeochemical effects. a. The seepage field equation is: Where n is porosity; s is saturation; ρ is fluid density; w is infiltration recharge; k is absolute permeability; k r is the relative permeability; μ is the fluid viscosity coefficient; P is the pressure; g is the gravitational acceleration; z is the position head; b. The solute transport control equation is derived based on the law of conservation of mass and is expressed as: Among them, C j is the j particle concentration; C i ′ is the concentration of the secondary species i particle corresponding to the j particle; υ ji is the stoichiometric number of all j and i particles related to the liquid phase reaction; N i is the total number of particles of minor species i; q is the groundwater velocity; D is the molecular diffusion tensor; Q j is the source and sink term of particle j; I m is the mineral reaction rate; jm ′ is the stoichiometric number of the j particle related to the mineral phase reaction; M is the total number of j particles involved in the mineral phase reaction; c. Mineral phase reaction rate I m Calculated by the law of conservation of mass; in, is the mineral volume fraction, is the molar volume of the mineral; d. According to the law of conservation of mass and energy of equilibrium reaction, the concentration of j particles of the same type is C j The corresponding particle concentration C of the secondary species i i Quantitative relationship of ′: Among them, K i is the reaction equilibrium constant of minor species i, γ j and γ i are the activity coefficients of particle j and its corresponding secondary species particle i, respectively.

6. A three-dimensional fracture network groundwater reactive solute transport simulation system, characterized in that: The system includes the following modules: Model parameter collection module, used to collect geological structure data and local hydrogeological parameters of fractured plots in the target study area; The three-dimensional fracture network construction module is used to perform modeling and parameterization processing on the collected geological structure data, set the lithological parameters of the formation, generate a discrete fracture network, and upscale the plane fracture grid to a three-dimensional cubic grid to map it into an equivalent continuous porous medium, and divide it into a three-dimensional fracture body grid to establish a geological structure model of the study area; The groundwater flow module is used to determine the initial water head conditions, boundary conditions, and permeability coefficients of fractures and fracture bodies of the model based on the established geological structure model and the collected hydrogeological parameters, and to solve the spatiotemporal distribution of groundwater levels in the study area to obtain the initial groundwater flow field; The solute transport module is used to set the physical and chemical parameters of pollutants and the source-sink phase conditions based on the initial groundwater flow field obtained, and determine the spatiotemporal distribution of solute concentration in fractured groundwater; The fracture reaction transport module is used to set relevant chemical reaction parameters based on the step solute transport model, and solve the results of reactive solute transport in fracture media driven by hydrodynamic and hydrogeochemical reactions according to the basic equations of mass and energy conservation, seepage field equations, and solute transport control equations to couple the water flow field and chemical field, and complete the numerical simulation of reactive solute transport in groundwater in a three-dimensional fracture network.

7. A three-dimensional fracture network groundwater reactive solute transport simulation system according to claim 6, characterized in that: The data include: (1) hydrogeological conditions, stratum lithology and thickness, burial depth, thickness and spatial distribution of aquifers and impermeable layers, recharge and drainage conditions of hydrogeological units, groundwater ion content, groundwater flow rate, and hydrogeochemical characteristics of the target study area; (2) fracture medium conditions, the overall number and spatial distribution of fractures, the occurrence, opening, development depth, roughness, physical and chemical properties of fracture fillings and fracture geological bodies, and spatial distribution of fracture network permeability.

8. A three-dimensional fracture network groundwater reactive solute transport simulation system according to claim 6, characterized in that: The three-dimensional fracture network model is constructed as follows: The characteristic parameters of the three-dimensional fracture network geological structure model are determined by using the discrete fracture network parameter control equation and the model upscaling grid generation control equation. (1) Discrete fracture network parameter control equation a. The direction of the normal vector of the crack surface is characterized by the von Mises-Fisher distribution function; in, is the crack average normal vector matrix; κ is the density coefficient; T represents the matrix transpose; I is the unit matrix; b. The crack radius R adopts a truncated power law distribution, defined by the following probability density function; Among them, R u and R l are the upper and lower bounds of the crack radius, respectively, and α is a constant coefficient; c. The crack opening φ is determined by the exponential function of the crack radius: φ=F·R γ (3) Among them, F and γ are parameter values ​​related to water permeability; d. Fracture permeability k is calculated based on fracture opening φ; (2) Model upscaling and meshing control equations a. The porosity n of the fractured medium is calculated based on the opening of the fractures in each grid; Where l is the length of the mesh edge, N is the total number of cracks passing through the mesh, and φ i represents the opening of the i-th crack in the target grid, which is calculated by formula (3); b. Fracture medium permeability k f The permeability of each fracture in the reference grid is assigned; Among them, k i The permeability of the i-th fracture in the target grid is calculated by formula (4).

9. A three-dimensional fracture network groundwater reactive solute transport simulation system according to claim 6, characterized in that: The pollutant parameters to be set include: groundwater background ion components and concentrations, pollution source ion components and concentrations, and source-sink boundary conditions; in step SS5, the chemical reaction parameters include: groundwater ion component reaction equations and reaction equilibrium constants, mineral components and their volume fractions, reaction rate constants, and mineral phase reaction specific surface areas.

10. A three-dimensional fracture network groundwater reactive solute transport simulation system according to claim 6, characterized in that: The groundwater reactive solute transport model is constructed as follows: The seepage field equation, solute transport control equation, mass and energy conservation equation are used to calculate the temporal and spatial evolution of the target pollutants in the model driven by hydrodynamic and hydrogeochemical effects. a. The seepage field equation is: Where n is porosity; s is saturation; ρ is fluid density; w is infiltration recharge; k is absolute permeability; k r is the relative permeability; μ is the fluid viscosity coefficient; P is the pressure; g is the gravitational acceleration; z is the position head; b. The solute transport control equation is derived based on the law of conservation of mass and is expressed as: Among them, C j is the j particle concentration; C i ′ is the concentration of the secondary species i particle corresponding to the j particle; υ ji is the stoichiometric number of all j and i particles related to the liquid phase reaction; N i is the total number of particles of minor species i; q is the groundwater velocity; D is the molecular diffusion tensor; Q j is the source and sink term of particle j; I m is the mineral reaction rate; jm ′ is the stoichiometric number of the j particle related to the mineral phase reaction; M is the total number of j particles involved in the mineral phase reaction; c. Mineral phase reaction rate I m Calculated by the law of conservation of mass; in, is the mineral volume fraction, is the molar volume of the mineral; d. According to the law of conservation of mass and energy of equilibrium reaction, the concentration of j particles of the same type is C j The corresponding particle concentration C of the secondary species i i Quantitative relationship of ′: Among them, K i is the reaction equilibrium constant of minor species i, γ j and γ i are the activity coefficients of particle j and its corresponding secondary species particle i, respectively.

Images