Mound aquifer well group coordinated dewatering simulation method and system based on digital twinning

Abstract

The application provides a butte aquifer well group coordinated dewatering simulation method and system based on digital twinning, relates to the technical field of digital twinning, and comprises the following steps: constructing an aquifer digital twinning model, performing virtual dewatering deduction, obtaining a water level space-time evolution field, identifying a hydraulic coupling region, according to a water level gradient distribution and a flow field deflection direction, performing coordinated regulation and control on well group pumping intensity, and determining an optimal pumping scheme. The application can effectively improve well group dewatering efficiency, realize overall water level uniform lowering of the aquifer, and reduce engineering cost.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of digital twinning, in particular to a simulation method and system for well group coordinated dewatering of a hill and mound aquifer based on digital twinning. BACKGROUND

[0002] In the fields of urban construction, mine exploitation, and water resource management, dewatering engineering of hill and mound aquifers is a common and important engineering measure. Hill and mound aquifers usually have complex geological structures and hydrological characteristics, and it is necessary to achieve precise dewatering by reasonably arranging well groups and scientifically pumping water. Traditional dewatering engineering of hill and mound aquifers mainly relies on empirical formulas and simplified models for design, and usually estimates the dewatering effect by simply superimposing single wells or multiple wells, which is difficult to accurately predict the dewatering process and effect under complex geological conditions in engineering practice.

[0003] With the development of computer technology and numerical simulation methods, groundwater numerical simulation technology is widely used in dewatering engineering design, which can more accurately describe the groundwater flow law. Current dewatering engineering simulation mainly uses numerical methods such as finite difference or finite element to predict the dewatering process and effect by solving the groundwater flow control equation. However, these methods still have some limitations and problems in practical application.

[0004] The defects and deficiencies of the prior art mainly exist in the following aspects: first, the traditional dewatering simulation method cannot fully consider the spatial heterogeneity and complex geological characteristics of the hill and mound aquifer, and the description of the stratum permeability distribution is too simplified, which cannot accurately reflect the groundwater flow law under actual geological conditions. Second, the existing technology does not accurately represent the interaction and hydraulic coupling relationship between well groups, often ignoring the superposition interference effect of the water level drawdown cone between adjacent well points, resulting in large errors in dewatering effect prediction. Finally, the traditional well group design method lacks a dynamic coordinated regulation mechanism for the pumping process, and cannot optimize and adjust the pumping intensity of each well point according to the real-time water level change, making it difficult to achieve the goal of uniform water level drop of the hill and mound aquifer as a whole, and easily causing problems of local over-dewatering or insufficient dewatering. SUMMARY

[0005] The embodiment of the present application provides a simulation method and system for well group coordinated dewatering of a hill and mound aquifer based on digital twinning, which can solve the problems in the prior art.

[0006] In a first aspect, the embodiment of the present application provides a simulation method for well group coordinated dewatering of a hill and mound aquifer based on digital twinning, comprising:

[0007] Obtaining geological parameters of the hillside aquifer, well group layout information and target precipitation depth requirements, constructing a digital twin model of the hillside aquifer, the digital twin model including stratum permeability distribution representation, point source pumping boundary conditions and hydraulic coupling relationship between well group and aquifer;

[0008] Based on the digital twin model, virtual precipitation deduction is carried out on the well group, and the spatio-temporal evolution field of the water level in the hillside aquifer is obtained by solving the seepage control equation set of three-dimensional unsteady groundwater in the hillside aquifer, and the pumping influence radius and water level drawdown funnel shape of each well point are extracted;

[0009] According to the superposition interference characteristics between the water level drawdown funnels of adjacent well points in the spatio-temporal evolution field of the water level, the hydraulic coupling region under the synergistic action of the well group is identified, and the water level gradient distribution and flow field deflection direction in the hydraulic coupling region are determined;

[0010] According to the water level gradient distribution and flow field deflection direction of the hydraulic coupling region, the pumping intensity of each well point in the well group is synergistically regulated, an evaluation index function that meets the target precipitation depth requirement and optimizes the overall water level drawdown uniformity of the hillside aquifer is generated, and the minimum of the evaluation index function is taken as the optimization goal, and the water level drawdown constraint condition is taken as the boundary limit to solve the well group synergistic pumping scheme.

[0011] Based on the digital twin model, virtual precipitation deduction is carried out on the well group, and the spatio-temporal evolution field of the water level in the hillside aquifer is obtained by solving the seepage control equation set of three-dimensional unsteady groundwater in the hillside aquifer, and the pumping influence radius and water level drawdown funnel shape of each well point are extracted including:

[0012] The stratum permeability distribution representation in the digital twin model is converted into a spatially discretized permeability coefficient field, and the well group layout information is converted into point source pumping boundary conditions;

[0013] Based on the permeability coefficient field and the point source pumping boundary conditions, the seepage control equation set of three-dimensional unsteady groundwater in the hillside aquifer is established, and a gravity potential correction term reflecting the hillside topographic relief characteristics is introduced into the seepage control equation set to represent the influence of topographic elevation on the flow direction and velocity distribution of groundwater;

[0014] The seepage control equation set is solved in the time and space domain to obtain the water level values of each spatial node in the hillside aquifer over a continuous time sequence, forming the spatio-temporal evolution field of the water level;

[0015] From the spatio-temporal evolution field of the water level, the spatial boundary around each well point that meets the preset drawdown criterion is extracted to determine the pumping influence radius of each well point;

[0016] The geometric features of the water level contour lines within the radius of influence of the pumping are extracted as the shape of the water level drawdown funnel. The shape of the water level drawdown funnel includes the offset of the funnel center position, the direction angle of the funnel's major axis and minor axis, and the asymmetry index of the funnel boundary.

[0017] Based on the permeability coefficient field and the point source pumping boundary conditions, a three-dimensional unsteady groundwater seepage control equation set for the hill aquifer is established. The seepage control equation set includes a gravitational potential energy correction term reflecting the hill topographic undulation characteristics, which includes:

[0018] A three-dimensional spatial discrete grid is obtained for the aquifer in the hilly area. The three-dimensional spatial discrete grid includes the grid node coordinates in the horizontal direction and the grid node elevations in the vertical direction.

[0019] For each grid node in the three-dimensional discrete grid, the permeability coefficient value at the corresponding location of the grid node is extracted from the permeability coefficient field, and the pumping intensity value at the corresponding location of the grid node is extracted from the point source pumping boundary condition;

[0020] The terrain elevation gradient of a grid node is calculated based on the grid node elevation. The terrain elevation gradient is obtained by calculating the ratio of the elevation difference between the grid node and its adjacent grid nodes to the horizontal distance.

[0021] The gravitational potential energy correction coefficient for the grid node is calculated based on the terrain elevation gradient of the grid node. The gravitational potential energy correction coefficient is obtained by multiplying the terrain elevation gradient by the gravitational acceleration constant and normalizing it by combining the aquifer thickness at the grid node.

[0022] Based on the permeability coefficient, pumping intensity, and gravitational potential energy correction coefficient of each grid node in the three-dimensional spatial discrete grid, a set of three-dimensional unsteady groundwater seepage control equations for the hill aquifer is constructed.

[0023] Based on the superimposed interference characteristics between the drawdown cones of adjacent well points in the spatiotemporal evolution field of water level, the hydraulic coupling region under the synergistic effect of the well group is identified, and the water level gradient distribution and flow field deflection direction within the hydraulic coupling region are determined, including:

[0024] Water level data at multiple time points at each spatial location are extracted from the spatiotemporal evolution field of water level, and the water level gradient vector at each spatial location is calculated;

[0025] Construct the hydraulic gradient vector field of the hill aquifer based on the water level gradient vector;

[0026] Calculate the divergence distribution of the hydraulic gradient vector field, construct the theoretical hydraulic gradient vector field under the condition of independent pumping of a single well, and calculate the theoretical divergence distribution of the theoretical hydraulic gradient vector field;

[0027] By comparing the divergence distribution with the theoretical divergence distribution, when the divergence distribution shows an abnormally positive value cluster or an abnormally negative value cluster and the region does not coincide with any well point, the region is determined to be a hydraulic gradient field disturbance region caused by the synergistic effect of the well group, and the hydraulic gradient field disturbance region is determined to be the hydraulic coupling region;

[0028] Within the hydraulic coupling region, the amplitude and direction distribution of the water level gradient vector are statistically analyzed to obtain the water level gradient distribution within the hydraulic coupling region;

[0029] Based on the water level gradient vector, a groundwater flow vector field is constructed within the hydraulic coupling region. By performing streamline tracing on the groundwater flow vector field, the deflection angle distribution of the groundwater streamlines relative to the radial flow direction of a single well is determined, and the dominant direction of the deflection angle distribution is determined as the flow field deflection direction.

[0030] Based on the water level gradient vector, a groundwater flow vector field is constructed within the hydraulic coupling region. By performing streamline tracing on the groundwater flow vector field, the deflection angle distribution of the groundwater streamlines relative to the radial flow direction of a single well is determined, and the dominant direction of the deflection angle distribution is determined as the flow field deflection direction, including:

[0031] Multiple initial seed points for streamline tracing are set in the groundwater flow vector field, and the initial seed points are uniformly distributed along the boundary of the hydraulic coupling region;

[0032] Starting from each of the initial seed points, streamline tracing is performed along the vector direction of the groundwater flow vector field to obtain the spatial distribution of groundwater streamlines within the hydraulic coupling region. The tangent direction vector at each point on the groundwater streamline is extracted, and the tangent direction vector represents the actual flow direction of groundwater at that point.

