One-dimensional soil column solute transport multi-parameter simultaneous inversion method based on finite difference method

Abstract

The application discloses a one-dimensional soil column solute transport multi-parameter synchronous inversion method based on a finite difference method, which comprises the following steps: obtaining the flow, water head and concentration-time data of a finite length one-dimensional soil column under the condition of instantaneous source injection, and performing non-dimensionalization on the concentration, space and time; establishing a mathematical model comprising a finite domain convection-dispersion control equation, a dynamic inlet boundary and an outlet zero concentration gradient boundary; solving the simulation concentration of each measuring point by using the finite difference method step by step; constructing a joint objective function according to the simulation concentration and the measured concentration of multiple tests, multiple measuring points and multiple time points; iteratively updating the to-be-inverted parameters under physical constraints until the convergence condition is met, and then outputting the porosity, pore water flow rate, longitudinal dispersion degree and derived parameters. The application is convenient for processing dynamic boundaries, parameter constraints and multiple test data, and improves the applicability and stability of parameter inversion.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation and parameter inversion technology of solute transport in soil and groundwater, and particularly to a method for synchronously inverting multiple hydraulic and solute transport parameters using finite difference calculations and multi-location concentration response data for instantaneous source injection tests of finite-length one-dimensional soil columns. Background Technology

[0002] One-dimensional soil column tracer tests are commonly used to study the convection, dispersion, and diffusion processes of solutes in porous media. The test typically involves applying a tracer solution of a certain concentration at the inlet of the soil column and recording the concentration changes over time at different locations or the outlet of the column. Subsequently, parameters such as porosity, flow velocity, and dispersion coefficient are identified using a convection-dispersion model.

[0003] For transient or attenuating sources formed by a finite-volume liquid supply container, the inlet concentration is not a constant value, but varies with the output of solute mass within the container. Treating the inlet as a constant concentration boundary or a pre-defined empirical function can easily isolate the inlet mass variation from the transport process within the soil column, and propagate boundary description errors to the parameter inversion results.

[0004] When considering dynamic inlet boundary, finite domain outlet boundary, multiple sets of experiments, multiple measurement point data and upper and lower limits of parameters, this application needs to be able to directly couple dynamic inlet mass balance with finite domain convection-dispersion equation, and be able to simultaneously identify the values ​​of multiple parameters using multiple sets of concentration data at multiple locations. Summary of the Invention

[0005] This invention aims to address the problems of insufficient coupling between the inlet boundary and the transport process within the soil column in the identification of instantaneous source solute transport parameters of existing one-dimensional soil columns, as well as the limited ability of single-point data to constrain multiple parameters. It provides a multi-parameter synchronous inversion method that uses the finite difference method as a solver and constructs a joint objective function using multiple sets of experiments and multiple measurement point observations.

[0006] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: A multi-parameter synchronous inversion method for solute transport in a one-dimensional soil column based on the finite difference method includes the following steps: 1) Construct a one-dimensional convection-dispersion model of the soil column under finite domain conditions. The model includes the dimensionless convection-dispersion governing equation, the initial concentration condition inside the soil column, the dynamic inlet boundary condition describing the change of solute mass in the supply water column, and the zero dispersion flux boundary condition at the outlet of the soil column. 2) Select the Peckley number Pe With dimensionless structural coefficient Using the core dimensionless parameter, the concentration, spatial, and temporal variables are dimensionless to obtain a dimensionless convection-diffusion mathematical model. 3) The spatial domain of the finite-domain one-dimensional convection-dispersion model is divided into multiple grids, and the dimensionless time domain is divided into multiple (two or more) time steps. The dimensionless governing equations and boundary conditions are discretized using the finite difference method; for multiple sets of... The number and dimensionless structure coefficients are solved step-by-step to obtain simulated concentration-time data at different monitoring locations, and a system is established. A library of standard curves relating number, dimensionless structure coefficient, and simulated concentration-time curves; 4) Obtain test data of a finite-length one-dimensional soil column under instantaneous source injection conditions. The test data shall include at least the length of the soil column test section, the radius of the soil column test section, the radius of the water column supplied above the inlet boundary, the initial height of the water column at the inlet, the permeability coefficient or hydraulic gradient and the measured flow rate, as well as the solute concentration-time data of multiple monitoring points set at different positions along the axial direction of the soil column. 5) Dimensionless processing of the monitoring data is performed. Dimensionless experimental concentration-time data from multiple monitoring locations are matched with simulated concentration-time data from the parameter curve library. A joint objective function is established based on the difference between the experimental and simulated concentrations. Corresponding parameters are selected from the breakthrough curve library and matched within the continuous parameter space to obtain the data corresponding to the experimental concentration-time data. Number and dimensionless structural coefficient ; 6) For any dimensionless structural coefficient, if the length of the soil column test section, the radius of the soil column test section, the radius of the supplied water column, and the initial water column height at the inlet are known, the porosity can be calculated. Then, the pore water velocity can be calculated using the porosity, the hydraulic conductivity or hydraulic gradient, and the measured flow rate. Calculate the longitudinal dispersion coefficient based on the number of pore water flow velocities; 7) Output The system outputs the experimental concentration-time curve, the optimal matching simulation curve, and the corresponding wiring error evaluation index, along with the dimensionless structural coefficient, porosity, average pore water flow velocity, and longitudinal dispersion coefficient.

