Method and system for simulating acidic water reaction migration under fracture-matrix upscaling model

Abstract

The invention discloses an acidic water reaction migration simulation method and system under a fracture-matrix upscaling model, and the model can carry out fracture parameter upscaling based on a locally encrypted unstructured grid, thereby constructing the fracture-matrix upscaling model. According to the method, the overall grid node operation load is reduced while the model simulation precision is improved, on the basis of an upscaling model, the hydrodynamic force-hydrochemical interaction of fractures and matrixes is considered, an acidic water reaction migration model in a fracture network is constructed, and the reaction migration process of acidic water characteristic components in the three-dimensional fracture network is researched. The method provided by the invention comprises the following steps: constructing a fracture-matrix parameter upscaling model with locally encrypted grids, considering the hydrodynamic-hydrochemical interaction of fractures and matrixes, constructing a reactive solute transport model based on the upscaling model, and carrying out acidic water reaction transport numerical simulation in a fracture-matrix system.

Description

Technical Field

[0001] The present invention relates to a method and system for simulating the reactive transport of acidic water under a fracture-matrix upscaling model, belonging to the field of earth science and engineering. Background Art

[0002] Acid mine drainage pollution is one of the main environmental problems in high-sulfur mining areas. Its pollutants continuously spread through the groundwater system, seriously threatening the regional water resource security and ecological balance. Existing studies have shown that constructing an accurate groundwater reactive transport model for fractured media is a fundamental technical bottleneck for pollution prevention and control. The core challenge stems from the dual complexity of the fracture system: on the one hand, fractured rock masses exhibit significant spatial heterogeneity and directional characteristics, making it difficult to effectively identify non-uniform flow fields and solute transport paths through conventional monitoring means; on the other hand, the hydrodynamics-hydrochemistry coupling mechanism involved in the migration of acidic water has not been systematically characterized in existing models.

[0003] The current mainstream simulation methods for fractured media mainly have the following technical defects: (1) Although the equivalent porous media model can reduce the modeling complexity through parameter equivalence processing, its continuum hypothesis is fundamentally in conflict with the discrete characteristics of fractured media. Especially when using a coarse grid to simplify the fracture structure, it will lead to the loss of geometric characteristics and connectivity information of key hydraulic channels; while using a fine grid to accurately depict fractures, it faces the technical bottleneck of exponential growth in the computational scale. (2) Although the discrete fracture network model can explicitly represent the geometric parameters and hydraulic characteristics of individual fractures, there is generally a problem of missing three-dimensional modeling dimensions in practical applications. Existing two-dimensional simplified models cannot accurately reflect the spatial topological structure of the fracture network and generally ignore two key elements: one is the water volume exchange and solute diffusion process between fractures and the matrix, and the other is the hydrogeochemical interaction between acidic components and surrounding rock minerals. This technical defect directly restricts the quantitative assessment of polluted sites and the optimization of treatment plans. Summary of the Invention

[0004] Object of the Invention: Aiming at the deficiencies of the above-mentioned existing technologies, the object of the present invention is to provide a method and system for simulating the reactive transport of acidic water under a fracture-matrix upscaling model. This method can perform fracture parameter upscaling based on locally refined unstructured grids, thereby constructing a fracture-matrix upscaling model, which improves the model simulation accuracy while reducing the overall grid node operation load. On the basis of the upscaling model, considering the hydrodynamics-hydrochemistry interaction between fractures and the matrix, a reactive transport model of acidic water in the fracture network is constructed to carry out numerical simulation of the reactive transport of acidic water in the fracture-matrix system.

[0005] Technical Solution: To solve the above technical problems, the present invention proposes a method for simulating the reactive transport of acidic water under a fracture-matrix upscaling model. The method includes the following steps:

[0006] Step SS1, determine the basic parameters required for the fracture-matrix upscaling model for local refined grid generation. The basic model parameters include: fracture distribution empirical function, number of fractures, fracture radius, fracture direction, fracture width, matrix permeability, matrix porosity, local refined grid encryption level, and basic grid cell length;

[0007] Step SS2, perform local refined grid meshing according to the grid setting parameters in Step SS1, and upscale the parameters of the corresponding fracture-matrix part of the grid to establish a fracture-matrix upscaled geological structure model;

[0008] Step SS3, based on the fracture-matrix upscaled geological structure model established in Step SS2, determine the groundwater parameters required for the model, including: initial head distribution, boundary conditions, and permeability after upscaling of the fracture-matrix parameters in Step SS2, solve for the fracture groundwater level and flow velocity distribution in the simulation area, and obtain the initial groundwater flow field;

[0009] Step SS4, based on the initial groundwater flow field obtained in Step SS3, determine the physicochemical properties, background concentration, and infiltration concentration of the acidic water pollution components to calculate the spatial distribution of solute components in the fracture-matrix of the simulation area;

[0010] Step SS5, based on the solute migration process obtained in Step SS4, determine the hydrogeochemical reaction parameters involved in the acidic water pollution migration process, including: chemical reaction formula and reaction rate, mineral volume fraction and specific surface area, couple the groundwater flow field - solute field - chemical field, establish an upscaled reactive solute transport model for acidic water in the fracture-matrix, and simulate the acidic water pollution migration process under the combined influence of fracture-matrix spatial heterogeneity and water-rock interaction.

