Unsaturated water flow and saturated water flow coupling estimation method and system

Abstract

The invention relates to an unsaturated water flow and saturated water flow coupling estimation method and system, and the method comprises the steps: building an underground water flow motion analysis model considering the influence of an unsaturated zone, solving the underground water flow motion analysis model, and obtaining an underground water head driven by the fluctuation of a river water level and a river-underground water interaction flux semi-analytical solution; collecting aquifer water level monitoring data of hydrogeological parameters to be estimated, and performing dimensionless processing on water level time sequences of a plurality of monitoring holes in the aquifer water level monitoring data to obtain dimensionless monitoring water level data; fitting wiring is carried out on the semi-analytic solution of the underground water head and the river-underground water interaction flux obtained by the analytic model and actually measured dimensionless monitoring water level data, a fitting curve meeting the set requirement is selected, and a parameter value corresponding to the fitting curve serves as a corresponding hydrogeological parameter estimated value of the aquifer; according to the method, the influence of unsaturated zone water flow is fully considered, so that the analytical model is closer to the actual situation.

Description

Technical Field

[0001] The present invention relates to the technical field of predicting the groundwater level in the riparian zone, and more specifically, to a method and system for coupling estimation of unsaturated flow and saturated flow. Background Art

[0002] Quantitatively predicting the changes in the groundwater level in the riparian zone is of great significance for water resources assessment and environmental protection. The traditional groundwater level prediction methods have complex modeling methods and high data requirements, and usually do not consider the influence of the unsaturated flow process. The unsaturated zone is the area below the ground surface and above the water table, where complex physical, chemical, and biological processes occur, with complex moisture states and migration mechanisms, which have an important impact on the water flow exchange process between surface water and groundwater in the riparian zone. Therefore, quantitatively studying the flow process in the unsaturated zone is crucial for water resources assessment and ecological environment protection.

[0003] To predict the water level changes in the unsaturated zone, simulating the coupled saturated-unsaturated flow process is widely adopted as a simple and effective method; this method is simpler and easier to use compared to the traditional groundwater level prediction methods in the riparian zone, and describes the hydrological process in the riparian zone more completely. However, the existing methods for simulating the coupled saturated-unsaturated flow process lack a more detailed description of the physical process of groundwater flow in the unsaturated zone, and the unsaturated zone has an undeniable influence on groundwater flow, thus affecting the accuracy of predicting the groundwater level in the riparian zone.

[0004] To solve this problem, a method and system for coupling estimation of unsaturated flow and saturated flow that can accurately predict the groundwater level in the riparian zone (including the unsaturated zone) are needed. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a method for coupling estimation of unsaturated flow and saturated flow in view of the above-mentioned defects of the prior art, and also provide a system for coupling estimation of unsaturated flow and saturated flow.

[0006] The technical solution adopted by the present invention to solve its technical problems is as follows:

[0007] Construct a method for coupling estimation of unsaturated flow and saturated flow, wherein the method includes the following steps:

[0008] Establish an analytical model of groundwater flow movement considering the influence of the unsaturated zone, and solve the analytical model of groundwater flow movement to obtain semi-analytical solutions of the groundwater head and river-groundwater interaction flux driven by river water level fluctuations;

[0009] Collect the water level monitoring data of the aquifer for the hydrogeological parameters to be estimated, and perform dimensionless processing on the water level time series of multiple monitoring wells in the aquifer water level monitoring data to obtain dimensionless monitored water level data;

[0010] Fit and match the semi-analytical solutions of the groundwater head and river-groundwater interaction flux obtained from the analytical model with the measured dimensionless monitored water level data, select the fitting curve that meets the set requirements, and take the parameter value corresponding to the fitting curve as the estimated value of the corresponding hydrogeological parameter of the aquifer.

[0011] In the unsaturated-saturated flow coupling estimation method described in the present invention, the establishment of an analytical model of groundwater flow movement considering the influence of the unsaturated zone and the solution of the analytical model of groundwater flow movement to obtain the semi-analytical solutions of the groundwater head and river-groundwater interaction flux driven by river water level fluctuations include:

[0012] For the part of groundwater flow driven by river water level changes, the one-dimensional Richards equation and the two-dimensional basic differential equation of groundwater flow are respectively selected to describe the water flow movement in the unsaturated zone and the saturated zone;

[0013] For the interface between the saturated zone and the unsaturated zone, the equality of water head and the equality of flux are used for constraint;

[0014] Using mathematical analysis methods, the semi-analytical solutions of the groundwater head and river-groundwater interaction flux are solved.

[0015] In the unsaturated-saturated flow coupling estimation method described in the present invention, the analytical model of groundwater flow movement is a one-dimensional unsaturated zone and two-dimensional saturated zone coupling profile model of a near-river unconfined aquifer, and the analytical model of groundwater flow movement satisfies the following conditions:

[0016] 1) Take the horizontal right direction as the positive direction, and the zero point is located at the intersection of the downward extension of the river bank zone and the impermeable interface;

[0017] 2) At the initial moment, the water levels in both the saturated and unsaturated zones are 0;

[0018] 3) The water flow movement in the unsaturated zone is one-dimensional vertical flow;

[0019] 4) The saturated zone is homogeneous, and the water flow movement in it is two-dimensional profile flow, and its hydraulic conductivity does not change with time;

[0020] 5) There are no other source-sink terms except river water recharge.

[0021] In the unsaturated-saturated flow coupling estimation method described in the present invention, the one-dimensional Richards equation and its boundary conditions of the water flow movement in the unsaturated zone are expressed by the formula:

[0022]

[0023] where K z is the permeability coefficient in the vertical direction, k(θ) is the relative hydraulic conductivity, u(z,t) is the total head at a certain position in the unsaturated zone, θ is the volumetric water content of the soil, C(θ) = dθ / dψ, C(θ) is the water capacity, ψ is the pressure head in the unsaturated zone, ξ is the position of the moving water table, and H u is the elevation at the ground surface;

[0024] The basic differential equation of two-dimensional groundwater flow and its initial and boundary conditions for the saturated zone water flow movement are expressed by the formula:

[0025]

[0026] h(x,z,0) = 0(2b)

[0027] h(x,z,t) = H(t), x = 0(2c)

[0028]

[0029] where K x is the permeability coefficient in the horizontal direction, h(x,z,t) is the total head at a certain position in the saturated zone, S s is the storage coefficient, L is the distance from the river to the watershed, -H s is the position of the bottom of the aquifer, and H(t) is the change in the river water level;

[0030] The control equation at the interface between the unsaturated zone and the saturated zone is described by the formula based on the equality of the head and flux at this position:

[0031] h(x,z,t) = u(z,t), z = ξ(3a)

[0032]

[0033] For the convenience of solving the analytical solution, linearization operations are performed on the control equations (1a), (1b), (2a), (2b), (2c), (2d), (2e), (3a), and (3b) to make the water table position at z = 0, and the Gardner-Kozeny water characteristic curve model is applied to the Richards equation. The simplified control equations are:

[0034]

[0035] h(x,z,0) = 0(4d)

[0036] h(x,z,t) = H(t), x = 0(4e)

[0037]

[0038] h(x,z,t) = u(z,t), z = 0(4h)

[0039]

[0040] In the formula, and are the zero - order approximations of the relative hydraulic conductivity k(θ0) and water capacity C(θ0) at the initial soil volumetric water content θ0 respectively. α c is the soil moisture retention index, α k is the relative hydraulic conductivity index, S y is the specific yield, d = ψ a - ψ k ψ a and ψ k are the pressure heads at the point where air begins to enter the saturated medium and the point where the relative hydraulic conductivity begins to equal 1 respectively;

[0041] For the convenience of subsequent solution, non - dimensional transformation is carried out on the governing equations (4a), (4b), (4c), (4d), (4e), (4f), (4g), (4h) and (4i). The definitions of non - dimensional variables can be adopted by the formula:

[0042]

[0043] Using equations (4a), (4b), (4c), (4d), (4e), (4f), (4g), (4h) and (4i), substituting equation (5) into them, the governing equations after non - dimensional transformation can be adopted by the formula:

[0044]

[0045] h D (x D ,z D ,0) = 0(6d)

[0046] h D (x D ,z D ,t D ) = H D (t D ), x D = 0(6e)

[0047]

[0048] h D (x D ,z D ,t D ) = u D (z D,t D ),z D = 0(6h)

[0049]

[0050] Based on the Laplace transform, h D (x D ,z D ,t D ) and u D (z D ,t D ) are transformed from the time domain to the Laplace domain and expressed in the following integral form:

[0051]

[0052]

[0053] In the formula, the variable superscript "-" represents the Laplace transform, p is the Laplace transform parameter, corresponding to time t;

[0054] To homogenize the saturation zone control equation, the variable substitution form can be adopted using the formula:

[0055]

[0056] Homogenize the control equation (6) and apply the Laplace transform to obtain a new control equation, expressed as:

[0057]

[0058] In the formula,

[0059] Using the method of integral transformation, and are transformed to eliminate the x term in the control equations (9a) and (9c). The transformation form adopted is:

[0060]

[0061]

[0062] In the formula, ψ(ω m ,x D ) is the kernel, which can be expressed as:

[0063]

[0064] Using equations (9e) and (9f), substitute equation (10c) to solve for the kernel ψ(ω m ,xD );

[0065]

[0066] In the formula,