[0007] The dimensionless relation in step 1) is: (2) In the formula: Let be the solute concentration at any point within the soil column. The initial solute concentration at the inlet boundary. Dimensionless concentration This refers to the axial position. The length of the soil column. For dimensionless distance, For time, The pore water flow velocity, For dimensionless time, The longitudinal dispersion coefficient is... For the Peckley number, This is the ratio of the pore water volume in the soil column to the initial solution volume in the supply container. Let the inner radius of the soil column be . Porosity The inner radius of the liquid supply container. This is the initial liquid level height in the supply container.

[0008] In step 2), the finite field model can achieve the coupled solution of the hydrodynamic field and the chemical field, where the length is... L Within a one-dimensional finite domain, the concentration satisfies the one-dimensional convection-dispersion equation: (3) After dimensionless transformation, we get: (4) Initial conditions of the chemical field: (5) After dimensionless transformation, we get: (6) Chemical field inlet boundary conditions: (7) In the formula: Solute concentration at the inlet boundary Because of the introduction If we represent the inlet boundary concentration, then we have a dimensionless inlet boundary concentration: (8) Then (6) after dimensionless transformation, we get: (9) In the formula: Dimensionless inlet boundary solute concentration In the finite field model, setting the outlet boundary as a constant concentration boundary, similar to the infinite field model, would result in a stable concentration gradient when the entire chemical field eventually stabilizes. This is unsuitable for models involving convection. Because under transient source excitation, the concentration gradient at the outlet remains zero until the excitation concentration reaches the outlet boundary. To accurately characterize the outlet boundary conditions of the finite field model, this paper chooses to ignore the dispersion effect at the outlet boundary and uses a zero-dispersion flux boundary as the outlet boundary, thus obtaining the chemical field outlet boundary conditions: (10) After dimensionless transformation, we obtain (11) The boundary value problem consisting of equations (2)-(10) contains two core dependent variables. and However, the boundary value problem only has one governing equation, equation (2), so it is impossible to solve it simultaneously. and To solve this boundary value problem, a governing equation must be added to obtain the analytical expression. Referring to the principles for constructing the water conservation equation at the well-aquifer interface in micro-water experiments, and based on the principles of water conservation and solute mass conservation at the water column-sand column interface at the inlet boundary, it is assumed that the mass of solute reduced in the water column above the one-dimensional sand column inlet boundary per unit time is equal to the mass of solute entering the sand column through convection-dispersion. Therefore, the solute mass conservation equation for the water column-sand column system and its initial conditions are established at the inlet boundary: (12) In the formula: Initial concentration (13) After dimensionless transformation, we get: (14) and its initial conditions: (15) Step 4) Using the finite difference method, the concentration at the inlet node in the program simultaneously represents the concentration in the supply container; therefore, it is initialized as follows: (16) In the formula: Initial concentration at the inlet node The remaining soil column nodes are initialized as follows: (17) In the formula: For the first Initial concentration at each node Divide the soil column axially into N parts The discrete equation for any internal grid node is: (18) in: Number the time steps. Number the spatial grid nodes. For dimensionless time steps, For dimensionless spatial step size, Indicates the first The time step, the first Dimensionless concentration of each spatial node Indicates the first The time step, the first Dimensionless concentration of each spatial node and They represent the first Nodes in each time step and Dimensionless concentration of adjacent nodes.

[0009] The dynamic inlet boundary conditions and the inlet concentration continuity conditions are implicitly discretized and written into the linear equation system at each time step. The program uses the backward Euler scheme to implicitly discretize the dynamic inlet equations and writes them into the first line of the tridiagonal equation system.

[0010] The discrete form corresponding to the entry node can be represented as: (19) in: and These represent the entry node at the [number]th [node]. and the Dimensionless concentration at each time step. This indicates the first internal node of the soil column adjacent to the entrance at the [missing information - likely a date or time value]. The dimensionless concentration at each time step, based on the inlet concentration continuity condition, means that the inlet node concentration also represents the dimensionless concentration in the supply container.

[0011] Discretize the zero-concentration gradient boundary at the outlet as follows: (20) in, For the concentration at the export node, The concentration is the concentration of the upstream adjacent node at the export point.

