Regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twinning

Abstract

The application discloses a regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twinning, relates to the technical field of hydrogeological numerical simulation and groundwater resource management, and comprises a multi-source data access module, a three-dimensional hydrogeological modeling module, a flow field simulation calculation module, a visual display module, a virtual-real deviation evaluation module and a hierarchical correction execution module. The hierarchical correction execution module is bidirectionally connected with the rest modules, corresponding correction strategies are executed according to deviation levels, and a closed-loop operation mechanism of monitoring input, simulation deduction, deviation evaluation, automatic correction and verification backflow updating is formed. The application adds the virtual-real deviation evaluation and hierarchical correction module, changes the operation mode of the traditional static model which relies on artificial periodic calibration, matches the correction strategies of corresponding depths according to deviation levels, avoids excessive correction under small deviation, can deal with complex situations of flow field mutation, and can long-term maintain the simulation accuracy of the flow field.

Description

Technical Field

[0001] This invention relates to the field of hydrogeological numerical simulation and groundwater resource management technology, specifically a regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins. Background Technology

[0002] Regional groundwater flow field simulation is a routine technical means in groundwater resource assessment, over-extraction area management, and geological environment risk prevention and control. With the popularization of digital technology, flow field simulation systems that integrate three-dimensional geological modeling and numerical calculation are being used more and more widely. They can intuitively display the spatiotemporal distribution and evolution of groundwater elements such as water level and flow velocity, providing data reference for daily groundwater management.

[0003] Existing groundwater numerical simulation systems mostly employ a static modeling approach, where the model structure and parameters remain largely unchanged after a one-time calibration. In actual operation, influenced by various factors such as groundwater extraction, climate fluctuations, and the long-term evolution of geological conditions, simulation accuracy gradually declines over time. This necessitates periodic manual parameter inversion and model calibration, which is not only time-consuming but also incurs high overall maintenance costs. Furthermore, existing simulation schemes incorporating digital twins are mostly in open-loop operation mode, primarily using monitoring data as simulation input. They cannot automatically adjust model parameters and structure based on simulation deviations, making it difficult to maintain long-term synchronization between the model and the actual groundwater system. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of existing technologies and provide a dynamic evolution simulation system for three-dimensional regional groundwater flow fields based on digital twins. This system, by adding a virtual-real deviation assessment module and a graded correction execution module to the groundwater flow field twin simulation system, constructs a closed-loop operation mechanism encompassing monitoring input, simulation deduction, deviation assessment, automatic correction, and verification backflow updates. This changes the traditional static model's reliance on periodic manual calibration. The system matches correction strategies to corresponding depths according to deviation levels: minor deviations adjust local hydraulic parameters, moderate deviations adjust source-sink boundaries and hydraulic connections, and severe deviations adjust the model topology and computational grid. This avoids over-correction when deviations are small and can handle abrupt changes in flow field structure. It eliminates the need for frequent manual full parameter inversion and can maintain the accuracy of flow field simulation over the long term.

[0005] To address the aforementioned technical problems, this invention provides the following technical solution: a regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins. This system includes a multi-source data access module, a three-dimensional hydrogeological modeling module, a flow field simulation calculation module, a visualization module, a virtual-real deviation assessment module, and a graded correction execution module. The multi-source data access module is used to access groundwater monitoring, hydrological and meteorological data, and water intake data from human activities. The accessed data covers various monitoring indicators such as water level, water temperature, rainfall, evaporation, and extraction volume. The data sampling frequency ranges from hourly to daily. After conversion to a unified spatiotemporal reference and preprocessing, the data is synchronously input to the flow field simulation calculation module and the virtual-real deviation evaluation module. The preprocessing process includes data cleaning, missing value interpolation, outlier removal, and time-series smoothing to ensure the continuity and reliability of the input data. The 3D hydrogeological modeling module is used to construct a virtual twin base model of the aquifer structure in the study area. The model is constructed based on borehole data, hydrogeological profiles and pumping test results in the study area. It divides different hydrogeological units and assigns corresponding parameters, providing the geometric boundaries and hydraulic parameters required for flow field simulation. The initial parameter values ​​are mainly derived from field pumping tests and regional hydrogeological survey results. For areas without test data, empirical values ​​of the same lithological strata are referenced. The flow field simulation calculation module is used to perform three-dimensional flow field dynamic evolution simulation calculation based on the twin substrate model. It uses the three-dimensional finite difference numerical method to solve the seepage control equation and outputs the water level, flow velocity, and hydraulic gradient flow field data at each time step. The default time step is 1 day, and the step size can be adaptively adjusted according to the flow field change rate. The virtual-real deviation assessment module is used to align the simulation results with the measured data of the same period on a cycle-by-cycle basis. By default, it performs deviation calculation once a day, and maps the simulation value and measured value of the monitoring point to each other in terms of spatial location and time step, calculates the virtual-real deviation in multiple dimensions, and divides it into three levels of deviation: mild, moderate and severe. The graded correction execution module is bidirectionally connected to the virtual-real deviation assessment module, the 3D hydrogeological modeling module, and the flow field simulation calculation module. It executes the corresponding correction strategy according to the deviation level. After the correction is completed, it is verified by back-calculation using the same measured data. The back-calculation verification is qualified when the comprehensive deviation index drops to the slight deviation range. After the verification is qualified, the updated model and parameters are sent back to the flow field simulation calculation module, forming a closed-loop operation mechanism consisting of monitoring input, simulation deduction, deviation assessment, automatic correction, and verification backflow update. This enables the model to continuously self-calibrate with the dynamic evolution of the actual groundwater system without frequent manual intervention.

