Hydrodynamic response analytical calculation method for layered heterogeneous cross-flow aquifer system

Abstract

The invention discloses an analytic calculation method for hydrodynamic response of a layered heterogeneous overflow aquifer system. The analytic calculation method comprises the following steps: constructing a hydrodynamic response control equation of a confined aquifer and an aquitard in the layered heterogeneous overflow aquifer system; deriving an aquitard over-flow expression represented by a dynamic water level steady-state component and a wave dynamic component of the confined aquifer; introducing a general water level variable to form a general ordinary differential equation set which simultaneously represents a steady-state component and a wave dynamic component; obtaining an analytical expression of the general water level of the confined aquifer; constructing a recursive relational expression of constant coefficient vectors in the general solution; and solving constant coefficient vectors of the lateral boundary partitions, substituting the constant coefficient vectors into the recursive relational expression to obtain an analytic solution constant coefficient vector of each horizontal partition, and obtaining a hydrodynamic response analytic solution of the layered heterogeneous cross-flow aquifer system. According to the method, the spatial distribution of the dynamic water level steady-state component, the wave dynamic component amplitude and the wave dynamic component phase difference of any confined aquifer can be quickly and accurately calculated.

Description

Technical Field

[0001] The present invention relates to the field of hydrogeology, and particularly to an analytical calculation method for hydrodynamic response of a layered heterogeneous leaky aquifer system. Background Art

[0002] The loose overburden layer in the plain area generally has a layered development characteristic, and at the same time, the interbedding of the aquifer and the semi-permeable layer is its typical binary structure characteristic. However, due to the complex and changeable sedimentary environment, the developed layered leaky aquifer system usually has significant horizontal heterogeneity. Clarifying the hydrodynamic response law in the layered heterogeneous leaky aquifer system plays a crucial role in groundwater resource management and engineering groundwater control. For the calculation of the hydrodynamic response of the layered heterogeneous leaky aquifer system, the existing analytical methods usually consider the cases of multi-layer homogeneous or single-layer heterogeneous. For the case of multi-layer and heterogeneous, there is still a lack of analytical calculation methods. To solve this technical bottleneck, the present invention proposes an analytical calculation method for the hydrodynamic response of a layered heterogeneous leaky aquifer system to meet the requirements of rapid and accurate calculation of the hydrodynamic response of the loose overburden layer in the plain area. Summary of the Invention

[0003] The purpose of the present invention is to provide an analytical calculation method for the hydrodynamic response of a layered heterogeneous leaky aquifer system. By constructing a mathematical model for the hydrodynamic response of a layered leaky aquifer system with an interbedding of aquifers and semi-permeable layers of any number of layers, and combining the method of separation of variables, matrix eigenvalue analysis method, and recurrence rule based on interface continuity conditions in the complex number domain, an analytical solution applicable to the hydrodynamic response of the layered heterogeneous leaky aquifer system is derived, so as to achieve rapid and accurate calculation of the spatial distributions of the steady-state component, amplitude of the fluctuation component, and phase difference of the fluctuation component of the dynamic water levels of multiple confined aquifers.

[0004] To achieve the above technical purpose, the technical solution adopted by the present invention is as follows:

[0005] The present invention discloses an analytical calculation method for the hydrodynamic response of a layered heterogeneous leaky aquifer system, and the method includes the following steps:

[0006] Construct the control equations for the hydrodynamic response of the confined aquifer and the semi-permeable layer in the layered heterogeneous leaky aquifer system;

[0007] Using the method of separation of variables in the complex number domain, decompose the dynamic water levels of the confined aquifer and the semi-permeable layer into steady state and wave states;

[0008] Taking the dynamic water level of the confined aquifer as the boundary condition for the vertical flow of the semi-permeable layer, derive the expression for the leakage flow of the semi-permeable layer characterized by the steady-state component and the wave state component of the dynamic water level of the confined aquifer;

[0009] Taking the steady-state component and the wave component of the piezometric aquifer's dynamic water level as the dependent variables respectively, establish an ordinary differential equation system; introduce a general water level variable into the ordinary differential equation system to form a general ordinary differential equation system that simultaneously represents the steady-state component and the wave component, and express it in matrix and vector forms;

[0010] Using the matrix eigenvalue analysis method, obtain the analytical expression of the general water level of the piezometric aquifer;

[0011] According to the water level and flow continuity conditions at the interface of the piezometric aquifer parameter partition, construct a recurrence relation for the constant coefficient vector in the general solution, and extend the recurrence relation to the lateral boundary of the system;

[0012] Obtain the constant coefficient vector of the lateral boundary partition, substitute it into the recurrence relation, calculate the constant coefficient vector of the analytical solution for each horizontal partition, and obtain the analytical solution of the hydrodynamic response of the layered heterogeneous leaky aquifer system;

[0013] Write the corresponding calculation program to calculate the spatial distributions of the steady-state component, the amplitude of the wave component, and the phase difference of the wave component of the piezometric aquifer's dynamic water level.

[0014] Furthermore, the process of constructing the hydrodynamic response control equations for the piezometric aquifer and the semi-pervious layer in the layered heterogeneous leaky aquifer system includes the following steps:

[0015] For a multi-level leaky aquifer system composed of N aquifers and N semi-pervious layers, number the aquifers sequentially from top to bottom, and the semi-pervious layer with the same number is above each aquifer; in the horizontal direction, the system is divided into M heterogeneous parameter partitions, numbered from left to right; within each horizontal partition, the hydrogeological parameters of the aquifer and the semi-pervious layer are considered homogeneous and isotropic, and have a uniform thickness;

[0016] Denote the hydraulic conductivity of the nth aquifer (n,m) located in horizontal partition m as T (n,m) [L 2 T -1 , the storage coefficient as S (n,m) [-], and the dynamic load efficiency coefficient as γ (n,m) [-]; denote the permeability coefficient of the nth semi-pervious layer (n,m) located in horizontal partition m as K′ (n,m) [LT -1 , the storage rate as S s ′ (n,m) [L -1 , and the thickness as b n [L]; L represents length, and T represents time;

[0017] Taking the right direction as the positive x-axis direction, denote the x coordinate of the interface between regions m and m + 1 as x m; For each semi-pervious layer, define a local vertical coordinate z n , with its positive direction upward, 0 ≤ z n ≤ b n ; The coordinate z n = 0 corresponds to the interface between aquifer n and semi-pervious layer n, and z n = b n corresponds to the interface between aquifer n - 1 and semi-pervious layer n; The governing equations for the horizontal Darcy flow in aquifer (n,m) and the vertical Darcy flow in semi-pervious layer (n,m) are respectively expressed as:

[0018]

[0019] where, h (n,m) = h (n,m) (x, t)[L] represents the water level at position x in aquifer (n,m) at time t, x m-1 ≤ x ≤ x m , h′ (n,m) = h′ (n,m) (x, z n , t)[L] represents the water level at horizontal position x and vertical position z n in semi-pervious layer (n,m) at time t; and represent the source / sink terms through the interface between the aquifer and the semi-pervious layer:

[0020]

[0021] Further, by using the method of separation of variables in the complex domain, the process of decomposing the dynamic water levels of the confined aquifer and the semi-pervious layer into steady-state and wave components includes the following steps:

[0022] In the vertical direction, the water level at the interface between the aquifer and the semi-pervious layer is continuous, and the expression is:

[0023]

[0024] In the horizontal direction, on the interface between aquifers (n,m) and (n,m + 1), the head and flux are both continuous, and the expressions are:

[0025]

[0026] where, and respectively represent the left and right sides of x m ;

[0027] Denote the periodic water level at the top boundary of semi-pervious layer (1,m) as h 0,m (t), and express it in complex form as:

[0028]

[0029] Among them, represents the average water level, denoted as h 0,m as the steady-state component; includes the amplitude and phase information of the water level, denoted as h 0,m as the fluctuating component, where \(i\) is the unit imaginary number and \(w\) is the frequency.

[0030] Furthermore, the process of deriving the expression for the leakage flow rate of the semi-pervious layer characterized by the steady-state component and the fluctuating component of the piezometric aquifer water level includes the following steps:

[0031] At the left and right boundaries of the system, the fluctuating water level or constant flow boundary conditions are:

[0032]

[0033] Among them, is the constant or known function corresponding to the left boundary of aquifer n, is the constant or known function corresponding to the right boundary of aquifer n; for the fluctuating water level boundary condition, set and let

[0034]

[0035] Among them, represents the steady-state component of the periodic water level on the left and right boundaries, represents the fluctuating component of the periodic water level on the left and right boundaries;

[0036] For the case where the left and right sides of aquifer n are constant flow boundaries, set and let

[0037]

[0038] Among them, represents the constant flow rate value at the left and right boundaries of aquifer n.

[0039] Furthermore, for the coupled ordinary differential system of the piezometric aquifer water level steady-state component, its matrix-vector form is expressed as:

[0040]

[0041] Among them, is a column vector containing elements, is a column vector with all elements being T mis an (N×N) diagonal matrix with diagonal elements being T (n,m) , E m is an (N×N) matrix with only one non-zero element located at E m (1,1) = K′ (1,m) / b 1 , D m is an (N×N) tridiagonal matrix with its non-zero elements expressed as:

[0042]

[0043] For the coupled ordinary differential system of the piezometric head fluctuation components, its matrix - vector form is expressed as:

[0044]

[0045] where, is a column vector containing elements, is a column vector with all elements being , S m and γ m are (N×N) diagonal matrices with diagonal elements being S (n,m) and γ (n,m) , G m is an (N×N) diagonal matrix with diagonal elements being:

[0046] G m (1,1) = f (1,m) +(g (1,m) -f (1,m) )γ′ (1,m) +(g (2,m) -f (2,m) )γ′ (2,m)

[0047] G m (n,n) = (g (n,m) -f (n,m) )γ′ (n,m) +(g (n+1,m) -f (n+1,m) )γ′ (n+1,m)

[0048] G m (N,N) = (g (N,m) -f (N,m) )γ′ (N,m)

[0049] F m is a tridiagonal matrix with diagonal elements being:

[0050] F m (1,1:2) = [g(1,m) +g (2,m) ,-f (2,m)

[0051] F m (n,n - 1:n + 1)=[-f (n,m) ,g (n,m) +g (n+1,m) ,-f (n+1,m)

[0052] F m (N,N - 1:N)=[-f (N,m) ,g (N,m) 。

[0053] Furthermore, a general water level variable is introduced into the ordinary differential equations where α = s represents the steady - state component and α = p represents the fluctuating component, forming a general ordinary differential equation system that simultaneously characterizes the steady - state component and the fluctuating component:

[0054]

[0055] Among them, for the steady - state component α = s, there is For the fluctuating component α = p, there is

[0056] Furthermore, the process of obtaining the analytical expression of the general water level in the confined aquifer by using the matrix eigenvalue analysis method includes the following steps:

[0057] The particular solution of the general ordinary differential equation system is calculated as:

[0058]

[0059] Solve the following eigenvalue problem to obtain the eigenvalues of the matrix and eigenvectors

[0060]

[0061] where I is the identity matrix; the eigenvectors are combined into an eigenvector matrix in the following way

[0062]

[0063] Then the general solution of the general ordinary differential equation system is:

[0064]

[0065] Among them, ​​Let \((2N\times1)\) be the constant coefficient column vector to be solved. is an \((N\times2N)\) matrix, and its non-zero elements are related to the eigenvalue and is expressed as

[0066]

[0067] Furthermore, according to the water level and flow continuity conditions at the interface of the confined aquifer parameter partition, the recurrence relation of the constant coefficient vector in the general solution is constructed as:

[0068]

[0069] where

[0070]

[0071] Furthermore, the process of obtaining the analytical solution of the hydrodynamic response of the layered heterogeneous leaky aquifer system includes the following steps:

[0072] Expand the recurrence relation to the left and right lateral boundaries of the system, and the relation of the constant coefficient column vectors corresponding to the heterogeneous parameter partitions 1 and M is:

[0073]

[0074] According to the lateral boundary conditions expressed in the fluctuating water level or constant flow boundary conditions, we get:

[0075]

[0076] where and are \((N\times1)\) column vectors, is an \((N\times2N)\) matrix, and the non-zero elements of the matrix and the column vector are expressed as:

[0077]

[0078] The following expression for the constant coefficient column vector is obtained:

[0079]

[0080] where is a column vector composed of respectively; Substitute the obtained into the recurrence relation of the constant coefficient vector in the general solution to obtain the constant coefficient column vector corresponding to each parameter partition. Then substitute into the general solution of the general ordinary differential equation system to obtain the closed solution of the aquifer water level ;

[0081] After obtaining the analytical solution of the wave dynamics, the amplitude and phase difference of the fluctuating water level are calculated by the following formulas:

[0082]

[0083] wherein, denotes the modulus of, respectively denote the imaginary and real parts of.

[0084] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0085] First, the analytical calculation method for the hydrodynamic response of the layered heterogeneous leaky aquifer system of the present invention can handle the layered leaky system with any number of aquifers compared with the existing analytical methods, and at the same time supports any number of parameter partitions in the horizontal extension direction, meeting the analytical calculation requirements under complex heterogeneous conditions, and has strong versatility.

[0086] Second, the analytical calculation method for the hydrodynamic response of the layered heterogeneous leaky aquifer system of the present invention is applicable to various hydrogeological conditions, and can be used for analyzing steady seepage analysis, and can also be used for analyzing the hydrodynamic response process under the action of a fluctuating water level boundary.

[0087] Third, for the analytical calculation method for the hydrodynamic response of the layered heterogeneous leaky aquifer system of the present invention, the lateral boundary can be uniformly set as a periodic water level or a constant flow boundary, or can be set as a periodic water level or a constant flow boundary for different aquifers respectively; the periodic boundary at the top can be uniformly set for all parameter partitions, or can be set for different parameter partitions respectively, and the boundary forms are diverse. Description of the Drawings

[0088] Figure 1 is a seepage concept diagram of a multi-level heterogeneous leaky aquifer system;

[0089] Figure 2 is a spatial distribution diagram of the steady-state component of the dynamic water level in a confined aquifer;

[0090] Figure 3 is a spatial distribution diagram of the amplitude of the fluctuating component of the dynamic water level in a confined aquifer;

[0091] Figure 4 is a spatial distribution diagram of the phase difference of the fluctuating component of the dynamic water level in a confined aquifer;

[0092] Figure 5 is a flow chart of an analytical calculation method for the hydrodynamic response of a layered heterogeneous leaky aquifer system disclosed by the present invention. Detailed Embodiments

[0093] The embodiments of the present invention will be further described in detail below with reference to the accompanying drawings.

[0094] See Figure 5 , the present invention discloses an analytical calculation method for the hydrodynamic response of a layered heterogeneous leaky aquifer system, and the main steps include:

[0095] 1. Establish a seepage control equation: Considering that the confined aquifer is a horizontal flow and the semi-pervious layer is a vertical flow, construct the hydrodynamic response control equations for the confined aquifer and the semi-pervious layer in the layered heterogeneous leaky aquifer system;

[0096] 2. Use the method of separation of variables in the complex domain: Considering the effect of the periodically fluctuating water level boundary, the water level in the leaky aquifer system also shows periodic fluctuation characteristics. Use the method of separation of variables in the complex domain to decompose the dynamic water levels of the confined aquifer and the semi-pervious layer into steady-state and wave components;

[0097] 3. Derive the analytical expression for the vertical flow of the semi-pervious layer: Taking the dynamic water level of the confined aquifer as the boundary condition for the vertical flow of the semi-pervious layer, derive the analytical expression for the dynamic water level of the semi-pervious layer characterized by the steady-state component and the wave component of the dynamic water level of the confined aquifer;

[0098] 4. Establish an expression for the leakage flow: Based on the analytical expression for the vertical flow, derive the expression for the leakage flow of the semi-pervious layer characterized by the steady-state component and the wave component of the dynamic water level of the confined aquifer;

[0099] 5. Construct a system of ordinary differential equations: Taking the steady-state component and the wave component of the dynamic water level of the confined aquifer as the dependent variables respectively, establish a system of ordinary differential equations, and express them in the form of matrix and vector operations;

[0100] 6. Construct a general system of ordinary differential equations: According to the systems of ordinary differential equations corresponding to the steady-state component and the wave component, introduce a general water level variable to form a general system of ordinary differential equations that can characterize both the steady-state component and the wave component;

[0101] 7. Eigenvalue analysis: Use the matrix eigenvalue analysis method to obtain the analytical expression for the general water level of the confined aquifer;

[0102] 8. Construct a recurrence relation: According to the water level and flow continuity conditions at the parameter partition interface of the confined aquifer, construct a recurrence relation for the constant coefficient vector in the general solution;

[0103] 9. Introduce boundary conditions for solution: Expand the recurrence relation to the lateral boundary of the system, and combine the given boundary conditions to obtain the constant coefficient vector for the lateral boundary partition;

[0104] 10. Analytical solution generation: Substitute the constant coefficient vector obtained from the lateral boundary conditions into the recurrence relation to obtain the analytical solution constant coefficient vector for each horizontal partition, thereby obtaining the analytical solution of the hydrodynamic response of the layered heterogeneous leaky aquifer system.

[0105] 11. Computational program development: Write the corresponding computational program according to the analytical solution to achieve fast and accurate calculation of the spatial distribution of the steady-state component of the dynamic water level, the amplitude of the dynamic component, and the phase difference of the dynamic component in any confined aquifer.

[0106] The analytical calculation method of the present invention has the following beneficial effects:

[0107] 1. Strong generality: Compared with the existing analytical methods, it can handle the layered leaky system with any number of aquifers, and at the same time supports any number of parameter partitions in the horizontal extension direction, meeting the analytical calculation requirements under complex heterogeneous conditions;

[0108] 2. High flexibility: Applicable to a variety of hydrogeological conditions, it can be used for analyzing steady seepage, and can also be used for analyzing the hydrodynamic response process under the action of fluctuating water level boundaries;

[0109] 3. Diverse boundary forms: The lateral boundary can be uniformly set as a periodic water level or a constant flow boundary, or can be set as a periodic water level or a constant flow boundary for different aquifers respectively; the periodic boundary at the top can be uniformly set for all parameter partitions or can be set for different parameter partitions respectively.

[0110] Consider a layered leaky aquifer system composed of N aquifers and n semi-pervious layers. The heterogeneity is reflected in that the system can be divided into multiple parameter partitions within a limited horizontal extension length. The aquifers are numbered sequentially from the top (n = 1) to the bottom (n = N), and above each aquifer is a semi-pervious layer with the same number (see Figure 1 ). In the horizontal direction, the system is divided into M heterogeneous parameter partitions, and the horizontal partitions are numbered from left to right, successively from m = 1 to m = M. Within each horizontal partition, the hydrogeological parameters of the aquifers and semi-pervious layers are considered homogeneous and isotropic and have uniform thicknesses. For the sake of convenience of expression, the nth aquifer located in partition m is simply referred to as aquifer (n, m). Similarly, the nth semi-pervious layer located in partition m is simply referred to as semi-pervious layer (n, m). Denote the hydraulic conductivity of aquifer (n, m) as T (n,m) [L 2 T -1 , the storage coefficient as S (n,m) [-], the dynamic load efficiency coefficient as γ (n,m) [-]; denote the permeability coefficient of semi-pervious layer (n, m) as K′ (n,m) [LT -1 , the storage rate as S s ′(n,m) [L -1 , the dynamic load efficiency coefficient is γ′ (n,m) [-], with a thickness of b n [L].

[0111] Taking Figure 1 the right direction shown as the positive x-axis direction, the x-coordinate of the interface between regions m and m + 1 is denoted as x m . For each semi-pervious layer, a local vertical coordinate z n (0 ≤ z n ≤ b n ) is defined, with its positive direction upward. The coordinate z n = 0 corresponds to the interface between aquifer n and semi-pervious layer n, while z n = b n corresponds to the interface between aquifer n - 1 and semi-pervious layer n. The governing equations for the horizontal Darcy flow in aquifer (n,m) and the vertical Darcy flow in semi-pervious layer (n,m) are respectively expressed as:

[0112]

[0113]

[0114] where h (n,m) = h (n,m) (x, t)[L](x m-1 ≤ x ≤ x m ) represents the water level at position x in aquifer (n,m) at time t, and h′ (n,m) = h′ (n,m) (x, z n , t)[L] represents the water level at horizontal position x and vertical position z n in semi-pervious layer (n,m) at time t; and represent the source-sink terms through the interface between the aquifer and the semi-pervious layer, that is, the leakage flux per unit area. Specifically, represents the flux flowing from aquifer (n,m) into the overlying semi-pervious layer (n,m), represents the flux flowing into the underlying semi-pervious layer (n + 1,m). and are respectively expressed as:

[0115]

[0116]

[0117] In the vertical direction, the water levels at the interface between the aquifer and the semi-pervious layer are continuous, and their expression is:

[0118]

[0119]

[0120] Horizontally, on the interface between aquifers (n,m) and (n,m+1), both the hydraulic head and the flux are continuous, and their expressions are:

[0121]

[0122]

[0123] where and represent the left and right sides of x m respectively.

[0124] The periodic water level at the top boundary of the semi-permeable layer (1,m) is denoted as h 0,m (t), and it is expressed in complex form as:

[0125]

[0126] where represents the average water level, i.e., the steady-state component of h 0,m =h 0,m (t), including the amplitude and phase information of the water level, is the fluctuating component of h 0,m (t), is the unit imaginary number, and w is the frequency.

[0127] At the left and right boundaries of the system, considering the fluctuating water level or the constant flow boundary condition, a mixed form is adopted for comprehensive expression:

[0128]

[0129]

[0130] where is the constant or known function corresponding to the left boundary of aquifer n, is the constant or known function corresponding to the right boundary of aquifer n. For the fluctuating water level boundary condition, it can be set that and let

[0131]

[0132]

[0133] where represents the steady-state components of the periodic water levels at the left and right boundaries, Indicates the fluctuating component of the periodic water level at the left and right boundaries.

[0134] For the case where the left and right sides of aquifer n are constant flow boundaries, it can be set And let

[0135]

[0136]

[0137] Among them, Indicates the constant flow values at the left and right boundaries of aquifer n.

[0138] Based on the method of separation of variables, the analytical solutions of the hydrodynamic responses of aquifer (n,m) and semi-pervious layer (n,m) in the complex domain are expressed as:

[0139]

[0140]

[0141] Among them, Respectively represent the steady-state component and the fluctuating component of the dynamic water level h (n,m) Of aquifer (n,m), Respectively represent the steady-state component and the fluctuating component of the dynamic water level h' (n,m) Of semi-pervious layer. It should be noted that Is a real number, while Is a complex number.

[0142] Substitute equations (6) and (10b) into equation (2), and after necessary rearrangement, the steady-state component and the fluctuating component of h' (n,m) Are respectively controlled by the following equations:

[0143]

[0144]

[0145] Among them,

[0146] Considering the water level continuity at the interface between the aquifer and the semi-pervious layer, the boundary conditions for the one-dimensional vertical seepage of semi-pervious layer (n,m) are:

[0147]

[0148]

[0149]

[0150]

[0151] Therefore, the analytical expressions of equations (11a) and (11b) are:

[0152]

[0153]

[0154] Substituting equation (13a) into equation (3a) and equation (13b) into equation (3b), the leakage discharge per unit area from the aquifer (n,m) to the semi - permeable layer (n,m) and (n + 1,m) is:

[0155]

[0156] where,

[0157] f (n,m) = K' (n,m) σ (n,m) csch(σ (n,m) b n ) (14b)

[0158] g (n,m) = K' (n,m) σ (n,m) coth(σ (n,m) b n ) (14c)

[0159] By substituting equations (14a) and (10b) into equation (1), a coupled ordinary differential equation (ODE) system for the steady - state component of the confined aquifer water level can be obtained, and its matrix - vector form is expressed as

[0160]

[0161] where, is a column vector containing elements, is a column vector with all elements being , T m is an (N×N) diagonal matrix with its diagonal elements being T (n,m) , E m is an (N×N) matrix with only one non - zero element located at E m (1,1)= K' (1,m) / b 1 , D m is an (N×N) tridiagonal matrix, and its non - zero elements are expressed as:

[0162]

[0163]

[0164]

[0165] Similarly, for the coupled ordinary differential equation (ODE) system of the fluctuation component, its matrix-vector form is expressed as:

[0166]

[0167] where, is a column vector containing elements, is a column vector with all elements being , S m and γ m are (N×N) diagonal matrices, and the diagonal elements are S (n,m) and γ (n,m) , respectively, G m is an (N×N) diagonal matrix, and the diagonal elements are:

[0168] G m (1,1) = f (1,m) + (g (1,m) - f (1,m) )γ′ (1,m) + (g (2,m) - f (2,m) )γ′ (2,m) (18a)

[0169] G m (n,n) = (g (n,m) - f (n,m) )γ′ (n,m) + (g (n+1,m) - f (n+1,m) )γ′ (n+1,m) (18b)

[0170] G m (N,N) = (g (N,m) - f (N,m) )γ′ (N,m) (18c)

[0171] F m is a tridiagonal matrix, and the diagonal elements are:

[0172] F m (1,1:2) = [g (1,m) + g (2,m) , - f (2,m) (19a)

[0173] F m (n,n - 1:n + 1) = [- f (n,m) , g (n,m) + g (n+1,m) , - f(n+1,m) (19b)

[0174] F m (N,N - 1:N)=[-f (N,m) ,g (N,m) (19c)

[0175] It should be noted that equations (15) and (17) have the same form. Therefore, whether it is the steady - state component or the fluctuation component, they follow the same solution process. Introduce the general water - level variable where α = s represents the steady - state component and α = p represents the fluctuation component. According to this general expression form, equations (15) and (17) can be unified as:

[0176]

[0177] Among them, for the steady - state component α = s, there is For the fluctuation component α = p, there is

[0178] It should be noted that is a constant. Therefore, the particular solution of equation (20) is

[0179]

[0180] Use the matrix eigenvalue analysis method to derive the general solution of the homogeneous part of equation (20). First, solve the following eigenvalue problem to obtain the eigenvalues of the matrix and the eigenvectors

[0181]

[0182] where I is the identity matrix. Combine the eigenvectors into the eigenvector matrix in the following way

[0183]

[0184] According to this equation, the general solution of equation (20) is:

[0185]

[0186] Among them, is a (2N×1) column vector of undetermined constant coefficients, is an (N×2N) matrix, whose non - zero elements are related to the eigenvalues and are expressed as

[0187]

[0188] Based on the water level and flow continuity conditions characterized by equations (5a) and (5b), the relational expressions for the constant coefficient column vectors can be obtained:

[0189]

[0190] where,

[0191]

[0192]

[0193]

[0194]

[0195] where, is an (N×2N) matrix, and its non-zero elements are expressed as:

[0196]

[0197] Using the recurrence rule, further from equation (26), we can get:

[0198]

[0199] where,

[0200]

[0201]

[0202] Extending the above recurrence relational expressions to the left and right lateral boundaries of the system, the relational expressions for the constant coefficient column vectors corresponding to the inhomogeneous parameter partitions 1 and M can be obtained as:

[0203]

[0204] Based on the lateral boundary conditions represented by equations (7a) and (7b), we can get

[0205]

[0206]

[0207] where, and are (N×1) column vectors, is an (N×2N) matrix. The non-zero elements of these matrices and column vectors are expressed as:

[0208]

[0209]

[0210]

[0211]

[0212]

[0213]

[0214] Combining equations (30), (31a) and (31b), we can obtain the following about the constant coefficient column vector The expression is:

[0215]

[0216] in, are column vectors, respectively. It is composed of and depends on the value of α (steady state or fluctuating state). Substituting into equation (28) we can obtain the constant coefficient column vector corresponding to each parameter partition: Then Substituting into equation (24) we can get the aquifer water level The closed solution of .

[0217] After obtaining the analytical solution of the wave dynamics, the amplitude and phase difference of the fluctuating water level can be calculated by the following formula:

[0218]

[0219]

[0220] in, express The model, Respectively The imaginary and real parts of .

[0221] According to the general analytical solution derived above, the corresponding calculation program code can be written to simultaneously obtain the spatial distribution of the steady-state component, the amplitude of the fluctuation component, and the phase difference of the fluctuation component of the dynamic water level of multiple aquifers. The following is a specific case to show the calculation results of the hydrodynamic response of the layered heterogeneous overflow aquifer system. To ensure generality, consider a layered heterogeneous overflow aquifer system (N=3) consisting of three confined aquifers and their adjacent weak permeable layers, where the horizontal extension length of the system is 150m and the thickness of the weak permeable layer is b. n =2m, left boundary coordinate x 0 =50m, right boundary coordinate x M =200m. Assume that at x1 =100m and x 2 = 150m, there are two heterogeneous parameter partition interfaces, so the confined aquifer and the weak permeable layer are divided into three parameter partitions in the horizontal direction (M = 3). In each partition, the hydrogeological parameters of the three confined aquifers and the three weak permeable layers are consistent. The specific parameters include: the hydraulic conductivity of the confined aquifer T (n,1) =50m 2 / d, T (n,2) =100m 2 / d, T (n,3) =150m 2 / d, water storage coefficient S (n,1) =10 -2 ,

[0222] S (n,2) =10 13 , S (n,3) =10 -2 , vertical permeability coefficient of weak permeable layer K′ (n,1) =10 -2 m / d, K′ (n,2) =10 -1 m / d, K′ (n,3) =10 -2 m / d, water storage rate S s ' (n,1) =10 -4 m -1 , S s ' (n,1) =5×10 -4 m -1 , S s ' (n,1) =10 -4 m -1 In terms of boundary condition setting, it is assumed that both sides of the overflow aquifer system are periodic water level boundaries with a period of 0.5d. The steady-state components of the boundary water level are The amplitude is 1m. In addition, it is assumed that the top plate of the weak permeable layer 1 has periodic water levels in the three partitions, and the period is also 0.5d. The steady-state components are The load effect of the top water level change on the hydrodynamic response of the aquifer and the weak permeable layer is not considered. According to the derived analytical solution, the spatial distribution of the steady-state component, the amplitude of the fluctuation component, and the phase difference of the fluctuation component of the dynamic water level of the three confined aquifers is calculated as follows: Figure 2 , Figure 3 , Figure 4 shown.

[0223] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of complete hardware embodiments, complete software embodiments, or embodiments in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code. The scheme in the embodiments of the present application can be implemented in various computer languages, for example, object-oriented programming language Java and literal scripting language JavaScript, etc.

[0224] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0225] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0226] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process in the computer or other programmable device. Figure 1 A process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0227] Although the preferred embodiments of the present application have been described, those skilled in the art may make other changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications falling within the scope of the present application.

[0228] Obviously, those skilled in the art can make various changes and modifications to the present application without departing from the spirit and scope of the present application. Thus, if these modifications and variations of the present application fall within the scope of the claims of the present application and their equivalents, the present application is also intended to include these modifications and variations.

Claims

1. A method for analytical calculation of the hydrodynamic response of a layered heterogeneous leaking aquifer system, characterized in that: The method comprises the following steps: Construct the governing equations of the hydrodynamic response of the confined aquifer and the aquitard in a layered heterogeneous leaking aquifer system; The dynamic water level of confined aquifer and aquitard is decomposed into steady state and fluctuating state by using the separation variable method in complex domain. Taking the dynamic water level of the confined aquifer as the boundary condition for the vertical flow in the aquitard, the expression of the aquitard overflow represented by the steady-state component and the fluctuating component of the dynamic water level of the confined aquifer is derived. Taking the steady-state component and the fluctuating component of the dynamic water level of the confined aquifer as dependent variables, a set of ordinary differential equations is established; a universal water level variable is introduced into the set of ordinary differential equations to form a set of universal ordinary differential equations that simultaneously characterizes the steady-state component and the fluctuating component, and is expressed in the form of matrix and vector; The analytical expression of the universal water level of confined aquifer is obtained by using matrix eigenvalue analysis method. According to the water level and flow continuity conditions of the confined aquifer parameter partition interface, the recursive relationship of the constant coefficient vector in the general solution is constructed and extended to the lateral boundary of the system. The constant coefficient vector of the lateral boundary partition is obtained and substituted into the recursive relationship to calculate the constant coefficient vector of the analytical solution for each horizontal partition, and the analytical solution of the hydrodynamic response of the layered heterogeneous leaking aquifer system is obtained; Write the corresponding calculation program to calculate the spatial distribution of the steady-state component of the dynamic water level of the confined aquifer, the amplitude of the wave-dynamic component, and the phase difference of the wave-dynamic component.

2. The analytical calculation method for hydrodynamic response of layered heterogeneous leaking aquifer system according to claim 1 is characterized in that: The process of constructing the governing equations for the hydrodynamic response of the confined aquifer and the aquitard in a layered heterogeneous leaking aquifer system includes the following steps: For a multi-level leaky aquifer system consisting of N aquifers and N aquitards, the aquifers are numbered from top to bottom, and above each aquifer is an aquitard with the same number; in the horizontal direction, the system is divided into M heterogeneous parameter partitions, and the horizontal partitions are numbered from left to right; in each horizontal partition, the hydrogeological parameters of the aquifer and the aquitard are considered to be homogeneous and isotropic, and have uniform thickness; The hydraulic conductivity of the nth aquifer (n, m) in the horizontal partition m is T (n，m) [L 2 T -1 ], the water storage coefficient is S (n，m) [-], dynamic load efficiency factor is γ (n，m) [-]; The permeability coefficient of the nth weak permeable layer (n, m) in the horizontal partition m is K′ (n，m) [LT -1 ], the water storage rate is S s ' (n，m) [L -1 ], thickness is b n [L]; L represents length, T represents time; With the right direction as the positive direction of the x-axis, the x-coordinate of the interface between regions m and m+1 is marked as x m ; For each aquitard, define a local vertical coordinate z n , its positive direction is upward, 0≤z n ≤ b n ; coordinate z n = 0 corresponds to the interface between aquifer n and aquitard n, z n =b n Corresponding to the interface between aquifer n-1 and aquitard n; the governing equations for the horizontal Darcy flow in aquifer (n, m) and the vertical Darcy flow in aquitard (n, m) are expressed as: Among them, h (n，m) =h (n，m) (x, t)[L] represents the water level at position x in the aquifer (n, m) at time t, x m-1 ≤x≤x m , h′ (n，m) =h′ (n，m) (x,z n , t) [L] represents the horizontal position x and vertical position z of the aquitard (n, m) at time t n The water level at and Represents the source and sink terms through the interface between the aquifer and the aquitard:

3. The analytical calculation method for hydrodynamic response of layered heterogeneous leaking aquifer system according to claim 2 is characterized in that: The process of decomposing the dynamic water level of the confined aquifer and the aquitard into steady state and fluctuating state by using the separation variable method in the complex domain includes the following steps: In the vertical direction, the water level at the interface between the aquifer and the aquitard is continuous, and the expression is: In the horizontal direction, at the interface between aquifers (n, m) and (n, m+1), the hydraulic head and flux are continuous, expressed as: in, and Respectively represent x m the left and right sides; The periodic water level at the top boundary of the aquitard (1, m) is denoted as h 0，m (t), and expressed in plural form as: in, represents the average water level, as h 0，m The steady-state component of Includes the amplitude and phase information of the water level, as h 0，m The fluctuation component of is a unit imaginary number, and w is the frequency.

4. The analytical calculation method for hydrodynamic response of layered heterogeneous leaking aquifer system according to claim 1 is characterized in that: The process of deriving the expression for the aquitard overflow represented by the steady-state component and the fluctuating component of the dynamic water level of the confined aquifer includes the following steps: At the left and right boundaries of the system, the fluctuating water level or constant flow boundary conditions are: in, is a constant or known function corresponding to the left boundary of aquifer n, is a constant or known function corresponding to the right boundary of aquifer n; for the fluctuating water level boundary condition, set And make in, represents the steady-state component of the periodic water level on the left and right boundaries, It represents the fluctuation component of the periodic water level at the left and right boundaries; For the case where the left and right sides of aquifer n are constant flow boundaries, set And make in, Represents the constant flow value at the left and right boundaries of aquifer n.

5. The analytical calculation method for hydrodynamic response of layered heterogeneous leaking aquifer system according to claim 1 is characterized in that: The coupled ordinary differential system of the steady-state component of the water level in the confined aquifer is expressed in matrix-vector form as follows: in, is included A column vector of elements, All elements are Column vector, T m is a (N×N) diagonal matrix whose diagonal elements are T (n，m) , E m is a (N×N) matrix with only one non-zero element, located in E m (1, 1) = K′ (1，m) / b1,D m is a (N×N) tridiagonal matrix whose non-zero elements are expressed as: For the coupled ordinary differential system of the water level fluctuation components of the confined aquifer, its matrix-vector form is expressed as: in, is included A column vector of elements, All elements are Column vector, S m and γ m It is a (N×N) diagonal matrix, and the diagonal elements are S (n，m) and γ (n，m) , G m is a (N×N) diagonal matrix with the following diagonal elements: G m (1, 1)=f (1，m) +(g (1，m) -f (1，m) )c′ (1，m) +(g (2，m) -f (2，m) )c′ (2，m) G m (n, n)=(g (n，m) -f (n，m) )c′ (n，m) +(g (n+1，m) -f (n+1，m) )c′ (n+1，m) G m (N, N)=(g (N，m) -f (N，m) )c′ (N，m) F m is a tridiagonal matrix with the following diagonal elements: F m (1,1:2)=[g (1，m) +g (2，m) ,-f (2，m) ] F m (n,n-1:n+1)=[-f (n，m) ,g (n，m) +g (n+1，m) ,-f (n+1，m) ] F m (N,N-1:N)=[-f (N，m) ,g (N，m) ]。 6. The analytical calculation method for hydrodynamic response of layered heterogeneous leaking aquifer system according to claim 1 is characterized in that: Introducing a universal water level variable into a system of ordinary differential equations Among them, when α=s, it represents the steady-state component, and when α=p, it represents the fluctuating component, forming a general set of ordinary differential equations that simultaneously characterize the steady-state component and the fluctuating component: Among them, for the steady-state component α=s, we have For the fluctuation component α=p, we have 7. The analytical calculation method for hydrodynamic response of layered heterogeneous leaking aquifer system according to claim 6 is characterized in that: The process of obtaining the analytical expression of the universal water level of the confined aquifer using the matrix eigenvalue analysis method includes the following steps: The special solution of the general ordinary differential equation system is calculated as: Solve the following eigenvalue problem and obtain the matrix The eigenvalue of and the eigenvector Where I is the identity matrix; the eigenvector The eigenvector matrix is ​​formed as follows Then the total solution of the general ordinary differential equation system is: in, (2N×1) is the constant coefficient column vector to be determined, is a (N×2N) matrix whose non-zero elements and eigenvalues Related to, expressed as 8. The analytical calculation method for hydrodynamic response of layered heterogeneous leaking aquifer system according to claim 1 is characterized in that: According to the water level and flow continuity conditions at the interface of the confined aquifer parameter partition, the recursive relationship of the constant coefficient vector in the general solution is constructed as follows: in, 9. The analytical calculation method for hydrodynamic response of layered heterogeneous leaking aquifer system according to claim 8 is characterized in that: The process of obtaining an analytical solution to the hydrodynamic response of a layered heterogeneous leaking aquifer system consists of the following steps: Extending the recursive relationship to the left and right lateral boundaries of the system, the relationship between the constant coefficient column vectors corresponding to the inhomogeneous parameter partitions 1 and M is obtained as follows: According to the lateral boundary conditions expressed in the fluctuating water level or constant flow boundary conditions, we obtain: in, and is a (N×1) column vector, is a (N×2N) matrix, and its non-zero elements of the matrix and column vector are represented as: We get the following constant coefficient column vector The expression is: in, are column vectors, respectively. to form; to Substitute the recursive relation of the constant coefficient vector in the general solution to obtain the constant coefficient column vector corresponding to each parameter partition Then Substituting the total solution of the general ordinary differential equation system into the aquifer water level The closed solution of After obtaining the analytical solution of the wave dynamics, the amplitude and phase difference of the fluctuating water level are calculated using the following formula: in, express The model, Respectively The imaginary and real parts of .

Images