One-dimensional aquifer hydraulic parameter joint inversion method for coupling steady flow field and transient disturbance

Abstract

The present application relates to a kind of one-dimensional aquifer hydraulic parameter joint inversion method of coupling steady flow field and transient disturbance.This method constructs a coupling seepage physical model containing natural steady background flow field and artificial transient disturbance simultaneously, and establishes the corresponding mathematical definite problem;Select dimensionless water storage rate and dimensionless time and other dimensionless factors to dimensionless the definite problem, analytically solve the definite problem after dimensionless using Laplace transform method, and obtain the analytical expression of head response with dimensionless water storage rate as core;Based on the solution, a series of theoretical standard curves are generated by programming calculation;The measured data of water head change with time collected by test is matched with standard curve cluster, and the permeability coefficient and water storage rate of aquifer, two key hydraulic parameters, are determined simultaneously, quickly and accurately by joint inversion.This application has the advantages of simple operation, high efficiency, good precision, rich information and the like.

Description

Technical Field

[0001] This invention belongs to the technical field of hydrogeological parameters and groundwater seepage analysis and measurement, specifically involving a joint inversion method for one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbances. This method establishes a frequency domain model through Laplace transform, obtains analytical and semi-analytical solutions by combining the residue theorem and Stehfest numerical inverse transform, and generates a standard curve family using Python programming to achieve simultaneous inversion of permeability coefficient and storage capacity based on experimental data. Background Technology

[0002] Permeability coefficient and storage capacity are two core parameters characterizing the hydraulic properties of aquifers. Accurate measurement of these parameters is fundamentally crucial for the sustainable development of groundwater resources, simulation of pollutant transport, prevention of land subsidence, and stability design and assessment of various geotechnical engineering projects (such as foundation reinforcement and pit dewatering). Currently, obtaining these parameters mainly relies on constant head and variable head tests conducted indoors or in the field. However, these methods have significant limitations. Traditional constant (variable) head tests typically only obtain the permeability coefficient and cannot reflect storage characteristics. While traditional micro-water tests can estimate both parameters simultaneously, their theoretical models are usually based on the assumption of an initial static flow field, ignoring the background hydraulic gradient that is prevalent under natural conditions, and cannot achieve effective coupling solutions for natural stable flow fields and artificially applied instantaneous disturbances. This challenges the applicability and accuracy of traditional methods in sites with significant background gradients. More importantly, micro-water tests rely solely on main well water level recovery data for storage capacity determination, failing to capture the early transient pressure propagation process at internal observation points, resulting in significant uncertainties and multiple solutions. Furthermore, in the field of geotechnical engineering, in-situ determination of the vertical permeability coefficient is a common requirement. Current dual-tube or other in-situ testing equipment is often complex to operate and time-consuming. Therefore, there is an urgent engineering need to develop a method that is easy to operate, fast to test, and can accurately reflect the vertical seepage characteristics.

[0003] This study aims to address the aforementioned theoretical and technical challenges. By establishing a novel physical and mathematical model coupling a stable background flow field with anthropogenic instantaneous disturbances, and developing a corresponding joint inversion method, simultaneous and accurate determination of permeability and storage capacity under real hydrogeological conditions is achieved. This method requires only a single experiment to jointly invert the two parameters, significantly improving the theoretical and technical level of parameter identification. Furthermore, due to its one-dimensional vertical-direct current field testing characteristics, it provides a simpler and more efficient alternative to traditional dual-pipe equipment for in-situ rapid determination of vertical permeability in geotechnical engineering, demonstrating significant application prospects in fields such as foundation reinforcement effect evaluation. Summary of the Invention

[0004] The purpose of this invention is to provide a one-dimensional aquifer hydraulic parameter joint inversion method that couples a stable flow field with instantaneous disturbances. This method overcomes the limitations of traditional constant-head or variable-head tests that cannot simultaneously consider both the natural background flow field and anthropogenic disturbances, and breaks through the technical bottleneck of existing methods' inaccurate determination of storage capacity. This invention aims to solve the technical problem of simultaneously and accurately determining permeability coefficient and storage capacity under natural hydraulic gradient conditions. By establishing an innovative coupling model and inversion method, two key hydrogeological parameters can be obtained jointly in a single test. This not only significantly improves the accuracy of parameter measurement but also provides a new technical means for the reliable determination of storage capacity, greatly enhancing the method's engineering applicability and theoretical value.

[0005] To solve the above-mentioned technical problems, the technical solution adopted by the present invention is as follows:

[0006] A method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a steady flow field and instantaneous disturbances includes the following steps:

[0007] S1: Establish a one-dimensional vertical seepage mathematical model. Based on Darcy's law and the law of conservation of mass, establish the mathematical boundary value problem of the coupled seepage system. This boundary value problem includes the head control equation, the initial condition characterizing the linearly distributed head, the lower boundary constant head condition, the upper boundary coupled disturbance condition, and a dynamic coupling equation characterizing the difference between the disturbance change rate and the boundary hydraulic gradient.

[0008] S2: The mathematical boundary value problem is made dimensionless by introducing a set of dimensionless factors, eliminating the original parameter μ. s ,K,L,r c ,r s The dimensional effects of H1, H2, and w0 were determined, and the Laplace transform method was used to obtain the dimensionless storage ratio μ. sD The analytical and semi-analytical solutions for the dimensionless head response are given as the key variables.