[0006] Furthermore, the virtual-real deviation assessment module calculates three deviation indicators: relative mean square error of water level, relative error of flow velocity, and hydraulic gradient deviation rate. These are then weighted and summed to obtain a comprehensive virtual-real deviation index, calculated using the following formula:

[0007] In the formula, To comprehensively measure the deviation between virtual and real values, This represents the relative mean square error of the water level. This is the relative error of the flow velocity. The hydraulic gradient deviation rate, , , The preset weights for the corresponding sub-items are determined based on the abundance of monitoring data and the flow field characterization ability of each indicator. Water level data is the most abundant and the characterization is the most intuitive, so it is usually given the highest weight. The sum of the three weights is 1, and the usual value ranges are 0.5 to 0.7, 0.2 to 0.3, and 0.1 to 0.2, respectively. The comprehensive virtual-real deviation index is divided into mild, moderate and severe deviations according to the preset three-level deviation thresholds. The three-level thresholds are determined based on the conventional simulation error statistics of the historical artificial calibration model of the study area. They can be flexibly adjusted according to the hydrogeological complexity and monitoring density of the study area. The conventional value range is 0.02 to 0.05 for mild threshold and 0.06 to 0.10 for moderate threshold.

[0008] Furthermore, the three-dimensional hydrogeological modeling module and the flow field simulation calculation module adopt a spatiotemporal multi-scale nested architecture, divided into three levels of mapping units: regional macro level, sub-regional meso level, and key point micro level. The regional macro level uses a kilometer-level basic grid to cover the entire study area, the sub-regional meso level uses a 500m-level medium-resolution grid, and the key point micro level uses a 200m-level high-precision densified grid for key areas such as over-extraction funnels or water sources. The virtual-real deviation evaluation module and the hierarchical correction execution module calculate the deviation and execute the correction in a hierarchical manner according to the three-level mapping units. The correction parameters and calibration results are transmitted bidirectionally between different levels of mapping units. The boundary conditions of the upper-level model are passed down to the lower-level sub-model, and the correction results of the lower-level model are aggregated and synchronized to the upper-level model, ensuring that the hydraulic parameters of the cross-scale flow field are matched and avoiding the problem of scale disconnection.

[0009] Furthermore, the hierarchical correction execution module incorporates a flow field singularity region identification unit. When a severe deviation is identified, the flow field singularity region identification unit identifies the spatial range of the drastic flow field change based on the characteristics of water level gradient and velocity abrupt change. The judgment criterion is that when the water level gradient of a grid cell exceeds twice the average water level gradient of the study area, it is classified as a flow field singularity region. The hierarchical correction execution module performs multi-level local densification on the three-dimensional computational grid in this range. The densification level is usually 2 to 3 levels, which improves the local grid resolution by 4 to 8 times. The computational grid in the smooth flow field region is merged and coarsened. While ensuring the accuracy of local calculations, the overall computational load is controlled, avoiding a significant decrease in computational efficiency due to global densification. The severe correction is completed simultaneously with the model topology reconstruction optimization.

