Runoff driving factor identification method integrating field pumping test and hydrodynamic simulation

Abstract

The invention belongs to the technical field of hydrogeology and water resources, and particularly relates to a runoff driving factor identification method fusing a field water pumping test and hydrodynamic simulation, which comprises the following steps: acquiring parameters such as nano-pore scale zeta potential, pH value and ionic strength in real time by an integrated in-situ sensing unit, constructing an electricity-flow coupling response function, and identifying a runoff driving factor of a water pumping test; dynamically correcting a critical Reynolds number threshold value of a non-Darcy flow in the regional hydrodynamic model; based on a bidirectional coupling iteration mechanism, the constitutive relation between the aquifer parameter field and the seepage flow is synchronously optimized, and the basic flow path, the supply intensity and the multi-source runoff driving factor contribution rate are accurately recognized. Physical self-consistent fusion of the microcosmic interface electrodynamic effect and macroscopic water flow simulation is achieved, and the accuracy and reliability of runoff cause analysis under the complex hydrogeological condition are remarkably improved.

Description

Technical Field

[0001] This invention belongs to the field of hydrogeology and water resources technology, specifically a method for identifying runoff driving factors that integrates field pumping tests and hydrodynamic simulation. Background Technology

[0002] In the fields of groundwater science and water resources management, accurately identifying the driving factors of runoff at the watershed scale, especially the contribution rate of baseflow to river discharge, is a core prerequisite for assessing aquifer sustainability, formulating water source protection strategies, and predicting extreme climate responses. Traditionally, this task relies on the collaborative analysis of field pumping tests and numerical hydrodynamic simulations: the former obtains aquifer parameters through measured drawdown and discharge relationships, while the latter constructs regional flow models based on Darcy's law to inversely determine the dynamic contributions of different recharge sources. In recent years, with the development of high-resolution grids and unstructured discretization techniques, advanced simulation platforms, represented by MODFLOW-USG, have significantly improved their ability to characterize complex geological structures. However, their theoretical foundation is still strictly based on the assumption of a macroscopic continuous medium, failing to fully incorporate the regulatory role of microscopic pore-scale physicochemical processes on macroscopic seepage behavior.

[0003] Specifically, existing hydrodynamic simulation methods generally treat aquifer pores as electrically neutral, static geometric channels, relying solely on empirical correlations of permeability coefficients based on particle size distribution or porosity, completely ignoring the dynamic evolution mechanism of the solid-liquid interface double layer within nanoscale pores (pore size less than 100 nm). However, numerous experimental observations show that in natural groundwater environments, even small fluctuations in pH can trigger protonation / deprotonation of functional groups on the pore surface, leading to drastic changes in the zeta potential (measured amplitudes often exceeding ±30 mV). This dynamic surface charge directly modulates the ionic atmosphere thickness and fluid slip boundary conditions within the pores, thus significantly altering the critical criterion for non-Darcy flow—the Reynolds number threshold. Under typical operating conditions, this critical value can increase from the traditionally accepted 10⁻⁻⁻⁶. 4 The abrupt jump to 10⁻² marks a sudden transition of seepage from the linear Darcy zone to the nonlinear inertial-dominated zone. This abrupt change is not a gradual parameter drift, but a phase-change response triggered by interfacial electrodynamics, causing unpredictable step-like distortions in the drawdown-time curves observed in pumping tests. Under such distortion, using the classical Theis or Jacob formulas to invert the permeability coefficient will introduce a systematic bias of up to 47%, severely distorting the true characterization of the aquifer's hydraulic conductivity.

[0004] Fundamentally, the mismatch between existing simulation systems and field measurement data does not stem from insufficient numerical algorithm accuracy or limited grid resolution, but rather from a structural blind spot in its underlying physical model: it cannot accommodate the nonlinear modulation of macroscopic seepage constitutive relations by nanopore electrodynamics. When simulation programs force the distorted pumping response curves to fit the Darcy framework, to compensate for the residuals between observation and simulation, the model is forced to "mask" the lack of physical mechanisms by irrationally adjusting regional hydraulic conductivity coefficients or source-sink term strengths. This ultimately leads to incorrect baseflow path identification, distorted recharge zoning, and even the invalidation of the entire runoff driving factor analysis conclusions. It is worth noting that this contradiction is inherently irreconcilable—even with higher-frequency field monitoring or more refined geological modeling, as long as a pore-scale electro-current coupling mechanism is not embedded in the governing equations, the simulation results will still systematically deviate from the actual hydrological process. Furthermore, since nanopores are widely present in clay minerals, organic cementation zones, and weathering fissure fillings, their charge sensitivity is particularly prominent under disturbance scenarios such as acid rain input, agricultural non-point source pollution, or seawater intrusion, causing the aforementioned distortion problems to be amplified exponentially in areas with dense human activity.