[0009] S3: Based on the analytical and semi-analytical solutions of the dimensionless head response, computer programming is used to generate solutions at different dimensionless storage saturations μ. sD A set of standard curves relating dimensionless head to dimensionless time for the given values;

[0010] S4: Construct a one-dimensional vertical seepage physical model. The lower boundary of the model maintains a constant head H2, while the upper boundary, based on the initial background head H1, is subject to an instantaneous head disturbance w(t), thus creating a coupled seepage system that simultaneously incorporates the natural stable background flow field and the anthropogenic instantaneous disturbance. In physical model tests or field tests, measure and record the measured data of the head change over time at observation points within the test section after the application of the instantaneous head disturbance.

[0011] S5: Analyze and process the experimental data, match the measured data curve with the standard curve, and determine the optimal dimensionless water storage ratio μ by finding the standard curve with the highest good of fit. sD Value and dimensionless time t D Values; based on the optimal matching value and the definition of the dimensionless factor, the permeability coefficient K and storage capacity μ of the aquifer are obtained by simultaneously solving the equations. s The specific value.

[0012] The above method overcomes the limitation of traditional constant head or variable head test methods, which can only obtain one parameter in a single test. It has the advantages of simple operation, high efficiency, good accuracy, and rich information.

[0013] The model construction and corresponding mathematical boundary value problem in step S1 above are described as follows:

[0014] The spatiotemporal evolution of water head in porous media is governed by the classical diffusion equation, which describes the flow and storage process of water under the influence of a pressure gradient:

[0015]

[0016] A key feature of this model is that it considers a pre-existing steady flow field to simulate the natural background flow field, which makes the initial conditions more realistic and more complex.

[0017]

[0018] The model bottom is connected to a regional aquifer or large body of water, maintaining a constant head boundary.

[0019] H(L,t)=H2,t≥0(1c)

[0020] The upper boundary condition organically combines natural background values ​​with anthropogenic disturbances:

[0021] H(0,t)=w(t)+H1,t≥0(1d)

[0022] This model innovatively proposes the aforementioned dynamic boundary condition equation, which strictly follows the law of mass conservation and, combined with Darcy's law, achieves quantitative coupling between the steady flow field and the dynamic disturbance at the upper boundary. The water conservation equation and corresponding initial conditions are as follows:

[0023]

[0024] w(0)=w0(1f)

[0025] This equation strictly adheres to the principle of mass conservation. The left side of the equation characterizes the rate of change of excess water volume caused by boundary disturbances; the right side quantifies the net flux contribution at the boundary, which physically represents the flux inherent in the natural background steady-state flow field subtracted from the current actual total flow. The difference between the two precisely characterizes the excess flow rate, which deviates from the background steady state and is entirely excited by transient disturbances. This flow rate directly drives the dynamic evolution of the boundary head. This equation plays a crucial coupling role in the mathematical model, not only clearly distinguishing between the natural background flow field and the transient head disturbance effect, but also establishing an organic coupling relationship between the two. This theoretical mechanism lays the core foundation of the method of this invention, enabling the joint inversion of permeability coefficient and storage capacity through a single excitation experiment.

[0026] In the formula: H and H(z,t) both represent the water head (L) at a certain point within the aquifer; μ s The water storage ratio (L) characterizes the water storage capacity of an aquifer. -1 K is the permeability coefficient, reflecting the permeability of porous media (LT). -1 ); z represents the vertical coordinate (L); L is the characteristic length of the seepage system, i.e., the seepage path (L) of the model; r c r is the radius (L) of the upper boundary water column; s Let be the radius of the aquifer sand column (L); t represents time (T). H1 and H2 represent the background constant heads (L) at the upper and lower boundaries of the model, respectively; w(t) is a time-varying perturbation term (L) applied based on the background head H1 at the upper boundary, with its initial value w0 being the initial amplitude (L) of the perturbation.

[0027] The units used in this application are expressed in international dimensions.

[0028] Furthermore, in S2 above, when the original boundary value problem is made dimensionless, a set of dimensionless factors are introduced to eliminate the influence of specific geometric scales and perturbation amplitudes, thereby obtaining a standardized mathematical model with universal applicability.

[0029] The dimensionless factor includes:

[0030]

[0031] By applying the dimensionless factor to equation -, the original boundary value problem is transformed into the following dimensionless boundary value problem:

[0032]