[0010] Furthermore, during the parameter optimization process for mild and moderate deviations, the graded correction execution module synchronously couples the effective stress calculation process. Based on the effective stress principle, a drop in water level causes a decrease in pore water pressure and an increase in formation effective stress, which in turn leads to formation compression and a decrease in porosity. The changes in formation effective stress and porosity are calculated based on the changes in water level in the flow field, and the three hydraulic parameters—permeability coefficient, specific yield, and release coefficient—are updated accordingly. The typical range for permeability coefficient in shallow loose layers is 0.1 m / d to 10 m / d, the typical range for specific yield is 0.05 to 0.25, and the typical range for release coefficient is 0.0. The values ​​range from 0.01 to 0.01, with specific values ​​adjusted based on the stratigraphic particle size distribution and density. The parameter optimization process satisfies both parameter value range constraints and spatial smoothness constraints. The parameter value range constraints ensure that the parameters of each unit are always limited to the empirical values ​​of the corresponding lithology, avoiding geologically unreasonable values. The spatial smoothness constraints require that the parameter difference between adjacent units of the same lithology does not exceed 15% of the average parameter value of adjacent units, ensuring that the parameter distribution conforms to the sedimentary pattern. Under slight deviation, the parameter adjustment range is smaller than that under moderate deviation. Under slight deviation, the parameter adjustment range does not exceed 10% of the initial value, and under moderate deviation, the parameter adjustment range does not exceed 20% of the initial value.

[0011] Furthermore, during the boundary adjustment process for moderate deviations, the graded correction execution module simultaneously incorporates two types of driving factors: natural hydrological data and human activity data. Natural hydrological data includes rainfall infiltration and surface water exchange, while human activity data includes groundwater extraction, artificial water transfer, and irrigation infiltration. These two types of data jointly drive the dynamic updating of the spatial distribution and intensity of the source-sink boundary. The data is converted into planar or point boundary conditions corresponding to the simulation grid through a spatial distribution method combining Thiessen polygons and land use type weighting. The recharge and discharge intensity of each grid unit is recalculated, and the interlayer overflow coefficient is adjusted to correct the hydraulic connection between aquifers. The interlayer overflow coefficient is calculated based on the permeability coefficient and thickness of the impermeable layer, with a typical value range of 10. -6 d -1 Up to 10 -4 d -1 .

[0012] Furthermore, the flow field simulation calculation module has a built-in multi-scenario parallel simulation unit, which can simultaneously launch multiple independent simulation instances based on the same twin base model. The number of simulation instances that can run simultaneously is determined according to the server's computing power configuration. Each simulation instance corresponds to a different mining scheme, water diversion scheme, or climate scenario. Each simulation instance shares real-time monitoring benchmark data, independently performs flow field evolution calculations, and synchronously outputs multiple sets of flow field evolution results, which facilitates the comparison of the effects of different control schemes and the selection of schemes without the need to repeatedly build models and set basic parameters.

[0013] Furthermore, the flow field simulation calculation module embeds a critical state early warning unit, which pre-sets flow field critical thresholds corresponding to two types of geological environmental risks: aquifer desiccation and ground subsidence triggering. The aquifer desiccation threshold corresponds to the elevation of the aquifer floor; when the groundwater level drops below the floor elevation, the aquifer is determined to be desiccated. The ground subsidence triggering threshold corresponds to the water level value corresponding to the critical effective stress of the stratum; when the water level drops and the effective stress exceeds the critical value, the risk of ground subsidence is triggered. The critical thresholds are determined based on the results of stratum compression tests and historical subsidence observation data in the study area. During the simulation, the flow field state and critical thresholds are compared in real time for each grid cell. After identifying the critical state point, the expected occurrence time, affected spatial range, and probability of occurrence of the corresponding state are output.

[0014] Furthermore, the multi-source data access module incorporates a physical constraint verification unit. Based on Darcy's law and the water conservation equation, it performs time-series rate of change verification and spatial hydraulic gradient verification on the accessed monitoring data. The time-series rate of change verification compares the rate of change of a single well with the average rate of change of adjacent monitoring wells. Data exceeding three times the average rate of change and not conforming to the seepage conduction law is marked as abnormal. The spatial hydraulic gradient verification compares the hydraulic gradient between adjacent wells with the empirical range of the corresponding lithology. Data exceeding the range is marked as abnormal. The module automatically identifies abnormal monitoring values ​​that do not conform to the hydrodynamic law and repairs them by combining the interpolation results of surrounding normal monitoring points before inputting them into the flow field simulation calculation module. The repair uses an inverse distance weighted interpolation method of 3 to 5 surrounding effective monitoring points to avoid abnormal data interfering with model calibration and simulation calculation.