[0033] Calculate the location of the well point closest to the groundwater flow line, and determine the theoretical direction vector of single-well radial flow at each point on the groundwater flow line. The theoretical direction vector of single-well radial flow is the radial unit vector pointing from each point on the groundwater flow line to the well point.

[0034] Calculate the angle between the tangential direction vector and the theoretical radial flow direction vector of the single well, define the angle as the deflection angle of the groundwater streamline relative to the radial flow direction of the single well, and statistically analyze the deflection angles of all groundwater streamlines within the hydraulic coupling region to form the deflection angle distribution;

[0035] The dominant direction of the deflection angle distribution is determined as the flow field deflection direction.

[0036] The evaluation index function for coordinating and controlling the pumping intensity of each well point in the well group to generate an evaluation index function that satisfies the target precipitation depth requirement and optimizes the overall uniformity of water level decline in the hill aquifer includes:

[0037] The water level gradient distribution is divided into multiple water level gradient grade intervals. The well point with the largest area-gradient weighted index in the water level gradient grade interval is determined as the key well point for gradient control. The spatial position deflection angle of each key well point for gradient control relative to the reference direction is calculated using the flow field deflection direction as the reference direction.

[0038] Based on the spatial location deflection angle, the key well points for gradient control are divided into a forward well point group and a reverse well point group. The well points in the forward well point group are located on the forward extension path of the flow field deflection direction, and the well points in the reverse well point group are located in the reverse hindrance region of the flow field deflection direction. For the forward well point group, a pumping intensity enhancement coefficient is set, and for the reverse well point group, a pumping intensity attenuation coefficient is set.

[0039] Based on the pumping intensity enhancement coefficient and the pumping intensity attenuation coefficient, the pumping intensity values ​​of each well point in the forward well point group and the reverse well point group are adjusted; based on the adjusted pumping intensity values ​​of each well point, the spatial distribution of water level contour lines within the hilly aquifer is extracted;

[0040] For the spatial distribution of the water level contour lines, the standard deviation of the curvature of the water level contour lines and the spatial uniformity of the density distribution of the water level contour lines are calculated. The standard deviation of curvature and the spatial uniformity are weighted and combined to construct an evaluation index function for the overall uniformity of water level decline in the hilly aquifer.

[0041] The method further includes:

[0042] The well group coordinated pumping scheme is simulated and verified in the digital twin model to obtain the deviation between the actual precipitation effect and the expected target. The deviation is then fed back to the digital twin model to update the correction coefficient of the formation permeability distribution characterization and hydraulic coupling relationship, thereby achieving dynamic alignment between the digital twin model and the real hydrological response characteristics of the hill aquifer.

[0043] A second aspect of the present invention provides a digital twin-based collaborative precipitation simulation system for well groups in hilly aquifers, comprising:

[0044] The first unit is used to obtain the geological parameters of the hill aquifer, well layout information, and target precipitation depth requirements, and to construct a digital twin model of the hill aquifer. The digital twin model includes a representation of the formation permeability distribution, point source pumping boundary conditions, and the hydraulic coupling relationship between the well group and the aquifer.

[0045] The second unit is used to perform virtual precipitation simulation of the well group based on the digital twin model. By solving the three-dimensional unsteady groundwater seepage control equations of the hill aquifer, the spatiotemporal evolution field of water level in the hill aquifer is obtained, and the pumping influence radius and water level drawdown funnel shape of each well point are extracted.

[0046] The third unit is used to identify the hydraulic coupling region under the synergistic effect of the well group based on the superimposed interference characteristics between the drawdown funnels of adjacent well points in the spatiotemporal evolution field of the water level, and to determine the water level gradient distribution and flow field deflection direction within the hydraulic coupling region;

[0047] The fourth unit is used to coordinate and regulate the pumping intensity of each well point in the well group based on the water level gradient distribution and flow field deflection direction in the hydraulic coupling region, generate an evaluation index function that meets the target drawdown depth requirement and optimizes the overall water level drop uniformity of the hill aquifer, and solve the well group coordinated pumping scheme with the minimization of the evaluation index function as the optimization objective and the drawdown constraint as the boundary condition.

[0048] A third aspect of the present invention provides an electronic device, comprising:

[0049] processor;

[0050] Memory used to store processor-executable instructions;

[0051] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0052] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0053] The beneficial effects of this application are as follows:

[0054] This invention achieves precise simulation and optimized control of the precipitation process of well groups by constructing a digital twin model of the hill aquifer. It can intuitively reflect the spatial distribution and dynamic evolution of the water level drawdown cone, thereby revealing the hydraulic coupling mechanism between well groups.

[0055] Based on digital twin technology, this invention can identify the hydraulic coupling region under the coordinated action of well groups, clarify the water level gradient distribution and flow field deflection direction, provide a scientific basis for optimizing the pumping intensity of well groups, and effectively avoid the blindness and uncertainty in traditional empirical methods.

[0056] This invention employs an optimization method that minimizes the evaluation index function. While meeting the target precipitation depth, it achieves optimal uniformity of the overall water level drop in the hilly aquifer. This effectively solves the common problems of local over-pumping or insufficient precipitation in traditional precipitation methods, improves precipitation efficiency, reduces energy consumption, and has significant economic and environmental benefits. Attached Figure Description

[0057] Figure 1 This is a flowchart illustrating the collaborative precipitation simulation method for well groups in hilly aquifers based on digital twins, according to an embodiment of the present invention.

[0058] Figure 2 This is a flowchart illustrating the implementation of the virtual precipitation simulation technology for well groups in hilly aquifers based on digital twin models, as described in this invention. Detailed Implementation

[0059] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0060] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0061] Figure 1 This is a flowchart illustrating the collaborative precipitation simulation method for well groups in hilly aquifers based on digital twins, as described in an embodiment of the present invention. Figure 1 As shown, the method includes:

[0062] Geological parameters, well layout information, and target precipitation depth requirements of the hill aquifer are obtained. A digital twin model of the hill aquifer is then constructed. The digital twin model includes a representation of the formation permeability distribution, point source pumping boundary conditions, and the hydraulic coupling relationship between the well group and the aquifer.

[0063] Based on the digital twin model, virtual precipitation simulation was performed on the well group. By solving the three-dimensional unsteady groundwater seepage control equations of the hill aquifer, the spatiotemporal evolution field of water level within the hill aquifer was obtained, and the pumping influence radius and water level drawdown cone shape of each well point were extracted.

[0064] Based on the superimposed interference characteristics between the drawdown funnels of adjacent well points in the spatiotemporal evolution field of water level, the hydraulic coupling region under the synergistic effect of the well group is identified, and the water level gradient distribution and flow field deflection direction within the hydraulic coupling region are determined.

[0065] Based on the water level gradient distribution and flow field deflection direction in the hydraulic coupling region, the pumping intensity of each well point in the well group is coordinated and controlled to generate an evaluation index function that meets the target drawdown depth requirement and optimizes the overall water level drop uniformity of the hill aquifer. The optimization objective is to minimize the evaluation index function, and the well group coordinated pumping scheme is solved with the drawdown constraint as the boundary condition.

[0066] In one optional implementation, virtual precipitation simulation is performed on the well group based on the digital twin model. The spatiotemporal evolution field of the water level within the hill aquifer is obtained by solving the three-dimensional unsteady groundwater seepage control equations of the hill aquifer, and the pumping influence radius and water level drawdown cone morphology of each well point are extracted, including:

[0067] The formation permeability distribution representation in the digital twin model is converted into a spatially discretized permeability coefficient field, and the well group layout information is converted into point source pumping boundary conditions;

[0068] Based on the permeability coefficient field and the point source pumping boundary conditions, a three-dimensional unsteady groundwater seepage control equation set for the hill aquifer is established. The seepage control equation set introduces a gravity potential energy correction term that reflects the topographic undulation characteristics of the hill to characterize the influence of topographic elevation difference on the direction and velocity distribution of groundwater flow.

[0069] Solving the seepage control equations in the spatiotemporal domain yields the water level values ​​of each spatial node within the hill aquifer over a continuous time series, thus forming the spatiotemporal evolution field of the water level.

[0070] Extract the spatial boundaries around each well point that satisfy the preset drop criterion from the spatiotemporal evolution field of the water level, and determine the pumping influence radius of each well point.

[0071] The geometric features of the water level contour lines within the radius of influence of the pumping are extracted as the shape of the water level drawdown funnel. The shape of the water level drawdown funnel includes the offset of the funnel center position, the direction angle of the funnel's major axis and minor axis, and the asymmetry index of the funnel boundary.

[0072] like Figure 2 As shown, the method includes:

[0073] The transformation of stratigraphic permeability distribution is based on geological parameter data from a digital twin model. The digital twin model stores the permeability distribution information of hill aquifers in the form of a three-dimensional geological volume, including the permeability coefficient values ​​and spatial distribution range of different geological units. The transformation process establishes a regular three-dimensional grid system, with the grid size set at 10 meters by 10 meters horizontally and 2 meters thick vertically. For each node in the grid, the geological unit to which the node belongs is determined through spatial positioning, and the permeability coefficient value of the corresponding geological unit is extracted from the digital twin model. When a grid node is located at the boundary of different geological units, a volume-weighted average method is used to calculate the comprehensive permeability coefficient, with the weight being the volume percentage of each geological unit within that grid unit.

[0074] The conversion of well cluster layout information transforms the well cluster configuration data in the digital twin model into the boundary condition format required for numerical computation. The well cluster layout information includes parameters such as the spatial coordinates, depth, pumping intensity, and pumping schedule for each well. The point source pumping boundary conditions establish each well as a point source in space, with the point source location corresponding to the well's coordinate position. The spatial distribution of pumping intensity adopts a radial allocation method around the wellbore, with grid nodes within a certain range from the well axis designated as pumping nodes. The numerical allocation of pumping intensity is based on the well's filter length and pumping layer, distributing the total pumping volume according to the distribution ratio of the filters in each grid layer.