[0033] Among them, z D Using dimensionless spatial coordinates, the physical domain [0,L] is mapped to the standard element [0,1]; t D For dimensionless time; w DH is a dimensionless change in boundary head. D H D (z D ,t D All are dimensionless heads. This dimensionless factor innovatively incorporates the background stable flow field head distribution, eliminating the influence of the initial linear head distribution and greatly simplifying the entire boundary value problem; μ sD The dimensionless water storage rate, which integrates the water storage capacity of the medium and the geometric scale effect, is a key parameter for the dynamic response of the control system.

[0034] This dimensionless boundary value problem, by introducing the above five dimensionless factors, transforms the original problem, which contained multiple characteristic physical parameters (μ... s ,K,L,r c ,r s The complex problem of H1, H2, and w0 is simplified to a problem involving only a single parameter μ. sD The standardization of control. This simplification not only makes subsequent theoretical solutions possible, but also greatly facilitates the generation and matching of standard curves in the parameter inversion process.

[0035] Furthermore, the analytical solution to the dimensionless boundary value problem in S3 is achieved through the Laplace transform method, specifically including the following process:

[0036] The dimensionless head diffusion equation with respect to the dimensionless time t D Performing a Laplace transform, we obtain the ordinary differential equations in the Laplace domain:

[0037]

[0038] The general solution to this equation can be expressed in hyperbolic function form:

[0039]

[0040] Using lower boundary conditions We can obtain A(s) = 0, so the general solution simplifies to:

[0041]

[0042] Using upper boundary conditions Establish relations in the Laplace domain:

[0043]

[0044] Thus, the explicit expression for the Laplace domain head distribution is obtained:

[0045]

[0046] Taking the Laplace transform of the dimensionless water conservation equation, we get:

[0047]

[0048] By taking the derivative, we can obtain:

[0049]

[0050] By combining the equations, we can obtain information about... The closing expression:

[0051]

[0052] Substituting the equation into the expression, we obtain the analytical expression for the dimensionless head in the Laplace domain:

[0053]

[0054] Through analysis The zeros in the denominator can determine the location of its poles, and the characteristic equation can be obtained through analysis:

[0055]

[0056] The transcendental equation has infinitely many positive roots in the interval (0,∞). Therefore, all the poles of the system in the Laplace domain are:

[0057]

[0058] These poles determine the decay modes and expansion forms of the time-domain solution, and are key to the construction of subsequent series solutions.

[0059] Continuing with the inverse Laplace transform using the residue theorem, first calculate the denominator function:

[0060]

[0061] At the pole s n The derivative at:

[0062]

[0063] By simplifying using trigonometric identities, we finally obtain the explicit form of the residue:

[0064]

[0065] Applying the residue theorem, the system response can be expressed as a series analytical solution in the time domain:

[0066]

[0067] In the formula: Both represent dimensionless head H DThe Laplace transform of; Both represent dimensionless boundary perturbations w D The Laplace transform of ; s is the complex frequency variable in the Laplace transform; A(s) and B(s) are the undetermined coefficient functions in the general solution, determined by the boundary conditions; D(s) is the denominator function; θ is the auxiliary variable introduced by solving the characteristic equation; θ n It is an eigenvalue, representing the nth positive root of the characteristic equation; s n yes At the poles of the complex plane; c n At the extreme point s n The residues calculated at point n represent the magnitude weights of the nth mode.

[0068] In this way, the originally complex coupled partial differential equation problem was transformed into an eigenvalue problem, and a complete analytical solution expression was given through series expansion. This result lays a solid theoretical foundation for subsequent experimental data matching and aquifer parameter inversion using standard curves.

[0069] Furthermore, in S3, as a replacement or supplement to the analytical solution method, the Laplace domain solution is obtained. and Then, the semi-analytical solution can be obtained directly using the Stehfest numerical inverse transform algorithm. The specific process is as follows:

[0070] For any time-domain function F(t) to be determined D Its Laplace transform is Then at a specific time t D The semi-analytical solution can be calculated by the following formula:

[0071]

[0072] Where the coefficient V i Determined by the following formula:

[0073]

[0074] In the formula: F(t) D Let be the time-domain function to be determined; Its Laplace transform; V i is the i-th weight coefficient; N is the total number of points in the Stehfest algorithm (an even number), which determines the calculation precision and the number of samples; i is the index variable in the summation loop, which takes the value of a continuous positive integer between 1 and N, and is used to traverse all sample points; k is the auxiliary index variable in the summation formula when calculating the weight coefficient, and the meanings of the other symbols are the same as above.

[0075] This solution method relies solely on the functional form of the Laplace domain solution, avoiding the complex process of solving the characteristic equation and its characteristic roots in analytical methods. It also eliminates the need for truncation and summation of infinite series, making it applicable to any situation where the Laplace domain solution can be obtained. The comprehensive application of this method enhances its universality and practicality.

[0076] Furthermore, the standard curve generation step in S3 is implemented through computer programming, specifically using Python for numerical calculation and visualization. The process includes:

[0077] First, the standard curve is plotted based on the time-domain solution expression obtained from the residue theorem, using the dimensionless storage ratio μ. sD Dimensionless distance z D Dimensionless time t D and characteristic roots θ n For key parameters, calculate the analytical solution in series form and sum them to obtain the boundary disturbance response w. D and head distribution H D The variation curve is then used to generate a standard curve family covering the predetermined parameter range. Simultaneously, the derived Laplace domain solution is directly used to convert the frequency domain solution to the time domain solution via the Stehfest numerical inverse transform algorithm. By adjusting the preset algorithm parameters, the boundary disturbance response w under different dimensionless parameter values ​​is calculated. D and head distribution H D The variation curve is then used to generate a corresponding standard curve family. Both methods include a data storage step, which stores the generated standard curve data (μ) sD z D t D Sequence and corresponding w D and H D The values ​​are stored as structured data files, and a standard curve chart is generated by calling a Python visualization library. By comparison, it can be found that the analytical solution obtained by the residue theorem and the semi-analytical solution obtained by the Stehfest algorithm fit the curves well under the same dimensionless parameters. This also verifies the applicability and accuracy of the two solution methods, providing a strong theoretical and numerical foundation for subsequent parameter inversion.

[0078] Furthermore, the physical model (experimental platform) used in S4 includes a sand column body, a water level control system, and a data acquisition system. The sand column body is a transparent cylindrical design, with the aquifer section filled with homogeneous fine sand. The water level control system uses a peristaltic pump to stably supply water and connects a constant pressure water tank to the bottom of the sand column body via pipeline to ensure stable and constant water head conditions. The data acquisition system includes a multi-channel data acquisition instrument, a wireless data transmitter, and high-precision pore water pressure sensors (range 0-100 kPa, accuracy 0.1% FS) arranged in layers along the sand column axis. It can synchronously and in real-time monitor and record continuous data on water head changes over time at different observation points during the experiment, and further analyze and process the data.

[0079] Furthermore, the analysis and processing of the collected data in S5 includes: a standard curve storage unit, which pre-stores the dimensionless parameter μ. sD The system includes a standard curve data set; a data processing unit that receives and preprocesses the head change data over time from the acquisition system; a parameter matching unit that compares and matches the measured data with the standard curve, and determines the best-fit parameters using an optimization algorithm; and a parameter calculation unit that calculates and outputs the permeability coefficient K and storage capacity μ based on the matching results. s The final result.

[0080] Furthermore, in step S5, [L] and [r] are recorded based on the device size and matching results. c ],[r s ],[μ sD ] and [t D The permeability coefficient K and storage capacity μ of an aquifer are calculated using the following formulas. s :

[0081]

[0082] Any techniques not mentioned in this invention are based on existing technologies.

[0083] Compared with the prior art, the present invention has the following advantages:

[0084] (1) A breakthrough in the simultaneous joint measurement of hydrogeological parameters has been achieved: the limitations of traditional tests that can only measure permeability coefficient and water storage rate are not accurately identified. The two key parameters of permeability coefficient and water storage rate can be obtained simultaneously through one test, solving the long-standing technical problem of the difficulty in accurately measuring water storage rate and significantly improving test efficiency and data quality.

[0085] (2) Improved measurement accuracy under complex hydrogeological conditions: By establishing a mathematical model that couples stable flow field and instantaneous disturbance, the comprehensive influence of natural background hydraulic gradient is effectively considered, eliminating the systematic error generated by traditional methods under the condition of natural flow field, and significantly improving the accuracy and engineering applicability of parameter measurement results.

[0086] (3) It provides a rigorous theoretical foundation and diversified solution path: the dimensionless analytical solution is obtained through rigorous mathematical derivation, and the Stehfest numerical inverse transformation method is introduced to form a dual technical route of mutual verification between analytical solution and semi-analytical solution, which ensures the theoretical rigor of the method.

[0087] (4) The method simplifies the test operation and expands the scope of engineering applications: It only requires recording the water level changes in the observation well through a high-frequency sensor, without the need for precise flow measurement, which greatly simplifies the test operation. It is particularly suitable for the need for in-situ rapid determination of vertical permeability coefficient in geotechnical engineering, and provides a more convenient and reliable testing method for engineering applications such as foundation reinforcement and drainage system performance evaluation. Attached Figure Description

[0088] Figure 1 This is a flowchart illustrating the operation of the method of the present invention;

[0089] Figure 2 This is a schematic diagram of the experimental setup;

[0090] Figure 3 Water level change curves at each observation point

[0091] Figure 4 A diagram showing the alignment of the standard curve with the observation point data; Detailed Implementation

[0092] The present invention will be further illustrated below with reference to the accompanying drawings and specific embodiments. It should be understood that these embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. After reading this invention, any modifications of the invention in various equivalent forms by those skilled in the art will fall within the scope defined by the appended claims.

[0093] like Figure 1 As shown, this invention provides a method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbance, comprising the following steps:

[0094] 1) Establish a one-dimensional vertical seepage mathematical model, and construct a boundary value problem based on Darcy's law and the law of conservation of mass, including the head control equation, linear initial conditions, constant head boundary conditions, and dynamic coupling boundary conditions; eliminate the dimensional influence of multiple original parameters by introducing a dimensionless factor, and use the Laplace transform method to obtain the solution with dimensionless storage ratio μ. sD These are the analytical and semi-analytical solutions for the key variables.

[0095] The spatiotemporal evolution of water head in porous media is governed by the classical diffusion equation, which describes the flow and storage process of water under the influence of a pressure gradient:

[0096]

[0097] A key feature of this model is that it considers a pre-existing steady flow field to simulate the natural background flow field, which makes the initial conditions more realistic and more complex.

[0098]

[0099] The model bottom is connected to a regional aquifer or large body of water, maintaining a constant head boundary.

[0100] H(L,t)=H2,t≥0(1c)

[0101] The upper boundary condition organically combines natural background values ​​with anthropogenic disturbances:

[0102] H(0,t)=w(t)+H1,t≥0(1d)

[0103] This model strictly follows the law of conservation of mass and, combined with Darcy's law, achieves quantitative coupling between the steady flow field and the dynamic disturbance at the upper boundary; the water conservation equation and the corresponding initial conditions are as follows:

[0104]

[0105] w(0)=w0(1f)

[0106] This equation strictly adheres to the principle of mass conservation. The left side of the equation characterizes the rate of change of excess water volume caused by boundary disturbances; the right side quantifies the net flux contribution at the boundary, which physically represents the flux inherent in the natural background steady-state flow field subtracted from the current actual total flow. The difference between the two precisely characterizes the excess flow rate, which deviates from the background steady state and is entirely excited by transient disturbances. This flow rate directly drives the dynamic evolution of the boundary head. This equation plays a crucial coupling role in the mathematical model, not only clearly distinguishing between the natural background flow field and the transient head disturbance effect, but also establishing an organic coupling relationship between the two. This theoretical mechanism lays the core foundation of the method of this invention, enabling the joint inversion of permeability coefficient and storage capacity through a single excitation experiment.

[0107] In the formula: H and H(z,t) both represent the water head (L) at a certain point within the aquifer; μ s The water storage ratio (L) characterizes the water storage capacity of an aquifer. -1 K is the permeability coefficient, reflecting the permeability of porous media (LT). -1); z represents the vertical coordinate (L); L is the characteristic length of the seepage system, i.e., the seepage path (L) of the model; r c r is the radius (L) of the upper boundary water column; s Let be the radius of the aquifer sand column (L); t represents time (T). H1 and H2 represent the background constant heads (L) at the upper and lower boundaries of the model, respectively; w(t) is the time-varying perturbation term (L) applied based on the background head H1 at the upper boundary, with its initial value w0 being the initial amplitude (L) of the perturbation. All units used are in international units.