[0011] Further, in Step SS2, the fracture-matrix upscaled geological structure model is constructed as follows:

[0012] a. After the grid scheme is determined, perform upscaling calculations for porosity and permeability tensors. The porosity φ of the fracture-matrix upscaling model is calculated as follows:

[0013]

[0014] In the formula, when the grid cell is the matrix, φ = φ M , when the grid cell is the fracture, φ = φ F ; V F is the total fracture volume intersecting within the model cell; V C is the model cell volume; N is the total number of fractures intersecting within the model cell; A f is the fracture area within the cell; b fis the fracture width;

[0015] b. The permeability tensor [k] of the fracture-matrix model upscaling model is expressed as:

[0016]

[0017]

[0018] In the formula, when the grid cell is the matrix, [k] = [k] M , where k xx , k yy , k zz are the permeability components of the matrix in the x, y, and z directions. When the grid cell is a fracture, [k] = [k] F ; φ f is the porosity of fracture f in the simulation unit; b f is the width of fracture f; is the coordinate transformation tensor; n1, n2, n3 are the orthogonal vectors of fracture f.

[0019] Furthermore, in step SS5, the acid water reactive solute transport model of the fracture-matrix upscaling is constructed as follows:

[0020] a. Groundwater flow equation

[0021]

[0022] In the formula, φ is the porosity; [k] is the permeability tensor; μ is the fluid viscosity coefficient; P is the pressure; g is the acceleration of gravity; ρ is the fluid density; w is the source-sink term; z is the position head;

[0023] b. Groundwater solute transport equation:

[0024]

[0025] In the formula, C j is the total concentration of the jth solute; V is the groundwater flow velocity; D is the hydrodynamic dispersion coefficient tensor; w j is the source-sink term of the jth solute; M is the total number of mineral phase reactions involved in the jth solute; ξ m is the stoichiometric coefficient of the jth solute in the reaction with minerals; R jm is the reaction rate of the jth solute participating in the mineral reaction;

[0026] c. Groundwater reaction system control equation:

[0027]

[0028] In the formula, C j is the total concentration of the jth solute; Xj is the free ion concentration of the j-th solute; N is the total number of the j-th solute participating in the liquid-phase reaction; ξ i is the stoichiometric number of the j-th solute in the i-th liquid-phase reaction; K i is the equilibrium constant of the i-th liquid-phase reaction; a i is the activity coefficient of the product of the i-th liquid-phase reaction; a j is the activity coefficient of the j-th solute in the i-th liquid-phase reaction;

[0029]

[0030] In the formula, R m is the mineral reaction rate; A m is the specific surface area of the mineral; k m (T, λ) is the effective reaction rate, T is the temperature influence parameter, and λ is the ion influence pre-parameter; K eq is the reaction equilibrium constant; Q is the ion activity product; α is the mineral saturation sensitivity parameter; β is the reaction affinity parameter; ρ m is the mineral density; M m is the molar mass of the mineral; η is the mineral volume fraction.

[0031] The present invention also provides an acidic water reaction and transport simulation system under a fracture-matrix upscaling model, and the system includes the following modules:

[0032] A model parameter input module, which is used to determine the basic parameters required for the fracture-matrix upscaling model for generating a locally refined grid. The basic model parameters include: a fracture distribution empirical function, the number of fractures, the fracture radius, the fracture direction, the fracture width, the matrix permeability, the matrix porosity, the locally refined grid encryption level, and the basic grid cell length;

[0033] A fracture-matrix upscaling module, which is used to perform local refined grid meshing according to the grid setting parameters in the model parameter input module, and upscale the parameters of the corresponding fracture-matrix part of the grid to establish a fracture-matrix upscaling geological structure model;

[0034] A groundwater dynamics calculation module, which is used to determine the groundwater parameters required for the model on the basis of the fracture-matrix upscaling geological structure model, including: the initial head distribution, the boundary conditions, and the permeability after upscaling the fracture-matrix parameters in step SS2, and solve the fracture groundwater level and flow velocity distribution in the simulation area to obtain the initial groundwater flow field;

[0035] A solute transport calculation module, which is used to determine the physicochemical properties, the background concentration, and the infiltration concentration of the acidic water pollution components on the basis of the obtained initial groundwater flow field, so as to calculate the spatial distribution of the solute components in the fracture-matrix of the simulation area;

[0036] The reactive solute transport module is used to determine the hydrogeochemical reaction parameters involved in the acidic water pollution migration process based on the obtained solute migration process, including: chemical reaction equations and reaction rates, mineral volume fractions and specific surface areas. It couples the groundwater flow field - solute field - chemical field, establishes a reactive solute transport model for acidic water with fracture - matrix upscaling, and simulates the acidic water pollution migration process under the combined influence of fracture - matrix spatial heterogeneity and water - rock interaction.

[0037] Furthermore, the fracture - matrix upscaled geological structure model is constructed as follows:

[0038] a. After determining the grid scheme, perform upscaling calculations of porosity and permeability tensors. The porosity φ of the fracture - matrix upscaled model is calculated as follows:

[0039]

[0040] In the formula, when the grid cell is the matrix, φ = φ M , and when the grid cell is the fracture, φ = φ F ; V F is the total fracture volume intersecting within the model cell; V C is the model cell volume; N is the total number of fractures intersecting within the model cell; A f is the fracture area within the cell; b f is the fracture width;

[0041] b. The permeability tensor [k] of the fracture - matrix model upscaled model is expressed as:

[0042]

[0043] In the formula, when the grid cell is the matrix, [k] = [k] M , where k xx , k yy , k zz are the permeability components of the matrix in the x, y, and z directions. When the grid cell is the fracture, [k] = [k] F ; φ f is the porosity of fracture f within the simulation cell; b f is the width of fracture f; is the coordinate transformation tensor; n1, n2, n3 are the orthogonal vectors of fracture f.

[0044] Furthermore, the reactive solute transport model for acidic water with fracture - matrix upscaling is constructed as follows:

[0045] a. Groundwater flow equation

[0046]

[0047] Wherein, φ is porosity; [k] is the permeability tensor; μ is the fluid viscosity coefficient; P is the pressure; g is the acceleration of gravity; ρ is the fluid density; w is the source-sink term; z is the position head;

[0048] b. Groundwater solute transport equation:

[0049]

[0050] Wherein, C j is the total concentration of the j-th solute; V is the groundwater flow velocity; D is the hydrodynamic dispersion coefficient tensor; w j is the source-sink term of the j-th solute; M is the total number of mineral phase reactions involved in the j-th solute; ξ m is the stoichiometric coefficient of the j-th solute in the reaction with the mineral; R jm is the reaction rate of the j-th solute participating in the mineral reaction;

[0051] c. Groundwater reaction system control equation:

[0052]

[0053] Wherein, C j is the total concentration of the j-th solute; X j is the free ion concentration of the j-th solute; N is the total number of liquid phase reactions participated in by the j-th solute; ξ i is the stoichiometric coefficient of the j-th solute in the i-th liquid phase reaction; K i is the equilibrium constant of the i-th liquid phase reaction; a i is the activity coefficient of the product of the i-th liquid phase reaction; a j is the activity coefficient of the j-th solute in the i-th liquid phase reaction;

[0054]

[0055] Wherein, R m is the mineral reaction rate; A m is the specific surface area of the mineral; k m (T, λ) is the effective reaction rate, T is the temperature influence parameter, λ is the ion influence pre-parameter; K eq is the reaction equilibrium constant; Q is the ion activity product; α is the mineral saturation sensitivity parameter; β is the reaction affinity parameter; ρ m is the mineral density; M m is the molar mass of the mineral; η is the mineral volume fraction.

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

[0057] In the fracture-matrix upscaling acidic water reactive solute transport model of the present invention, due to the local refinement of the grid for the fracture-matrix upscaling model, it is possible to improve the simulation accuracy of the model while reducing the overall grid node operation load, further quantitatively calculate the excretion intensity of characteristic components in the acidic water pollution site of the fracture medium, provide a quantitative tool for the migration and evolution process of acidic water under the hydrodynamic-hydrochemical interaction between fractures and the matrix, and provide scientific and reasonable decision-making support for the prevention and control of acidic water pollution in the fracture medium.

[0058] The simulation technology provided by the present invention can quantitatively calculate the excretion intensity of characteristic components in the acidic water pollution site of the fracture medium, accurately depict the migration and evolution process of acidic water under the hydrodynamic-hydrochemical interaction between fractures and the matrix, and provide scientific and reasonable decision-making support for the prevention and control of acidic water pollution in the fracture medium.

[0059] In view of the above defects, the present invention proposes an acidic water reactive transport simulation method and system under a fracture-matrix upscaling model. This technology can generate a fracture-matrix upscaling model with locally refined grids, reduce the overall grid node operation load while improving the simulation accuracy of the model, consider the influence of the three-dimensional fracture space structure and the hydrodynamic-hydrochemical interaction between fractures and the matrix on the reactive transport of groundwater, quantitatively calculate the excretion intensity of characteristic components in the acidic water pollution site of the fracture-matrix system, accurately depict the migration and evolution process of acidic water, and provide scientific and reasonable decision-making support for the prevention and control of acidic water pollution in the fracture medium. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] Figure 1 is a flow chart of the acidic water reactive transport simulation method under the fracture-matrix upscaling model;

[0061] Figure 2 is a conceptual model diagram of the reactive transport of acidic water in the fracture-matrix body;

[0062] Figure 3 is a diagram of the locally refined fracture-matrix upscaling model;

[0063] Figure 4 is for the simulation of the spatial distribution and typical cross-section distribution of H + concentration in the fracture-matrix body in the 10th year;

[0064] Figure 5 is for the spatial distribution of Ca 2+ , tracer concentration and the concentration difference between the two at the typical cross-section of the fracture-matrix at different simulation times;

[0065] Figure 6 is for the average concentration curve diagrams of H + and Ca 2+ and their corresponding tracers at the downstream boundary. Detailed implementation manners