[0067] The newly obtained control equation after integral transformation is:

[0068]

[0069] The general solutions of equations (11a) and (11c) are respectively:

[0070]

[0071] In the formula, J n and Y n are the Bessel functions of the first kind and the second kind of order n, respectively;

[0072] Using equations (11b), (11d), (11e) and (11f), determine the expressions of C1(ω m ), C2(ω m ), C3(ω m ) and C4(ω m ):

[0073]

[0074] In the formula, the expressions of the introduced intermediate variables are:

[0075]

[0076] P = α kD J n (B) + γnJ n (B) - γBJ n+1 (B) (13g)

[0077] Q = α kD Y n (B) + γnY n (B) - γBY n+1 (B) (13h)

[0078]

[0079] The water head at a specific position in the saturated zone and the water head at a specific position in the unsaturated zone in the Laplace domain adopt the formula: Adopt the formula:

[0080]

[0081] The input signal of the aquifer system where the groundwater flow is located is the river water fluctuation on the left side of the aquifer, and the output signal is the interaction flux between the aquifer and the river; in the time domain, the interaction flux Q(t) (L 2 T -1 ) is expressed by the formula:

[0082]

[0083] Based on the dimensionless transformation of Equation (5), the dimensionless form of the interaction flux is expressed as:

[0084]

[0085] The expression of the interaction flux in the Laplace domain is given by the formula:

[0086]

[0087] Based on the result of Equation (14a), the Laplace transform form of the interaction flux is expressed as:

[0088]

[0089] For the coupled estimation method of unsaturated flow and saturated flow described in the present invention, wherein, the groundwater head obtained from the analytical model and the semi-analytical solution of the river-groundwater interaction flux are fitted with the measured dimensionless monitoring water level data, and the fitting curve that meets the set requirements is selected and the parameter values corresponding to the fitting curve are used as the estimated values of the corresponding hydrogeological parameters of the aquifer, including:

[0090] Draw an image of the water level change over time of the measured dimensionless monitoring water level data;

[0091] Use the semi-analytical solutions of the groundwater head and the river-groundwater interaction flux to gradually adjust the various hydrogeological parameters involved for fitting until the analytical solution can fit the measured dimensionless monitoring water level data;

[0092] The parameter values corresponding to the analytical solution are used as the estimated values of the corresponding hydrogeological parameters of the aquifer.

[0093] A coupled estimation system for unsaturated flow and saturated flow, used to implement the coupled estimation method of unsaturated flow and saturated flow as described above, wherein the system includes a model construction unit, a monitoring data acquisition unit, and a fitting unit;

[0094] The model construction unit is used to establish an analytical model of groundwater flow movement considering the influence of the unsaturated zone, and solve the analytical model of groundwater flow movement to obtain the semi-analytical solutions of the groundwater head driven by river water level fluctuations and the river-groundwater interaction flux;

[0095] The monitoring data acquisition unit is used to collect the water level monitoring data of the aquifer for the hydrogeological parameters to be estimated, and perform dimensionless processing on the water level time series of multiple monitoring wells in the aquifer water level monitoring data to obtain dimensionless monitored water level data;

[0096] The fitting and matching unit is used to fit and match the semi-analytical solutions of the groundwater head and river-groundwater interaction flux obtained from the analytical model with the measured dimensionless monitored water level data, select the fitting curve that meets the set requirements, and use the parameter values corresponding to the fitting curve as the estimated values of the corresponding hydrogeological parameters of the aquifer.

[0097] The beneficial effects of the present invention are as follows: Based on the method of saturated-unsaturated zone coupled flow simulation, the present invention derives the response relationship of the near-river unconfined aquifer system to river water fluctuations, namely head changes and interaction fluxes, according to the existing experimental observation data, and performs wiring with the water level changes of the measured data. The hydrogeological parameters of the corresponding aquifer (including the unsaturated zone) are determined by using the optimal fitting and wiring, connecting the measured results with the theoretical results, and fully exploiting the utilization value of the existing data;

[0098] The present invention fully considers the influence of unsaturated zone flow, making the analytical model closer to the actual situation. The present invention is reasonable and reliable, the derivation is rigorous, with a certain degree of innovation, and the operation of the present invention is simple and easy to apply, providing a new method and idea for the prediction of riverbank water levels. Description of the Drawings

[0099] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will further illustrate the present invention in conjunction with the drawings and embodiments. The drawings in the following description are only partial embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts:

[0100] Figure 1 It is a flowchart of the coupled estimation method for unsaturated flow and saturated flow in the preferred embodiment of the present invention;

[0101] Figure 2 It is a sectional conceptual model diagram of saturated-unsaturated flow driven by river water fluctuations;

[0102] Figure 3 It is a schematic diagram of the water level change and interaction flux change of the observation points of the analytical solution under different amplitude conditions;