[0005] Therefore, the present invention provides a method for identifying runoff driving factors that integrates field pumping tests and hydrodynamic simulation. Summary of the Invention

[0006] To achieve the above-mentioned objectives, this invention provides a method for identifying runoff driving factors that integrates field pumping tests and hydrodynamic simulations. Its core lies in sensing the dynamic surface potential of nanopores in situ and reconstructing the critical Reynolds number threshold in the non-Darcy flow seepage control equation in real time, thereby establishing physically self-consistent seepage boundary conditions in the regional hydrodynamic model and ultimately achieving high-fidelity identification of baseflow paths, recharge intensity, and multi-source runoff driving factors.

[0007] The method of the present invention includes the following steps:

[0008] First, several pumping test well groups were set up in the target watershed. Each pumping well was equipped with an integrated in-situ sensing unit, which integrates a microscale fluid sampling probe, a solid-liquid interface potential detection module and a water level-flow synchronous recording device.

[0009] Secondly, based on the time series data of groundwater pH, ionic strength, colloidal concentration and zeta potential of pore walls obtained by the in-situ sensing unit, a nanopore-scale electro-current coupling response function is constructed.

[0010] Furthermore, the response function is embedded into the non-Darcy flow criterion module of the regional hydrodynamic simulation platform to dynamically update the critical Reynolds number threshold corresponding to each computational grid cell;

[0011] Finally, based on the modified seepage constitutive relation, a two-way coupled iteration of pumping response inversion and regional flow simulation is performed to output a spatial distribution map of aquifer parameter field and baseflow contribution rate that conforms to physical reality.

[0012] Furthermore, the integrated in-situ sensing unit comprises a waterproof encapsulation shell, a multi-channel microfluidic chip, a double-layer potential sensor array, a piezoresistive level gauge, and an electromagnetic flow meter. The multi-channel microfluidic chip features a biomimetic nanopore network with a pore size distribution ranging from 2 nanometers to 100 nanometers, simulating the typical pore structure of clay minerals and organic cementation bands in natural aqueous media. The double-layer potential sensor array uses a gold-indium tin oxide composite electrode pair configuration, employing AC impedance spectroscopy to analyze the protonation state of functional groups on the pore walls in real time, and outputting a phase difference signal that monotonically maps to the zeta potential. The piezoresistive level gauge and electromagnetic flow meter are connected to a central data acquisition unit via an RS485 bus, ensuring strict alignment of the timestamps for water level drawdown and pumping flow rate, with a time synchronization accuracy better than 10 milliseconds.

[0013] As a preferred embodiment of the present invention, the construction process of the nanopore-scale electro-current coupling response function is as follows: using the measured ζ-potential in the microfluidic chip as the input variable, combined with the synchronously acquired fluid conductivity and dielectric constant, the effective slip length expression is derived using the Gouy-Chapman-Stern double-layer theory; then, this slip length is substituted into the wall boundary conditions of the Navier-Stokes equation, and a modified Darcy-Forchheimer type seepage equation is obtained through dimensionless processing, in which the inertia term coefficient is exponentially correlated with the ζ-potential; finally, the nonlinear relationship curves between seepage velocity and pressure gradient under different ζ-potentials are numerically solved based on this equation, and the Reynolds number is extracted to be 10⁻. 4 The critical velocity value corresponding to the time is defined as the non-Darcy flow initiation criterion under the current operating conditions.

[0014] Furthermore, the regional hydrodynamic simulation platform adopts an unstructured grid discrete architecture, with each grid cell bound to an independent set of seepage constitutive parameters. This set of parameters includes the traditional hydraulic conductivity coefficient, storage coefficient, and the critical Reynolds number threshold dynamically generated by the aforementioned response function. During the simulation initialization phase, the platform reads the initial zeta potential data uploaded by the in-situ sensing units of all pumping wells and calls the pre-stored electro-current coupling response function library to complete the initial assignment of the critical Reynolds number field for the entire region. During the pumping test, the platform receives the updated measured zeta potential values ​​of each well point at a fixed time step and triggers the recalculation operation of the critical threshold of the local grid cells. The recalculation range is dynamically expanded according to the chemical diffusion radius of pore water to ensure that the spatial propagation effect of potential disturbance is accurately captured.

[0015] As a key technical feature of this invention, the bidirectional coupling iterative mechanism comprises two nested loops: the outer loop performs pumping response fitting, using the regional flow field output by the inner loop as a constraint, adjusting the local permeability coefficient around the pumping well until the sum of squared residuals between the simulated drawdown curve and the measured curve converges to a preset tolerance; the inner loop, based on the updated permeability coefficient field, resolves the global nonlinear seepage control equations and simultaneously refreshes the critical Reynolds number thresholds for each grid cell; when the actual Reynolds number in a grid cell exceeds its current threshold, the system automatically activates the Forchheimer inertia term; otherwise, it maintains the linear Darcy flow mode; this switching logic is controlled by a state machine, and the state transition condition is the real-time comparison result between the Reynolds number and the critical threshold, and each transition triggers a Jacobian matrix reconstruction to ensure the numerical stability of the nonlinear equation solution.