[0066] The present invention will be further described below with reference to the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solutions of the present invention and cannot be used to limit the protection scope of the present invention.

[0067] As Figure 1 shown, the present invention provides a simulation method for the reactive transport of acidic water under a fracture-matrix upscaling model. The method includes the following steps:

[0068] Step SS1: Determine the basic parameters required for the fracture-matrix upscaling model for generating a locally refined grid. The basic parameters of the model include: the empirical function of fracture distribution, the number of fractures, the radius of fractures, the direction of fractures, the width of fractures, the matrix permeability, the matrix porosity, the encryption level of the locally refined grid, and the length of the basic grid unit.

[0069] Step SS2: Perform a local refined grid meshing according to the grid setting parameters in Step SS1, and upscale the parameters of the corresponding fracture-matrix part of the grid to establish a fracture-matrix upscaled geological structure model.

[0070] Step SS3: On the basis of the fracture-matrix upscaled geological structure model established in Step SS2, determine the groundwater parameters required for the model, including: the initial head distribution, the boundary conditions, and the permeability after upscaling the fracture-matrix parameters in Step SS2, solve for the fracture groundwater level and flow velocity distribution in the simulation area to obtain the initial groundwater flow field.

[0071] Step SS4: On the basis of the initial groundwater flow field obtained in Step SS3, determine the physicochemical properties, background concentration, and infiltration concentration of the acidic water pollution components to calculate the spatial distribution of solute components in the fracture-matrix of the simulation area.

[0072] Step SS5: On the basis of the solute migration process obtained in Step SS4, determine the hydrogeochemical reaction parameters involved in the acidic water pollution migration process, including: the chemical reaction formula and reaction rate, the mineral volume fraction and specific surface area, couple the groundwater flow field - solute field - chemical field, establish an acidic water reactive solute transport model for the fracture-matrix upscaling, and simulate the acidic water pollution migration process under the combined influence of the fracture-matrix spatial heterogeneity and water-rock interaction.

[0073] In Step SS2, the fracture-matrix upscaled geological structure model is constructed as follows:

[0074] a. After the grid scheme is determined, perform the upscaling calculation of porosity and permeability tensors. The porosity φ of the fracture-matrix upscaling model is calculated as follows:

[0075]

[0076] In the formula, when the grid cell is the matrix, φ = φ M , and when the grid cell is a fracture, φ = φ F ; V F is the total fracture volume intersecting within the model unit; V C is the model unit volume; N is the total number of fractures intersecting within the model unit; A f is the fracture area within the unit; b f is the fracture width;

[0077] b. The permeability tensor [k] of the fracture-matrix model upscaling model is expressed as:

[0078]

[0079] In the formula, when the grid cell is the matrix, [k] = [k] M , where k xx , k yy , k zz are the permeability components of the matrix in the x, y, and z directions. When the grid cell is a fracture, [k] = [k] F ; φ f is the porosity of fracture f within the simulation unit; b f is the width of fracture f; is the coordinate transformation tensor; n1, n2, n3 are the orthogonal vectors of fracture f.

[0080] In step SS5, the acid water reactive solute transport model for the fracture-matrix upscaling is constructed as follows:

[0081] a. Groundwater flow equation

[0082]

[0083] In the formula, φ is the porosity; [k] is the permeability tensor; μ is the fluid viscosity; P is the pressure; g is the acceleration due to gravity; ρ is the fluid density; w is the source-sink term; z is the position head;

[0084] b. Groundwater solute transport equation:

[0085]

[0086] In the formula, C j is the total concentration of the jth solute; V is the groundwater flow velocity; D is the hydrodynamic dispersion coefficient tensor; w j is the source-sink term of the jth solute; M is the total number of mineral phase reactions involved in the jth solute; ξ m is the stoichiometric coefficient of the jth solute in the reaction with the mineral; R jm is the reaction rate of the jth solute participating in the mineral reaction;

[0087] c. Governing equations of groundwater reaction system:

[0088]

[0089] In the formula, C j is the total concentration of the j-th solute; X j is the free ion concentration of the j-th solute; N is the total number of the j-th solute participating in the liquid-phase reaction; ξ i is the stoichiometric number of the j-th solute in the i-th liquid-phase reaction; K i is the equilibrium constant of the i-th liquid-phase reaction; a i is the activity coefficient of the product of the i-th liquid-phase reaction; a j is the activity coefficient of the j-th solute in the i-th liquid-phase reaction;

[0090]

[0091] In the formula, R m is the mineral reaction rate; A m is the specific surface area of the mineral; k m (T, λ) is the effective reaction rate, T is the temperature influence parameter, and λ is the ion influence pre-parameter; K eq is the reaction equilibrium constant; Q is the ion activity product; α is the mineral saturation sensitivity parameter; β is the reaction affinity parameter; ρ m is the mineral density; M m is the molar mass of the mineral; η is the mineral volume fraction.

[0092] The present invention also provides an acidic water reaction and transport simulation system under a fracture-matrix upscaling model, and the system includes the following modules:

[0093] Model parameter input module, which is used to determine the basic parameters required for the fracture-matrix upscaling model of local refined grid generation. The basic model parameters include: fracture distribution empirical function, number of fractures, fracture radius, fracture direction, fracture width, matrix permeability, matrix porosity, local refined grid encryption level, basic grid unit length;

[0094] Fracture-matrix upscaling module, which is used to perform local refined grid meshing according to the grid setting parameters in the model parameter input module, and upscale the parameters of the corresponding fracture-matrix part of the grid to establish a fracture-matrix upscaling geological structure model;

[0095] The groundwater dynamics calculation module is used to determine the groundwater parameters required for the model based on the fracture-matrix upscaled geological structure model, including: the initial head distribution, boundary conditions, and the permeability after upscaling of the fracture-matrix parameters in step SS2, solve for the fracture groundwater level and velocity distribution in the simulation area, and obtain the initial groundwater flow field;

[0096] The solute transport calculation module is used to determine the physicochemical properties, background concentration, and infiltration concentration of the acidic water pollution components based on the obtained initial groundwater flow field, and calculate the spatial distribution of solute components in the fracture-matrix in the simulation area;

[0097] The reactive solute transport module is used to determine the hydrogeochemical reaction parameters involved in the acidic water pollution migration process based on the obtained solute transport process, including: chemical reaction equations and reaction rates, mineral volume fractions and specific surface areas, couple the groundwater flow field - solute field - chemical field, establish an acidic water reactive solute transport model for fracture-matrix upscaling, and simulate the acidic water pollution migration process under the combined influence of fracture-matrix spatial heterogeneity and water-rock interaction.

[0098] The fracture-matrix upscaled geological structure model is constructed as follows:

[0099] a. After determining the grid scheme, perform upscaling calculations of porosity and permeability tensors. The porosity φ of the fracture-matrix upscaled model is calculated as follows:

[0100]

[0101] In the formula, when the grid cell is the matrix, φ = φ M , when the grid cell is the fracture, φ = φ F ; V F is the total fracture volume intersecting within the model cell; V C is the model cell volume; N is the total number of fractures intersecting within the model cell; A f is the fracture area within the cell; b f is the fracture width;

[0102] b. The permeability tensor [k] of the fracture-matrix model upscaled model is expressed as:

[0103]

[0104] In the formula, when the grid cell is the matrix, [k] = [k] M , where k xx , k yy , k zz are the permeability components of the matrix in the x, y, and z directions. When the grid cell is the fracture, [k] = [k] F ; φ fis the porosity of the fracture f within the simulation unit; b f is the width of the fracture f; is the coordinate transformation tensor; n1, n2, n3 are the orthogonal vectors of the fracture f.

[0105] The upscaling acidic water reactive solute transport model for the fracture-matrix is constructed as follows:

[0106] a. Groundwater flow equation

[0107]

[0108] In the formula, φ is the porosity; [k] is the permeability tensor; μ is the fluid viscosity coefficient; P is the pressure; g is the acceleration due to gravity; ρ is the fluid density; w is the source-sink term; z is the position head;

[0109] b. Groundwater solute transport equation:

[0110]

[0111] In the formula, C j is the total concentration of the jth solute; V is the groundwater flow velocity; D is the hydrodynamic dispersion coefficient tensor; w j is the source-sink term of the jth solute; M is the total number of mineral phase reactions involved in the jth solute; ξ m is the stoichiometric coefficient of the jth solute in the reaction with the mineral; R jm is the reaction rate of the jth solute participating in the mineral reaction;

[0112] c. Groundwater reaction system control equation:

[0113]

[0114] In the formula, C j is the total concentration of the jth solute; X j is the free ion concentration of the jth solute; N is the total number of liquid phase reactions participated in by the jth solute; ξ i is the stoichiometric coefficient of the jth solute in the ith liquid phase reaction; K i is the equilibrium constant of the ith liquid phase reaction; a i is the activity coefficient of the product of the ith liquid phase reaction; a j is the activity coefficient of the jth solute in the ith liquid phase reaction;

[0115]

[0116] In the formula, R m is the mineral reaction rate; A m is the specific surface area of the mineral; k m(T, λ) is the effective reaction rate, T is the temperature influence parameter, and λ is the ion influence pre-parameter; K eq is the reaction equilibrium constant; Q is the ion activity product; α is the mineral saturation sensitivity parameter; β is the reaction affinity parameter; ρ m is the mineral density; M m is the mineral molar mass; η is the mineral volume fraction.

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

[0118] 1) Conceptual model

[0119] The example model is a cube region with a side length of 20 m. It is assumed that the matrix is a homogeneous isotropic medium, and its permeability, porosity, and tortuosity are 1.0×10 -15 m 2 , 0.01, 0.5 respectively, and the longitudinal, transverse, and vertical dispersivities are set to 1.0 m, 0.001 m, and 0.001 m respectively. The fracture radius ranges from 2 to 10 m, and the direction is randomly distributed. The fracture width has a positive correlation with the fracture radius. The conceptual model is as Figure 2 shown. The YZ plane of the coordinate is set as the given head boundary, the direction of groundwater flow is consistent with the X-axis direction, and the other surfaces are all set as impermeable boundaries or zero-flux boundaries. The head difference between the upstream and downstream boundaries is 10 m. The locally refined fracture-matrix upscaling model generated by the example is as Figure 3 shown. This simulation mainly conducts the simulation of groundwater flow and solute transport in the fractured confined aquifer. The pH of the background groundwater in the model is 8, which reaches equilibrium with the calcite aquifer, and acidic water (pH = 4) infiltrates from the upstream boundary.

