Inversion calculation method for elevation and parameters of air inlet surface under driving of periodic water head

Abstract

The invention discloses an air inlet surface elevation and parameter inversion calculation method under periodic water head driving. The method comprises the following steps: S1, extracting a specific yield mu ha at a parameter air inlet value surface, a saturated permeability coefficient Kx, an infiltration supply rate W at the air inlet value surface, an air inlet value water head ha, an initial air inlet value surface elevation Zha0 and a starting moment t0; s2, the parameters are substituted into a formula, and the elevation Zha of the air inlet value face is worked out. The method relates to the technical field of simulation and analysis of groundwater flow dynamics under periodic water head driving. In order to overcome the limitation of a classical model of which a free surface is positioned on a diving surface in describing the groundwater dynamic state under the driving of the periodic water head, the invention constructs a groundwater dynamic model of which the free surface is positioned on an air inlet value surface under the driving of the periodic water head; and the method is used for describing the groundwater dynamic state under the driving of the periodic water head and carrying out aquifer parameter inversion.

Description

Technical Field

[0001] The present invention relates to the technical field of simulating and analyzing the dynamic state of groundwater flow driven by a periodic water head, and specifically relates to a method for inverse calculation of the elevation and parameters of an air intake surface under the drive of a periodic water head. Background Technique

[0002] One-dimensional saturated flow models and wet-dry cycle models for locating the free surface of the phreatic water surface are often used as important tools for simulating and analyzing the dynamic state of groundwater flow driven by a periodic water head. However, these classical models using the free surface to locate the phreatic water surface will have three limitations: 1) inconsistent assumptions in the models; 2) the seepage rate and the velocity of the free surface movement in the model calculation results are always zero at the same time; 3) the hydrodynamic mechanism differences between the free surface movement and infiltration recharge are ignored.

[0003] The inconsistent assumptions are manifested in that these classical models all use both the sudden interface assumption and the non-sudden interface assumption at the same time. The sudden interface assumption does not consider the existence of moisture in the sediment above the free surface. For example, in classical models, the effective porosity n (equal to the difference between the saturated water content and the residual water content) is used as a parameter of the model to characterize the drainage capacity of the free surface, and this parameter is derived from the sudden interface assumption. The non-sudden interface assumption considers the existence of moisture above the free surface. In classical models, the equivalent saturated zone height is adopted, which is defined as the sum of the saturated zone height and the equivalent saturated height of the unsaturated zone. This means that the concept of the equivalent saturated zone height uses the non-sudden interface assumption. Using mutually contradictory sudden and non-sudden interfaces in the same equation is not advisable in a strict physical sense.

[0004] According to these classical models, in the calculation of the groundwater dynamics in a 1D soil column under the drive of a periodic water head, it is always possible to deduce that the seepage rate and the free surface movement velocity are zero at the same time. In fact, the experimental observation results show that the seepage rate and the free surface movement velocity are not zero at the same time. They are only zero at the same time in some special cases, such as when the seepage reaches a steady state.

[0005] Classical models ignore the hydrodynamic mechanism differences between the free surface movement and infiltration recharge. They confuse the process of recharging or draining the unconfined aquifer when the free surface rises or falls with the process of water flowing from the unsaturated zone above the free surface downward through the free surface into the aquifer or the process of water flowing from the aquifer upward through the free surface into the unsaturated zone above the free surface.

[0006] To overcome the limitations of classical models in describing the groundwater dynamics under the drive of a periodic water head, there is an urgent need for a method for inverse calculation of the elevation and parameters of an air intake surface under the drive of a periodic water head. Summary of the Invention

[0007] To achieve the above object, the present invention is realized by the following technical solutions: A method for inverse calculation of the elevation and parameters of the intake surface under periodic head driving, comprising the following steps:

[0008] S1. Estimate the initial values of the parameters and set the upper and lower limits. The parameters include: the specific yield μ at the intake value surface ha , the saturated hydraulic conductivity K z , the infiltration recharge rate W at the intake value surface, the intake value head h a , the initial elevation z of the intake value surface ha0 and the starting time t0.

[0009] S2. Substitute the above initial parameter values into the calculation formula for the elevation z of the groundwater intake value surface under periodic head driving with the free surface located at the intake value surface to calculate the elevation z of the intake value surface ha , and the formula is:

[0010]

[0011] In the formula, A is the amplitude of the periodic head, ω is the frequency of the periodic head; t is the time;

[0012] S3. Substitute the above calculation results into the total head formula with the free surface located at the intake value surface to calculate the total head H, and the formula is:

[0013]

[0014] S4. Input the observed value of the total head H at the observation point.

[0015] S5. Calculate the objective function value f through the objective function formula, and the formula is:

[0016]

[0017] In the formula, H cij, H tij are respectively the calculated value and the observed value of the total head at the i th -th observation point at the j th -th moment.

[0018] S6. Judge whether the termination condition is satisfied according to the calculated objective function value f;

[0019] S7. If the termination condition is satisfied, output the optimized parameters: μ ha , K z , W, h a , z ha0 , t0, and obtain the inversion result.

[0020] S8. Substitute the parameters output in step S7 into the calculation formula for the elevation of the groundwater intake value surface under the periodic water head with the free surface located at the intake value surface, and calculate the elevation Zha of the intake value surface.

[0021] S9. Substitute the elevation of the intake value surface calculated in step S8 into the total water head formula with the free surface located at the intake value surface, predict the total water head profile under the periodic water head, and obtain the forward prediction result.

[0022] Preferably, the calculation formula for the elevation of the groundwater intake value surface under the periodic water head with the free surface located at the intake value surface in step S2 can also calculate the elevation z of the intake value surface through the following approximate formula ha approximate solution, the formula is:

[0023]

[0024] Preferably, in step S6, if the termination condition is not met, it is necessary to regenerate the parameter values μ ha 、K z 、W, h a 、z ha0 、t0, and compare the regenerated parameter values with the parameter upper and lower limits to determine whether they meet the parameter upper and lower limits. If not, regenerate the parameter values; if so, substitute the regenerated parameter values into step S2 and repeat the operations from step S2 to step S6.

[0025] Preferably, the formula for the parameters to meet the parameter upper and lower limits is:

[0026] μ hamin <μ ha <μ hamax

[0027] K zmin <K z <K zmax

[0028] W min <W<W max

[0029] h amin <h a <h amax

[0030] z ha0min <z ha0 <z ha0max

[0031] t 0min <t0<t 0max

[0032] K z >-W

[0033] In the formula, the parameter μ ha , K z , W, h a , z ha0 , the subscripts min and max of t0 are the parameters μ ha , K z , W, h a , z ha0 , the upper and lower limit values of t0; K z > -W restricts the saturated hydraulic conductivity K z and the relative magnitude of the infiltration recharge rate W at the intake value surface.

[0034] The present invention provides a method for inverse calculation of the elevation of the intake surface and parameters under periodic head driving. To overcome the limitations of the classical model in describing the groundwater dynamics under periodic head driving, the present invention gives the calculation formula for the elevation of the groundwater intake value surface with the free surface located at the intake value surface and the total head formula with the free surface located at the intake value surface, and uses them to describe the groundwater dynamics under periodic head driving and conduct inverse aquifer parameter calculation. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 is the flow block diagram of the present invention;

[0036] Figure 2 is the comparison diagram of the calculated total head value and the observed value at the observation point of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0037] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0038] First embodiment, please refer to Figure 1 , the present invention provides a technical solution: a method for inverse calculation of the elevation of the intake surface and parameters under periodic head driving, including the following steps:

[0039] S1. Estimate the initial values of the parameters and set the upper and lower limits. The parameters include: the specific yield μ at the intake value surface ha , the saturated hydraulic conductivity K z , the infiltration recharge rate W at the intake value surface, the intake value head h a , the initial elevation z of the intake value surface ha0 and the starting time t0.

[0040] S2. Substitute the above initial parameter values into the calculation formula for the elevation z of the groundwater intake value surface under the periodic head drive with the free surface located at the intake value surface to calculate the elevation of the intake value surface z ha , and the formula is:

[0041]

[0042] In the formula, A is the amplitude of the periodic head, ω is the frequency of the periodic head; t is the time;

[0043] S3. Substitute the above calculation results into the total head formula with the free surface located at the intake value surface to calculate the total head H, and the formula is:

[0044]

[0045] S4. Input the observed value H of the total head at the observation point.

[0046] S5. Calculate the objective function value f through the objective function formula, and the formula is:

[0047]

[0048] In the formula, H cij , H tij are respectively the calculated value and the observed value of the total head at the i th observation point at the j th moment.

[0049] S6. Judge whether the termination condition is satisfied according to the calculated objective function value f;

[0050] S7. If f does not satisfy the termination condition, re-extract the parameter values μha, Kx, W, ha, Zha0, t0, and substitute the extracted parameter values into the upper and lower limit formulas

[0051] μ hamin <μ ha <μ hamax

[0052] K zmin <K z <K zmax

[0053] W min <W<W max

[0054] h amin <h a <h amax

[0055] z ha0min <z ha0 <z ha0max

[0056] t 0min <t0<t 0max

[0057] K z > -W;

[0058] When determining whether the upper and lower limits are satisfied, if not, re-extract the function value. If satisfied, substitute the extracted parameter value into Step 2 and repeat the operations from Step 2 to Step 6;

[0059] S8. If the termination condition is satisfied, output the optimized parameters: μ ha , K z , W, h a , z ha0 , t0, to obtain the inversion result.

[0060] S9. Substitute the parameters output in Step S6 into the calculation formula for the elevation of the groundwater intake value surface under the periodic water head drive with the free surface located at the intake value surface, and calculate the elevation Zha of the intake value surface;

[0061] S10. Substitute the elevation of the intake value surface calculated in Step S8 into the total water head formula with the free surface located at the intake value surface, predict the total water head profile under the periodic water head drive, and obtain the forward prediction result.

[0062] For the second embodiment, please refer to Figure 1 - Figure 2 , the present invention provides a technical solution: A method for inverse calculation of the elevation and parameters of the intake surface under periodic water head drive, including the following steps: S1. After measurement, the saturated permeability coefficient of a certain sand layer is 5×10-5 m / s, and the intake value ha is -0.08 m. The initial total water head in the sand layer is 1 m, and a periodic water head of 0.2×sin(1 / 5000×2π×t)+1 m acts at the bottom of the sand layer. Observation points are set at 0.5 m and 0.7 m above the bottom of the sand layer, and the total water heads at the observation points are obtained as shown in Table 1.

[0063] Table 1 Total water heads at observation points 0.5 m and 0.7 m from the bottom of the sand layer

[0064]

[0065]

[0066]

[0067]

[0068]

[0069]

[0070]

[0071]

[0072]

[0073]

[0074]

[0075]

[0076]

[0077]

[0078]

[0079] S2. Estimate the initial values of the parameters and set the upper and lower limits. The parameters include: the specific yield μ at the intake value surface ha , the saturated hydraulic conductivity K z , the infiltration recharge rate W at the intake value surface, the intake value head h a , the elevation z of the initial intake value surface ha0 and the starting time t0, as shown in Table 2.

[0080] Table 2 Initial values of the parameters to be inverted and their upper and lower limits