[0016] Furthermore, the generation of the baseflow contribution rate spatial distribution map relies on the source-sink term decomposition module. Based on the finally converged flow field, this module performs mass conservation integrals on each river confluence node to separate the instantaneous recharge fluxes formed by precipitation infiltration, surface recharge, and deep confined water overflow. Simultaneously, it uses the particle tracking algorithm to integrate along the streamlines in reverse to identify the upstream aquifer units corresponding to each recharge flux and to calculate their proportion in the total baseflow. In this process, all streamline calculations adopt the fourth-order Runge-Kutta scheme and use a dynamically updated seepage velocity field as input to ensure that the streamline bending effect in the non-Darcy region is accurately reflected.

[0017] As another important component of the present invention, the biomimetic nanopore network of the microfluidic chip is prepared on a single-crystal silicon substrate by focused ion beam etching. The orientation and branch topology of the pores are generated by fractal dimension matching based on the CT scan image of the target aquifer core. The inner wall of the pores is treated with oxygen plasma and then grafted with carboxyl or amino functional groups to simulate the acid-base active sites on the surface of natural minerals. Before the pumping test is started, the system injects a prepared solution with the same composition as the groundwater on site into the microfluidic chip and applies the same pressure gradient as the pumping well, so that the flow state inside the chip is similar to the hydrodynamics of the field pore environment.

[0018] Furthermore, the phase difference signal of the double-layer potential sensor array is input to the embedded signal processor after analog-to-digital conversion. The processor runs a sliding window Fourier transform algorithm to extract the phase offset between the AC excitation signal and the response current in real time. The phase offset is then fed into a pre-trained neural network mapper, which uses a large amount of calibration experimental data as a training set to establish a phase-ζ potential nonlinear mapping relationship. The output of the neural network is the effective ζ potential value under the current pore environment, and its update frequency is no less than once per second to ensure that the dynamic changes in potential are fully sampled.

[0019] As a system-level innovation of this invention, a closed-loop feedback link based on power line carrier communication is established between the regional hydrodynamic simulation platform and the field pumping test device. This link allows the simulation platform to actively send adjustment commands to the corresponding pumping wells when it detects that the seepage state in a certain area is about to cross the non-Darcy critical point, temporarily changing the pumping rate to avoid data distortion in the strong nonlinear region. At the same time, the in-situ sensing unit at the pumping well end can also automatically trigger the platform's global parameter refresh process when it detects that the zeta potential change exceeds the preset safety threshold, preventing local electrochemical disturbances from causing overall model instability.

[0020] Furthermore, the non-Darcy flow criterion module has a built-in hysteresis compensation mechanism to handle the non-instantaneous relationship between zeta potential changes and seepage response. This mechanism constructs a first-order inertial link transfer function based on measured historical data, and outputs the current zeta potential input as an effective control signal after being weighted by a time constant. The value of the time constant is determined based on pore water chemical equilibrium kinetic experiments and is dynamically adjusted with temperature and ion type to ensure that the time delay characteristics of the electro-current coupling process are accurately modeled.

[0021] As an engineering guarantee for the implementation of this invention, all in-situ sensing units are deployed in dedicated monitoring wells. These wells employ a double-casing structure, with the inner tube used for pumping operations and the outer tube annularly filled with highly permeable quartz sand and fitted with a water-stop plug to isolate vertical interference between different aquifers. The in-situ sensing units are fixed to the bottom of the inner tube via flange interfaces, with their microfluidic chip inlets facing the main flow direction to ensure that the sampled fluid represents the actual seepage path. Power is supplied by an underground supercapacitor energy storage module, which is periodically charged by a ground-based photovoltaic array via armored cables, supporting continuous operation for no less than 30 calendar days.

[0022] Furthermore, after the regional flow simulation is completed, the system automatically executes the sensitivity analysis module. This module freezes each driving factor (such as precipitation intensity, evaporation rate, and artificial extraction) in sequence and reruns the simplified model to quantify the independent contribution of each factor to the baseflow output. All sensitivity tests are conducted under corrected non-Darcy flow boundary conditions to avoid attribution bias caused by the lack of physical mechanisms. The final output runoff driving factor ranking table is arranged in descending order of contribution rate and the key aquifer unit numbers involved in the action path of each factor are marked.

[0023] The beneficial effects of this invention are as follows:

[0024] This invention presents a method for identifying runoff driving factors that integrates field pumping tests and hydrodynamic simulations. By sensing the dynamics of surface potential in nanopores in situ and correcting the non-Darcy flow critical criterion in real time, it avoids systematic misinterpretations of pumping response curves caused by neglecting interfacial electrodynamic effects in traditional hydrodynamic models. It enables the reconstruction of seepage constitutive relations that conform to microscopic physical mechanisms at the regional scale, ensuring that the aquifer parameter inversion results truly reflect the water transport capacity of the geological medium. It achieves accurate characterization of nonlinear seepage zones during baseflow path identification, preventing misalignment of recharge zones due to flow regime misjudgment. It ensures that the analysis of runoff driving factors remains physically consistent under charge-sensitive scenarios such as acid rain, agricultural pollution, or seawater intrusion. Furthermore, through the fluid dynamic similarity design of microfluidic chips and field hydrological conditions, it effectively bridges the scale gap between laboratory and field scales. Attached Figure Description

[0025] The invention will now be further described with reference to the accompanying drawings.

[0026] Figure 1 This is a flowchart of the method for identifying runoff driving factors that integrates field pumping tests and hydrodynamic simulation as described in this invention.

[0027] Figure 2 This is a structural block diagram of the integrated in-situ sensing unit in this invention;

[0028] Figure 3 This is a system block diagram of the non-Darcy flow criterion module and the bidirectional coupling iterative mechanism in the regional hydrodynamic simulation platform of this invention. Detailed Implementation

[0029] To make the technical means, creative features, objectives and effects of this invention easier to understand, the invention will be further described below in conjunction with specific embodiments.

[0030] like Figure 1 - Figure 3As shown in the embodiment of the present invention, a method for identifying runoff driving factors that integrates field pumping tests and hydrodynamic simulation includes deploying several pumping test well groups within a target watershed, with each well equipped with an integrated in-situ sensing unit. This integrated in-situ sensing unit consists of a waterproof enclosure, a multi-channel microfluidic chip, a double-layer potential sensor array, a piezoresistive level gauge, and an electromagnetic flowmeter. The waterproof enclosure is made of 316L stainless steel, achieving an overall sealing rating of IP68, allowing it to withstand long-term corrosion and high pressure in groundwater environments. The multi-channel microfluidic chip contains a biomimetic nanopore network with a pore size distribution ranging from 2 nanometers to 100 nanometers, used to simulate the typical pore structure of clay minerals and organic cementation bands in natural aquifers. This biomimetic nanoporous network was fabricated on a single-crystal silicon substrate using focused ion beam etching. The pore orientation and branch topology were generated based on fractal dimension matching of CT scan images of the target aquifer core. The inner walls of the pores were treated with oxygen plasma and then grafted with carboxyl or amino functional groups to simulate the acid-base active sites on the surface of natural minerals. Before the pumping test was initiated, a solution with the same composition as the groundwater in the field was injected into the microfluidic chip, and the same pressure gradient as the pumping well was applied, making the internal flow pattern of the chip similar to the hydrodynamics of the field pore environment.

[0031] The double-layer potential sensor array employs a gold-indium tin oxide composite electrode pair. It uses AC impedance spectroscopy to analyze the protonation state of functional groups on the pore wall in real time and outputs a phase difference signal that monotonically maps to the zeta potential. This phase difference signal is converted from analog to digital and input to an embedded signal processor. The processor runs a sliding window Fourier transform algorithm to extract the phase shift between the AC excitation signal and the response current in real time. This phase shift is then fed into a pre-trained neural network mapper, which uses a large amount of calibration experimental data as a training set to establish a nonlinear phase-zeta potential mapping relationship. The neural network output is the effective zeta potential value under the current pore environment, with an update frequency of no less than once per second, ensuring that dynamic potential changes are fully sampled. A piezoresistive level gauge and an electromagnetic flowmeter are connected to a central data acquisition unit via an RS485 bus, ensuring strict alignment of the timestamps for water level drawdown and pumping flow rate, with a time synchronization accuracy better than 10 milliseconds.

[0032] Based on the time-series data of groundwater pH, ionic strength, colloid concentration, and zeta potential at the pore wall obtained by the integrated in-situ sensing unit, a nanopore-scale electro-current coupling response function is constructed. The construction process of this response function is as follows: using the measured zeta potential in the microfluidic chip as the input variable, combined with the synchronously acquired fluid conductivity and dielectric constant, the expression for the effective slip length is derived using the Gouy-Chapman-Stern double-layer theory:

[0033]

[0034] in, Indicates the effective slip length (unit: m). is the relative permittivity of the fluid (dimensionless). The vacuum permittivity (unit: F / m, with values ​​of...) F / m), The zeta potential at the pore wall (unit: V). The viscosity is the fluid dynamic viscosity (unit: Pa·s). The slip velocity is given in m / s. This slip length is then substituted into the wall boundary conditions of the Navier-Stokes equations, and the modified Darcy-Forchheimer type seepage equation is obtained through dimensionless processing:

[0035]

[0036] in, Pressure (unit: Pa). The viscosity is the fluid dynamic viscosity (unit: Pa·s). Permeability (unit: m²). This is the seepage velocity vector (unit: m / s). Fluid density (unit: kg / m³). The inertial drag coefficient (unit: 1 / m) is expressed as follows:

[0037]

[0038] in, The reference inertia coefficient (unit: 1 / m). The constant is an empirical attenuation constant (unit: V⁻¹), both determined through microfluidic chip calibration experiments. Finally, based on this equation, the nonlinear relationship curves between seepage velocity and pressure gradient at different zeta potentials are numerically solved, and the Reynolds number is extracted. The critical velocity value corresponding to the specified time is defined as the non-Darcy flow initiation criterion under the current operating condition. Reynolds number. The definition of is:

[0039]

[0040] in, The characteristic pore diameter (unit: m) is taken as the average hydraulic diameter of the microfluidic chip channels.

[0041] The regional hydrodynamic simulation platform employs an unstructured grid discrete architecture. Each grid cell is bound to an independent set of seepage constitutive parameters, including the traditional hydraulic conductivity coefficient, storage coefficient, and a critical Reynolds number threshold dynamically generated by the aforementioned response function. During the simulation initialization phase, the platform reads the initial zeta potential data uploaded by the in-situ sensing units of all pumping wells and calls the pre-stored electro-current coupled response function library to complete the initial assignment of the critical Reynolds number field across the entire region. During the pumping test, the platform receives updated measured zeta potential values ​​from each well point at fixed time steps (e.g., 10 seconds) and triggers a recalculation of the critical threshold for local grid cells. The recalculation range dynamically expands based on the pore water chemical diffusion radius. Estimated by the following formula:

[0042]

[0043] in, is the ion diffusion coefficient (unit: m² / s). The time since the ζ-potential change occurred (in seconds). This mechanism ensures that the spatial propagation effect of potential perturbations is accurately captured.

[0044] The bidirectional coupling iterative mechanism comprises two nested loops. The outer loop performs pumping response fitting, using the regional flow field output by the inner loop as a constraint, adjusting the local permeability coefficient around the pumping well until the sum of squared residuals between the simulated drawdown curve and the measured curve converges to a preset tolerance (e.g., ...). (m²). The inner loop, based on the updated permeability coefficient field, resolves the global nonlinear seepage control equations and simultaneously refreshes the critical Reynolds number thresholds for each grid cell. When the actual Reynolds number in a grid cell exceeds its current threshold, the system automatically activates the Forchheimer inertia term; otherwise, it maintains the linear Darcy flow mode. This switching logic is controlled by a state machine, with the state transition condition being the real-time comparison between the Reynolds number and the critical threshold. Each transition triggers a Jacobian matrix reconstruction to ensure the numerical stability of the nonlinear equations. Jacobian matrix. elements Defined as:

[0045]

[0046] in, For the first The residuals of the governing equations For the first An unknown variable (such as pressure or velocity component).

[0047] The generation of the spatial distribution map of baseflow contribution rate relies on the source-sink term decomposition module. Based on the finally converged flow field, this module performs mass conservation integrals over each river confluence node to separate the instantaneous recharge fluxes formed by precipitation infiltration, surface recharge, and deep confined water overflow. Simultaneously, using the particle tracking algorithm, it integrates backwards along streamlines to identify the upstream aquifer units corresponding to each recharge flux and calculates their proportion in the total baseflow. Throughout this process, all streamline calculations employ the fourth-order Runge-Kutta scheme, using a dynamically updated seepage velocity field as input to ensure accurate reflection of streamline tortuosity effects outside the Darcy zone. The iterative formula for the fourth-order Runge-Kutta scheme is:

[0048]

[0049]

[0050]

[0051]

[0052]

[0053] in, For the first Step particle position (unit: m). For the corresponding time (unit: seconds). The time step (in seconds). Let be the velocity field function.

[0054] All in-situ sensing units are deployed in dedicated monitoring wells with a double-casing structure. The inner casing is used for pumping operations, while the outer casing is annularly filled with highly permeable quartz sand and fitted with a stop plug to isolate vertical interference between different aquifers. The in-situ sensing units are fixed to the bottom of the inner casing via flange interfaces, with their microfluidic chip inlets facing the main flow direction to ensure that the sampled fluid represents the actual seepage path. Power is supplied by a downhole supercapacitor energy storage module, periodically charged by a ground-based photovoltaic array via armored cables, supporting continuous operation for at least 30 calendar days.