[0120] Table 1 Concentrations of main components in the example of the invention (mol / L)

[0121]

[0122] Table 2 Chemical reaction parameters in the example of the invention

[0123]

[0124] The upstream boundary of the solute transport model is set as a fixed concentration boundary, and the given acidic water pollution concentration and the background concentration of fracture groundwater are shown in Table 1. The main simulated components include H + , Ca 2+ , HCO3 2- and tracers, and 1 kind of calcite mineral with an initial volume fraction of 1×10 -5 . The main reaction is the dissolution-precipitation reaction of carbonate rocks, and the chemical reaction parameters are shown in Table 2.

[0125] 2) Determine the control equations

[0126] (1) Governing equations of the fracture-matrix upscaling geological structure model:

[0127] a. After determining the grid scheme, perform upscaling calculations of porosity and permeability tensors. The porosity φ of the fracture-matrix upscaling model is calculated as follows:

[0128]

[0129] In the formula, when the grid cell is the matrix, φ = φ M , when the grid cell is the fracture, φ = φ F ; V F is the total fracture volume intersecting within the model cell; V C is the model cell volume; N is the total number of fractures intersecting within the model cell; A f is the fracture area within the cell; b f is the fracture width;

[0130] b. The permeability tensor [k] of the fracture-matrix model upscaling model is expressed as:

[0131]

[0132] In the formula, when the grid cell is the matrix, [k] = [k] M , where k xx , k yy , k zz are the permeability components of the matrix in the x, y, and z directions. When the grid cell is the fracture, [k] = [k] F ; φ f is the porosity of fracture f within the simulation cell; b f is the width of fracture f; is the coordinate transformation tensor; n1, n2, n3 are the orthogonal vectors of fracture f.

[0133] (2) Governing equations of the reactive solute transport model for acidic water in the fracture-matrix upscaling:

[0134] a. Groundwater flow equation:

[0135]

[0136] In the formula, φ is the porosity; [k] is the permeability tensor; μ is the fluid viscosity; P is the pressure; g is the acceleration due to gravity; ρ is the fluid density; w is the source-sink term; z is the position head;

[0137] b. Groundwater solute transport equation:

[0138]

[0139] In the formula, Cj is the total concentration of the j-th solute; V is the groundwater flow velocity; D is the hydrodynamic dispersion coefficient tensor; w j is the source-sink term of the j-th solute; M is the total number of mineral phase reactions involving the j-th solute; ξ m is the stoichiometric coefficient of the j-th solute in the reaction with minerals; R jm is the reaction rate of the j-th solute participating in mineral reactions;

[0140] c. Governing equations for the groundwater reaction system:

[0141]

[0142] In the formula, C j is the total concentration of the j-th solute; X j is the free ion concentration of the j-th solute; N is the total number of liquid phase reactions involving the j-th solute; ξ i is the stoichiometric coefficient of the j-th solute in the i-th liquid phase reaction; K i is the equilibrium constant of the i-th liquid phase reaction; a i is the activity coefficient of the product in the i-th liquid phase reaction; a j is the activity coefficient of the j-th solute in the i-th liquid phase reaction;

[0143]

[0144] In the formula, R m is the mineral reaction rate; A m is the specific surface area of the mineral; k m (T,λ) is the effective reaction rate, T is the temperature influence parameter, λ is the ion influence pre-parameter; K eq is the reaction equilibrium constant; Q is the ion activity product; α is the mineral saturation sensitivity parameter; β is the reaction affinity parameter; ρ m is the mineral density; M m is the molar mass of the mineral; η is the mineral volume fraction.

[0145] 3) Analysis of numerical simulation results of the reactive solute transport model for acidic water in fracture-matrix upscaling

[0146] As Figure 4 shown, in the 10th year of solute migration, the H + concentration spatial distribution map shows that acidic water mainly migrates in the preferential flow channels formed by the fracture network. At the same time, due to the reaction with carbonate rocks, it diffuses to the surrounding matrix. At the inflow boundary, the H + concentration diffuses and migrates into the matrix by about 3m. However, at this time, the H+ in the fracture channel has migrated to the outflow boundary, indicating that the preferential flow formed by fractures controls the solute migration range.

[0147] AsFigure 5 As shown, since acidic water reacts with calcite to form Ca 2+ , and its concentration value is greater than that of the tracer. It can be seen from the typical cross-section of the fracture-matrix continuum that a relatively high concentration of Ca 2+ is formed around the fracture channel, and the Ca 2+ concentration in the matrix adjacent to the fracture diffuses continuously with time. However, the concentration distribution of the tracer does not show a significant increase around the fracture, and its concentration value gradually decreases with time.