[0012] Combining the entrance boundary equations, internal node equations, and exit boundary equations forms a tridiagonal linear system of equations: (twenty one) in: This is the tridiagonal coefficient matrix corresponding to the current time step; Let the concentration vector at the new time step be the one to be determined. Composed of the concentration from the previous time step. Due to the use of variable time steps, the concentration at each step... The dimensionless time increments to the th Time step used at each dimensionless moment Because the values ​​are different, the program recalculates the tridiagonal matrix coefficients at each time step. The program uses the Thomas algorithm, also known as the tridiagonal matrix chasing method, to solve for the concentration vector at the current time step. The solution process consists of: forward elimination of the tridiagonal matrix; forward substitution of the right-hand side; back substitution from the exit node to the entry node; and obtaining the concentration vector at all spatial nodes. The dimensionless concentration at time step [time]. After solving, the current concentration is used as a known condition for the next time step.

[0013] In the parameter inversion process, the wiring method is used to extract the optimal model parameters from the experimental data. For each set of parameters to be evaluated (including the Peckley number), and dimensionless structure coefficient ), execute the following process: 1) Based on the current The value is used to calculate the corresponding porosity. ; 2) Utilize the obtained porosity The average pore water flow velocity was further calculated. ; 3) Record the actual time intervals during the experiment. Convert to dimensionless time under this set of parameters This is to align with the time reference of the simulation curve; 4) Interpolate the simulated curves so that they correspond to the position of each sampling point of the experimental data on the dimensionless time axis; 5) Finally, at the same time point, the concentration values ​​measured in the experiment are compared point by point with the concentration values ​​in the experimental concentration-time data and parameter curve library.

[0014] Construct a joint objective function that can comprehensively evaluate the goodness of fit across all monitoring locations and all time points: (twenty two) (twenty three) in, Represents the total number of monitoring points. For the first The number of valid sampling points at each monitoring location For the first Dimensionless coordinates of each monitoring location, This is the actual elapsed time from the moment of injection to the sampling point. For candidate structure coefficients The average pore water flow velocity calculated based on the corresponding porosity. Number the monitoring location. For the first The sampling point number at each monitoring location, and These are the weighting coefficients assigned to each data point, used to adjust the importance of different data points in the fit. To test dimensionless concentration, To simulate dimensionless concentration, This represents the maximum absolute value of the dimensionless concentration tested. Based on this, the final parameters are determined, and the error between the experimental curve and all candidate curves in the parameter curve library is calculated using the above objective function. From these, several (two or more) groups with the smallest error are selected. and These are initially selected as candidate results. Then, the group that performs best among these candidate results—that is, the group with the smallest NRMSE value under the same valid experimental data, weight settings, and parameter constraints—is used as the initial guess. A local fine-grained search is then performed in the continuous parameter space, and through further wiring optimization, the optimal parameters are finally obtained. and .

[0015] The optimal dimensionless parameters are obtained through the wiring method. and Then, based on the set parameters of the experimental setup, it is possible to further reverse-engineer the results. , , Three hydrogeological parameters. The specific conversion process is as follows: 1) Back-calculation of porosity Using the finally identified structural coefficients Based on the geometric dimensions of the test column, the porosity is calculated using the following formula: (twenty four) in, and These represent the characteristic radius and characteristic height (related to the solute or the solid phase), respectively. Let be the radius of the cylinder. Let be the length of the column.

[0016] 2) Back-calculation of average pore water velocity Both methods are acceptable, and the appropriate method can be selected based on the actual measurement conditions. If the hydraulic conductivity coefficient is known and hydraulic gradient Then, using Darcy's law, the value is converted to the pore size: (25) in, This is the hydraulic conductivity coefficient, expressed in m / s. The dimensionless hydraulic gradient; If the measured flow rate is known Then simply divide the flow rate by the pore cross-sectional area: (26) 3) Inverse calculation of longitudinal dispersion coefficient After obtaining and the optimal Peckley number Then, the longitudinal dispersion coefficient is obtained by inversely solving the definition of the Peckley number: (27) Any techniques not mentioned in this invention are based on existing technologies.

[0017] Beneficial effects: Compared with the prior art, the present invention has the following advantages: 1) This invention directly uses the finite difference method in the time domain to uniformly discretize the one-dimensional convection-diffusion control equation, dynamic inlet boundary conditions, and outlet zero concentration gradient boundary conditions, without the need for Laplace domain frequency domain solution derivation and numerical inverse transformation. This shortens the calculation chain from the control equation to the time domain concentration response curve, and facilitates unified control of the discretization format, grid size, and time step.

[0018] 2) This invention adopts a variable time step approach, using denser time nodes during the rapid concentration change phase and sparser time nodes during the slow concentration decay phase, and reconstructs the tridiagonal coefficient matrix according to the current time step, thereby taking into account the computational resolution and computational load of different stages of the instantaneous source response process.

[0019] 3) This invention uses the backward Euler scheme to implicitly discretize the time term, the central difference for the dispersion term, and the first-order upwind difference for the convection term. It also uses the Thomas algorithm to solve the tridiagonal linear equation system, which helps to suppress non-physical numerical oscillations that may occur under large Péclet numbers, while reducing the matrix solution overhead for a single time step.