[0075] The seepage control equations are established based on the fundamental control principle of three-dimensional groundwater seepage. The core of the equations is the mass conservation equation, which expresses the balance between the input, output, and storage changes of groundwater at each grid node. The storage change term is calculated by multiplying the time derivative of the water level at the grid node by the storage coefficient and the grid cell volume. The seepage term includes horizontal and vertical water exchange. The horizontal seepage flow is calculated by multiplying the water level difference between adjacent nodes, the permeability coefficient, and the cross-sectional area of ​​the water passage. The vertical seepage flow is calculated by multiplying the water level difference between adjacent nodes, the vertical permeability coefficient, and the horizontal cross-sectional area of ​​the water passage.

[0076] The introduction of the gravitational potential energy correction term considers the impact of hilly terrain undulations on groundwater flow. The correction term is calculated based on the product of the terrain elevation gradient and gravitational acceleration, reflecting the change in gravitational potential energy caused by terrain tilt. For each node in the grid, the terrain elevation gradient at that node is calculated, including gradient components in the east-west and north-south directions. The terrain elevation gradient is calculated by dividing the elevation difference between adjacent nodes by the distance between nodes. The gravitational potential energy correction term is added as an additional driving force to the seepage control equation. The correction coefficient is calculated by dividing the product of the terrain elevation gradient and gravitational acceleration by the reference permeability coefficient, ensuring dimensional consistency between the correction term and other seepage terms.

[0077] The spatiotemporal domain solution employs an implicit difference scheme using the finite difference method. Time domain discretization divides continuous time into equal-length time steps, with each time step set to one hour. Spatial domain discretization is based on an established three-dimensional mesh system, where each mesh node corresponds to an unknown in the equation system. Within each time step, the seepage control equations constitute a linear algebraic equation system with a sparse, symmetric, positive definite coefficient matrix. The equation system is solved using a preconditional conjugate gradient iteration method, with the precondition matrix employing either diagonal preconditions or incomplete Cholsky decomposition preconditions. The iteration termination condition is that the L2 norm of the residual vector is less than a set convergence threshold.

[0078] The spatiotemporal evolution field of water level is formed by collecting the solution results of all time steps. For each spatial node in the grid, its water level value at each time step is recorded, forming the water level time series of that node. The water level time series of all spatial nodes are combined to form a complete spatiotemporal evolution field of water level. The data structure of the evolution field adopts a four-dimensional array form, containing three spatial dimensions and one temporal dimension. To facilitate subsequent analysis, data interpolation processing is performed on the spatiotemporal evolution field of water level to achieve continuous representation in both spatial and temporal dimensions.

[0079] The radius of influence of pumping is determined based on the spatial distribution characteristics of water level drop. A preset drop criterion is set at a drop of 0.1 meters, reflecting the critical condition under which pumping significantly impacts the surrounding groundwater level. For each well, water level distribution data under stable pumping conditions are extracted from the spatiotemporal evolution field of the water level. The water level drop at each spatial node is calculated, i.e., the difference between the initial and current water levels at that node. A contour map of the water level drop is plotted using spatial interpolation, identifying contour lines where the drop is equal to 0.1 meters. The boundary of the area enclosed by these contour lines is the pumping influence range, and the radius of influence is defined as the radius of the equivalent circle of this area, calculated by dividing the area of ​​the influence region by the square root of pi.

[0080] The extraction of the water level drawdown funnel shape is based on the geometric analysis of water level contour lines within the pumping influence radius. Multiple water level contour lines are selected for morphological feature analysis. The contour line values ​​are selected sequentially outwards at 0.5-meter intervals, starting from the lowest water level near the well point. Geometric shape identification is performed on each contour line, and the optimal ellipse approximation is determined using an ellipse fitting method. The ellipse fitting employs the least squares method, determining the ellipse parameters by minimizing the sum of the squared distances from each point on the contour line to the ellipse.

[0081] The offset of the funnel center position is calculated by comparing the positional difference between the theoretical funnel center and the actual funnel center. The theoretical funnel center is the coordinate position of the well point, while the actual funnel center is determined by the center of the ellipse obtained from the contour ellipse fitting result. A representative contour line is selected where the water level drop is half the maximum drop at the well point. The distance between the center of this contour ellipse and the well point position is calculated; this distance is the offset of the funnel center position. The direction of the offset is represented by the azimuth angle of the ellipse center relative to the well point position.

[0082] The direction angles of the major and minor axes of the funnel are determined based on the geometric parameters of the contour ellipse fitting results. The direction angle of the major axis of the ellipse is defined as the angle between the major axis and true north, and the direction angle of the minor axis is the angle between the minor axis and true north. The direction angles are calculated directly through the rotation angle parameter of the ellipse, which reflects the inclination of the ellipse's principal axis relative to the coordinate axes. Statistical analysis is performed on the ellipse fitting results of multiple contour lines to calculate the mean and standard deviation of the direction angles of the major and minor axes, thus obtaining the dominant directional characteristics of the falling funnel.

[0083] The asymmetry index of the funnel boundary is calculated based on the symmetry analysis of the contour lines. The asymmetry index is obtained by comparing the differences in the radius length of the contour lines in different directions. The contour lines are divided into eight fan-shaped regions around the center of the funnel, each with an angle of 45 degrees. The average distance from the contour lines to the center within each fan-shaped region is calculated, and the standard deviation of the average distance in the eight directions is statistically analyzed. The asymmetry index is defined as the ratio of the standard deviation to the mean of the average distance; the larger this ratio, the more asymmetric the funnel shape. A circular funnel with good symmetry has an asymmetry index close to zero, while highly asymmetric elliptical or irregularly shaped funnels have a larger asymmetry index.

[0084] The temporal evolution analysis of funnel morphology parameters is achieved by tracking changes in funnel morphology at different times. Multiple time points are selected during the pumping process, and the funnel morphology parameters are calculated at each time point. Time series analysis is used to identify the evolution trend of the funnel morphology, including the rate of expansion of funnel size, the stability of shape changes, and the consistency of directional characteristics. The stability of the funnel morphology is judged based on the analysis of the rate of change of morphology parameters between consecutive time steps; when the rate of change is less than a preset threshold, the funnel morphology is considered to have reached a stable state.

[0085] In one optional implementation, a three-dimensional unsteady groundwater seepage control equation set for the hill aquifer is established based on the permeability coefficient field and the point source pumping boundary conditions. The seepage control equation set includes a gravity potential energy correction term reflecting the hill topographic undulations, comprising:

[0086] A three-dimensional spatial discrete grid is obtained for the aquifer in the hilly area. The three-dimensional spatial discrete grid includes the grid node coordinates in the horizontal direction and the grid node elevations in the vertical direction.

[0087] For each grid node in the three-dimensional discrete grid, the permeability coefficient value at the corresponding location of the grid node is extracted from the permeability coefficient field, and the pumping intensity value at the corresponding location of the grid node is extracted from the point source pumping boundary condition;

[0088] The terrain elevation gradient of a grid node is calculated based on the grid node elevation. The terrain elevation gradient is obtained by calculating the ratio of the elevation difference between the grid node and its adjacent grid nodes to the horizontal distance.

[0089] The gravitational potential energy correction coefficient for the grid node is calculated based on the terrain elevation gradient of the grid node. The gravitational potential energy correction coefficient is obtained by multiplying the terrain elevation gradient by the gravitational acceleration constant and normalizing it by combining the aquifer thickness at the grid node.

[0090] Based on the permeability coefficient, pumping intensity, and gravitational potential energy correction coefficient of each grid node in the three-dimensional spatial discrete grid, a set of three-dimensional unsteady groundwater seepage control equations for the hill aquifer is constructed.

[0091] The three-dimensional discrete mesh was obtained based on the geometric boundaries and geological structure characteristics of the hilly aquifer. A two-dimensional horizontal mesh covering the entire study area was established according to the horizontal range of the aquifer. The mesh cells adopted a square structure with a side length of ten meters. The node coordinates of the horizontal mesh were determined by the minimum bounding rectangle of the area boundary. Starting from the southwest corner, the east-west and north-south coordinates of each node were generated sequentially according to the mesh spacing. In the vertical direction, based on the elevation data of the top and bottom plates of the aquifer, the vertical space was divided into multiple layers of equal thickness, each with a thickness of two meters. The elevation of the vertical mesh nodes was calculated by subtracting the accumulated layer thickness from the elevation of the top plate of the aquifer, ensuring that the mesh nodes covered the entire thickness range of the aquifer.

[0092] The extraction of permeability coefficient values ​​at grid nodes is achieved using a three-dimensional interpolation method. The permeability coefficient field is stored as a three-dimensional data volume, containing the horizontal and vertical permeability coefficient values ​​for each spatial location. For any node in the three-dimensional grid, its spatial location is determined using its three-dimensional coordinates within the permeability coefficient field. When the grid node's location completely coincides with a sampling point in the permeability coefficient field, the corresponding permeability coefficient value is directly extracted. When the grid node's location lies between sampling points in the permeability coefficient field, a trilinear interpolation method is used to calculate the permeability coefficient value at that location. The trilinear interpolation process involves bilinear interpolation in the horizontal plane, followed by linear interpolation in the vertical direction, ultimately yielding the permeability coefficient value at the grid node.