[0148] As Figure 6 shown in a, due to the reaction of acidic water with calcite, the concentration curves of H + all show lower average concentration values at the outflow boundary during the simulation period compared to the tracer. Due to the reaction of acidic water, Ca 2+ is continuously released by the reaction. As can be seen from Figure 6 b, the change in Ca 2+ concentration first increases slightly and then changes steadily, and finally increases rapidly, indicating that the dynamic change of the reactive solute concentration under the action of fracture preferential flow and acidic water reaction is strongly non-linear. For the tracer concentration corresponding to Ca 2+ , it shows a slow decreasing trend, which is mainly because the Ca 2+ tracer with a lower concentration flowing in at the boundary continuously dilutes the tracer concentration in the fracture preferential channel, resulting in a slow decreasing trend of the average concentration at the outflow boundary, indicating that the fracture preferential flow channel controls the flux of solute outflow downstream.

[0149] The above are only the preferred embodiments of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the technical principle of the present invention, several improvements and deformations can be made, and these improvements and deformations should also be regarded as the protection scope of the present invention.

Claims

1. A method for simulating the reaction and migration of acidic water under a fracture-matrix upscaling model, characterized in that: The method comprises the following steps: Step SS1, determining the basic parameters required for the fracture-matrix upscaling model generated by the local refined grid, the basic model parameters including: fracture distribution empirical function, fracture number, fracture radius, fracture direction, fracture width, matrix permeability, matrix porosity, local refined grid encryption level, and basic grid unit length; Step SS2, performing local fine grid division according to the grid setting parameters in step SS1, and scaling the corresponding fracture-matrix part grid parameters to establish a fracture-matrix scaled geological structure model; Step SS3, based on the fracture-matrix upscaled geological structure model built in step SS2, determining the groundwater parameters required by the model, including: initial head distribution, boundary conditions, permeability after the fracture-matrix parameters are upscaled in step SS2, solving the fracture groundwater level and flow velocity distribution in the simulated area to obtain the initial groundwater flow field; Step SS4, based on the initial groundwater flow field obtained in step SS3, determining the physicochemical properties, background concentration, and infiltration concentration of the acidic water pollution components to calculate the spatial distribution of the solute components in the fracture-matrix of the simulation area; Step SS5, based on the solute migration process obtained in step SS4, determine the hydrogeochemical reaction parameters involved in the acidic water pollution migration process, including: chemical reaction formula and reaction rate, mineral volume fraction and specific surface area, couple the groundwater flow field-solute field-chemical field, establish a fracture-matrix up-scale acidic water reactive solute migration model, and simulate the acidic water pollution migration process under the combined influence of fracture-matrix spatial heterogeneity and water-rock interaction.

2. The method for simulating the reaction and migration of acidic water under a fracture-matrix upscaling model according to claim 1, characterized in that: In step SS2, the fracture-matrix upscaling geological structure model is constructed as follows: a. After the grid scheme is determined, the porosity and permeability are upscaled. The porosity φ of the fracture-matrix upscaling model is calculated as follows: In the formula, when the grid unit is the matrix, φ = φ M ,φ M is the matrix unit porosity. When the grid unit is a crack, φ = φ F ,φ F is the fracture unit porosity; V F is the total fracture volume intersecting within the model unit; V C is the volume of the model unit; N is the total number of fractures intersecting within the model unit; A f is the crack area within the unit; b f is the crack width; b. The permeability tensor [k] of the fracture-matrix model upscaling model is expressed as: In the formula, when the grid unit is the matrix, [k] = [k] M , [k] M is the matrix unit permeability tensor, where k xx , k yy , k zz is the permeability component of the matrix in the x, y, and z directions. When the grid unit is a fracture, [k] = [k] F , [k] F is the fracture unit permeability tensor; φ f is the porosity of the fracture f in the simulation unit; b f is the width of the crack f; is the coordinate transformation tensor; n1, n2, n3 are the orthogonal vectors of the crack f.

3. The method for simulating the reaction and migration of acidic water under a fracture-matrix upscaling model according to claim 1, characterized in that: In step SS5, the fracture-matrix up-scaling acidic water reactive solute transport model is constructed as follows: a. Groundwater flow equation Where φ is porosity; [k] is the permeability tensor; μ is the fluid viscosity; P is pressure; g is the acceleration of gravity; ρ is the fluid density; w is the source and sink term; z is the position head; b. Groundwater solute transport equation In the formula, C j is the total concentration of the jth solute; V is the groundwater velocity; D is the hydrodynamic diffusion coefficient tensor; w j is the source and sink term of the j-th solute; M is the total number of mineral phase reactions involved in the j-th solute; ξ m is the stoichiometric coefficient of the reaction between the jth solute and the mineral; R jm is the rate at which the jth solute participates in the mineral reaction; c. Governing equations of groundwater reaction system In the formula, C j is the total concentration of the jth solute; X j is the free ion concentration of the j-th solute; N1 is the total number of the j-th solute participating in the liquid phase reaction; ξ i is the stoichiometric number of the jth solute in the i-th liquid phase reaction; K i is the equilibrium constant of the i-th liquid phase reaction; a i is the activity coefficient of the i-th liquid phase reaction product; a j is the activity coefficient of the jth solute in the i-th liquid phase reaction; In the formula, R m is the mineral reaction rate; A m is the specific surface area of ​​the mineral; k m (T,λ) is the effective reaction rate, T is the temperature effect parameter, λ is the ion effect pre-parameter; K eq is the reaction equilibrium constant; Q is the ion activity product; α is the mineral saturation sensitive parameter; β is the reaction affinity parameter; ρ m is the mineral density; M m is the mineral molar mass; η is the mineral volume fraction.