[0103] Figure 4 It is a schematic diagram of the water level change and interaction flux change of the observation points of the analytical solution under different period conditions;

[0104] Figure 5 It is the best fitting wiring diagram of the water level changes of different monitoring wells and the analytical solution;

[0105] Figure 6 It is a schematic block diagram of a non - saturated flow and saturated flow coupling estimation system according to a preferred embodiment of the present invention. Specific embodiments

[0106] In order to make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be described clearly and completely below. Obviously, the described embodiments are partial 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 without creative efforts shall fall within the protection scope of the present invention.

[0107] The non - saturated flow and saturated flow coupling estimation method according to a preferred embodiment of the present invention is as Figure 1 shown, and referring to Figures 2 - 5 simultaneously, the method includes the following steps:

[0108] S01: Establish an analytical model of groundwater flow movement considering the influence of the unsaturated zone, and solve the analytical model of groundwater flow movement to obtain semi - analytical solutions of the groundwater head and river - groundwater interaction flux driven by river - water level fluctuations;

[0109] S02: Collect the monitoring data of the water level of the aquifer for the hydrogeological parameters to be estimated, and perform dimensionless processing on the water - level time series of multiple monitoring wells in the monitoring data of the aquifer water level to obtain dimensionless monitored water - level data;

[0110] S03: Fit and match the semi - analytical solutions of the groundwater head and river - groundwater interaction flux obtained from the analytical model with the measured dimensionless monitored water - level data, select the fitting curve that meets the set requirements, and take the parameter values corresponding to the fitting curve as the estimated values of the corresponding hydrogeological parameters of the aquifer;

[0111] Based on the method of saturated - unsaturated zone coupled flow simulation, the present invention deduces the response relationship of the near - river unconfined aquifer system to river - water fluctuations, that is, the head change and interaction flux, according to the existing experimental observation data, and matches it with the water - level change of the measured data. The optimal fitting and matching are used to determine the hydrogeological parameters of the corresponding aquifer (including the unsaturated zone), connecting the measured results with the theoretical results, and fully exploiting the utilization value of the existing data;

[0112] The present invention fully considers the influence of the unsaturated - zone flow, making the analytical model closer to the actual situation. The present invention is reasonable and reliable, with rigorous derivation, having a certain degree of innovation. Moreover, the present invention is simple to operate and easy to apply, providing a new method and idea for the prediction of the water level in the riparian zone.

[0113] Preferably, an analytical model of groundwater flow considering the influence of the unsaturated zone is established, and the semi-analytical solutions of the groundwater head and river-groundwater interaction flux driven by river water level fluctuations obtained by solving the analytical model of groundwater flow include:

[0114] For the part of groundwater flow driven by river water level changes, the one-dimensional Richards equation and the two-dimensional basic differential equation of groundwater flow are respectively selected to describe the water flow in the unsaturated zone and the saturated zone;

[0115] For the interface between the saturated zone and the unsaturated zone, equal head and equal flux are used for constraint;

[0116] Using mathematical analytical methods, the semi-analytical solutions of the groundwater head and river-groundwater interaction flux are obtained.

[0117] Preferably, the analytical model of groundwater flow is a one-dimensional unsaturated zone and two-dimensional saturated zone coupling profile model for the unconfined aquifer near the river, and the analytical model of groundwater flow satisfies the following conditions:

[0118] 1) The horizontal right direction is taken as the positive direction, and the zero point is located at the intersection of the downward extension of the riverbank zone and the impermeable interface;

[0119] 2) The water levels in the saturated and unsaturated zones are both 0 at the initial moment;

[0120] 3) The water flow in the unsaturated zone is one-dimensional vertical flow;

[0121] 4) The saturated zone is homogeneous, and the water flow in it is two-dimensional profile flow, and its hydraulic conductivity does not change with time;

[0122] 5) There are no other source-sink terms except river water recharge.

[0123] The one-dimensional Richards equation and its boundary conditions for water flow in the unsaturated zone are expressed by the formula:

[0124]

[0125] In the formula, K z is the hydraulic conductivity in the z direction (vertical direction) (unit: LT -1 ), k(θ) is the relative hydraulic conductivity (-), u(z,t) is the total head at a certain position in the unsaturated zone (unit: L), θ is the soil volumetric water content (-), C(θ) = dθ / dψ is the water capacity (unit: L -1 ), ψ is the pressure head in the unsaturated zone (unit: L), ξ is the position of the moving water table (unit: L), H u is the elevation at the ground surface (the top of the unsaturated zone) (unit: L);

[0126] The basic differential equation of two-dimensional groundwater flow for the movement of water in the saturated zone and its initial and boundary conditions are expressed by the formulas:

[0127]