[0055] A closed-loop feedback link based on power line carrier communication is established between the regional hydrodynamic simulation platform and the field pumping test device. This link allows the simulation platform to proactively send adjustment commands to the corresponding pumping wells when it detects that the seepage state in a certain area is about to cross the non-Darcy critical point, temporarily changing the pumping rate to avoid data distortion in the strongly nonlinear region. Simultaneously, the in-situ sensing unit at the pumping well end can automatically trigger the platform's global parameter refresh process when it detects a sudden change in zeta potential exceeding a preset safety threshold (e.g., ±50mV), preventing local electrochemical disturbances from causing overall model instability.

[0056] The non-Darcy flow criterion module incorporates a hysteresis compensation mechanism to handle the non-instantaneous relationship between zeta potential changes and seepage response. This mechanism constructs a first-order inertial element transfer function based on measured historical data.

[0057]

[0058] in, For gain (dimensionless). The time constant (unit: s). This is a Laplace variable. The current zeta potential input is weighted by a time constant, and the output is the effective control signal. The value of the time constant is determined based on pore water chemical equilibrium kinetic experiments and is dynamically adjusted with temperature and ion species. For example, in a Na⁺-dominant system... s, in the Ca²⁺-dominant system s, for every 10°C increase in temperature, Shortened by 15%.

[0059] After the regional flow simulation is completed, the system automatically executes the sensitivity analysis module. This module freezes each driving factor (such as precipitation intensity, evaporation rate, and artificial extraction) sequentially and reruns the simplified model to quantify the independent contribution of each factor to the baseflow output. All sensitivity tests are conducted under corrected non-Darcy flow boundary conditions to avoid attribution bias caused by missing physical mechanisms. The final output runoff driving factor ranking table is arranged in descending order of contribution rate and indicates the key aquifer unit numbers involved in the action path of each factor.

[0060] The technical effects of the present invention are further illustrated below through three specific embodiments and three comparative examples.

[0061] Example 1: A field experiment was conducted in the frontal region of an alluvial fan in the North China Plain. Five pumping wells were deployed, spaced 500 m apart and 30 m deep, penetrating a Holocene silty sand aquifer. An integrated in-situ sensing unit uploaded zeta potential, pH, conductivity, and flow rate data every 5 seconds. The initial groundwater pH was 7.2, ionic strength was 0.015 mol / L, and zeta potential was −32 mV. The pumping rate was 0.5 L / s for 48 hours. The regional hydrodynamic model used a triangular unstructured mesh with an average cell size of 100 m and a total of 12,000 cells. The electro-current coupled response function library contained 100 pre-calibrated curves. The bidirectional coupled iteration converged to the 7th outer cycle, with a residual sum of squares of 8.7 × 10−5 m². Ultimately, the baseflow was identified as mainly originating from the northwest precipitation infiltration zone, contributing 68.3%, followed by the southern irrigation and recharge zone (22.1%) and deep confined water overflow (9.6%).

[0062] Example 2: An experiment was conducted in a coastal aquifer in the Yangtze River Delta affected by seawater intrusion. Four pumping wells, 25 m deep, were installed in the brackish water transition zone. The initial zeta potential was −18 mV (due to double-layer compression caused by increased Cl⁻ and Na⁺ concentrations). During pumping, the zeta potential was monitored to decrease from −18 mV to −8 mV within 6 hours, indicating a significant decay of interfacial charge. The model dynamically updated the critical Reynolds number field, finding that the non-Darcy flow region expanded from the initial 3.2% to 12.7%. Inversion results showed that the traditional Darcy model overestimated the permeability coefficient by approximately 23%, while the parameter field obtained by the method of this invention had an error of less than 5% compared to the core permeability test results. In the baseflow contribution rate, density flow driven by seawater intrusion accounted for 15.4%, far higher than the 3.1% estimated by the traditional method.

[0063] Example 3 involved experiments conducted in a fractured-porous medium zone of the southwestern red bed. The lithology consisted mainly of mudstone interbedded with thin sandstone layers, and the pores were predominantly nanoscale clay pores. Carboxyl groups were grafted onto the inner walls of the microfluidic chip pores to simulate the surface of montmorillonite. The initial zeta potential was −45 mV, slowly increasing to −38 mV as pumping progressed (due to Ca²⁺ desorption). The model employed a hysteresis compensation mechanism with a time constant of 280 s. Particle tracking revealed that 78% of the baseflow originated from residual rainwater during the rainy season, exhibiting a distinct non-Darcy characteristic when flowing through the low-permeability mudstone matrix. Ignoring this effect would cause the streamlines to be incorrectly biased towards the high-permeability fracture zone, resulting in a recharge zone shift of up to 1.2 km.