[0093] The extraction of pumping intensity values ​​for grid nodes is based on the spatial distribution characteristics of point source boundary conditions. Point source pumping boundary conditions are defined as a set of discrete point sources, each containing parameters such as pumping location coordinates, pumping intensity, and pumping influence range. For each node in the 3D grid, the 3D Euclidean distance from that node to all point sources is calculated. When the distance between a grid node and a point source is less than the influence radius of that point source, the pumping intensity of that point source is assigned to that grid node. When a grid node is simultaneously within the influence range of multiple point sources, a distance-inverse weighted method is used to calculate the comprehensive pumping intensity, which is the sum of the products of the pumping intensity of each point source and the reciprocal of its distance to the grid node, divided by the sum of the reciprocals of those distances. When a grid node is not within the influence range of any point source, its pumping intensity value is set to zero.

[0094] The calculation of terrain elevation gradient is based on the spatial geometric relationship of grid nodes. For any node in the 3D grid, its four adjacent nodes in the horizontal plane are selected: the east, west, south, and north adjacent nodes. The elevation difference between this node and its east-west adjacent node is calculated and divided by the horizontal distance between the two nodes to obtain the east-west elevation gradient component. Similarly, the elevation difference between this node and its north-west adjacent node is calculated and divided by the horizontal distance between the two nodes to obtain the north-south elevation gradient component. The east-west and north-south elevation gradient components are vector-synthesized, i.e., the square root of the sum of the squares of the two components, to obtain the total terrain elevation gradient of the node. For nodes located at the grid boundary, the elevation gradient is calculated using the one-sided difference method, that is, only the elevation information of the inner adjacent nodes is used.

[0095] The calculation of the gravitational potential energy correction coefficient combines the terrain elevation gradient and physical parameters. The terrain elevation gradient value of each grid node is multiplied by the gravitational acceleration constant, which is taken as 9.8 m / s². The product represents the gravitational potential energy gradient at that node due to terrain inclination. For normalization, the aquifer thickness information at that grid node needs to be obtained. The aquifer thickness is calculated by the difference between the elevation of the top and bottom plates of the aquifer at that node. The gravitational potential energy gradient is divided by the aquifer thickness to obtain the gravitational potential energy correction intensity per unit thickness. To ensure dimensional consistency of the correction coefficient, the gravitational potential energy correction intensity per unit thickness is then divided by the reference permeability coefficient value. The reference permeability coefficient is selected as the average permeability coefficient within the study area, finally yielding the dimensionless gravitational potential energy correction coefficient.

[0096] The seepage control equations are constructed based on a three-dimensional extension of the principle of mass conservation and Darcy's law. For each node in the three-dimensional mesh, a water balance equation is established for that node. The left side of the water balance equation represents the rate of change of water storage at that node, calculated by multiplying the rate of change of water level at that node by the storage coefficient and the mesh cell volume. The right side of the equation contains three parts: a horizontal seepage term, a vertical seepage term, and a source-sink term.

[0097] The calculation of the horizontal seepage term considers the water exchange in the east-west and north-south directions. The east-west seepage flow is calculated by dividing the water level difference between the node and its eastern neighbor by the grid spacing, then multiplying by the horizontal permeability coefficient and the cross-sectional area. The cross-sectional area is the product of the grid spacing and the aquifer thickness. The calculation method for the north-south seepage flow is the same as that for the east-west direction. The total horizontal seepage term is the algebraic sum of the east-west and north-south seepage flows.

[0098] The calculation of the vertical seepage term reflects the water exchange between upper and lower grid nodes. The vertical seepage rate is calculated by dividing the water level difference between the node and the adjacent node in the upper layer by the vertical grid spacing, and then multiplying by the vertical permeability coefficient and the horizontal flow area. The horizontal flow area is the horizontal projected area of ​​the grid cell, i.e., the square of the grid spacing. For multi-layer aquifer structures, it is necessary to calculate the seepage exchange with the upper and lower layer nodes separately.

[0099] The calculation of the source-sink term comprises two components: pumping intensity and gravitational potential energy correction. The pumping intensity term directly uses the pumping intensity value of the grid node; a positive value indicates pumping, and a negative value indicates recharge. The gravitational potential energy correction term is calculated by multiplying the gravitational potential energy correction coefficient by the water level gradient. The calculation method for the water level gradient is the same as that for the topographic elevation gradient, but water level data is used instead of elevation data. The gravitational potential energy correction term reflects the impact of topographic relief on groundwater flow, exhibiting a significant effect in areas with steep terrain slopes.

[0100] Boundary conditions are handled using appropriate mathematical expressions for different types of boundaries. For constant water level boundaries, the water level at the boundary nodes is directly set to a given value, and the boundary equation simplifies to water level equal to the boundary water level value. For constant flow boundaries, the seepage flow at the boundary is set to a given value, and the boundary equation is expressed as seepage flow perpendicular to the boundary direction equal to the boundary flow value. For free boundaries, an iterative solution method is used to determine the boundary location and boundary conditions, and the water level at the boundary nodes is equal to the topographic elevation at that location.

[0101] The initial conditions are set based on the initial water level distribution of the aquifer. For a stable initial state, the initial water level of each grid node is obtained by interpolation using static water level observation data. The interpolation method uses the inverse distance weighting method, with the weight being the reciprocal of the distance from the observation point to the grid node. For an unstable initial state, the initial water level distribution is determined based on historical pumping records and water level change trends. The initial water level setting must satisfy water balance constraints to ensure that the seepage flow of each node is balanced under the initial state.

[0102] The numerical solution of the equations employs an implicit scheme using the finite difference method, discretizing the time domain with a time step of one hour. Within each time step, the seepage control equations constitute a linear system with a sparse coefficient matrix. The linear system is solved using the preconditioned conjugate gradient method, with the precondition matrix employing incomplete LU decomposition. The convergence criterion for iterative solutions is that the relative residual is less than 10 to the power of -6. For nonlinear problems, an outer iteration method is used, updating the nonlinear coefficients in each outer iteration until the solution converges.

[0103] In one optional implementation, based on the superimposed interference characteristics between the drawdown funnels of adjacent well points in the spatiotemporal evolution field of the water level, the hydraulic coupling region under the synergistic effect of the well group is identified, and the distribution of the water level gradient and the direction of flow field deflection within the hydraulic coupling region are determined, including:

[0104] Water level data at multiple time points at each spatial location are extracted from the spatiotemporal evolution field of water level, and the water level gradient vector at each spatial location is calculated;

[0105] Construct the hydraulic gradient vector field of the hill aquifer based on the water level gradient vector;

[0106] Calculate the divergence distribution of the hydraulic gradient vector field, construct the theoretical hydraulic gradient vector field under the condition of independent pumping of a single well, and calculate the theoretical divergence distribution of the theoretical hydraulic gradient vector field;

[0107] By comparing the divergence distribution with the theoretical divergence distribution, when the divergence distribution shows an abnormally positive value cluster or an abnormally negative value cluster and the region does not coincide with any well point, the region is determined to be a hydraulic gradient field disturbance region caused by the synergistic effect of the well group, and the hydraulic gradient field disturbance region is determined to be the hydraulic coupling region;

[0108] Within the hydraulic coupling region, the amplitude and direction distribution of the water level gradient vector are statistically analyzed to obtain the water level gradient distribution within the hydraulic coupling region;

[0109] Based on the water level gradient vector, a groundwater flow vector field is constructed within the hydraulic coupling region. By performing streamline tracing on the groundwater flow vector field, the deflection angle distribution of the groundwater streamlines relative to the radial flow direction of a single well is determined, and the dominant direction of the deflection angle distribution is determined as the flow field deflection direction.

[0110] The extraction of water level spatiotemporal evolution field data adopted a multi-time-node sampling method. Water level data was extracted at six-hour intervals starting from the moment pumping began in the well group, continuing for 72 hours. Spatially, a 5m x 5m regular grid system was used to discretize the hilly aquifer. For each grid node, the water level value at each time node was recorded, forming a spatiotemporal water level data matrix.

[0111] The water level gradient vector is calculated using the spatial difference method. For each grid node, the east-west gradient component is obtained by dividing the water level difference between the node's eastern and western adjacent nodes by twice the grid spacing. The north-south gradient component is obtained by dividing the water level difference between the node's northern and southern adjacent nodes by twice the grid spacing. The gradient components in both directions are combined to form the water level gradient vector for that node.

[0112] The hydraulic gradient vector field is constructed by spatially organizing the water level gradient vectors of all grid nodes. The hydraulic gradient vector field is represented using a two-dimensional vector field data structure, with each grid node corresponding to a vector element. The vector field is visualized using arrow graphics; the direction of the arrow indicates the gradient direction, and the length of the arrow indicates the gradient magnitude.

[0113] The divergence distribution is calculated using a numerical differential method. For each grid node in the hydraulic gradient vector field, the divergence value of the vector field at that node is calculated. The divergence is obtained by adding the east-west derivative of the east-west gradient component and the north-south derivative of the north-south gradient component at that node. The derivative is calculated using the central difference method, that is, the difference between corresponding components of adjacent nodes is divided by the grid spacing.

[0114] The theoretical hydraulic gradient vector field is constructed based on a single-well independent pumping theoretical model. For each well point, the theoretical drawdown at each grid node is calculated using the well function method, based on its pumping intensity and aquifer parameters. The calculation of the theoretical drawdown considers the distance from the well point to the grid node, pumping time, and the aquifer's hydraulic conductivity and storage coefficient. The theoretical drawdowns generated by all well points at each grid node are superimposed to obtain the theoretical water level distribution. Based on the theoretical water level distribution, the theoretical hydraulic gradient vector field is constructed using the same calculation method as the actual water level gradient.