[0128] h(x,z,0) = 0(2b)

[0129] h(x,z,t) = H(t), x = 0(2c)

[0130]

[0131] Where K x is the hydraulic conductivity in the x-direction (horizontal direction) (unit: LT -1 ), h(x,z,t) is the total head at a certain position in the saturated zone (unit: L), S s is the specific storage (unit: L -1 ), L is the distance from the river to the watershed (unit: L), -H s is the position of the bottom of the aquifer (unit: L), H(t) is the change in the river water level (unit: L);

[0132] The governing equations describing the interface between the unsaturated zone and the saturated zone are based on the equality of head and flux at this interface and are expressed by the formulas:

[0133] h(x,z,t) = u(z,t), z = ξ(3a)

[0134]

[0135] For the convenience of solving the analytical solution, linearization operations are performed on the governing equations (1a), (1b), (2a), (2b), (2c), (2d), (2e), (3a) and (3b) to make the position of the water table at z = 0, and the Gardner-Kozeny water characteristic curve model is applied to the Richards equation. The simplified governing equations are:

[0136]

[0137] h(x,z,0) = 0(4d)

[0138] h(x,z,t) = H(t), x = 0(4e)

[0139]

[0140] h(x,z,t) = u(z,t), z = 0(4h)

[0141]

[0142] Where and are the zero-order approximations of the relative hydraulic conductivity k(θ0) and water capacity C(θ0) at the initial soil volumetric water content θ0, respectively. α c is the soil moisture retention index (reflecting the water holding capacity of the soil) (unit: L -1 ), α k is the relative hydraulic conductivity index (reflecting the permeability of the soil) (unit: L -1 ), S y is the specific yield (-), d = ψ a - ψ k , ψ a and ψ k are the pressure heads at the points where air begins to enter the saturated medium and where the relative hydraulic conductivity begins to equal 1, respectively;

[0143] For the convenience of subsequent solution, dimensionless transformations are performed on the governing equations (4a), (4b), (4c), (4d), (4e), (4f), (4g), (4h), and (4i). The dimensionless variables can be defined using the formula:

[0144]

[0145]

[0146] Using the equations (4a), (4b), (4c), (4d), (4e), (4f), (4g), (4h), and (4i), substituting equation (5), the governing equations after dimensionless transformation can be expressed by the formula:

[0147]

[0148] h D (x D , z D , 0) = 0 (6d)

[0149] h D (x D , z D , t D ) = H D (t D ), x D = 0 (6e)

[0150]

[0151] h D (x D , z D , t D ) = u D (z D , t D),z D = 0(6h)

[0152]

[0153] Based on the Laplace transform, transform h D (x D ,z D ,t D ) and u D (z D ,t D ) from the time domain to the Laplace domain, and represent them in the following integral form:

[0154]

[0155] In the formula, the variable superscript "-" represents the Laplace transform, p is the Laplace transform parameter, and it corresponds to the time t;

[0156] In order to homogenize the saturation zone control equation, the variable substitution form can adopt the formula:

[0157]

[0158] Homogenize the control equation (6), and use the Laplace transform to obtain a new control equation, which is expressed as:

[0159]

[0160] In the formula,

[0161] Using the method of integral transform, transform and to eliminate the x term in the control equations (9a) and (9c). The transformation form adopted is:

[0162]

[0163] In the formula, ψ(ω m ,x D ) is the kernel, can be expressed as:

[0164]

[0165] Using equations (9e) and (9f), substitute equation (10c) to solve for the kernel ψ(ω m ,x D );

[0166]

[0167] In the formula,

[0168] The newly obtained governing equations after integral transformation are as follows:

[0169]

[0170] The general solutions of equations (11a) and (11c) are respectively:

[0171]

[0172] where J n and Y n are the Bessel functions of the first and second kind of order n respectively;

[0173] Using equations (11b), (11d), (11e) and (11f), the expressions of C1(ω m ), C2(ω m ), C3(ω m ) and C4(ω m ) are determined:

[0174]

[0175] where the expressions of the introduced intermediate variables are:

[0176]

[0177] P = α kD J n (B) + γnJ n (B) - γBJ n+1 (B) (13g)

[0178] Q = α kD Y n (B) + γnY n (B) - γBY n+1 (B) (13h)

[0179]

[0180] The water head at a specific position in the saturated zone in the Laplace domain and the water head at a specific position in the unsaturated zone Adopt the formula:

[0181]

[0182] The input signal of the aquifer system where the groundwater flow is located is the river water fluctuation on the left side of the aquifer, and the output signal is the interaction flux between the aquifer and the river; in the time domain, the interaction flux Q(t) (L 2 T-1 ) Using the formula:

[0183]

[0184] Based on the dimensionless transformation of Equation (5), the dimensionless form of the interaction flux is expressed as:

[0185]

[0186] The expression of the interaction flux in the Laplace domain adopts the formula:

[0187]

[0188] Based on the result of Equation (14a), the Laplace transform form of the interaction flux is expressed as:

[0189]

[0190] Preferably, the hydrogeological data can be obtained by laboratory simulation using the water level data of the monitoring wells in the sand tank experiment. When simulating, time series data with more monitoring wells than the set number are acquired; the sand tank experiment simulates the water level changes in the aquifer under two conditions of periodic river water level fluctuations and segmented river water fluctuations respectively.

[0191] Preferably, the groundwater head obtained from the analytical model and the semi-analytical solution of the river-groundwater interaction flux are fitted with the measured dimensionless monitoring water level data. The fitting curve that meets the set requirements is selected, and the parameter values corresponding to the fitting curve are used as the estimated values of the corresponding hydrogeological parameters of the aquifer, including:

[0192] Draw the image of the water level change over time of the measured dimensionless monitoring water level data;

[0193] Use the semi-analytical solution of the groundwater head and the river-groundwater interaction flux to gradually adjust the various hydrogeological parameters involved for fitting and wiring until the analytical solution can fit the measured dimensionless monitoring water level data;

[0194] The parameter values corresponding to the analytical solution are used as the estimated values of the corresponding hydrogeological parameters of the aquifer.

[0195] As Figures 2 - 5 shown, it is further explained as follows:

[0196] As Figure 2 shown, in the present invention, a groundwater flow model in an unconfined aquifer driven by river water level fluctuations is established. In this model, the left boundary is the river boundary, the right boundary is the impermeable boundary, with a distance L between the left and right. The interface between the unsaturated zone and the saturated zone is the water table. The water flow in the unsaturated zone only flows in the vertical direction, and the water flow in the saturated zone flows in a two-dimensional plane. There is no other source-sink term except river recharge;

[0197] As Figure 3 shown, the present method can be used to predict the change of water head at any observation position and the change of the interaction flux between the river and the aquifer under different amplitude fluctuations of the river water. Figure 3 It is a time series change diagram of one of the observation points.

[0198] As Figure 4 shown, the present method can be used to predict the change of water head at any observation position and the change of the interaction flux between the river and the aquifer under different periodic fluctuations of the river water. Figure 4 It is a time series change diagram of one of the observation points.

[0199] The above steps will be explained below in combination with a specific simulation example of the water level data of the monitoring points in the sand tank experiment.

[0200] Specific implementation case data:

[0201] As Figure 5 shown, using the same set of parameters, the theoretical analytical expressions of the water levels at different observation points are compared and matched with the monitoring data at the same time. When the optimal matching combination appears, the estimated values of each hydrogeological parameter are determined as: d = 0.04m, α c = 5m -1 , α k = 60m -1 , S y = 0.2, S s = 6×10 -3 , K x = 14m / d, K z = 3m / d.

[0202] In addition, due to the difficulty in obtaining hydrological data, the water level simulation data selected by the present method is only for illustrative purposes and is not used for limitation.

[0203] A coupled estimation system for unsaturated flow and saturated flow, which is used to implement the coupled estimation method of unsaturated flow and saturated flow as described above. As Figure 6 shown, the system includes a model construction unit 10, a monitoring data acquisition unit 11, and a fitting and matching unit 12;

[0204] The model construction unit 10 is used to establish an analytical model of groundwater flow movement considering the influence of the unsaturated zone, and solve the analytical model of groundwater flow movement to obtain the semi-analytical solutions of the groundwater head and the river-groundwater interaction flux driven by the river water level fluctuation;

[0205] The monitoring data acquisition unit 11 is used to collect the aquifer water level monitoring data of the hydrogeological parameters to be estimated, and perform dimensionless processing on the water level time series of multiple monitoring wells in the aquifer water level monitoring data to obtain dimensionless monitoring water level data;

[0206] The fitting and matching unit 12 is used to fit and match the semi-analytical solutions of the groundwater head and river-groundwater interaction flux obtained from the analytical model with the measured dimensionless monitoring water level data, select the fitting curve that meets the set requirements, and use the parameter values corresponding to the fitting curve as the estimated values of the corresponding hydrogeological parameters of the aquifer.

[0207] The present invention starts from the Richards equation describing the groundwater flow in the unsaturated zone, and establishes a coupled mathematical model for the saturated-unsaturated zone flow to better reflect the interaction process and mechanism of water between the river and groundwater. By using the method of Laplace transform to obtain the semi-analytical expressions of the head and interaction flux in the Laplace domain, the groundwater level changes including the unsaturated zone can be effectively predicted, which is of great significance for the development of groundwater level prediction work in the riparian zone. Through this method, the hydrogeological process in the riparian zone can be understood more accurately, providing reliable data support for water resources assessment and environmental protection.