[0064] Comparative Example 1 used the same site and pumping scheme as Example 1, but disabled the dynamic update function of the zeta potential and fixed the critical Reynolds number to 10⁻⁴ / 10⁻⁴. The inverted permeability field showed false high-value bands, which did not match the geological profile. The baseflow path was misjudged as being distributed along the eastern paleochannel, and the tracer in the actual monitoring wells was not detected. The contribution rate error reached ±31%.

[0065] Comparative Example 2 uses the traditional Darcy flow model in the scenario of Example 2, without introducing the Forchheimer term. Although the drawdown curve is fitted (residual sum of squares 9.1 × 10⁻⁵ m²), the inverted permeability coefficient is 27% higher than the true value, and the nonlinear decay of the flow rate in the later stages of pumping cannot be explained. The contribution of seawater intrusion in the baseflow is completely ignored.

[0066] In Comparative Example 3, the microfluidic chip used smooth straight channels (without fractal structure or functional group modification), causing the measured zeta potential to deviate from the actual interface state by approximately 15 mV. The resulting response function failed, and the critical Reynolds number update lagged by 6 hours. The particle tracking streamline deviated significantly from the actual tracer path, with a maximum offset distance of 1.8 km.

[0067] The table below summarizes the key performance indicators of each embodiment and the comparative example:

[0068] Case Number Scene type Enable dynamic potential updates? Does it use a biomimetic microfluidic chip? Permeability coefficient inversion error Baseflow path location error (km) Non-Darcy area identification accuracy Example 1 Alluvial fan yes yes <5% 0.12 94.7% Example 2 Seawater intrusion yes yes <5% 0.18 91.3% Example 3 Red layer medium yes yes <6% 0.21 89.5% Comparative Example 1 Alluvial fan no yes 23% 1.05 42.1% Comparative Example 2 Seawater intrusion no yes 27% 1.37 38.6% Comparative Example 3 Red layer medium yes no 19% 1.80 53.2%

[0069] The above data indicate that high-precision aquifer parameter inversion and baseflow path identification can only be achieved when zeta potential dynamic sensing, biomimetic microfluidic chips, and non-Darcy flow adaptive criteria are simultaneously enabled. This invention establishes a cross-scale physically consistent seepage modeling paradigm by coupling interfacial electrodynamic behavior at the nanopore scale with regional hydrodynamic processes, providing a reliable technical path for the quantitative analysis of runoff driving factors under complex hydrogeological conditions.

[0070] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for identifying runoff driving factors that integrates field pumping tests and hydrodynamic simulation, characterized in that, Includes the following steps: Several pumping test well groups were set up in the target watershed. Each pumping well was equipped with an integrated in-situ sensing unit. The integrated in-situ sensing unit was used to simultaneously acquire time-series data of groundwater pH value, ionic strength, colloidal concentration and zeta potential of pore wall, and record water level drawdown and pumping flow rate. Based on the time-series data, a nanopore-scale electro-current coupling response function is constructed, which characterizes the dynamic mapping relationship between the zeta potential and the non-Darcy critical Reynolds number threshold. The electro-current coupling response function is embedded into the non-Darcy flow criterion module of the regional hydrodynamic simulation platform, and a critical Reynolds number threshold is dynamically assigned to each computational grid cell. The pumping response inversion and regional flow simulation are coupled in a two-way iterative process: the outer circulation uses the measured drawdown curve as a constraint to adjust the local permeability coefficient until the simulation residual converges; the inner circulation solves the global nonlinear seepage control equation based on the updated permeability coefficient field, and switches between Darcy flow or Forchheimer non-Darcy flow mode in real time according to the comparison between the actual Reynolds number of each grid cell and its dynamic critical threshold. After iterative convergence, source-sink term decomposition and particle tracking are performed based on the final flow field to generate a spatial distribution map of the baseflow contribution rate and a ranking table of multi-source runoff driving factors.

2. The method for identifying runoff driving factors integrating field pumping tests and hydrodynamic simulation as described in claim 1, characterized in that, The integrated in-situ sensing unit includes a waterproof encapsulation shell, a multi-channel microfluidic chip, a double-layer potential sensor array, a piezoresistive level gauge, and an electromagnetic flow meter. The multi-channel microfluidic chip features a biomimetic nanopore network with a pore size distribution of 2–100 nanometers. The pore topology is generated based on fractal dimension matching of CT images of the target aquifer core, and carboxyl or amino functional groups are grafted onto the inner walls of the pores to simulate the acid-base active sites on the surface of natural minerals. The double-layer potential sensor array uses a gold-indium tin oxide composite electrode pair and outputs a phase difference signal that is monotonically mapped to the zeta potential using AC impedance spectroscopy. The piezoresistive level gauge and the electromagnetic flow meter are connected to a central data acquisition unit via an RS485 bus, with a time synchronization accuracy better than 10 milliseconds.