[0020] 4) This invention constructs a weighted NRMSE joint objective function that integrates multiple monitoring locations and multiple sampling times. Data weights can be set according to measurement variance, number of effective samples, or concentration magnitude to avoid excessive influence of a single measuring point, a large number of sampling points, or a measuring point with a large concentration amplitude on the inversion results, thereby improving the consistency of multi-point experimental data utilization. Attached Figure Description

[0021] Figure 1 A schematic diagram of a one-dimensional solute transport model under transient source excitation; Figure 2 This is a flowchart illustrating the multi-parameter synchronous inversion method for solute transport in a one-dimensional soil column based on the finite difference method, as an example of the present invention. Figure 3 For different Pe Comparison of dimensionless concentration-time response curves at different dimensionless spatial locations under the influence of the number and dimensionless structure coefficient; Figure 4 Examples of the present invention - Matching diagram of fitted curve and standard curve library; Figure 5 This is a comparison chart of the NRMSE of five advantage curves in an example of the present invention; Figure 6 This is a comparison diagram of the Stehfest semi-analytical solution and the finite difference solution in the example of this invention; Detailed Implementation

[0022] The following examples are for illustrative purposes only and are not intended to limit the scope of protection of the invention. Unless otherwise specified, all physical quantities are expressed in a consistent unit system.

[0023] A multi-parameter synchronous inversion method for solute transport in a one-dimensional soil column based on the finite difference method is proposed. The method acquires flow rate, head, and concentration-time data at multiple locations for a finite-length one-dimensional soil column under instantaneous source injection conditions, and performs dimensionless processing of concentration, spatial, and temporal parameters. A mathematical model is established, including finite-domain convection-dispersion control equations, dynamic inlet boundary, and outlet zero-concentration gradient boundary. The finite difference method is used to solve for the simulated concentration at each measuring point step-by-step. A joint objective function is constructed based on the simulated and measured concentrations from multiple sets of experiments, multiple measuring points, and multiple time points. The parameters to be inverted are iteratively updated under physical constraints until the convergence condition is met, and the porosity, pore water velocity, and longitudinal dispersion are output. The invention is further described below with reference to the accompanying drawings and embodiments.

[0024] Example 1: This example provides a multi-parameter synchronous inversion method for solute transport in a one-dimensional soil column based on the finite difference method. A sodium chloride transient source tracer experiment was conducted using a vertical one-dimensional sand column to obtain concentration-time data at different axial positions, as shown below. Figure 1 A schematic diagram of a one-dimensional solute transport model under transient source excitation is shown. Based on the convection-dispersion governing equations, using boundary control conditions as the framework, and the finite difference method as the solution method, a series of solutions under different... Pe Number and sum β A standard curve library is provided under certain conditions. After obtaining the standard library, the corresponding parameters of any experimental dataset can be obtained through the matching method. By substituting the experimental setup data, a series of hydrogeological parameters can be further obtained, such as... Figure 2 The flowchart of the multi-parameter synchronous inversion method for solute transport in a one-dimensional soil column based on the finite difference method is shown, including the following steps: 1) Construct a convection-dispersion model of the soil column under finite domain conditions. The model includes dimensionless convection-dispersion governing equations, initial concentration conditions inside the soil column, dynamic inlet boundary conditions describing the mass change of solute in the supply water column, and zero dispersion flux boundary conditions at the outlet of the soil column. 2) Select appropriate dimensionless factors to make key variables such as solute concentration and time dimensionless, and obtain a dimensionless convection-dispersion mathematical model; 3) The spatial domain of the finite-domain one-dimensional convection-dispersion model is divided into multiple grids, and the dimensionless time domain is divided into multiple time steps. The finite difference method is used to discretize the governing equations and boundary conditions; for multiple sets of... The number and dimensionless structure coefficients are solved step-by-step to obtain simulated concentration-time data at different monitoring locations, and a system is established. A library of standard curves relating number, dimensionless structure coefficient, and simulated concentration-time curves; 4) Obtain test data of a finite-length one-dimensional soil column under instantaneous source injection conditions. The test data shall include at least the length of the soil column test section, the radius of the soil column test section, the radius of the water column supplied above the inlet boundary, the initial height of the water column at the inlet, the permeability coefficient or hydraulic gradient and the measured flow rate, as well as the solute concentration-time data of multiple monitoring points set at different positions along the axial direction of the soil column. 5) Dimensionless processing of the monitoring data is performed. Dimensionless experimental concentration-time data from multiple monitoring locations are matched with simulated concentration-time data from the parameter curve library. A joint objective function is established based on the difference between the experimental and simulated concentrations. Corresponding parameters are selected from the breakthrough curve library and matched within the continuous parameter space to obtain the data corresponding to the experimental concentration-time data. Number and dimensionless structural coefficient ; 6) For any dimensionless structural coefficient, if the length of the soil column test section, the radius of the soil column test section, the radius of the supplied water column, and the initial water column height at the inlet are known, the porosity can be calculated. Then, the pore water velocity can be calculated using the porosity, the hydraulic conductivity or hydraulic gradient, and the measured flow rate. Calculate the longitudinal dispersion coefficient based on the number of pore water flow velocities; 7) Output The system outputs the experimental concentration-time curve, the optimal matching simulation curve, and the corresponding wiring error evaluation index, along with the dimensionless structural coefficient, porosity, average pore water flow velocity, and longitudinal dispersion coefficient.