[0208] It should be understood that for those of ordinary skill in the art, improvements or transformations can be made according to the above description, and all such improvements and transformations shall fall within the protection scope of the appended claims of the present invention.

Claims

1. A method for coupling estimation of unsaturated water flow and saturated water flow, characterized in that, The method includes the following steps: Establish an analytical model of groundwater flow considering the influence of the unsaturated zone, and solve the analytical model of groundwater flow to obtain semi-analytical solutions of groundwater head and river-groundwater interaction flux driven by river water level fluctuations; Collect the water level monitoring data of the aquifer for the hydrogeological parameters to be estimated, and perform dimensionless processing on the water level time series of multiple monitoring wells in the aquifer water level monitoring data to obtain dimensionless monitored water level data; Fit and match the semi-analytical solutions of groundwater head and river-groundwater interaction flux obtained from the analytical model with the measured dimensionless monitored water level data, select the fitting curve that meets the set requirements, and take the parameter values corresponding to the fitting curve as the estimated values of the corresponding hydrogeological parameters of the aquifer.

2. The unsaturated flow and saturated flow coupling estimation method according to claim 1, characterized in that The establishment of an analytical model of groundwater flow considering the influence of the unsaturated zone and the solution of the analytical model of groundwater flow to obtain semi-analytical solutions of groundwater head and river-groundwater interaction flux driven by river water level fluctuations include: For the part of groundwater flow driven by river water level changes, the one-dimensional Richards equation and the two-dimensional basic differential equation of groundwater flow are respectively selected to describe the water flow movement in the unsaturated zone and the saturated zone; For the interface between the saturated zone and the unsaturated zone, the equality of water head and flux is used for constraint; Using mathematical analytical methods, semi-analytical solutions of groundwater head and river-groundwater interaction flux are solved.

3. The unsaturated-saturated flow coupling estimation method according to claim 2, wherein The analytical model of groundwater flow is a one-dimensional unsaturated zone and two-dimensional saturated zone coupling profile model for a near-river unconfined aquifer, and the analytical model of groundwater flow satisfies the following conditions: 1) Take the horizontal right direction as the positive direction, and the zero point is located at the intersection of the extension of the riverbank zone downward and the impermeable interface; 2) At the initial moment, the water levels in the saturated and unsaturated zones are both 0; 3) The water flow movement in the unsaturated zone is one-dimensional vertical flow; 4) The saturated zone is homogeneous, and the water flow movement in it is two-dimensional profile flow, and its hydraulic conductivity does not change with time; 5) There is no other source-sink term except river water recharge.