[0108] When the original boundary value problem is made dimensionless, a set of dimensionless factors are introduced to eliminate the influence of specific geometric scales and perturbation amplitudes, thereby obtaining a standardized mathematical model with universal applicability.

[0109] The dimensionless factor includes:

[0110]

[0111] By applying the dimensionless factor to equation -, the original boundary value problem is transformed into the following dimensionless boundary value problem:

[0112]

[0113] Among them, z D Using dimensionless spatial coordinates, the physical domain [0,L] is mapped to the standard element [0,1]; t D For dimensionless time; w D H is a dimensionless change in boundary head. D H D (z D ,t D All are dimensionless heads. This dimensionless factor innovatively incorporates the background stable flow field head distribution, eliminating the influence of the initial linear head distribution and greatly simplifying the entire boundary value problem; μ sD The dimensionless water storage rate, which integrates the water storage capacity of the medium and the geometric scale effect, is a key parameter for the dynamic response of the control system.

[0114] This dimensionless boundary value problem, by introducing the above five dimensionless factors, transforms the original problem, which contained multiple characteristic physical parameters (μ... s ,K,L,r c ,r s The complex problem of H1, H2, and w0 is simplified to a problem involving only a single parameter μ. sD The standardization of control. This simplification not only makes subsequent theoretical solutions possible, but also greatly facilitates the generation and matching of standard curves in the parameter inversion process.

[0115] The dimensionless head diffusion equation with respect to the dimensionless time t D Performing a Laplace transform, we obtain the ordinary differential equations in the Laplace domain:

[0116]

[0117] The general solution to this equation can be expressed in hyperbolic function form:

[0118]

[0119] Using lower boundary conditions We can obtain A(s) = 0, so the general solution simplifies to:

[0120]

[0121] Using upper boundary conditions Establish relations in the Laplace domain:

[0122]

[0123] Thus, the explicit expression for the Laplace domain head distribution is obtained:

[0124]

[0125] Taking the Laplace transform of the dimensionless water conservation equation, we get:

[0126]

[0127] By taking the derivative, we can obtain:

[0128]

[0129] By combining the equations, we can obtain information about... The closing expression:

[0130]

[0131] Substituting the equation into the expression, we obtain the analytical expression for the dimensionless head in the Laplace domain:

[0132]

[0133] Through analysis The zeros in the denominator can determine the location of its poles, and the characteristic equation can be obtained through analysis:

[0134]

[0135] The transcendental equation has infinitely many positive roots in the interval (0,∞). Therefore, all the poles of the system in the Laplace domain are:

[0136]

[0137] These poles determine the decay modes and expansion forms of the time-domain solution, and are key to the construction of subsequent series solutions.

[0138] Continuing with the inverse Laplace transform using the residue theorem, first calculate the denominator function:

[0139]

[0140] At the pole s n The derivative at:

[0141]

[0142] By simplifying using trigonometric identities, we finally obtain the explicit form of the residue:

[0143]

[0144] Applying the residue theorem, the system response can be expressed as a series analytical solution in the time domain:

[0145]

[0146] In the formula: Both represent dimensionless head H D The Laplace transform of; Both represent dimensionless boundary perturbations w D The Laplace transform of ; s is the complex frequency variable in the Laplace transform; A(s) and B(s) are the undetermined coefficient functions in the general solution, determined by the boundary conditions; D(s) is the denominator function; θ is the auxiliary variable introduced by solving the characteristic equation; θ n It is an eigenvalue, representing the nth positive root of the characteristic equation; s n yes At the poles of the complex plane; c n At the extreme point s n The residues calculated at point n represent the magnitude weights of the nth mode.

[0147] In this way, the originally complex coupled partial differential equation problem was transformed into an eigenvalue problem, and a complete analytical solution expression was given through series expansion. This result lays a solid theoretical foundation for subsequent experimental data matching and aquifer parameter inversion using standard curves.

[0148] In obtaining the Laplace domain solution and Then, the semi-analytical solution can be obtained directly using the Stehfest numerical inverse transform algorithm. For any time-domain function F(t) to be solved... D Its Laplace transform is Then at a specific time t D The semi-analytical solution can be calculated by the following formula:

[0149]

[0150] Where the coefficient V i Determined by the following formula:

[0151]

[0152] In the formula: F(t) D Let be the time-domain function to be determined; Its Laplace transform; V i is the i-th weight coefficient; N is the total number of points in the Stehfest algorithm (an even number), which determines the calculation precision and the number of samples; i is the index variable in the summation loop, which takes the value of a continuous positive integer between 1 and N, and is used to traverse all sample points; k is the auxiliary index variable in the summation formula when calculating the weight coefficient, and the meanings of the other symbols are the same as above.

[0153] This solution method relies solely on the functional form of the Laplace domain solution, avoiding the complex process of solving the characteristic equation and its characteristic roots in analytical methods. It also eliminates the need for truncation and summation of infinite series, making it applicable to any situation where the Laplace domain solution can be obtained. The comprehensive application of this method enhances its universality and practicality.

[0154] 2) The standard curve was generated using computer programming, specifically employing Python for numerical calculations and visualization. First, the standard curve was plotted based on the time-domain solution expression obtained from the residue theorem, using the dimensionless storage capacity μ. sD Dimensionless distance z D Dimensionless time t D and characteristic roots θ n For key parameters, calculate the analytical solution in series form and sum them to obtain the boundary disturbance response w. D and head distribution H D The variation curve is then used to generate a standard curve family covering the predetermined parameter range. Simultaneously, the derived Laplace domain solution is directly used to convert the frequency domain solution to the time domain solution via the Stehfest numerical inverse transform algorithm. By adjusting the preset algorithm parameters, the boundary disturbance response w under different dimensionless parameter values ​​is calculated. D and head distribution H D The variation curve is then used to generate a corresponding standard curve family. Both methods include a data storage step, which stores the generated standard curve data (μ) sD z D t D Sequence and corresponding w D and HD The values ​​are stored as structured data files, and a standard curve chart is generated by calling a Python visualization library. By comparison, it can be found that the analytical solution obtained by using the residue theorem and the semi-analytical solution obtained by the Stehfest algorithm fit the curves well under the same dimensionless parameters. This also verifies the applicability and accuracy of the two solution methods, providing a strong theoretical and numerical foundation for subsequent parameter inversion.

[0155] 3) The experimental platform includes a sand column body, a water level control system, and a data acquisition system. The sand column body adopts a transparent cylindrical design, with the aquifer section in the middle filled with homogeneous fine sand. The water level control system uses a peristaltic pump to stably supply water and connects a constant pressure water tank to the bottom of the sand column body through a pipeline to ensure stable and constant water head conditions. The data acquisition system includes a multi-channel data acquisition instrument, a wireless data transmitter, and high-precision pore water pressure sensors (range 0-100 kPa, accuracy 0.1% FS) arranged in layers along the sand column axis. It can synchronously monitor and record the continuous data of water head changes over time at different observation points during the experiment, and further analyze and process the data.