[0025] The dimensionless relation in step 1) is: (28) In the formula: Let be the solute concentration at any point within the soil column. The initial solute concentration at the inlet boundary. Dimensionless concentration This refers to the axial position. The length of the soil column. For dimensionless distance, For time, The pore water flow velocity, For dimensionless time, The longitudinal dispersion coefficient is... For the Peckley number, This is the ratio of the pore water volume in the soil column to the initial solution volume in the supply container. Let the inner radius of the soil column be . Porosity The inner radius of the liquid supply container. This is the initial liquid level height in the supply container.

[0026] In step 2), the finite field model can achieve the coupled solution of the hydrodynamic field and the chemical field, where the length is... L Within a one-dimensional finite domain, the concentration satisfies the one-dimensional convection-dispersion equation: (29) After dimensionless transformation, we get: (30) Initial conditions of the chemical field: (31) After dimensionless transformation, we get: (32) Chemical field inlet boundary conditions: (33) In the formula: Solute concentration at the inlet boundary Because of the introduction If we represent the inlet boundary concentration, then we have a dimensionless inlet boundary concentration: (34) Then (6) after dimensionless transformation, we get: (35) In the formula: Dimensionless inlet boundary solute concentration In the finite field model, setting the outlet boundary as a constant concentration boundary, similar to the infinite field model, would result in a stable concentration gradient when the entire chemical field eventually stabilizes. This is unsuitable for models involving convection. Because under transient source excitation, the concentration gradient at the outlet remains zero until the excitation concentration reaches the outlet boundary. To accurately characterize the outlet boundary conditions of the finite field model, this paper chooses to ignore the dispersion effect at the outlet boundary and uses a zero-dispersion flux boundary as the outlet boundary, thus obtaining the chemical field outlet boundary conditions: (36) After dimensionless transformation, we obtain (37) The boundary value problem consisting of equations (2)-(10) contains two core dependent variables. and However, the boundary value problem only has one governing equation, equation (2), so it is impossible to solve it simultaneously. and To solve this boundary value problem, a governing equation must be added to obtain the analytical expression. Referring to the principles for constructing the water conservation equation at the well-aquifer interface in micro-water experiments, and based on the principles of water conservation and solute mass conservation at the water column-sand column interface at the inlet boundary, it is assumed that the mass of solute reduced in the water column above the one-dimensional sand column inlet boundary per unit time is equal to the mass of solute entering the sand column through convection-dispersion. Therefore, the solute mass conservation equation for the water column-sand column system and its initial conditions are established at the inlet boundary: (38) In the formula: Initial concentration (39) After dimensionless transformation, we get: (40) and its initial conditions: (41) Step 4) Using the finite difference method, the concentration at the inlet node in the program simultaneously represents the concentration in the supply container; therefore, it is initialized as follows: (42) In the formula: Initial concentration at the inlet node The remaining soil column nodes are initialized as follows: (43) In the formula: For the first Initial concentration at each node Divide the soil column axially into N parts The discrete equation for any internal grid node is: (44) in: Number the time steps. Number the spatial grid nodes. For dimensionless time steps, For dimensionless spatial step size, Indicates the first The time step, the first Dimensionless concentration of each spatial node Indicates the first The time step, the first Dimensionless concentration of each spatial node and They represent the first Nodes in each time step and Dimensionless concentration of adjacent nodes.

[0027] The dynamic inlet boundary conditions and the inlet concentration continuity conditions are implicitly discretized and written into the linear equation system at each time step. The program uses the backward Euler scheme to implicitly discretize the dynamic inlet equations and writes them into the first line of the tridiagonal equation system.

[0028] The discrete form corresponding to the entry node can be represented as: (45) in: and These represent the entry node at the [number]th [node]. and the Dimensionless concentration at each time step. This indicates the first internal node of the soil column adjacent to the entrance at the [missing information - likely a date or time value]. The dimensionless concentration at each time step, based on the inlet concentration continuity condition, means that the inlet node concentration also represents the dimensionless concentration in the supply container.

[0029] Discretize the zero-concentration gradient boundary at the outlet as follows: (46) in, For the concentration at the export node, The concentration is at the upstream adjacent node of the export.