4. The unsaturated-saturated flow coupling estimation method according to claim 2, characterized in that The one-dimensional Richards equation and its boundary conditions for water flow movement in the unsaturated zone are expressed by the formula: where K z is the permeability coefficient in the vertical direction, k(θ) is the relative hydraulic conductivity, u(z,t) is the total head at a certain position in the unsaturated zone, θ is the volumetric water content of the soil, C(θ) = dθ / dψ, C(θ) is the water capacity, ψ is the pressure head in the unsaturated zone, ξ is the position of the moving water table, H u is the elevation at the ground surface; The two-dimensional basic differential equation of groundwater flow and its initial and boundary conditions for water flow movement in the saturated zone are expressed by the formula: h(x,z,0) = 0(2b) h(x,z,t) = H(t), x = 0(2c) where K x is the permeability coefficient in the horizontal direction, h(x, z, t) is the total head at a certain position in the saturated zone, S s is the storage rate, L is the distance from the river to the watershed, -H s is the position of the bottom of the aquifer, and H(t) is the change in the river water level; Describe the control equation at the interface between the unsaturated zone and the saturated zone. According to the condition of the equality of water head and flux between the unsaturated zone and the saturated zone, the formula is used: h(x,z,t) = u(z,t), z = ξ(3a) For the convenience of solving the analytical solution, linearization operations are performed on the control equations (1a), (1b), (2a), (2b), (2c), (2d), (2e), (3a) and (3b) to make the position of the water table at z = 0, and the Gardner-Kozeny water characteristic curve model is applied to the Richards equation. The simplified control equations are: h(x,z,0) = 0(4d) h(x,z,t) = H(t), x = 0(4e) h(x,z,t) = u(z,t), z = 0(4h) In the formula, and are the zero-order approximations of the relative hydraulic conductivity k(θ0) and the water capacity C(θ0) at the soil volumetric water content θ0 at rest, respectively. α c is the soil moisture retention index, and α k is the relative hydraulic conductivity index. S y is the specific yield. d = ψ a - ψ k , where ψ a and ψ k are the pressure heads at the point where air starts to enter the saturated medium and the point where the relative hydraulic conductivity starts to equal 1, respectively. For the convenience of subsequent solution, non - dimensional transformation is performed on the governing equations (4a), (4b), (4c), (4d), (4e), (4f), (4g), (4h) and (4i). The definition of non - dimensional variables can adopt the formula: α kD = α k L, α cD = α c L, γ = α kD -α cD , (5) Using equations (4a), (4b), (4c), (4d), (4e), (4f), (4g), (4h) and (4i), substituting equation (5) into them, the governing equations after non - dimensional transformation can adopt the formula: h D (x D ,z D ,0) = 0(6d) h D (x D ,z D ,t D ) = H D (t D ),x D = 0(6e) h D (x D ,z D ,t D ) = u D (z D ,t D ),z D = 0(6h) Based on the Laplace transform, transform h D (x D ,z D ,t D ) and u D (z D ,t D ) from the time domain to the Laplace domain, and represent them in the following integral form: In the formula, the variable top - mark "-" represents the Laplace transform, p is the Laplace transform parameter, which corresponds to time t; In order to homogenize the governing equation of the saturated zone, the form of variable substitution adopted can be expressed by the formula: Performing homogenization on the governing equation (6) and applying the Laplace transform, a new governing equation can be obtained, expressed as: In the formula, Using the method of integral transformation, for and perform the transformation to eliminate the x-term in the governing equations (9a) and (9c). The transformation form adopted is: where ψ(ω m , x D ) is the kernel, which can be expressed as: Using equations (9e) and (9f), substitute equation (10c) and solve for the kernel ψ(ω m ,x D ); In the formula, The newly obtained governing equation after integral transformation is: The general solutions of equations (11a) and (11c) are respectively: wherein, J n and Y n are the Bessel functions of the first and second kind of order n, respectively; Using equations (11b), (11d), (11e) and (11f), determine the expressions for C1(ω m ), C2(ω m ), C3(ω m ) and C4(ω m ): In the formula, the expressions of each introduced intermediate variable are: P = α kD J n (B) + γnJ n (B) - γBJ n+1 (B) (13g) Q = α kD Y n (B) + γnY n (B) - γBY n+1 (B) (13h) Head at a specific position in the saturated zone in the Laplace domain and head at a specific position in the unsaturated zone Adopt the formula: The input signal of the aquifer system where the groundwater flow is located is the river water fluctuation on the left side of the aquifer, and the output signal is the interaction flux between the aquifer and the river; in the time domain, the interaction flux Q(t) (L 2 T -1 ) is calculated using the formula: Based on the non - dimensional transformation of equation (5), the non - dimensional form of the interaction flux is expressed as: The expression of the interaction flux in the Laplace domain adopts the formula: Based on the result of equation (14a), the Laplace transform form of the interaction flux is expressed as:

5. The unsaturated-saturated flow coupling estimation method according to claim 1, wherein Fitting and matching the semi - analytical solutions of the groundwater head and river - groundwater interaction flux obtained from the analytical model with the measured non - dimensionalized monitoring water - level data, selecting the fitting curve that meets the set requirements and taking the parameter values corresponding to the fitting curve as the estimated values of the corresponding hydrogeological parameters of the aquifer, including: Drawing the image of the water - level change with time of the measured non - dimensionalized monitoring water - level data; Using the semi - analytical solutions of the groundwater head and river - groundwater interaction flux to gradually adjust each involved hydrogeological parameter for fitting and matching until the analytical solution can fit the measured non - dimensionalized monitoring water - level data; The parameter values corresponding to the analytical solution are used as the estimated values of the corresponding hydrogeological parameters of the aquifer.

6. A non-saturated flow and saturated flow coupling estimation system for implementing the non-saturated flow and saturated flow coupling estimation method according to any one of claims 1-5, characterized in that The system includes a model construction unit, a monitoring data acquisition unit and a fitting and matching unit; The model construction unit is used to establish an analytical model of groundwater flow movement considering the influence of the unsaturated zone, solve the analytical model of groundwater flow movement to obtain the semi - analytical solutions of the groundwater head and river - groundwater interaction flux driven by river - water - level fluctuations; The monitoring data acquisition unit is used to collect the water - level monitoring data of the aquifer for the hydrogeological parameters to be estimated, perform non - dimensionalization processing on the water - level time series of multiple monitoring wells in the aquifer water - level monitoring data to obtain non - dimensionalized monitoring water - level data; The fitting and matching unit is used to fit and match the semi - analytical solutions of the groundwater head and river - groundwater interaction flux obtained from the analytical model with the measured non - dimensionalized monitoring water - level data, select the fitting curve that meets the set requirements and take the parameter values corresponding to the fitting curve as the estimated values of the corresponding hydrogeological parameters of the aquifer.

Images