[0156] 4) The collected data is analyzed and processed, and its configuration includes: a standard curve storage unit, which pre-stores the dimensionless parameter μ. sD The system includes a standard curve data set; a data processing unit that receives and preprocesses the head change data over time from the acquisition system; a parameter matching unit that compares and matches the measured data with the standard curve, and determines the best-fit parameters using an optimization algorithm; and a parameter calculation unit that calculates and outputs the permeability coefficient K and storage capacity μ based on the matching results. s The final result.

[0157] 5) Record [L] and [r] according to the device dimensions and matching results. c ],[r s ],[μ sD ] and [t D The permeability coefficient K and storage capacity μ of an aquifer are calculated using the following formulas. s :

[0158]

[0159] To verify the actual effect of the method of the present invention, the above scheme is applied in this embodiment, and the specific experimental process is as follows:

[0160] 1. Assembly of the test apparatus and filling of the test medium

[0161] An acrylic cylinder is used as the main body of the one-dimensional percolator (e.g., Figure 2As shown, the inner diameter of the plexiglass cylinder is 0.3m and its height is 2.5m. The bottom of the plexiglass cylinder is connected to a constant pressure water tank (the constant pressure water tank has an outlet at the bottom) via a pipe; the top of the plexiglass cylinder is supplied with water stably via a peristaltic pump. Standard quartz sand (particle size of about 0.045mm) that has been dried and sieved is used as the test medium and is uniformly filled into the percolator using a layered compaction method. The filling height is 1.5m, and the density of each layer is controlled to be consistent during the filling process.

[0162] 2. Deployment of the data acquisition system and saturation of the test medium

[0163] With the top surface of the sand layer as the zero point (z = 0), and the vertical downward direction as the positive direction, four monitoring points are arranged along the z-axis: z = 0 to -0.1m (water column zone, monitoring inlet water pressure, corresponding to...). Figure 2 At z0, the water pressure sensor is used to monitor changes in water pressure at the inlet. The specific installation location of the sensor is not strictly limited, as long as the entire sensor is within the water column. The sensor is located in the upper water column of the seepage meter and along the vertical direction along the side wall of the seepage meter (0.3m, 0.75m, and 1.35m from the origin, corresponding to...). Figure 2 z in 1-3 High-precision pore water pressure sensors (range 0-100 kPa, accuracy 0.1% FS) were installed on each sample. All sensors were connected to a computer via a data acquisition unit, with a sampling frequency set to 0.1-5 s / time. The sample was saturated using a circulating water injection method until water continuously flowed from the top outlet and the readings of each sensor stabilized, indicating that the sample was fully saturated.

[0164] 3. Formal commencement of the experiment and recording of experimental data

[0165] The peristaltic pump flow rate and outlet height were adjusted to maintain stable head levels at the upper and lower boundaries. After seepage stabilized, the initial head distribution at each measuring point was recorded to confirm a linear distribution. A momentary disturbance was introduced at t=0 to raise the head, creating the upper boundary disturbance condition. The data acquisition system was started in advance to continuously record head changes at each pressure sensor until the head changes at each measuring point stabilized. The ambient temperature was kept constant at 20±1℃ throughout the experiment.

[0166] 4. Data processing and analysis, and parameter calculation

[0167] The collected raw head data is converted into a dimensionless head H by subtracting the initial linear distribution background value. D The data was processed using Python programming, matching the head time history data of each measuring point with the standard curve generated in advance through analytical solutions. Figure 3 ), to obtain the optimal matching parameters. Record [L], [r] c ],[r s ],[μ sD] and [t D The permeability coefficient K and storage capacity μ of an aquifer are calculated using the formula. s .

[0168] Experiments were conducted using the above-described apparatus and methods, and the results of a typical set of confirmatory experiments are as follows:

[0169]

[0170] The baseline value for the permeability coefficient in this experiment was 3.20 × 10⁻⁶. -4 m / s, the benchmark value for water storage rate is 1.50×10 -3 m -1 The permeability coefficient measured using the method of this invention is 3.00 × 10⁻⁶. -4 m / s, water storage rate is 1.53×10 -3 m -1 The results are in good agreement with the baseline values ​​obtained by traditional testing methods, with a relative error of 6.25% for the permeability coefficient and 2% for the storage rate. This not only verifies the reliability of the method of the present invention, but more importantly, it shows that while maintaining the measurement accuracy of traditional methods, the present method successfully solves the technical problem that traditional constant head and variable head test methods cannot simultaneously measure the storage rate.

Claims

1. A method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a steady flow field and instantaneous disturbance, characterized in that, Includes the following steps: S1: Establish a one-dimensional vertical seepage mathematical model. Based on Darcy's law and the law of conservation of mass, establish the mathematical boundary value problem of the coupled seepage system. This boundary value problem includes the head control equation, the initial condition characterizing the linearly distributed head, the lower boundary constant head condition, the upper boundary coupled disturbance condition, and a dynamic coupling equation characterizing the difference between the disturbance change rate and the boundary hydraulic gradient. S2: The mathematical boundary value problem is made dimensionless by introducing a set of dimensionless factors, eliminating the original parameters. μ s , K , L , r c , r s , H 1, H 2 and w The dimensionless influence of 0 was considered, and the Laplace transform method was used to obtain the dimensionless storage ratio. μ sD The dimensionless head response analytical and semi-analytical solutions are given for the key variables; among them, μ s The water storage rate characterizes the water storage capacity of an aquifer. K Permeability coefficient, which reflects the permeability of porous media; L This is the characteristic length of the seepage system, i.e., the seepage path of the model; r c The radius of the upper boundary water column; r s The radius of the aquifer sand column; H 1 and H 2 represents the background constant head at the upper and lower boundaries of the model, respectively; w 0 represents the background head at the upper boundary. H The initial amplitude of the time-varying disturbance term applied on the basis of 1; S3: Based on the analytical and semi-analytical solutions of the dimensionless head response, computer programming is used to generate solutions at different dimensionless storage rates. μ sD A set of standard curves relating dimensionless head to dimensionless time for the given values; S4: Construct a one-dimensional vertical seepage physical model, with a constant head maintained at the lower boundary of the model. H 2. The upper boundary is at the initial background water head. H Based on 1, apply an instantaneous head disturbance. w ( t This constructs a coupled seepage system that simultaneously includes a natural stable background flow field and an anthropogenic instantaneous disturbance; in physical model tests or field tests, the measured data of the change in head over time at the observation points within the test section after the instantaneous head disturbance is applied are measured and recorded. S5: Analyze and process the experimental data, match the measured data curve with the standard curve, and determine the optimal dimensionless water storage rate by finding the standard curve with the highest good of fit. μ sD Value and dimensionless time t D value; Based on the optimal matching value and the definition of the dimensionless factor, the permeability coefficient of the aquifer can be obtained by solving a series of simultaneous equations. K and water storage rate μ s The specific value.

2. The method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbance as described in claim 1, characterized in that, The process of establishing the one-dimensional vertical seepage mathematical model and the corresponding mathematical boundary value problem in step S1 is as follows: The spatiotemporal evolution of water head in porous media is governed by the classical diffusion equation, which describes the flow and storage process of water under the influence of a pressure gradient: (1a) The initial conditions are: (1b) The model bottom is connected to a regional aquifer or large body of water, maintaining a constant head boundary. (1c) The upper boundary condition organically combines natural background values ​​with anthropogenic disturbances: (1d) The water conservation equation and the corresponding initial conditions are as follows: (1e) (1f) In the formula: H , H ( z , t () All of these represent the water head at a specific point within the aquifer; z Represents the vertical coordinate; t Indicates time; w ( t (This refers to the background water head at the upper boundary.) H The time-varying perturbation term applied on the basis of 1, its initial value w 0 represents the initial amplitude of the disturbance.

3. The method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbance as described in claim 2, characterized in that, In S2, when making the mathematical boundary value problem dimensionless, a set of dimensionless factors are introduced; The dimensionless factor includes: (2a) By performing a dimensionless transformation on equations (1a)-(1f) using the aforementioned dimensionless factor, the original boundary value problem is transformed into the following dimensionless boundary value problem: (2b) in, z D For dimensionless spatial coordinates, the physical domain [0, L Mapped to standard cell [0, 1]; t D Time is dimensionless; w D This is a dimensionless change in boundary head. H D , H D ( z D , t D All are dimensionless heads; μ sD It is a dimensionless water storage rate.

4. The method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbance as described in claim 3, characterized in that, In S3, the analytical solution to the dimensionless boundary value problem is achieved through the Laplace transform method, which specifically includes the following process: The dimensionless head diffusion equation with respect to dimensionless time t D Performing a Laplace transform, we obtain the ordinary differential equations in the Laplace domain: (3a) The general solution to this equation can be expressed in hyperbolic function form: (3b) Using lower boundary conditions , can be obtained Therefore, the general solution simplifies to: (3c) Using upper boundary conditions Establish relations in the Laplace domain: (3d) Thus, the explicit expression for the Laplace domain head distribution is obtained: (3e) Taking the Laplace transform of the dimensionless water conservation equation, we get: (4a) By taking the derivative, we can obtain: (4b) Combining equations (4a) and (4b), we can obtain the following about The closing expression: (4c) Substituting equation (4c) into equation (3e) yields the analytical expression for the dimensionless head in the Laplace domain: (4d) Through analysis The zeros in the denominator determine the location of its poles, and the characteristic equation can be obtained through analysis: (5a) The characteristic equation is in Since there are infinitely many positive roots within the interval, we can deduce that all the poles of the system in the Laplace domain are: (5b) Continuing with the inverse Laplace transform using the residue theorem, we first calculate the denominator function: (6a) At the pole The derivative at: (6b) By simplifying using trigonometric identities, we finally obtain the explicit form of the residue: (6c) Applying the residue theorem, the system response can be expressed as a series analytical solution in the time domain: (7a) (7b) In the formula: Both represent dimensionless head. The Laplace transform of; Both represent dimensionless boundary perturbations. The Laplace transform of; It is the complex frequency variable in the Laplace transform; The undetermined coefficients in the general solution are determined by the boundary conditions; The denominator function; These are auxiliary variables introduced when solving the characteristic equation; It is an eigenvalue, representing the eigenvalue of the characteristic equation. n A positive root; yes At the poles of the complex plane; At the extreme point The residue calculated at point represents the number of points. n The amplitude weights of each mode.

5. The method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbance as described in claim 4, characterized in that, S3 serves as a replacement or supplement to the analytical solution method, in obtaining the Laplace domain solution. and Then, the semi-analytical solution was obtained directly using the Stehfest numerical inverse transform algorithm. The specific process is as follows: For any time-domain function to be determined Its Laplace transform is Then at a specific moment The semi-analytical solution can be calculated by the following formula: (8a) Where the coefficient Determined by the following formula: (8b) In the formula: Let be the time-domain function to be determined; Its Laplace transform; For the first i Each weighting coefficient; The total number of points in the Stehfest algorithm, taken as an even number, determines the computational precision and the number of samples. The index variable in the summation loop takes values ​​from 1 to... A series of consecutive positive integers between 1 and 2, used to iterate through all sampling points; This is an auxiliary index variable used in the summation formula when calculating the weighting coefficients.

6. The method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbance as described in claim 5, characterized in that, The standard curve generation step in S3 is implemented through computer programming, specifically using Python for numerical calculation and visualization. The process includes: First, the standard curve is plotted based on the time-domain solution expression obtained from the residue theorem, using the dimensionless storage ratio. μ sD Dimensionless distance z D Dimensionless time t D and characteristic roots θ n Calculate the analytical solution in series form and sum them to obtain the boundary disturbance response. w D and water head distribution H D The variation curves are then used to generate a standard curve family covering the predetermined parameter range; simultaneously, the derived Laplace domain solution is directly used to convert the frequency domain solution into the time domain solution through the Stehfest numerical inverse transform algorithm; and then, by adjusting the preset algorithm parameters, the boundary disturbance response under different dimensionless parameter values ​​is calculated. w D and water head distribution H D The variation curve is then used to generate a corresponding set of standard curves.

7. The method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbance as described in claim 5, characterized in that, The physical model used in S4 includes a sand column body, a water level control system, and a data acquisition system. The sand column body is a transparent cylindrical structure with a homogeneous fine sand filling the middle aquifer section. The water level control system uses a peristaltic pump to stably supply water and connects a constant pressure water tank to the bottom of the sand column body through a pipeline to ensure stable and constant water head conditions. The data acquisition system includes a multi-channel data acquisition instrument, a wireless data transmitter, and high-precision pore water pressure sensors arranged in layers along the height of the sand column body. It can synchronously monitor and record the continuous data of water head changes over time at different observation points during the experiment, and analyze and process the data.

8. The method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbance as described in claim 7, characterized in that, The S5 analyzes and processes the collected data, and its configuration includes: a standard curve storage unit with pre-stored dimensionless parameters. μ sD The system includes a standard curve data set; a data processing unit that receives and preprocesses the head change data over time from the acquisition system; a parameter matching unit that compares and matches the measured data with the standard curve, and determines the best-fit parameters using an optimization algorithm; and a parameter calculation unit that calculates and outputs the permeability coefficient based on the matching results. K and water storage rate μ s The final result.

9. The method for joint inversion of one-dimensional aquifer hydraulic parameters coupled with a stable flow field and instantaneous disturbance as described in claim 8, characterized in that, S5 records the data based on device dimensions and matching results. The permeability coefficient of an aquifer is calculated using the following formula. K With water storage rate μ s : (9a) (9b)

Images