[0030] Combining the entrance boundary equations, internal node equations, and exit boundary equations forms a tridiagonal linear system of equations: (47) in: This is the tridiagonal coefficient matrix corresponding to the current time step; Let the concentration vector at the new time step be the one to be determined. Composed of the concentration from the previous time step. Due to the use of variable time steps, the concentration at each step... The dimensionless time increments to the th Time step used at each dimensionless moment Because the values ​​are different, the program recalculates the tridiagonal matrix coefficients at each time step. The program uses the Thomas algorithm, also known as the tridiagonal matrix chasing method, to solve for the concentration vector at the current time step. The solution process consists of: forward elimination of the tridiagonal matrix; forward substitution of the right-hand side; back substitution from the exit node to the entry node; and obtaining the concentration vector at all spatial nodes. The dimensionless concentration at time step [time]. After solving, the current concentration is used as a known condition for the next time step.

[0031] By performing multiple rounds of iterative calculations on the tridiagonal matrix, the results are obtained from different groups. Pe Number and sum β Under these conditions, there will be corresponding penetration curves. Nine representative sets of data are selected and listed, such as... Figure 3 different Pe The figure shows a comparison of dimensionless concentration-time response curves at different dimensionless spatial locations under the number and dimensionless structure coefficient.

[0032] After designing and fabricating experimental equipment, filling soil samples, and conducting basic hydrogeological tests, the basic hydrogeological parameters required for inversion can be obtained. An inner radius of [missing information] is used. Effective fill length is L A one-dimensional soil column, the main body of which is filled with homogeneous quartz sand, is connected above the entrance of the soil column with an inner radius of... The initial liquid level height is The supply containers are shown in Table 1. Four concentration monitoring points were set along the axial direction of the soil column, and the sampling time and concentration at each monitoring point were recorded. Deionized water was continuously added from the supply section to form a controllable seepage field, and the liquid level in the supply container, the water head at both ends of the soil column, and the outflow rate were recorded simultaneously. Repeated timestamps were merged, missing values ​​were marked without forced interpolation, and data that deviated significantly from the repeated measurement range were reviewed or removed. In the indoor solute transport experiment, a series of scattered data points could be obtained by observing the probe software interface on the computer connected to the probe. These data points were fitted with a standard curve library, and the optimal fitting curve could be found by comparing them one by one. Figure 4 As shown

[0033] Table 1. Device parameters and their hydrogeological parameters

[0034] In the parameter inversion process, the wiring method is used to extract the optimal model parameters from the experimental data. For each set of parameters to be evaluated (including the Peckley number), and dimensionless structure coefficient ), execute the following process: 1) Based on the current The value is used to calculate the corresponding porosity. ; 2) Utilize the obtained porosity The average pore water flow velocity was further calculated. ; 3) Record the actual time intervals during the experiment. Convert to dimensionless time under this set of parameters This is to align with the time reference of the simulation curve; 4) Interpolate the simulated curves so that they correspond to the position of each sampling point of the experimental data on the dimensionless time axis; 5) Finally, at the same time point, the concentration values ​​measured in the experiment are compared point by point with the concentration values ​​in the experimental concentration-time data and parameter curve library.

[0035] Construct a joint objective function that can comprehensively evaluate the goodness of fit across all monitoring locations and all time points: (48) (49) in, Represents the total number of monitoring points. For the first The number of valid sampling points at each monitoring location For the first Dimensionless coordinates of each monitoring location, This is the actual elapsed time from the moment of injection to the sampling point. For candidate structure coefficients The average pore water flow velocity calculated based on the corresponding porosity. Number the monitoring location. For the first The sampling point number at each monitoring location, and These are the weighting coefficients assigned to each data point, used to adjust the importance of different data points in the fit. To test dimensionless concentration, To simulate dimensionless concentration, This represents the maximum absolute value of the dimensionless concentration tested. Weight The objective function can be set based on the reciprocal of the measurement variance, or it can be normalized first to the residuals of each measuring point according to its effective sample size or concentration order of magnitude, to avoid a single measuring point with a large number of sampling points or high concentration dominating the objective function. Head and flow rate data can be used as additional residual terms. Based on this, the final parameters are determined, and the error between the experimental curve and all candidate curves in the parameter curve library is calculated using the above objective function. From these, several sets with the smallest errors are selected. and like Figure 5 As shown, these serve as preliminary candidate results. Then, the best-performing set of these candidate results is used as the initial guess, and a local fine-grained search is performed in the continuous parameter space. Through further wiring optimization, the optimal parameters are finally obtained. and .

[0036] The optimal dimensionless parameters are obtained through the wiring method. and Then, based on the set parameters of the experimental setup, it is possible to further reverse-engineer the results. , , Three hydrogeological parameters. The specific conversion process is as follows: 1) Back-calculation of porosity Using the finally identified structural coefficients Based on the geometric dimensions of the test column, the porosity is calculated using the following formula: (50) in, and These represent the characteristic radius and characteristic height (related to the solute or the solid phase), respectively. Let be the radius of the cylinder. Let be the length of the column.

[0037] Given calculations: Therefore The formula simplifies to:

[0038] 2) Back-calculation of average pore water velocity Both methods are acceptable, and the appropriate method can be selected based on the actual measurement conditions. If the hydraulic conductivity coefficient is known and hydraulic gradient Then, using Darcy's law, the value is converted to the pore size: (51) Calculate: Pore water flow velocity

[0039] 3) Inverse calculation of longitudinal dispersion coefficient After obtaining and the optimal Peckley number Then, the longitudinal dispersion coefficient is obtained by inversely solving the definition of the Peckley number: (51) Given calculations:

[0040] Output , , , and At the same time, the measured curves, root mean square error, and coefficient of determination of each measuring point are obtained.

[0041] Comparative Example 1

[0042] Under the same solute transport control equation, initial conditions, and boundary conditions, a new set of databases is obtained based on the semi-analytical solution of the Laplace transform and Stehfest numerical inverse transform. After multiple wiring, the optimal fitting curve that best matches the experimental data in Example 1 can be obtained. The values ​​of each item of this curve are compared with the curve obtained by the finite difference decomposition to obtain Table 2.

[0043] Table 2 Comparison of Stehfest semi-analytical solution and finite difference decomposition

[0044] As shown in the table above, the fitted Pe and β numbers are very close, and the joint objective function test reveals that the difference between them is small. Taking the fitting time of the two methods, it is clear that the solution process using the finite difference method is significantly faster than the other group. Correspondingly, in the entire multiple optimization fitting process, the total fitting time of the finite difference method is also much shorter than that of the Stehfest numerical inverse transform method.

Claims

1. A multi-parameter synchronous inversion method for solute transport in a one-dimensional soil column based on the finite difference method, characterized in that, Includes the following steps: 1) Construct a one-dimensional convection-dispersion model of the soil column under finite domain conditions. The model includes the convection-dispersion control equation inside the soil column, the initial concentration condition inside the soil column, the dynamic inlet boundary condition describing the change of solute mass in the supply water column, and the zero dispersion flux boundary condition at the outlet of the soil column. 2) Select the Peckley number Pe With dimensionless structural coefficient Using the core dimensionless parameter, the concentration, spatial, and temporal variables are dimensionless to obtain a dimensionless convection-diffusion mathematical model. 3) The spatial domain of the finite-domain one-dimensional convection-dispersion model is divided into multiple grids, and the dimensionless time domain is divided into multiple time steps. The dimensionless governing equations and boundary conditions are discretized using the finite difference method; for multiple sets of... The number and dimensionless structure coefficients are solved step-by-step to obtain simulated concentration-time data at different monitoring locations, and a system is established. A library of standard curves relating number, dimensionless structure coefficient, and simulated concentration-time curves; 4) Obtain test data of a finite-length one-dimensional soil column under instantaneous source injection conditions. The test data shall include at least the length of the soil column test section, the radius of the soil column test section, the radius of the water column supplied above the inlet boundary, the initial height of the water column at the inlet, the permeability coefficient or hydraulic gradient and the measured flow rate, as well as the solute concentration-time data of multiple monitoring points set at different positions along the axial direction of the soil column. 5) Dimensionless processing of the monitoring data is performed. Dimensionless experimental concentration-time data from multiple monitoring locations are matched with simulated concentration-time data from the parameter curve library. A joint objective function is established based on the difference between the experimental and simulated concentrations. Corresponding parameters are selected from the breakthrough curve library and matched within the continuous parameter space to obtain the data corresponding to the experimental concentration-time data. Number and dimensionless structural coefficient ; 6) For any dimensionless structural coefficient, if the length of the soil column test section, the radius of the soil column test section, the radius of the supplied water column, and the initial water column height at the inlet are known, the porosity can be calculated. Then, the pore water velocity can be calculated using the porosity, the hydraulic conductivity or hydraulic gradient, and the measured flow rate. Calculate the longitudinal dispersion coefficient based on the number of pore water flow velocities; 7) Output The system outputs the experimental concentration-time curve, the optimal matching simulation curve, and the corresponding wiring error evaluation index, along with the dimensionless structural coefficient, porosity, average pore water flow velocity, and longitudinal dispersion coefficient.