[0081] <![CDATA[μ ha > <![CDATA[K z (m / s)]]> W (m / s) <![CDATA[h a (m)]]> <![CDATA[z ha0 (m)]]> <![CDATA[t0(s)]]> Initial value 0.01 <![CDATA[3.5×10 -5 > 0.0001 -0.07 0.8 10 Lower limit value 0 <![CDATA[3×10 -9 > -0.001 -0.15 0.7 -500 Upper limit value 0.13 <![CDATA[1×10 -4 > 0.001 0.01 1.3 500

[0082] S3. Substitute the initial values of the parameters in Table 2 above into the calculation formula for the elevation z of the groundwater intake value surface under the periodic head drive with the free surface located at the intake value surface to calculate the elevation z of the intake value surface ha , and the formula is:

[0083]

[0084] In the formula, A is the amplitude of the periodic head, ω is the frequency of the periodic head; t is the time;

[0085] S4. Substitute the above calculation results into the total head formula with the free surface located at the intake value surface to calculate the total head H, and the formula is:

[0086]

[0087] S5. Input the observed value of the total head H at the observation point, Table 1.

[0088] S6. Calculate the objective function value f through the objective function formula, and the formula is:

[0089]

[0090] where H c0.5j , H t0.5j , H c0.7j , H t0.7j are respectively the calculated values and the observed values of the total head at 0.5 m and 0.7 m above the elevation of the observation point at time j th .

[0091] S7. Judge whether the termination condition is satisfied according to the calculated objective function value f; if not, regenerate the parameter values, and substitute the regenerated parameter values into step S3 and repeat the operations from step S3 to step S7.

[0092] S8. If the termination condition is satisfied, output the optimized parameters: μ ha , K z , W, h a , z ha0 , t0, and obtain the inversion result. After iterative calculation, the objective function value is 0.036, and the parameter inversion result is shown in Table 3.

[0093] Table 3 Parameter Inversion Results

[0094] <![CDATA[μ ha > <![CDATA[K z (m / s)]]> W (m / s) <![CDATA[h a (m)]]> <![CDATA[z ha0 (m)]]> <![CDATA[t0(s)]]> Lower limit value 0.07 <![CDATA[2.3×10 -5 > <![CDATA[1.7×10 -7 > -0.99 1.09 8.4

[0095] Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art and related fields without creative efforts shall fall within the protection scope of the present invention. The structures, devices, and operation methods not specifically described and explained in the present invention shall be implemented by conventional means in the art without special instructions and limitations.

Claims

1. A method for inverse calculation of the elevation and parameters of the air intake surface driven by periodic water head, characterized in that, It includes the following steps: S1. Estimate the initial values of the parameters and set the upper and lower limits. The parameters include: the specific yield μ at the intake value surface ha , the saturated hydraulic conductivity K z , the infiltration recharge rate W at the intake value surface, the intake value head h a , the initial elevation z of the intake value surface ha0 and the starting time t0; S2. Substitute the above initial parameter values into the calculation formula for the elevation z of the groundwater intake value surface under the periodic head drive with the free surface positioned at the intake value surface to calculate the elevation of the intake value surface ha , and the formula is: where A is the amplitude of the periodic head, ω is the frequency of the periodic head; t is time; S3. Substitute the above calculation results into the total head formula with the free surface located at the intake value surface to calculate the total head H. The formula is: S4. Input the observed value of the total head H at the observation point; S5. Calculate the objective function value f through the objective function formula as follows: In the formula, H cij, H tij are respectively the calculated value and the observed value of the total head at the i th observation point at the j th moment. S6. Determine whether the termination condition is satisfied according to the calculated objective function value f;; S7. If the termination condition is satisfied, output the optimized parameters: μ ha , K z , W, h a , z ha0 , t0, and obtain the inversion result; S8. Substitute the parameters output in step S7 into the calculation formula for the elevation z of the groundwater intake value surface under the periodic water head drive with the free surface located at the intake value surface, and calculate the elevation z of the intake value surface. ha ; S9. Substitute the elevation of the intake value surface calculated in step S8 into the total head formula with the free surface located at the intake value surface to predict the total head profile under periodic head drive and obtain the forward prediction result.

2. A method for inverse calculation of the elevation and parameters of the air intake surface driven by a periodic water head, as claimed in claim 1, wherein: In the step S2, the calculation formula for the elevation of the groundwater intake value surface driven by the periodic water head with the free surface located at the intake value surface can also calculate the elevation z of the intake value surface through the following approximate formula ha The approximate solution of which is given by the formula:

3. A method for inverse calculation of the elevation and parameters of the air intake surface under periodic water head drive according to claim 1, characterized in that: In the step S6, if the termination condition is not satisfied, it is necessary to regenerate the parameter values of μ ha , K z , W, h a , z ha0 , t0, and compare the regenerated parameter values with the upper and lower limits of the parameters to determine whether the upper and lower limits of the parameters are satisfied. If not, regenerate the parameter values; If satisfied, substitute the newly generated parameter values into step S2 and repeat the operations from step S2 to step S6.

4. A method for inverse calculation of the elevation and parameters of the air intake surface driven by a periodic water head according to claim 3, characterized in that: The formula for the parameters to satisfy the upper and lower limits of the parameters is: μ hamin <μ ha <μ hamax K zmin <K z <K zmax W min <W<W max h amin <h a <h amax z ha0min <z ha0 <z ha0max t 0min <t0<t 0max K z > -W where the parameter μ ha , K z , W, h a , z ha0 , the lower and upper limits of t0 are the parameters μ ha , K z , W, h a , z ha0 , respectively; K z > -W restricts the relative magnitude of the saturated hydraulic conductivity K z and the infiltration recharge rate W at the intake value surface.

Images