3. The method for identifying runoff driving factors by integrating field pumping tests and hydrodynamic simulation as described in claim 1, characterized in that, The process of constructing the nanoporous scale current-coupled response function includes: Using measured zeta potential, fluid conductivity, and dielectric constant as inputs, the expression for the effective slip length is derived based on the Gouy-Chapman-Stern double-layer theory. ; Substituting the effective slip length into the wall boundary conditions of the Navier-Stokes equation, and then dimensionlessly converting it, yields the modified Darcy-Forchheimer type seepage equation. Among them, the inertial drag coefficient , Determined by microfluidic calibration experiments; Numerical solution of the nonlinear relationship between seepage velocity and pressure gradient under different zeta potentials, and extraction of Reynolds number. equal The critical velocity corresponding to the time is used as the criterion for the initiation of non-Darcy flow under the current operating conditions.

4. The method for identifying runoff driving factors by integrating field pumping tests and hydrodynamic simulation as described in claim 1, characterized in that, The regional hydrodynamic simulation platform employs an unstructured grid discrete architecture, with each grid cell bound to an independent set of seepage constitutive parameters, including the hydraulic conductivity coefficient, storage coefficient, and dynamic critical Reynolds number threshold. During pumping tests, the platform receives updated zeta potential values ​​from each well point at fixed time steps and triggers recalculation of the critical threshold for the local grid. The recalculation range is based on the pore water chemical diffusion radius. Dynamic expansion, where D is the ion diffusion coefficient and t is the time since the change in self-potential occurred.

5. The method for identifying runoff driving factors integrating field pumping tests and hydrodynamic simulation as described in claim 1, characterized in that, In the bidirectional coupling iterative mechanism, when the actual Reynolds number of a certain grid cell exceeds its dynamic critical threshold, the system activates the Forchheimer inertial term; otherwise, it maintains the Darcy flow mode. The flow state switching is controlled by a state machine, and each state transition triggers Jacobian matrix reconstruction to ensure the numerical stability of the solution of the nonlinear equations.

6. The method for identifying runoff driving factors integrating field pumping tests and hydrodynamic simulation as described in claim 1, characterized in that, The generation of the spatial distribution map of the baseflow contribution rate includes: performing mass conservation integrals on the converging flow field at the river confluence nodes to separate the instantaneous recharge fluxes of precipitation infiltration, surface recharge, and deep confined water overflow; and using a fourth-order Runge-Kutta scheme to perform particle tracking along streamlines to identify the upstream aquifer units corresponding to each recharge flux and to calculate their proportion in the total baseflow; wherein the streamline calculation uses a dynamically updated seepage velocity field as input to accurately reflect the streamline bending effect in the non-Darcy region.

7. The method for identifying runoff driving factors integrating field pumping tests and hydrodynamic simulation as described in claim 2, characterized in that, The phase difference signal output by the double-layer potential sensor array is converted from analog to digital, and then the phase offset is extracted by the sliding window Fourier transform performed by the embedded signal processor. The phase offset is then input into a pre-trained neural network mapper. The mapper establishes a phase-ζ potential nonlinear mapping relationship based on calibration experimental data and outputs an effective ζ potential value with an update frequency of not less than once per second.

8. The method for identifying runoff driving factors integrating field pumping tests and hydrodynamic simulation as described in claim 1, characterized in that, A closed-loop feedback link based on power line carrier communication is established between the regional hydrodynamic simulation platform and the field pumping test device: when the platform predicts that the seepage state in a certain area is about to cross the non-Darcy critical point, it sends a command to the corresponding pumping well to temporarily adjust the pumping rate; when the in-situ sensing unit detects that the ζ potential change exceeds the preset safety threshold, it automatically triggers the platform's global parameter refresh process.

9. The method for identifying runoff driving factors integrating field pumping tests and hydrodynamic simulation as described in claim 1, characterized in that, The non-Darcy flow criterion module incorporates a hysteresis compensation mechanism, which is based on the transfer function of a first-order inertial element. The ζ potential input is time-weighted to output an effective control signal; the time constant τ is determined based on the chemical equilibrium kinetics experiment of pore water and is dynamically adjusted with temperature and ion type.

10. The method for identifying runoff driving factors integrating field pumping tests and hydrodynamic simulation as described in claim 1, characterized in that, After the regional flow simulation is completed, the system performs sensitivity analysis: the driving factors such as precipitation intensity, evaporation rate or artificial extraction are frozen in sequence, and the simplified model is rerun to quantify the independent contribution of each factor to the baseflow output; all sensitivity tests are conducted under corrected non-Darcy flow boundary conditions, and the final output is a ranking table of runoff driving factors arranged in descending order of contribution rate, and the key aquifer unit numbers involved in the action path of each factor are marked.

Images