2. The method for simultaneous inversion of solute transport in a one-dimensional soil column based on finite difference forward modeling according to claim 1, characterized in that, In step 3), the finite difference method is used. The concentration at the inlet node in the program simultaneously represents the concentration in the supply container; therefore, it is initialized as follows: (2) In the formula: The initial concentration at the inlet node; The remaining soil column nodes are initialized as follows: (3) In the formula: For the first Initial concentration at each node Divide the soil column axially into N parts; The discrete equation for any internal grid node is: (4) in: Number the time steps. Number the spatial grid nodes. For dimensionless time steps, The dimensionless spatial step size, Indicates the first The time step, the first Dimensionless concentration of each spatial node Indicates the first The time step, the first Dimensionless concentration of each spatial node and They represent the first Nodes in each time step and Dimensionless concentration of adjacent nodes; implicitly discretize the dynamic inlet boundary conditions and the inlet concentration continuity conditions and write them into the linear equation system at each time step. The program uses the backward Euler scheme to implicitly discretize the dynamic inlet equations and writes them into the first line of the tridiagonal equation system. The discrete form corresponding to the entry node can be represented as: (5) in: and These represent the entry node at the [number]th [node]. and the Dimensionless concentration at each time step; This indicates the first internal node of the soil column adjacent to the entrance at the [missing information - likely a date or time value]. The dimensionless concentration at each time step, based on the continuous inlet concentration condition, means that the inlet node concentration also represents the dimensionless concentration in the supply container. Discretize the zero-concentration gradient boundary at the outlet as follows: (6) in, For the concentration at the export node, The concentration of the upstream adjacent node at the export point; Combining the entrance boundary equations, internal node equations, and exit boundary equations forms a tridiagonal linear system of equations: (7) in: This is the tridiagonal coefficient matrix corresponding to the current time step; Let the concentration vector at the new time step be the one to be determined. Composed of the concentration from the previous time step; due to the use of variable time steps, the concentration of each step is... The dimensionless time increments to the th Time step used at each dimensionless moment The process differs, therefore the program recalculates the tridiagonal matrix coefficients at each time step. The program employs the Thomas algorithm, i.e., the tridiagonal matrix chasing method, to solve for the concentration vector at the current time step. The solution process involves: forward elimination of the tridiagonal matrix; forward substitution of the right-hand side; and back substitution from the exit node to the entry node; resulting in the concentration vector at all spatial nodes. The dimensionless concentration at time step 1; after solving, the current concentration is used as a known condition for the next time step.

3. The method according to claim 1 or 2, characterized in that, In step 5, the wiring method is used to extract the optimal model parameters from the experimental data. For each set of parameters to be evaluated, the following procedure is executed: 1) Based on the current The value is used to calculate the corresponding porosity. ; 2) Utilize the obtained porosity The average pore water flow velocity was further calculated. ; 3) Record the actual time intervals during the experiment. Convert to dimensionless time under this set of parameters This is to align with the time reference of the simulation curve; 4) Interpolate the simulated curves so that they correspond to the position of each sampling point of the experimental data on the dimensionless time axis; 5) At the same time point, compare the dimensionless concentration measured in the experiment with the dimensionless concentration in the experimental concentration-time data and the parameter curve library point by point.

4. The method according to claim 3, characterized in that, In step 5, a joint objective function is constructed that can comprehensively evaluate the goodness of fit for all monitoring locations and all time points: (8) (9) in, Represents the total number of monitoring points. For the first The number of valid sampling points at each monitoring location For the first Dimensionless coordinates of each monitoring location, This is the actual elapsed time from the moment of injection to the sampling point. For candidate structure coefficients The average pore water flow velocity calculated based on the corresponding porosity. Number the monitoring location. For the first The sampling point number at each monitoring location, and These are the weighting coefficients assigned to each data point, used to adjust the importance of different data points in the fit. To test dimensionless concentration, To simulate dimensionless concentration, This represents the maximum absolute value of the dimensionless concentration tested.

5. The method according to claim 4, characterized in that, In step 5, the specific method for determining the parameters is as follows: calculate the error between the experimental curve and all candidate parameter combinations in the standard curve library using the joint objective function, and then select the Peckley numbers with the smallest errors. and dimensionless structure coefficient The first set of candidates is used as the initial candidate; then, the best-performing set from the initial candidate is used as the initial guess, and a local fine-grained search is performed in the continuous parameter space. Through wiring optimization, the optimal Pecley number is finally obtained. and dimensionless structure coefficient Among them, the best performance refers to the minimum value of the joint objective function NRMSE under the same valid experimental data, weight settings and parameter constraints.

6. The method according to claim 5, characterized in that, In step 6), the porosity is further calculated back based on the set parameters of the experimental setup. Average pore water velocity Longitudinal dispersion coefficient The specific conversion process for the three hydrogeological parameters is as follows: 1) Back-calculation of porosity Using the finally identified structural coefficients Based on the geometric dimensions of the test column, the porosity is calculated using the following formula: (10) in, and These represent the characteristic radius and characteristic height (related to the solute or the solid phase), respectively. Let be the radius of the cylinder. The length of the column; 2) Back-calculation of average pore water velocity Two approaches are used: If the hydraulic conductivity coefficient is known and hydraulic gradient Then, using Darcy's law, the value is converted to the pore size: (11) in, This is the hydraulic conductivity coefficient, expressed in m / s. The dimensionless hydraulic gradient; If the measured flow rate is known Then simply divide the flow rate by the pore cross-sectional area: (12) 3) Inverse calculation of longitudinal dispersion coefficient After obtaining and the optimal Peckley number Then, the longitudinal dispersion coefficient is obtained by inversely solving the definition of the Peckley number: (13)。

Images