[0115] The calculation process for the theoretical divergence distribution is the same as that for the actual divergence distribution. For each grid node in the theoretical hydraulic gradient vector field, the theoretical divergence value at that node is calculated. The theoretical divergence distribution reflects the theoretical characteristics of the hydraulic gradient field under independent pumping conditions of a single well.

[0116] Divergence distribution comparative analysis is achieved through difference calculation and anomaly identification. The difference between the actual and theoretical divergence distributions at each grid node is calculated to form a divergence difference distribution map. Anomaly identification employs a statistical threshold method, identifying regions where the absolute value of the divergence difference exceeds twice the standard deviation as anomalous areas. Aggregates of positive anomalous values ​​indicate abnormal divergence of the actual hydraulic gradient field in that region, while aggregates of negative anomalous values ​​indicate abnormal convergence.

[0117] The determination of the hydraulic gradient field disturbance zone is based on spatial analysis of the anomalous region. Spatial connectivity analysis is performed on the identified clusters of positive and negative anomalous values, merging spatially adjacent anomalous grid nodes into the same disturbance zone. The boundary of the disturbance zone is determined using the contour line method, with contour line values ​​set to twice the standard deviation of the divergence difference. For each disturbance zone, it is checked whether it coincides with any well point; the criterion for coincidence is that the shortest distance from the boundary of the disturbance zone to the well point is less than ten meters.

[0118] The determination of the hydraulic coupling region involves spatially merging all hydraulic gradient field disturbance regions that meet the conditions. Disturbance regions with a spatial distance of less than fifty meters are merged into a single hydraulic coupling region. The boundary of the hydraulic coupling region is determined using the convex hull algorithm to ensure the smoothness and integrity of the region boundary. Hydraulic coupling regions with an area less than one thousand square meters are considered pseudo-regions caused by calculation errors and are discarded.

[0119] The statistical analysis of water level gradient distribution was performed on all grid nodes within the hydraulic coupling region. The amplitude distribution statistics included the minimum, maximum, average, standard deviation, and quantile distribution of the gradient amplitude. The amplitude distribution was represented using a histogram, dividing the gradient amplitude range into twenty equal intervals and counting the number of nodes within each interval. The directional distribution statistics divided the gradient direction angles into thirty-six sector intervals at ten-degree intervals, and counted the number of nodes within each directional interval.

[0120] The groundwater flow vector field is constructed based on Darcy's law vector transformation. The water level gradient vector of each grid node within the hydraulic coupling region is multiplied by a negative one to obtain the groundwater flow direction vector at that node. The magnitude of the flow vector is determined by multiplying the gradient magnitude by the aquifer permeability coefficient.

[0121] Streamline tracing was implemented using a numerical integration method. Streamline starting points were established at the upstream boundary of the hydraulic coupling region, with a spacing of ten meters between starting points. Starting from each starting point, streamline tracing was performed according to the direction of the groundwater flow vector at that point. The tracing step size was set to one meter, and the tracing terminated when the streamline reached the region boundary or entered a region of extremely low flow velocity.

[0122] The deflection angle distribution is calculated based on the comparison between the streamline tangent direction and the radial flow direction. For each coordinate point on a streamline, the tangent vector and the radial vector pointing to the nearest well point are calculated. The tangent vector is obtained through the directional difference between adjacent coordinate points, and the radial vector is determined by the positional relationship between the coordinate point and the nearest well point. The deflection angle is the angle between the two direction vectors and is calculated using the vector dot product method.

[0123] The dominant direction is determined based on peak value analysis of the frequency distribution of deflection angles. The angle interval with the highest frequency is identified in the frequency distribution histogram, and the center angle of this interval is the dominant deflection angle. The precise calculation of the dominant direction uses a weighted average method, with the frequency of each interval as the weight to calculate the weighted average deflection angle. The flow field deflection direction is characterized by comparing the dominant deflection angle with a reference direction, which is the natural groundwater flow direction under conditions of no pumping disturbance.

[0124] In one optional implementation, a groundwater flow vector field within the hydraulic coupling region is constructed based on the water level gradient vector. By performing streamline tracing on the groundwater flow vector field, the deflection angle distribution of the groundwater streamlines relative to the radial flow direction of a single well is determined, and the dominant direction of the deflection angle distribution is determined as the flow field deflection direction, including:

[0125] Multiple initial seed points for streamline tracing are set in the groundwater flow vector field, and the initial seed points are uniformly distributed along the boundary of the hydraulic coupling region;

[0126] Starting from each of the initial seed points, streamline tracing is performed along the vector direction of the groundwater flow vector field to obtain the spatial distribution of groundwater streamlines within the hydraulic coupling region. The tangent direction vector at each point on the groundwater streamline is extracted, and the tangent direction vector represents the actual flow direction of groundwater at that point.

[0127] Calculate the location of the well point closest to the groundwater flow line, and determine the theoretical direction vector of single-well radial flow at each point on the groundwater flow line. The theoretical direction vector of single-well radial flow is the radial unit vector pointing from each point on the groundwater flow line to the well point.

[0128] Calculate the angle between the tangential direction vector and the theoretical radial flow direction vector of the single well, define the angle as the deflection angle of the groundwater streamline relative to the radial flow direction of the single well, and statistically analyze the deflection angles of all groundwater streamlines within the hydraulic coupling region to form the deflection angle distribution;

[0129] The dominant direction of the deflection angle distribution is determined as the flow field deflection direction.

[0130] The construction of the groundwater flow vector field is based on the spatial distribution characteristics of the water level gradient vector. After obtaining the water level gradient vectors of each grid node within the hydraulic coupling region, these vectors need to be converted into groundwater flow vectors. The conversion process follows the basic principle of Darcy's law, where the direction of groundwater flow is opposite to the direction of the water level gradient. Specifically, the water level gradient vector of each grid node is multiplied by a negative one to obtain the groundwater flow direction vector for that node. The magnitude of the groundwater flow vector is determined by multiplying the water level gradient value by the aquifer permeability coefficient. For anisotropic aquifers, the permeability coefficients in the horizontal and vertical directions need to be considered separately. The horizontal component of the water level gradient vector is multiplied by the horizontal permeability coefficient, and the vertical component is multiplied by the vertical permeability coefficient to form a complete groundwater flow vector field.

[0131] The initial seed points are set using a uniform distribution strategy along the upstream boundary of the hydraulic coupling region. A seed point distribution line is established with a spacing of five meters between points to ensure appropriate density for streamline tracing. The coordinates of the seed points are determined through a parameterized representation of the boundary line, which is parameterized according to the cumulative arc length, with a seed point placed every five meters. For irregularly shaped boundary lines, a piecewise linear approximation method is used, decomposing the boundary line into multiple straight line segments. Seed points are placed on each segment at equal intervals. The total number of seed points is determined based on the total length of the boundary line to ensure sufficient coverage of the entire upstream boundary.

[0132] Streamline tracing is implemented using a fourth-order Runge-Kutta numerical integration method. Starting from each initial seed point, it advances along the direction of the groundwater flow vector at that point. The integration step size is set to one meter to ensure a balance between tracing accuracy and computational efficiency. In each integration step, the groundwater flow vector at the current location is calculated, and the direction and position of the next step are determined based on the vector direction. Streamline tracing terminates under three conditions: the streamline reaches the downstream boundary of the hydraulic coupling region, the streamline enters an area with extremely low groundwater velocity, or the tracing length exceeds the preset maximum tracing distance. The maximum tracing distance is set to twice the diagonal length of the hydraulic coupling region to prevent streamline tracing from entering an infinite loop.

[0133] The spatial distribution of groundwater streamlines is obtained by recording all coordinate points during streamline tracing. For each streamline, coordinate points along its path are recorded according to the integration step size, forming a discretized representation of the streamline. Streamline coordinate points are stored using a linked list data structure, supporting dynamic addition of coordinate points. For areas with excessively high streamline density, adaptive sparsity processing is employed to retain key streamline characteristics while reducing computational burden. Streamline quality is verified by calculating the smoothness and continuity of the streamlines, removing abnormal streamlines with obvious jumps or reversals.

[0134] The extraction of the tangent direction vector is based on the geometric analysis of discrete points on the streamline. For each coordinate point on the streamline, one adjacent point before and one point after it are selected to form a three-point combination. The tangent direction vector at that coordinate point is obtained by calculating the direction vector from the preceding point to the following point. The tangent direction vector is calculated using the central difference method, i.e., subtracting the coordinates of the preceding point from the coordinates of the following point, and then dividing by the distance between the two points, yields the unit tangent direction vector. For the starting and ending points of the streamline, the tangent direction vector is calculated using forward difference and backward difference methods respectively. The direction angle of the tangent direction vector is calculated using the arctangent function, with an angle range of 0 to 360 degrees.

[0135] The location of the nearest wellpoint is determined through spatial distance calculation. For each coordinate point on the streamline, the Euclidean distance from that point to all wellpoints within the hydraulic coupling region is calculated. The distance calculation employs the square root method of the sum of squares of the two-dimensional coordinate differences: the square of the east-west coordinate difference between the coordinate point and the wellpoint, plus the square of the north-south coordinate difference, is taken as the square root. The wellpoint with the smallest distance is selected as the nearest wellpoint to that coordinate point. To improve computational efficiency, a spatial indexing method is used to pre-establish a spatial distribution index of the wellpoints, and a quadtree data structure is used to accelerate the nearest neighbor search process.

[0136] The determination of the theoretical direction vector for radial flow in a single well is based on geometric analysis. For each coordinate point on the streamline, the direction vector pointing from that point to the nearest well point is calculated. The direction vector is calculated by subtracting the coordinate point's coordinates from the coordinate point's coordinates, resulting in the vector pointing from the coordinate point to the well point. This vector is then normalized by dividing each component of the vector by its magnitude, yielding the unit direction vector. The theoretical direction vector for radial flow in a single well represents the theoretical flow direction of groundwater from that coordinate point to the well point under single-well pumping conditions, and this direction always points towards the center of the well point.