4. A system for simulating the reaction and migration of acidic water under a fracture-matrix upscaling model, characterized in that: The system includes the following modules: A model parameter input module is used to determine the basic parameters required for the fracture-matrix upscaling model generated by the local refined grid, wherein the basic model parameters include: fracture distribution empirical function, fracture number, fracture radius, fracture direction, fracture width, matrix permeability, matrix porosity, local refined grid encryption level, and basic grid unit length; The fracture-matrix upscaling module is used to perform local fine grid division according to the grid setting parameters in the model parameter input module, and to perform parameter upscaling on the corresponding fracture-matrix part grid to establish a fracture-matrix upscaling geological structure model; The groundwater dynamics calculation module is used to determine the groundwater parameters required by the model based on the fracture-matrix upscaling geological structure model, including: initial water head distribution, boundary conditions, permeability after the fracture-matrix parameters are upscaled in step SS2, solve the fracture groundwater level and flow velocity distribution in the simulation area, and obtain the initial groundwater flow field; The solute transport calculation module is used to determine the physicochemical properties, background concentration, and infiltration concentration of acidic water pollution components based on the initial groundwater flow field obtained, so as to calculate the spatial distribution of solute components in the fracture-matrix of the simulation area; The reactive solute transport module is used to determine the hydrogeochemical reaction parameters involved in the acidic water pollution migration process based on the obtained solute migration process, including: chemical reaction formula and reaction rate, mineral volume fraction and specific surface area, couple the groundwater flow field-solute field-chemical field, establish a fracture-matrix up-scale acidic water reactive solute migration model, and simulate the acidic water pollution migration process under the combined influence of fracture-matrix spatial heterogeneity and water-rock interaction.

5. The acidic water reaction migration simulation system under the fracture-matrix upscaling model according to claim 4, characterized in that: The fracture-matrix upscaling geological structure model is constructed as follows: a. After the grid scheme is determined, the porosity and permeability tensors are upscaled. The porosity φ of the fracture-matrix upscaling model is calculated as follows: In the formula, when the grid unit is the matrix, φ = φ M ,φ M is the matrix unit porosity. When the grid unit is a crack, φ = φ F ,φ F is the fracture unit porosity; V F is the total fracture volume intersecting within the model unit; V C is the volume of the model unit; N is the total number of fractures intersecting within the model unit; A f is the crack area within the unit; b f is the crack width; b. The permeability tensor [k] of the fracture-matrix model upscaling model is expressed as: In the formula, when the grid unit is the matrix, [k] = [k] M , [k] M is the matrix unit permeability tensor, where k xx , k yy , k zz is the permeability component of the matrix in the x, y, and z directions. When the grid unit is a fracture, [k] = [k] F , [k] F is the fracture unit permeability tensor; φ f is the porosity of the fracture f in the simulation unit; b f is the width of the crack f; is the coordinate transformation tensor; n1, n2, n3 are the orthogonal vectors of the crack f.

6. The acidic water reaction migration simulation system under the fracture-matrix upscaling model according to claim 4, characterized in that: The fracture-matrix up-scale model of acidic water reactive solute transport is constructed as follows: a. Groundwater flow equation: Where φ is porosity; [k] is the permeability tensor; μ is the fluid viscosity; P is pressure; g is the acceleration of gravity; ρ is the fluid density; w is the source and sink term; z is the position head; b. Groundwater solute transport equation: In the formula, C j is the total concentration of the jth solute; V is the groundwater velocity; D is the hydrodynamic diffusion coefficient tensor; w j is the source and sink term of the j-th solute; M is the total number of mineral phase reactions involved in the j-th solute; ξ m is the stoichiometric coefficient of the reaction between the jth solute and the mineral; R jm is the rate at which the jth solute participates in the mineral reaction; c. Groundwater reaction system control equation: In the formula, C j is the total concentration of the jth solute; X j is the free ion concentration of the j-th solute; N is the total number of the j-th solute participating in the liquid phase reaction; ξ i is the stoichiometric number of the jth solute in the i-th liquid phase reaction; K i is the equilibrium constant of the i-th liquid phase reaction; a i is the activity coefficient of the i-th liquid phase reaction product; a j is the activity coefficient of the jth solute in the i-th liquid phase reaction; In the formula, R m is the mineral reaction rate; A m is the specific surface area of ​​the mineral; k m (T,λ) is the effective reaction rate, T is the temperature effect parameter, λ is the ion effect pre-parameter; K eq is the reaction equilibrium constant; Q is the ion activity product; α is the mineral saturation sensitive parameter; β is the reaction affinity parameter; ρ m is the mineral density; M m is the mineral molar mass; η is the mineral volume fraction.

Images