[0015] Furthermore, the system also includes a ground settlement simulation module, which is bidirectionally data-coupled with the flow field simulation calculation module. The water level change data output by the flow field simulation calculation module serves as the driving input for the ground settlement simulation module, calculating the compression amount and ground settlement rate of each stratum. The stratum compressibility coefficient is derived from indoor geotechnical test results or regionally accepted empirical values, with the conventional range for the compressibility coefficient of the shallow loose layer being 0.1 MPa. -1 up to 0.5MPa -1 The ground compaction parameters calculated by the ground settlement simulation module are fed back to update the aquifer hydraulic parameters of the flow field simulation module, realizing the bidirectional coupling calculation of flow field evolution and ground settlement, so that the flow field simulation results and ground settlement calculation results can be mutually verified and updated synchronously.

[0016] Compared with existing technologies, this digital twin-based regional groundwater three-dimensional flow field dynamic evolution simulation system has the following advantages: I. This invention adds a virtual-to-real deviation assessment module and a graded correction execution module to the groundwater flow field twin simulation system, forming a closed-loop operation mechanism of monitoring input, simulation deduction, deviation assessment, automatic correction, verification backflow update, and changing the traditional static model operation mode that relies on manual periodic calibration. The system can match the correction strategy corresponding to the depth according to the deviation level between simulation results and measured data. Mild deviations adjust local hydraulic parameters, moderate deviations adjust source and sink boundaries and hydraulic connections, and severe deviations adjust the model topology and computational grid. It will not overcorrect when the deviation is small, and can also cope with the complex situation of sudden changes in flow field structure. It does not require frequent manual full parameter inversion and can maintain the accuracy of flow field simulation for a long time.

[0017] Second, this invention adopts a spatiotemporal multi-scale nested twin modeling architecture, combined with supporting technologies such as effective stress coupling calculation, adaptive grid adjustment in singular regions of the flow field, and multi-scenario parallel simulation. It simultaneously considers the computational efficiency of large-area flow field simulation and the simulation accuracy of key areas. The system can reproduce the evolution process of formation parameters caused by water level changes, and supports the simultaneous simulation of multiple mining and water transfer schemes. It can also be connected to the ground subsidence module for collaborative calculation, matching the actual use needs of daily management and risk analysis of regional groundwater.

[0018] Other advantages, objectives and features of the invention will be set forth in part in the description which follows, and in part will be apparent to those skilled in the art from the following examination or study, or may be learned from the practice of the invention. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are merely some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without any creative effort.

[0020] Figure 1 This is a schematic diagram of the overall operation process of the system of the present invention; Figure 2 This is a schematic diagram of the hierarchical correction execution process of the present invention; Figure 3 This is a schematic diagram of the virtual-real deviation calculation and classification process of the present invention. Detailed Implementation

[0021] To further illustrate the technical means and effects of the present invention in achieving its intended purpose, the following detailed description of the specific implementation methods, structures, features, and effects of the present invention, in conjunction with the accompanying drawings and preferred embodiments, is provided below.

[0022] This paper describes the specific implementation of a digital twin-based three-dimensional dynamic evolution simulation system for regional groundwater flow field in a shallow groundwater over-extraction area of ​​a plain. The study area covers approximately 4800 km². 2 The area contains shallow porous aquifers, which have long been affected by agricultural irrigation and industrial mining, resulting in a significant groundwater drawdown cone. Flow field changes are highly susceptible to human activity. The system is deployed on a cloud computing server, connecting to existing monitoring station data in the study area, enabling dynamic flow field evolution simulation and autonomous model iteration. Figure 1 As shown, the system adopts a closed-loop operating architecture, and the data flow and interaction relationships between the modules are as follows.