[0137] The deflection angle is calculated through vector angle analysis. For each coordinate point on the streamline, the angle between the tangent vector at that point and the theoretical radial flow direction vector of the single well is calculated. The angle calculation uses the vector dot product method, where the dot product of two unit vectors equals the cosine of the angle. The radian value of the angle is calculated using the inverse cosine function and then converted to an angle value. The deflection angle is defined as the clockwise deflection angle of the tangent vector relative to the radial flow direction vector, ranging from zero to 360 degrees. When the deflection angle is zero, it indicates that the groundwater flow direction is consistent with the radial flow direction of the single well; when the deflection angle is 90 degrees, it indicates that the groundwater flow direction is perpendicular to the radial flow direction.

[0138] The distribution of deflection angles was statistically analyzed using frequency analysis. The angle range from 0 to 360 degrees was divided into 36 equally spaced intervals, each spanning 10 degrees. The deflection angles at all coordinate points on all streamlines were statistically analyzed, and the frequency of occurrence of each angle interval was calculated. A frequency distribution histogram of deflection angles was generated, with the horizontal axis representing the angle interval and the vertical axis representing the frequency of occurrence of the deflection angle within that interval. The frequency distribution was smoothed using a moving average method to reduce the impact of statistical noise.

[0139] The dominant direction is determined based on peak value analysis of the frequency distribution. In the histogram of deflection angle frequency distribution, the angle interval with the highest frequency is identified, and the center angle of this interval is the dominant direction of the deflection angle distribution. Determining the dominant direction also requires considering the frequency distribution of adjacent intervals, using a weighted average method to calculate a more accurate dominant direction angle. The weights are set based on the frequency value of each interval; the higher the frequency of the interval, the greater the weight. The dominant direction angle is calculated by summing the products of the center angle of each interval and its corresponding weight, divided by the total sum of the weights.

[0140] The direction of flow field deflection is characterized by comparing the dominant direction angle with the reference direction. The natural groundwater flow direction within the hydraulic coupling region is selected as the reference direction, which is determined by analyzing the groundwater flow field under conditions without pumping interference. The direction of flow field deflection is defined as the difference between the dominant direction angle and the reference direction angle. A positive difference indicates a clockwise deflection of the groundwater flow field relative to the natural flow direction; a negative difference indicates a counterclockwise deflection. The magnitude of the deflection angle reflects the intensity of the disturbance to the groundwater flow field caused by well pumping; a larger deflection angle indicates a stronger disturbance.

[0141] In one optional implementation, the pumping intensity of each well point in the well group is coordinated and controlled to generate an evaluation index function that satisfies the target precipitation depth requirement and optimizes the overall uniformity of water level decline in the hilly aquifer, including:

[0142] The water level gradient distribution is divided into multiple water level gradient grade intervals. The well point with the largest area-gradient weighted index in the water level gradient grade interval is determined as the key well point for gradient control. The spatial position deflection angle of each key well point for gradient control relative to the reference direction is calculated using the flow field deflection direction as the reference direction.

[0143] Based on the spatial location deflection angle, the key well points for gradient control are divided into a forward well point group and a reverse well point group. The well points in the forward well point group are located on the forward extension path of the flow field deflection direction, and the well points in the reverse well point group are located in the reverse hindrance region of the flow field deflection direction. For the forward well point group, a pumping intensity enhancement coefficient is set, and for the reverse well point group, a pumping intensity attenuation coefficient is set.

[0144] Based on the pumping intensity enhancement coefficient and the pumping intensity attenuation coefficient, the pumping intensity values ​​of each well point in the forward well point group and the reverse well point group are adjusted; based on the adjusted pumping intensity values ​​of each well point, the spatial distribution of water level contour lines within the hilly aquifer is extracted;

[0145] For the spatial distribution of the water level contour lines, the standard deviation of the curvature of the water level contour lines and the spatial uniformity of the density distribution of the water level contour lines are calculated. The standard deviation of curvature and the spatial uniformity are weighted and combined to construct an evaluation index function for the overall uniformity of water level decline in the hilly aquifer.

[0146] The implementation process of well cluster coordinated regulation begins with the precise analysis of water level gradient distribution. After obtaining the spatiotemporal evolution field of water level in the hilly aquifer, water level distribution data under stable conditions is selected, and a regular grid system is established to spatially discretize the study area. The grid size is set to 5 meters by 5 meters to ensure a balance between computational accuracy and efficiency. For each grid node, its water level gradient components in the east-west and north-south directions are calculated. Specifically, the water level difference between this node and its adjacent eastern node is divided by the grid spacing, and the water level difference between this node and its adjacent northern node is divided by the grid spacing. The gradient components in the two directions are then vector-synthesized to obtain the total water level gradient value of the node.

[0147] The water level gradient classification intervals were determined using statistical analysis. Water level gradient values ​​from all grid nodes within the study area were collected, sorted by value, and the minimum, 20th percentile, 40th percentile, 60th percentile, 80th percentile, and maximum values ​​were calculated. Based on this, five gradient classification intervals were established: the very low gradient interval covers the minimum to the 20th percentile; the low gradient interval covers the 20th to 40th percentile; the medium gradient interval covers the 40th to 60th percentile; the high gradient interval covers the 60th to 80th percentile; and the very high gradient interval covers the 80th percentile to the maximum value.

[0148] The area-gradient weighted index is calculated individually for each well point. A circular influence zone is established with the well point location as the center and half the radius of the well point's pumping influence as the radius. This circular zone is overlaid with the grid system to identify grid cells completely within the circular zone and those partially within it. For grid cells completely within the circular zone, their contribution area is the entire grid area; for grid cells partially within the circular zone, their specific area within the circular zone is determined through geometric calculations. The contribution area of ​​each grid cell is multiplied by its water level gradient value, and the product results of all relevant grid cells are summed to obtain the area-gradient weighted index for that well point.

[0149] The key well points for gradient control are determined through a hierarchical comparison method. Within each gradient level interval, the area-gradient weighted index values ​​of all well points are compared, and the well point with the largest value is selected as the key well point for gradient control in that interval. If no well points are distributed within a certain gradient level interval, no key well point is set for that interval. Using this method, a maximum of five key well points for gradient control can be obtained, each representing the control focus at different gradient levels.

[0150] The direction of flow field deflection is determined based on streamline tracing analysis. Uniformly distributed streamline starting points are set at the upstream boundary of the hydraulic coupling region, with a spacing of ten meters between them. Starting from each point, the streamline path is traced along the hydraulic gradient direction until the streamline reaches the downstream boundary of the hydraulic coupling region or flows out of the study area. The path coordinates of each streamline within the hydraulic coupling region are recorded, and the average flow direction angle of the streamline is calculated. A weighted average of all streamline average flow direction angles is then calculated, with the weight being the path length of each streamline within the hydraulic coupling region, to obtain the overall dominant flow direction angle. This dominant flow direction angle is the reference direction.

[0151] The calculation of the spatial position deflection angle is based on coordinate transformation. A rotating coordinate system is established, using the reference direction as the true north direction of the new coordinate system. The original coordinates of each key wellpoint for gradient control are transformed into the new coordinate system through coordinate rotation. The angle between the wellpoint's position vector in the new coordinate system and the true north direction vector is calculated; this angle is the spatial position deflection angle. The angle is calculated using the arctangent function method, obtaining the angle value from the wellpoint's horizontal and vertical coordinates in the new coordinate system. The angle range is from -180 degrees to +180 degrees.

[0152] The division between forward and reverse well point groups is based on the numerical range of spatial location deflection angles. Key gradient control well points with spatial location deflection angles ranging from -45° to +45° are classified as forward well point groups. These well points are located on the forward extension path of the reference direction and can effectively guide groundwater flow in the expected direction. Key gradient control well points with spatial location deflection angles ranging from -180° to -135° and from +135° to +180° are classified as reverse well point groups. These well points are located in the reverse hindrance region of the reference direction and have a hindrance effect on the groundwater flow field.

[0153] The pumping intensity enhancement coefficient is set using a piecewise linear function method. For well points in the forward well point group, the enhancement coefficient is determined based on the absolute value of their spatial deflection angle. When the absolute value of the deflection angle is zero degrees, the enhancement coefficient is set to 2.0; when the absolute value of the deflection angle is forty-five degrees, the enhancement coefficient is set to 1.2. Within the range of zero to forty-five degrees, the enhancement coefficient varies linearly, and the specific value is calculated through linear interpolation. The interpolation calculation method is that the enhancement coefficient equals 2.0 minus the result of multiplying the absolute value of the deflection angle by 0.0178.

[0154] The pumping intensity attenuation coefficient is also set using a piecewise linear function method. For well points in the reverse well point group, the attenuation coefficient is determined based on the degree of deviation of their spatial position angle from 180 degrees or -180 degrees. When the deviation angle is ±180 degrees, the attenuation coefficient is set to 0.5; when the deviation angle is ±135 degrees, the attenuation coefficient is set to 0.9. Within the corresponding angle range, the attenuation coefficient changes linearly, and the specific value is calculated through linear interpolation.

[0155] The pumping intensity of wellpoints is adjusted through a coefficient product method. For each wellpoint in the forward wellpoint group, its original pumping intensity value is multiplied by the corresponding enhancement coefficient to obtain the adjusted pumping intensity. For each wellpoint in the reverse wellpoint group, its original pumping intensity value is multiplied by the corresponding attenuation coefficient to obtain the adjusted pumping intensity. For gradient control key wellpoints that do not belong to either the forward or reverse wellpoint groups, as well as all non-key wellpoints, their original pumping intensity values ​​are kept unchanged.