[0023] The multi-source data access module is integrated with real-time monitoring data from over 120 water level monitoring wells, 15 rain gauge weather stations, and 32 key water intakes within the study area. It also integrates monthly irrigation area statistics and manual water diversion scheduling data. The accessed data covers various monitoring indicators, including water level, water temperature, rainfall, evaporation, and extraction volume. The sampling frequency for different data types varies from hourly to daily, and after integration, the data is uniformly converted to the same spatiotemporal reference. The module has a built-in physical constraint verification unit that performs two levels of rationality verification on the accessed data based on Darcy's law and the water conservation equation. At the single-point time-series verification level, the daily water level change of a single monitoring well is first calculated, and then compared with the average water level change of three adjacent monitoring wells in the same layer. If the rate of change of a single well exceeds three times the average rate and does not conform to the groundwater seepage and conduction law, it is marked as an outlier. At the regional equilibrium verification level, on the one hand, the hydraulic gradient is obtained by calculating the ratio of water level difference to distance between adjacent monitoring wells. This gradient is then compared with the empirical range of hydraulic gradients corresponding to the aquifer lithology, and points exceeding the range are marked as outliers. On the other hand, the equilibrium relationship between the total recharge and total discharge in the study area is calculated. If the difference between the two exceeds a reasonable range, the monitoring data for the corresponding time period is simultaneously reviewed. Identified outliers are corrected using inverse distance weighted interpolation results from 3 to 5 surrounding valid monitoring points to avoid interference with model calibration. In addition to outlier correction, the preprocessing includes data cleaning, missing measurement interpolation, and temporal smoothing. Data with jumps or discontinuities are repaired and smoothed to ensure the continuity and reliability of the input data. After verification, all data are unified to the same spatiotemporal reference. A continuous monitoring dataset matching the simulation grid and simulation step size is generated through spatiotemporal interpolation and simultaneously input into the flow field simulation calculation module and the virtual-real deviation evaluation module.

[0024] The 3D hydrogeological modeling module, based on data from 86 boreholes, 7 hydrogeological profiles, and pumping tests in the study area, constructs a 3D structured model including shallow aquifers and underlying aquitards. Initial hydraulic parameters such as permeability coefficient, specific yield, release coefficient, and porosity are assigned to each hydrogeological unit. The initial parameter values ​​are primarily derived from field pumping tests and regional hydrogeological surveys; for areas without test data, empirical values ​​from strata of similar lithology are used. Specifically, the permeability coefficient of the shallow loose layer typically ranges from 0.1 m / d to 10 m / d, the specific yield from 0.05 to 0.25, the release coefficient from 0.001 to 0.01, and the porosity from 0.3 to 0.45, with specific values ​​adjusted according to the stratigraphic particle size distribution and density. The model adopts a spatiotemporal multi-scale nested architecture, with three levels of mapping units from top to bottom. The macroscopic layer covers the entire study area, using a 1000m×1000m base grid to analyze the overall evolution trend of the flow field and regional water balance. The mesoscopic layer, divided into six sub-regions according to hydrogeological units, uses a 500m×500m medium-resolution grid for zoning deviation assessment and control analysis. The microscopic layer, targeting the central area of ​​the over-extraction funnel, employs a high-precision 200m×200m grid for detailed simulation of local flow field characteristics. A bidirectional data interface is established between the three levels of units, allowing boundary conditions from the upper-level model to be passed down to the lower-level sub-models, and allowing corrections from the lower-level model to be aggregated and synchronized back to the upper-level model, ensuring the matching of hydraulic parameters across different scales of the flow field.

[0025] The virtual-to-real deviation assessment module performs deviation calculations once a day by default. Figure 3 As shown, the simulation results for the corresponding time period were spatiotemporally aligned with the measured data of the same period, and three deviation indicators were calculated respectively. The relative mean square error of water level was calculated by first calculating the mean square error between the simulated and measured water levels at all monitoring points, and then dividing it by the average measured water level of the same period to obtain the relative error. The relative error of flow velocity was calculated by selecting three key hydrogeological sections in the study area, calculating the absolute value of the difference between the average simulated and measured flow velocities of the sections, and then dividing it by the measured flow velocity to obtain the relative error. The deviation rate of hydraulic gradient was calculated by selecting 20 pairs of adjacent monitoring wells, calculating the absolute value of the difference between the simulated and measured hydraulic gradient values ​​between the wells, and then dividing it by the measured hydraulic gradient to obtain the deviation rate.

[0026] The three deviation indicators are weighted and summed to obtain the comprehensive virtual-real deviation index, calculated using the following formula:

[0027] In the formula, To comprehensively measure the deviation between virtual and real values, This represents the relative mean square error of the water level. This is the relative error of the flow velocity. The hydraulic gradient deviation rate, , , The preset weights are assigned to the corresponding sub-items. In this embodiment, the weights of the three items are set to 0.6, 0.25, and 0.15, respectively. The weight values ​​are determined based on the abundance of monitoring data and the flow field characterization ability of each indicator. Water level data is the most abundant and has the most intuitive characterization, so it is given the highest weight.

[0028] The comprehensive virtual-real deviation index is used to classify deviation levels according to a preset three-level threshold: less than 0.03 indicates slight deviation, greater than or equal to 0.03 and less than 0.08 indicates moderate deviation, and greater than or equal to 0.08 indicates severe deviation. The three-level thresholds are determined based on the conventional simulation error statistics of the historical artificial calibration model of the study area and can be flexibly adjusted according to the hydrogeological complexity and monitoring density of the study area. Deviation calculation is performed separately for each of the three mapping units, and the deviation level and spatial distribution of each unit are output as the basis for graded correction.

[0029] The graded correction execution module executes the corresponding correction strategy based on the deviation assessment results, such as... Figure 2 As shown, after the correction is completed, the measured data from the previous cycle is used for back-calculation and verification. A passing result is considered acceptable if the comprehensive deviation index drops below 0.03. After successful verification, the updated model and parameters are sent back to the flow field simulation calculation module, completing one round of closed-loop iteration. If the deviation still does not meet the standard after a single round of correction, the correction intensity is automatically increased by one level and the process is repeated.

[0030] When a slight deviation is identified, it indicates that the overall model structure and boundary conditions are reasonable, with only minor deviations in local parameters. In this case, only local hydraulic parameter optimization is performed. Guided by the spatial distribution of the deviation, local hydrogeological units with deviations exceeding the threshold are identified, and a constrained particle swarm optimization algorithm is used to adjust the permeability coefficient and specific yield parameters. Two constraints are set during the parameter optimization process: first, a parameter value range constraint, ensuring that the parameters of each unit are always limited to the empirical values ​​for the corresponding lithology to avoid parameter values ​​that do not conform to geological laws; second, a spatial smoothing constraint, ensuring that the parameter difference between adjacent units of the same lithology does not exceed 15% of the average parameter value of adjacent units, guaranteeing that the parameter distribution conforms to sedimentary laws. Under slight deviation, the parameter adjustment amplitude does not exceed 10% of the initial value. The correction process is simultaneously coupled with effective stress calculation. Based on the principle of effective stress, a drop in water level will cause a decrease in pore water pressure and an increase in formation effective stress, leading to formation compression and a decrease in porosity. The changes in formation effective stress and porosity are calculated based on the water level change, and then the change in permeability is obtained through empirical seepage formulas. The hydraulic parameters are updated accordingly to restore the aquifer parameter evolution characteristics caused by water level fluctuations. The formation compressibility coefficient was derived from indoor geotechnical tests in the study area. When no test data was available, regionally accepted empirical values ​​were used. The conventional range for the compressibility coefficient of shallow loose layers was 0.1 MPa. -1 up to 0.5MPa -1 .

[0031] When a deviation is determined to be moderate, it indicates that adjusting parameters alone is insufficient to eliminate the deviation, and the source of the deviation is often a change in boundary conditions or alteration of the interlayer hydraulic connection. In this case, parameter optimization and boundary adjustment are performed simultaneously. The parameter optimization range is expanded to the entire corresponding sub-region, and the adjustment range is widened to 20% of the initial value, while still retaining the value range constraints and spatial smoothing constraints. Boundary adjustment incorporates two types of driving factors: natural hydrological data and human activity data. Natural hydrological data includes rainfall infiltration and surface water exchange, while human activity data includes groundwater extraction, artificial water transfer, and irrigation infiltration. Both types of data are converted into areal or point boundary conditions corresponding one-to-one with the simulation grid using a Thiessen polygon combined with land use type weighted spatial distribution. The recharge and discharge intensity of each grid cell is recalculated, and the interlayer overflow coefficient is adjusted to correct the hydraulic connection between aquifers. The interlayer overflow coefficient is calculated based on the permeability coefficient and thickness of the impermeable layer, with a typical value range of 10. -6 d -1 Up to 10 -4 d -1 .

[0032] When a severe deviation is identified, it indicates a significant discrepancy between the model structure and actual geological conditions, potentially caused by factors such as aquifer desiccation and changes in water-conducting channels. At this point, three operations are performed simultaneously: First, singularity zone identification and mesh adjustment are performed. The criterion is that if the water level gradient of a mesh cell exceeds twice the average water level gradient of the study area, that cell is classified as a singularity zone. The 3D computational mesh in this region is then locally refined at levels 2 to 3. Mesh in areas with gentle flow is merged and coarsened to improve local accuracy while controlling the overall computational load. Next, model topology correction is performed, verifying the spatial distribution boundaries of the aquifer, updating the desiccation zone range and water-conducting channel distribution, and adjusting the model topology to match long-term changes in geological conditions. Finally, global parameter inversion is performed. Based on the updated mesh and topology, a geologically constrained variational adjoint method is used for global parameter inversion. During the inversion process, parameter value ranges and spatial smoothing constraints are applied simultaneously to complete the model's depth calibration. Severe correction involves a large amount of computation and is triggered on a monthly basis by default. If the comprehensive deviation index is still in the moderate or above range after two consecutive weekly moderate corrections, severe correction will be automatically triggered.