[0156] The extraction of the spatial distribution of water level contour lines is based on recalculating the water level field using the updated wellpoint pumping intensity. The adjusted pumping intensity of each wellpoint is used as the new boundary condition to re-establish and solve the three-dimensional unsteady groundwater seepage control equations. Steady-state water level distribution data from the solution results are selected, and contour lines at different water level elevations are extracted using a contour line tracing algorithm. The extraction interval for contour lines is set to 0.5 meters; starting from the highest water level, one contour line is extracted every 0.5 meters until the lowest water level. Contour line tracing employs a linear interpolation method, determining the precise location of the contour lines through interpolation between adjacent grid nodes.

[0157] The standard deviation of the curvature of water level contour lines is calculated using a geometric analysis method. Each contour line is sampled at equal intervals (two meters apart) to obtain the coordinates of discrete points on the contour lines. For each discrete point, one sampling point before and one before it are selected, forming a three-point combination. An arc is determined using these three points, and the curvature value of this arc is calculated; the curvature value is equal to the reciprocal of the arc's radius. The curvature values ​​of all sampling points on all contour lines are collected, and the mean and standard deviation of these curvature values ​​are calculated. The standard deviation of curvature reflects the irregularity of the contour line shape; a larger value indicates a more irregular shape.

[0158] The spatial uniformity of the water level contour density distribution was assessed using a grid statistical method. The study area was divided into regular statistical grids of 100 meters by 100 meters, and the number of contour lines crossing each grid was counted. The number of contour lines was determined using a line segment intersection method, checking whether a contour line intersects the boundary of the statistical grid; if an intersection occurs, it is counted in the grid's contour count. The mean and standard deviation of the contour line counts for all statistical grids were calculated, and the spatial uniformity was equal to the mean divided by the standard deviation. A higher spatial uniformity value indicates a more even distribution of contour lines in space.

[0159] The evaluation index function is constructed using a weighted combination method. The reciprocal of the standard deviation of curvature is taken and normalized. The normalization method is to subtract the minimum value among all candidate schemes from the index and then divide by the difference between the maximum and minimum values. Spatial uniformity is directly normalized in the same way. The reciprocal of the normalized standard deviation of curvature and the normalized spatial uniformity are then weighted and summed according to a weight ratio of 0.6 to 0.4 to obtain the final evaluation index function value. The value of this function ranges from zero to one; a larger value indicates better uniformity of the overall water level drop in the hilly aquifer.

[0160] In one optional implementation, the method further includes:

[0161] The well group coordinated pumping scheme is simulated and verified in the digital twin model to obtain the deviation between the actual precipitation effect and the expected target. The deviation is then fed back to the digital twin model to update the correction coefficient of the formation permeability distribution characterization and hydraulic coupling relationship, thereby achieving dynamic alignment between the digital twin model and the real hydrological response characteristics of the hill aquifer.

[0162] When simulating and verifying the well group coordinated pumping scheme, the established scheme is imported into a digital twin model. This scheme includes key parameters such as the specific location coordinates of multiple wells, pumping flow rate, and pumping sequence. For example, in a hilly aquifer control project in a mining area, 10 pumping wells were installed, each 150 meters deep and 300 millimeters in diameter. The pumping flow rate of each well was set between 50 and 120 cubic meters per hour, with continuous 24-hour operation and maintenance every 7 days.

[0163] After importing the pumping scheme, simulation calculations were initiated in the digital twin model to simulate the dynamic changes in groundwater level over time. During the simulation, the system calculated the changes in hydraulic head at each grid node based on pre-established hydrogeological parameter distributions and boundary conditions. The simulation time step was set to 1 day, with a total simulation period of 180 days to fully reflect the changing trends in the scope and extent of the pumping's impact.

[0164] After the simulation is completed, the water level change data of key monitoring points are automatically extracted and compared with the expected precipitation target. The expected precipitation target is usually determined by engineering safety requirements, such as requiring the water level of the aquifer above the working face to drop to 30 meters below the mining elevation during mining operations. By comparing the actual precipitation effect with the expected target, the deviation is calculated. The deviation is calculated using the absolute difference between the actual water level and the target water level. For example, if at monitoring point A, the simulated water level after 180 days is -25 meters (relative to the mining elevation), while the expected target is -30 meters, then the deviation at that point is 5 meters.

[0165] The system collects deviation data from all key monitoring points to create a deviation distribution map. Based on spatial distribution characteristics, the system automatically identifies abnormal areas based on the deviation data for different regions. For example, in the northern part of a mining area, the deviation is generally between 3 and 7 meters, significantly higher than the 0-2 meters in other areas, indicating a large error in the geological parameters of that region.

[0166] An adaptive correction mechanism for model parameters is initiated in response to detected deviations. This mechanism first analyzes the spatial distribution characteristics of the deviations to identify parameter regions requiring focused correction. For example, for the larger deviations in the aforementioned northern region, the system will focus on correcting the permeability coefficient in that region. The correction process employs an inversion algorithm, performing multiple iterative calculations to adjust the correction coefficients representing the formation permeability distribution and its hydraulic coupling relationship until the error between the model output and the actual observation data is reduced to an acceptable range (generally controlled within 2 meters).

[0167] During parameter calibration, the sensitivity of each parameter is evaluated, and parameters with a greater impact on the model results are adjusted first. For example, in a hilly aquifer model, the calibration range for the permeability coefficient is set to 0.5-2.0 times the original value, and the calibration range for the storage coefficient is 0.8-1.2 times the original value. Through multiple rounds of iterative calculations, the system determines that the permeability coefficient in the northern region should be reduced to 0.7 times the original value, while the storage coefficient remains unchanged, so that the model calculation results are closest to the actual observation data.

[0168] After calibration, the simulation calculation is re-executed using the updated parameters to verify the calibration effect. If the deviation still exceeds the acceptable range (e.g., the deviation at some points is greater than 3 meters), parameter calibration continues. If the deviation at all monitoring points is within the acceptable range, the model calibration is considered successful, and the dynamic alignment of the digital twin model with the actual hydrological response characteristics of the hill aquifer is completed.

[0169] To ensure long-term operational reliability, actual monitoring data is compared with model predictions periodically (e.g., every 30 days), and model parameters are continuously updated. This dynamic correction mechanism can adapt to changes in formation parameters over time, such as formation compaction effects caused by long-term pumping and seasonal rainfall variations.

[0170] In a practical application case at a mining area, the initial digital twin model predicted an average water level drop of 28 meters after 180 days, while actual monitoring data showed an average drop of only 23 meters, a deviation of 5 meters. After two rounds of parameter correction, the model recalculated an average water level drop of 22.5 meters, reducing the deviation from the actual monitoring data to within 0.5 meters, meeting the engineering accuracy requirements. The adjusted model successfully predicted the water level change trend for the following 90 days, with an average prediction error controlled within 1.2 meters, providing a reliable basis for further optimization of the well group operation plan.

[0171] Through the aforementioned dynamic alignment process, the digital twin model can continuously absorb information contained in the actual monitoring data, improve its own prediction accuracy, provide a scientific basis for the continuous optimization and adjustment of the well group collaborative pumping scheme, and ultimately achieve precise control of groundwater in the hilly aquifer.

[0172] A second aspect of the present invention provides a digital twin-based collaborative precipitation simulation system for well groups in hilly aquifers, comprising:

[0173] The first unit is used to obtain the geological parameters of the hill aquifer, well layout information, and target precipitation depth requirements, and to construct a digital twin model of the hill aquifer. The digital twin model includes a representation of the formation permeability distribution, point source pumping boundary conditions, and the hydraulic coupling relationship between the well group and the aquifer.

[0174] The second unit is used to perform virtual precipitation simulation of the well group based on the digital twin model. By solving the three-dimensional unsteady groundwater seepage control equations of the hill aquifer, the spatiotemporal evolution field of water level in the hill aquifer is obtained, and the pumping influence radius and water level drawdown funnel shape of each well point are extracted.

[0175] The third unit is used to identify the hydraulic coupling region under the synergistic effect of the well group based on the superimposed interference characteristics between the drawdown funnels of adjacent well points in the spatiotemporal evolution field of the water level, and to determine the water level gradient distribution and flow field deflection direction within the hydraulic coupling region;

[0176] The fourth unit is used to coordinate and regulate the pumping intensity of each well point in the well group based on the water level gradient distribution and flow field deflection direction in the hydraulic coupling region, generate an evaluation index function that meets the target drawdown depth requirement and optimizes the overall water level drop uniformity of the hill aquifer, and solve the well group coordinated pumping scheme with the minimization of the evaluation index function as the optimization objective and the drawdown constraint as the boundary condition.

[0177] A third aspect of the present invention provides an electronic device, comprising:

[0178] processor;

[0179] Memory used to store processor-executable instructions;