[0033] The flow field simulation module is built on a three-dimensional finite difference numerical kernel, with a default time step of 1 day, which can be adaptively adjusted according to the flow field change rate. The module incorporates a multi-scenario parallel simulation unit, which can simultaneously launch multiple independent simulation instances based on the same twin-based model. The number of simulation instances that can run simultaneously depends on the server's computing power configuration. Under normal configuration, it supports parallel calculation of multiple scenarios. Each instance can independently set mining schemes, water diversion schemes, or climate scenarios, share real-time monitoring benchmark data, independently complete flow field calculations, and synchronously output multiple sets of flow field evolution results for comparison and selection of control schemes. The module also embeds a critical state early warning unit, which pre-sets the flow field critical thresholds corresponding to two types of geological environmental risks: aquifer dewatering threshold corresponding to the aquifer floor elevation, and ground subsidence trigger threshold corresponding to the water level value corresponding to the critical effective stress of the stratum. During the simulation, the flow field state is compared with the critical threshold in real time for each grid cell. After identifying the critical state point, the expected occurrence time, affected spatial range, and probability of occurrence of the corresponding state are output.

[0034] The system also includes a land subsidence simulation module, which is bidirectionally data-coupled with the flow field simulation module. The unit-by-unit water level change data output by the flow field simulation module serves as the driving input for the land subsidence simulation module, calculating the compression and land subsidence rate of each stratum. The formation porosity change parameters calculated by the land subsidence simulation module are fed back to update the aquifer hydraulic parameters of the flow field simulation module, achieving bidirectional coupled calculation of flow field evolution and land subsidence.

[0035] The visualization module constructs hydrogeological scenes based on a 3D rendering engine, supports the display of aquifer layers, arbitrary profile cutting, and dynamic display of flow field streamlines. It can overlay deviation levels, correction records, and early warning information on the 3D model and supports interactive querying of historical flow fields and simulation results.

[0036] This embodiment achieves autonomous model iteration through a closed-loop self-evolution mechanism. The average error of the water level simulation, after one year of continuous operation, remains stable within 2.5%. No manual full-scale model calibration is required, and the calibration cycle is shortened from the traditional 3-6 months to a single day, reducing the workload of manual calibration and maintenance investment. This system is not only suitable for over-extraction areas of groundwater in loose plains, but can also be adapted to different scenarios such as karst development areas, bedrock fissure areas, and centralized water supply sources by adjusting parameters and model structure, meeting the groundwater flow field simulation and management needs of different regions.

[0037] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins, characterized in that, The system includes a multi-source data access module, a 3D hydrogeological modeling module, a flow field simulation calculation module, a visualization module, a virtual-real deviation assessment module, and a graded correction execution module. The multi-source data access module is used to access groundwater monitoring, hydrological and meteorological data, and water intake data from human activities at the physical end. After preprocessing, the data is synchronously input into the flow field simulation calculation module and the virtual-real deviation evaluation module. The 3D hydrogeological modeling module is used to construct a virtual twin base model of the aquifer structure in the study area, providing the geometric boundaries and hydraulic parameters required for flow field simulation; The flow field simulation calculation module is used to perform three-dimensional flow field dynamic evolution simulation calculations based on the twin-based model, and outputs the water level, flow velocity, and hydraulic gradient flow field data at each time step; The virtual-real deviation assessment module is used to align the simulation results with the measured data of the same period on a cycle-by-cycle basis, calculate the virtual-real deviation in multiple dimensions, and classify the deviation into three levels: mild, moderate and severe. The graded correction execution module is bidirectionally connected to the virtual-real deviation assessment module, the three-dimensional hydrogeological modeling module, and the flow field simulation calculation module. It executes the corresponding correction strategy according to the deviation level. After the correction is completed, it is verified by back-calculation using the same measured data. After the verification is qualified, the updated model and parameters are sent back to the flow field simulation calculation module, forming a closed-loop operation mechanism consisting of monitoring input, simulation deduction, deviation assessment, automatic correction, verification backflow update.

2. The regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins according to claim 1, characterized in that, The virtual-real deviation assessment module calculates three deviation indicators: relative mean square error of water level, relative error of flow velocity, and hydraulic gradient deviation rate. A comprehensive virtual-real deviation index is obtained by weighted summation, and the calculation formula is as follows: ； In the formula, To comprehensively measure the deviation between virtual and real values, This represents the relative mean square error of the water level. This is the relative error of the flow velocity. The hydraulic gradient deviation rate, , , Preset weights for the corresponding sub-items; The comprehensive virtual-real deviation index is divided into mild, moderate and severe deviations according to the preset three-level deviation threshold.

3. The regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins according to claim 1, characterized in that, The three-dimensional hydrogeological modeling module and flow field simulation calculation module adopt a spatiotemporal multi-scale nested architecture, which is divided into three-level mapping units: regional macro layer, sub-region meso layer, and key point micro layer. The virtual-real deviation evaluation module and the hierarchical correction execution module calculate the deviation and execute the correction in a hierarchical manner according to the three-level mapping units. The correction parameters and calibration results are transmitted bidirectionally between different levels of mapping units.

4. The regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins according to claim 1, characterized in that, The graded correction execution module has a built-in flow field singularity region identification unit. When a severe deviation is determined, the flow field singularity region identification unit identifies the spatial range of the drastic flow field change based on the characteristics of water level gradient and flow velocity change. The graded correction execution module performs multi-level local densification of the three-dimensional computational grid in this range, merges and coarsens the computational grid in the smooth flow field region, and completes the severe correction in sync with the model topology reconstruction optimization.

5. The regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins according to claim 1, characterized in that, The graded correction execution module simultaneously couples the effective stress calculation process during the parameter optimization process for mild and moderate deviations. It calculates the changes in formation effective stress and porosity based on the changes in flow field water level, and updates the three types of hydraulic parameters, namely permeability coefficient, specific yield, and release coefficient. The parameter optimization process satisfies the parameter value range constraint and spatial smoothness constraint, and the parameter adjustment range is smaller under mild deviation than under moderate deviation.

6. The regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins according to claim 1, characterized in that, During the boundary adjustment process for moderate deviation, the graded correction execution module simultaneously incorporates two types of driving factors: natural hydrological data and human activity data. Natural hydrological data includes rainfall infiltration and surface water exchange, while human activity data includes groundwater extraction, artificial water transfer, and irrigation infiltration. These two types of data jointly drive the dynamic update of the spatial distribution and intensity of the source-sink boundary.

7. The regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins according to claim 1, characterized in that, The flow field simulation calculation module has a built-in multi-scenario parallel simulation unit, which can simultaneously launch multiple independent simulation instances based on the same twin base model. Each simulation instance corresponds to a different mining scheme, water diversion scheme, or climate scenario. Each simulation instance shares real-time monitoring benchmark data, independently performs flow field evolution calculations, and synchronously outputs multiple sets of flow field evolution results.

8. The regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins according to claim 1, characterized in that, The flow field simulation calculation module has an embedded critical state early warning unit, which pre-sets the flow field critical thresholds corresponding to two types of geological environmental risks: aquifer desiccation and ground subsidence. During the simulation, the flow field state is compared with the critical threshold in real time for each grid cell. After identifying the critical state point, the expected occurrence time, affected spatial range and occurrence probability of the corresponding state are output.

9. The regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins according to claim 1, characterized in that, The multi-source data access module has a built-in physical constraint verification unit. Based on Darcy's law and the water conservation equation, it performs time-series change rate verification and spatial hydraulic gradient verification on the accessed monitoring data. It automatically identifies abnormal monitoring values ​​that do not conform to the hydrodynamic law, and repairs them by combining the interpolation results of surrounding normal monitoring points before inputting them into the flow field simulation calculation module.

10. The regional groundwater three-dimensional flow field dynamic evolution simulation system based on digital twins according to claim 1, characterized in that, The system also includes a ground settlement simulation module, which is bidirectionally data coupled with the flow field simulation calculation module. The water level change data output by the flow field simulation calculation module serves as the driving input for the ground settlement simulation module, and the ground compaction parameters calculated by the ground settlement simulation module are fed back to update the aquifer hydraulic parameters of the flow field simulation calculation module.

Images