[0180] The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0181] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0182] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0183] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A simulation method for coordinated precipitation in a well group of hilly aquifers based on digital twins, characterized in that, include: Geological parameters, well layout information, and target precipitation depth requirements of the hill aquifer are obtained. A digital twin model of the hill aquifer is then constructed. The digital twin model includes a representation of the formation permeability distribution, point source pumping boundary conditions, and the hydraulic coupling relationship between the well group and the aquifer. Based on the digital twin model, virtual precipitation simulation was performed on the well group. By solving the three-dimensional unsteady groundwater seepage control equations of the hill aquifer, the spatiotemporal evolution field of water level within the hill aquifer was obtained, and the pumping influence radius and water level drawdown cone shape of each well point were extracted. Based on the superimposed interference characteristics between the drawdown cones of adjacent well points in the spatiotemporal evolution field of water level, the hydraulic coupling region under the synergistic effect of the well group is identified, and the water level gradient distribution and flow field deflection direction within the hydraulic coupling region are determined, including: Water level data at multiple time points at each spatial location are extracted from the spatiotemporal evolution field of water level, and the water level gradient vector at each spatial location is calculated; Construct the hydraulic gradient vector field of the hill aquifer based on the water level gradient vector; Calculate the divergence distribution of the hydraulic gradient vector field, construct the theoretical hydraulic gradient vector field under the condition of independent pumping of a single well, and calculate the theoretical divergence distribution of the theoretical hydraulic gradient vector field; By comparing the divergence distribution with the theoretical divergence distribution, when the divergence distribution shows an abnormally positive value cluster or an abnormally negative value cluster and the region does not coincide with any well point, the region is determined to be a hydraulic gradient field disturbance region caused by the synergistic effect of the well group, and the hydraulic gradient field disturbance region is determined to be the hydraulic coupling region; Within the hydraulic coupling region, the amplitude and direction distribution of the water level gradient vector are statistically analyzed to obtain the water level gradient distribution within the hydraulic coupling region; Based on the water level gradient vector, a groundwater flow vector field is constructed within the hydraulic coupling region. By performing streamline tracing on the groundwater flow vector field, the deflection angle distribution of the groundwater streamlines relative to the radial flow direction of a single well is determined, and the dominant direction of the deflection angle distribution is determined as the flow field deflection direction. Based on the water level gradient distribution and flow field deflection direction in the hydraulic coupling region, the pumping intensity of each well point in the well group is coordinated and controlled to generate an evaluation index function that meets the target drawdown depth requirement and optimizes the overall water level drop uniformity of the hill aquifer. The optimization objective is to minimize the evaluation index function, and the well group coordinated pumping scheme is solved with the drawdown constraint as the boundary condition.

2. The method according to claim 1, characterized in that, Based on the digital twin model, virtual precipitation simulation was performed on the well group. By solving the three-dimensional unsteady groundwater seepage control equations of the hill aquifer, the spatiotemporal evolution field of the water level within the hill aquifer was obtained. The pumping influence radius and water level drawdown cone morphology of each well point were extracted, including: The formation permeability distribution representation in the digital twin model is converted into a spatially discretized permeability coefficient field, and the well group layout information is converted into point source pumping boundary conditions; Based on the permeability coefficient field and the point source pumping boundary conditions, a three-dimensional unsteady groundwater seepage control equation set for the hill aquifer is established. The seepage control equation set introduces a gravity potential energy correction term that reflects the topographic undulation characteristics of the hill to characterize the influence of topographic elevation difference on the direction and velocity distribution of groundwater flow. Solving the seepage control equations in the spatiotemporal domain yields the water level values ​​of each spatial node within the hill aquifer over a continuous time series, thus forming the spatiotemporal evolution field of the water level. Extract the spatial boundaries around each well point that satisfy the preset drop criterion from the spatiotemporal evolution field of the water level, and determine the pumping influence radius of each well point. The geometric features of the water level contour lines within the radius of influence of the pumping are extracted as the shape of the water level drawdown funnel. The shape of the water level drawdown funnel includes the offset of the funnel center position, the direction angle of the funnel's major axis and minor axis, and the asymmetry index of the funnel boundary.

3. The method according to claim 2, characterized in that, Based on the permeability coefficient field and the point source pumping boundary conditions, a three-dimensional unsteady groundwater seepage control equation set for the hill aquifer is established. The seepage control equation set includes a gravitational potential energy correction term reflecting the hill topographic undulation characteristics, which includes: A three-dimensional spatial discrete grid is obtained for the aquifer in the hilly area. The three-dimensional spatial discrete grid includes the grid node coordinates in the horizontal direction and the grid node elevations in the vertical direction. For each grid node in the three-dimensional discrete grid, the permeability coefficient value at the corresponding location of the grid node is extracted from the permeability coefficient field, and the pumping intensity value at the corresponding location of the grid node is extracted from the point source pumping boundary condition; The terrain elevation gradient of a grid node is calculated based on the grid node elevation. The terrain elevation gradient is obtained by calculating the ratio of the elevation difference between the grid node and its adjacent grid nodes to the horizontal distance. The gravitational potential energy correction coefficient for the grid node is calculated based on the terrain elevation gradient of the grid node. The gravitational potential energy correction coefficient is obtained by multiplying the terrain elevation gradient by the gravitational acceleration constant and normalizing it by combining the aquifer thickness at the grid node. Based on the permeability coefficient, pumping intensity, and gravitational potential energy correction coefficient of each grid node in the three-dimensional spatial discrete grid, a set of three-dimensional unsteady groundwater seepage control equations for the hill aquifer is constructed.

4. The method according to claim 1, characterized in that, Based on the water level gradient vector, a groundwater flow vector field is constructed within the hydraulic coupling region. By performing streamline tracing on the groundwater flow vector field, the deflection angle distribution of the groundwater streamlines relative to the radial flow direction of a single well is determined, and the dominant direction of the deflection angle distribution is determined as the flow field deflection direction, including: Multiple initial seed points for streamline tracing are set in the groundwater flow vector field, and the initial seed points are uniformly distributed along the boundary of the hydraulic coupling region; Starting from each of the initial seed points, streamline tracing is performed along the vector direction of the groundwater flow vector field to obtain the spatial distribution of groundwater streamlines within the hydraulic coupling region. The tangent direction vector at each point on the groundwater streamline is extracted, and the tangent direction vector represents the actual flow direction of groundwater at that point. Calculate the location of the well point closest to the groundwater flow line, and determine the theoretical direction vector of single-well radial flow at each point on the groundwater flow line. The theoretical direction vector of single-well radial flow is the radial unit vector pointing from each point on the groundwater flow line to the well point. Calculate the angle between the tangential direction vector and the theoretical radial flow direction vector of the single well, define the angle as the deflection angle of the groundwater streamline relative to the radial flow direction of the single well, and statistically analyze the deflection angles of all groundwater streamlines within the hydraulic coupling region to form the deflection angle distribution; The dominant direction of the deflection angle distribution is determined as the flow field deflection direction.

5. The method according to claim 1, characterized in that, The evaluation index function for coordinating and controlling the pumping intensity of each well point in the well group to generate an evaluation index function that satisfies the target precipitation depth requirement and optimizes the overall uniformity of water level decline in the hill aquifer includes: The water level gradient distribution is divided into multiple water level gradient grade intervals. The well point with the largest area-gradient weighted index in the water level gradient grade interval is determined as the key well point for gradient control. The spatial position deflection angle of each key well point for gradient control relative to the reference direction is calculated using the flow field deflection direction as the reference direction. Based on the spatial location deflection angle, the key well points for gradient control are divided into a forward well point group and a reverse well point group. The well points in the forward well point group are located on the forward extension path of the flow field deflection direction, and the well points in the reverse well point group are located in the reverse hindrance region of the flow field deflection direction. For the forward well point group, a pumping intensity enhancement coefficient is set, and for the reverse well point group, a pumping intensity attenuation coefficient is set. Based on the pumping intensity enhancement coefficient and the pumping intensity attenuation coefficient, the pumping intensity values ​​of each well point in the forward well point group and the reverse well point group are adjusted; based on the adjusted pumping intensity values ​​of each well point, the spatial distribution of water level contour lines within the hilly aquifer is extracted; For the spatial distribution of the water level contour lines, the standard deviation of the curvature of the water level contour lines and the spatial uniformity of the density distribution of the water level contour lines are calculated. The standard deviation of curvature and the spatial uniformity are weighted and combined to construct an evaluation index function for the overall uniformity of water level decline in the hilly aquifer.

6. The method according to claim 1, characterized in that, The method further includes: The well group coordinated pumping scheme is simulated and verified in the digital twin model to obtain the deviation between the actual precipitation effect and the expected target. The deviation is then fed back to the digital twin model to update the correction coefficient of the formation permeability distribution characterization and hydraulic coupling relationship, thereby achieving dynamic alignment between the digital twin model and the real hydrological response characteristics of the hill aquifer.

7. A digital twin-based collaborative precipitation simulation system for well groups in hilly aquifers, used to implement the method as described in any one of claims 1-6, characterized in that, include: The first unit is used to obtain the geological parameters of the hill aquifer, well layout information, and target precipitation depth requirements, and to construct a digital twin model of the hill aquifer. The digital twin model includes a representation of the formation permeability distribution, point source pumping boundary conditions, and the hydraulic coupling relationship between the well group and the aquifer. The second unit is used to perform virtual precipitation simulation of the well group based on the digital twin model. By solving the three-dimensional unsteady groundwater seepage control equations of the hill aquifer, the spatiotemporal evolution field of water level in the hill aquifer is obtained, and the pumping influence radius and water level drawdown funnel shape of each well point are extracted. The third unit is used to identify the hydraulic coupling region under the synergistic effect of the well group based on the superimposed interference characteristics between the drawdown funnels of adjacent well points in the spatiotemporal evolution field of the water level, and to determine the water level gradient distribution and flow field deflection direction within the hydraulic coupling region; The fourth unit is used to coordinate and regulate the pumping intensity of each well point in the well group based on the water level gradient distribution and flow field deflection direction in the hydraulic coupling region, generate an evaluation index function that meets the target drawdown depth requirement and optimizes the overall water level drop uniformity of the hill aquifer, and solve the well group coordinated pumping scheme with the minimization of the evaluation index function as the optimization objective and the drawdown constraint as the boundary condition.

8. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 6.

9. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 6.

Images