Analytic method for seepage calculation of multi-level heterogeneous cross-flow aquifer system

Abstract

The invention discloses an analytical method for seepage calculation of a multi-level heterogeneous overflow aquifer system, which comprises the following steps of: establishing a stable seepage control equation of a multi-level aquifer and an adjacent aquitard system, and converting the stable seepage control equation into an ordinary differential equation set in a matrix and vector form; solving the ordinary differential equation set through a matrix eigenvalue analysis method to obtain a general solution of the water level of the confined aquifer; constructing a recursive relational expression of a constant coefficient vector in a general solution on the basis of the continuity conditions of the water level and flux of a heterogeneous parameter partition interface; expanding the recursive relational expression to a lateral boundary of the system, and obtaining a global analytical solution in combination with a lateral boundary condition; and writing a calculation program of the global analytic solution, and calculating the actual underground water seepage problem. According to the method, rapid and high-precision calculation of multi-layer water level spatial distribution can be realized.

Description

Technical Field

[0001] The invention relates to the technical field of hydrogeology, and in particular to an analytical method for seepage calculation of a multi-level heterogeneous overflow aquifer system. Background Art

[0002] The seepage characteristics of groundwater play a vital role in groundwater resource management, environmental protection and engineering groundwater control. For the seepage calculation of multi-level overflow aquifer systems, the existing analytical methods usually assume that the system is homogeneous or weakly heterogeneous. However, for complex systems in which both the aquifer and the aquitard have significant heterogeneity, traditional analytical methods are difficult to meet the requirements of high-precision calculations. Although numerical simulation methods have advantages in flexibility, their computational costs are high and it is difficult to meet the needs of high-efficiency calculations. Therefore, for multi-level and strongly heterogeneous overflow aquifer systems, there is still a lack of efficient and accurate analytical calculation methods. In order to solve this technical bottleneck, the present invention proposes a seepage calculation method for a multi-level heterogeneous overflow aquifer system based on an analytical solution to meet the dual requirements of calculation accuracy and efficiency in practical applications. Summary of the invention

[0003] The purpose of the present invention is to provide an analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system. By constructing a multi-level system mathematical model consisting of an aquifer and a weakly permeable layer, and combining matrix eigenvalue operations with interface continuity conditions, an analytical solution for seepage suitable for a multi-level, strongly heterogeneous leaking aquifer system is derived, thereby achieving fast and high-precision calculation of the spatial distribution of multi-layer water levels.

[0004] In order to achieve the above technical objectives, the technical solution adopted by the present invention is:

[0005] An analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system, the method comprising the following steps:

[0006] Establish the stable seepage control equation of the multi-level aquifer and adjacent weak permeable layer system and transform it into a set of ordinary differential equations in matrix and vector form;

[0007] The ordinary differential equations are solved by matrix eigenvalue analysis method to obtain the general solution of the water level in the confined aquifer.

[0008] Based on the continuity conditions of water level and flux at the interface of heterogeneous parameter partitions, a recursive relation of the constant coefficient vector in the general solution is constructed;

[0009] The recursive relationship is extended to the lateral boundary of the system, and combined with the lateral boundary conditions, a global analytical solution is obtained;

[0010] Write a calculation program for the global analytical solution and use it to calculate actual groundwater seepage problems.

[0011] Furthermore, the process of establishing the stable seepage control equation of the multi-level aquifer and the adjacent weak permeable layer system and converting it into a set of ordinary differential equations in matrix and vector form includes the following steps:

[0012] Construct the seepage control equations of the multi-layer leaky aquifer system, including the horizontal seepage of the confined aquifer and the vertical seepage of the weak permeable layer;

[0013] Taking the water level of the confined aquifer as the boundary condition for vertical seepage in the aquitard, the analytical expression for vertical seepage is derived.

[0014] Based on the analytical expression of vertical seepage, the expression of aquitard overflow represented by the water level of the confined aquifer is derived.

[0015] Taking the water level of the confined aquifer as the dependent variable, a set of ordinary differential equations is established and expressed in the form of matrix and vector operations.

[0016] Furthermore, the process of constructing the seepage control equation of the multi-level leaking aquifer system includes the following steps:

[0017] 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;

[0018] The conductivity of the nth aquifer (n, m) in the horizontal partition m is T (n,m) [L 2 T -1 ], the permeability coefficient of the nth weak permeable layer (n,m) in the horizontal partition m is K′ (n,m) [LT -1 ], thickness is b n [L], L stands for length, T stands for time;

[0019] 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 nCorresponding to the interface between aquifer n-1 and aquitard n; the governing equations for horizontal Darcy flow in aquifer (n, m) and vertical Darcy flow in aquitard (n, m) are expressed as:

[0020]

[0021] Among them, h (n,m) =h (n,m) (x)[L] represents the water level in the aquifer (n,m), x m-1 ≤x≤x m , h′ (n,m) =h′ (n,m) (x,z n )[L] represents the position z in the aquitard (n,m) n The water level at and They represent the source and sink terms through the interface between the aquifer and the aquitard respectively:

[0022]

[0023] Furthermore, the process of converting the steady seepage control equation into a set of ordinary differential equations in the form of matrix and vector includes the following steps:

[0024] In the vertical direction, the hydraulic head at the interface between the aquifer and the aquitard is continuous, and the expressions are:

[0025]

[0026]

[0027] In the horizontal direction, at the interface between the aquifers (n,m) and (n,m+1), the hydraulic head and flux are continuous, and the expressions are:

[0028]

[0029] in, and Respectively represent x m the left and right sides;

[0030] The water level at the top boundary of the aquitard (1, m) is recorded as At the left and right boundaries of the system, considering the constant water level boundary conditions, we have:

[0031]

[0032] in, and Respectively represent the constant water level of the left and right boundaries, n = 1, 2, ..., N;

[0033] The boundary condition for vertical one-dimensional seepage in a weak permeable layer (n, m) is:

[0034] h′ (n,m) (x,0)=h (n,m) (x)(7a);

[0035] h′ (n,m) (x,b n )=h (n-1,m) (x)(7b);

[0036] The general solution of equation (2) is:

[0037]

[0038] Substituting formula (8) into formulas (3a) and (3b), we obtain:

[0039]

[0040] By substituting formula (9) into formula (1), the coupled ordinary differential system of the stable water level of the confined aquifer is obtained, which is expressed in matrix-vector form as follows:

[0041]

[0042] in, It contains h (n,m) 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) / b 1 , D m is a (N×N) tridiagonal matrix whose non-zero elements are expressed as:

[0043]

[0044] Furthermore, the process of solving the ordinary differential equations by matrix eigenvalue analysis method includes the following steps:

[0045] make Then formula (10) becomes

[0046]

[0047] The particular solution of equation (12) is:

[0048]

[0049] The matrix eigenvalue analysis method is used to derive the general solution of the homogeneous part of equation (12); specifically, the following eigenvalue problem is solved to obtain the matrix ξ m The eigenvalue of and the eigenvector

[0050]

[0051] Where I is the identity matrix; the eigenvector The eigenvector matrix V is composed as follows m :

[0052]

[0053] The overall solution of equation (12) is:

[0054]

[0055] in, (2N×1) is the constant coefficient column vector to be determined, U m (x) is a (N×2N) matrix whose non-zero elements and eigenvalues Related, expressed as:

[0056]

[0057] Furthermore, the process of constructing the recursive relation of the constant coefficient vector in the general solution includes the following steps:

[0058] According to the continuity conditions represented by formulas (5a) and (5b), the constant coefficient column vector of adjacent inhomogeneous parameter partitions is obtained: The relationship is:

[0059]

[0060] in,

[0061] Ω m =Φ m+1 (x m ) -1 Φ m (x m )(19a);

[0062]

[0063]

[0064] in, is a (N×2N) matrix, whose non-zero elements are represented as:

[0065]

[0066] Using the recursive rule, we can further obtain the following equation (18):

[0067]

[0068] in,

[0069]

[0070] Furthermore, the recursive relationship is extended to the lateral boundary of the system, and combined with the lateral boundary conditions, the process of obtaining the global analytical solution includes the following steps:

[0071] The recursive relationship is extended to the lateral boundary of the system, and the constant coefficient vector of the lateral boundary partition is obtained in combination with the given boundary conditions; the constant coefficient vector obtained from the lateral boundary conditions is substituted into the recursive relationship to obtain the analytical solution constant coefficient vector of each horizontal partition, and the analytical solution for seepage calculation of multi-level heterogeneous leaking aquifer system is constructed.

[0072] Furthermore, the process of obtaining a global analytical solution includes the following steps:

[0073] 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:

[0074]

[0075] According to the lateral constant water level boundary conditions expressed by formulas (6a) and (6b), we obtain:

[0076]

[0077] in, and is a column vector composed of the constant water levels on the left and right boundaries, and is a (N×2N) matrix, whose non-zero elements are represented as:

[0078]

[0079] Combining formula (22), formula (23a) and formula (23b), we get the following about the constant coefficient column vector The expression is:

[0080]

[0081] What you want Substituting into formula (20) we can get the constant coefficient column vector corresponding to each parameter partition: Then Substitute into formula (16) to obtain the aquifer water level The closed solution of .

[0082] Compared with the prior art, the present invention has the following beneficial effects:

[0083] First, compared with the existing analytical methods, the analytical method for seepage calculation of a multi-level heterogeneous overflow aquifer system of the present invention can handle groundwater seepage systems composed of any number of aquifers, and at the same time support arbitrary parameter partitioning in the horizontal extension direction, meeting the analytical calculation requirements under complex heterogeneous conditions and having strong adaptability.

[0084] Second, compared with numerical simulation methods, the analytical method for seepage calculation of a multi-level heterogeneous overflow aquifer system of the present invention does not require grid division, has higher calculation accuracy, and significantly improves efficiency.

[0085] Third, the analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system of the present invention is applicable to a variety of hydrogeological conditions. Through parameter adjustment, it can meet the seepage calculation requirements of homogeneous aquifers, impermeable weak permeable layers, single-layer leaking aquifers, etc., and has high flexibility.

[0086] Fourthly, the analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system of the present invention can achieve high-precision approximate calculation by increasing the number of leaking aquifers and adjusting the thickness and permeability coefficient of the weak permeable layer in the case of vertical heterogeneity inside a single aquifer, and has strong approximate calculation capability. BRIEF DESCRIPTION OF THE DRAWINGS

[0087] Figure 1 This is a conceptual diagram of seepage in a multi-level heterogeneous overflow aquifer system of the present invention;

[0088] Figure 2 It is the spatial distribution map of water level in confined aquifer;

[0089] Figure 3 The present invention is a flow chart of an analytical method for seepage calculation of a multi-level heterogeneous overflow aquifer system. DETAILED DESCRIPTION

[0090] The embodiments of the present invention are further described in detail below in conjunction with the accompanying drawings.

[0091] See also Figure 3 The present invention discloses an analytical method for calculating seepage in a multi-level heterogeneous leaking aquifer system, characterized in that the method comprises the following steps:

[0092] Establish the stable seepage control equation of the multi-level aquifer and adjacent weak permeable layer system and transform it into a set of ordinary differential equations in matrix and vector form;

[0093] The ordinary differential equations are solved by matrix eigenvalue analysis method to obtain the general solution of the water level in the confined aquifer.

[0094] Based on the continuity conditions of water level and flux at the interface of heterogeneous parameter partitions, a recursive relation of the constant coefficient vector in the general solution is constructed;

[0095] The recursive relationship is extended to the lateral boundary of the system, and combined with the lateral boundary conditions, a global analytical solution is obtained;

[0096] Write a calculation program for the global analytical solution and use it to calculate actual groundwater seepage problems.

[0097] The specific steps of the analytical method of the present invention are described below.

[0098] The first step is to establish the stable seepage control equation of the multi-level aquifer and adjacent weak permeable layer system and transform it into a set of ordinary differential equations in matrix and vector form. Specifically, it includes:

[0099] (1) Establishing the seepage control equation: Establishing the seepage control equation of the multi-level leaky aquifer system, including the horizontal seepage of the confined aquifer and the vertical seepage of the weak permeable layer;

[0100] (2) Derive the analytical expression of vertical seepage: Taking the water level of the confined aquifer as the boundary condition for the vertical seepage of the weak permeable layer, derive its analytical expression;

[0101] (3) Establish overflow expression: Based on the analytical expression of vertical seepage, the overflow expression of the aquitard represented by the water level of the confined aquifer is derived;

[0102] (4) Constructing a set of ordinary differential equations: Taking the water level of the confined aquifer as the dependent variable, establish a set of ordinary differential equations and express them in the form of matrix and vector operations.

[0103] Consider a multi-layer leaky aquifer system consisting of N aquifers and N aquitards, with a finite horizontal extension. The aquifers are numbered from top (n = 1) to bottom (n = N), and above each aquifer is a quitard 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, m=1 to m=M. In each horizontal partition, the hydrogeological parameters of the aquifer and the weak permeability layer are considered to be homogeneous and isotropic, and have uniform thickness. For the convenience of expression, the nth aquifer located in the partition m is referred to as aquifer (n, m). Similarly, the nth weak permeability layer located in the partition m is referred to as weak permeability layer (n, m). The hydraulic conductivity of the aquifer (n, m) is denoted by T (n,m) [L 2 T -1 ], the permeability coefficient of the aquitard (n,m) is K′(n,m) [LT -1 ], thickness is b n [L].

[0104] by Figure 1 The right direction is the positive direction of the x-axis, and 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 (0≤z n ≤ b n ), with its positive direction pointing upward. n = 0 corresponds to the interface between aquifer n and aquitard n, while z n =b n Corresponding to the interface between aquifer n-1 and aquitard n. The governing equations for horizontal Darcy flow in aquifer (n, m) and vertical Darcy flow in aquitard (n, m) are respectively expressed as:

[0105]

[0106] Among them, h (n,m) =h (n,m) (x)[L](x m-1 ≤x≤x m ) represents the water level in the aquifer (n, m), h′ (n,m) =h′ (n,m) (x,z n )[L] represents the position z in the aquitard (n,m) n The water level at and represents the source and sink term through the interface between the aquifer and the aquitard, that is, the overflow per unit area. Specifically, represents the flux from the aquifer (n,m) into the overlying aquitard (n,m), represents the flux into the underlying aquitard (n+1,m). and They are expressed as:

[0107]

[0108] In the vertical direction, the water level at the interface between the aquifer and the aquitard is continuous, and its expression is:

[0109]

[0110] In the horizontal direction, at the interface between the aquifers (n,m) and (n,m+1), the hydraulic head and flux are continuous, and their expressions are:

[0111]

[0112] in, and Respectively represent x m on the left and right sides.

[0113] The water level at the top boundary of the aquitard (1, m) is recorded as At the left and right boundaries of the system, the constant water level boundary condition is considered, and its expression is:

[0114]

[0115] in, and Represent the constant water levels at the left and right boundaries respectively.

[0116] Considering the water level continuity at the interface between the aquifer and the aquitard, the boundary condition for vertical one-dimensional seepage in the aquitard (n, m) is:

[0117] h′ (n,m) (x,0)=h (n,m) (x)(7a)

[0118] h′ (n,m) (x,b n )=h (n-1,m) (x)(7b)

[0119] The general solution of equation (2) is

[0120]

[0121] Substituting formula (8) into formula (3a) and formula (3b) yields

[0122]

[0123] By substituting formula (9) into formula (1), we can obtain the coupled ordinary differential equation (ODE) system of the stable water level of the confined aquifer, which is expressed in matrix-vector form as follows:

[0124]

[0125] in, It contains h (n,m) 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) / b 1 , Dm is a (N×N) tridiagonal matrix whose non-zero elements are expressed as:

[0126]

[0127] make Then formula (10) becomes

[0128]

[0129] It should be noted that is a constant, so the particular solution of equation (12) is

[0130]

[0131] In the second step, the matrix eigenvalue analysis method is used to obtain the general solution of the water level in the confined aquifer.

[0132] The matrix eigenvalue analysis method is used to derive the general solution of the homogeneous part of equation (12). First, solve the following eigenvalue problem to obtain the matrix ξ m The eigenvalue of and the eigenvector

[0133]

[0134] Where I is the identity matrix. The eigenvector matrix V is composed as follows m :

[0135]

[0136] According to this formula, the total solution of equation (12) is:

[0137]

[0138] in, (2N×1) is the constant coefficient column vector to be determined, U m (x) is a (N×2N) matrix whose non-zero elements and eigenvalues Related to, expressed as

[0139]

[0140] According to the continuity conditions represented by formulas (5a) and (5b), the constant coefficient column vector of adjacent inhomogeneous parameter partitions can be obtained: The relationship is:

[0141]

[0142] in,

[0143] Ω m =Φ m+1 (x m ) -1 Φ m (x m )(19a)

[0144]

[0145] in, is a (N×2N) matrix, whose non-zero elements are represented as:

[0146]

[0147] In the third step, the recursive relation of the constant coefficient vector in the analytical solution is constructed according to the water level and flow continuity conditions at the aquifer parameter partition interface.

[0148] Specifically, using the recursive rule, according to formula (18), we can further obtain

[0149]

[0150] in,

[0151]

[0152] The fourth step is to extend the recursive relationship to the lateral boundary of the system, and combine the given boundary conditions to obtain the constant coefficient vector of the lateral boundary partition; substitute the constant coefficient vector obtained from the lateral boundary conditions into the recursive relationship to obtain the analytical solution constant coefficient vector of each horizontal partition, thereby constructing the analytical solution for the seepage calculation of the multi-level heterogeneous overflow aquifer system.

[0153] By extending the above recursive relationship to the left and right lateral boundaries of the system, we can obtain the relationship between the constant coefficient column vectors corresponding to the inhomogeneous parameter partitions 1 and M as follows:

[0154]

[0155] According to the lateral constant water level boundary conditions expressed by equations (6a) and (6b), we can obtain

[0156]

[0157] in, and is a column vector composed of the constant water levels on the left and right boundaries, and is a (N×2N) matrix, whose non-zero elements are represented as:

[0158]

[0159]

[0160] Combining formula (22), formula (23a) and formula (23b), we can get the following about the constant coefficient column vector The expression is:

[0161]

[0162] What you want Substituting into formula (20) we can get the constant coefficient column vector corresponding to each parameter partition: Then Substituting into formula (16), we can get the aquifer water level The closed solution of .

[0163] The fifth step is to write a seepage calculation program based on the analytical solution to realize the spatial distribution calculation of the water level in any aquifer.

[0164] According to the above analytical solution, the corresponding calculation program code can be written to obtain the spatial distribution of water levels in multiple aquifers at the same time. The following is a specific case to show the seepage calculation results of a multi-level heterogeneous leaking aquifer system. To ensure generality, consider a multi-level leaking 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 x 1 =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, vertical permeability coefficient of aquitard K′ (n,1) =10 - 2 m / d, K′ (n,2) =10 -1 m / d, K′ (n,3) =10 -2 m / d (n=1,2,3). In terms of boundary condition setting, it is assumed that both sides of the overflow aquifer system are constant water level boundaries, and the water levels are In addition, it is assumed that the water level of the top plate of the weak permeable layer 1 in the three partitions is

[0165] According to the above parameter settings, the spatial distribution of water levels in the three confined aquifers is calculated as follows: Figure 2 As shown in the figure. It can be seen that in the multi-level heterogeneous leakage aquifer system, the water level distribution of the three confined aquifers has significant inter-layer differences and zoning characteristics. Specifically, among the three aquifers, the water level of confined aquifer 1 is always the highest, mainly because it has the best conditions for receiving vertical leakage recharge at the constant water level boundary of the top plate. On the contrary, the conditions for receiving vertical leakage recharge of confined aquifers 2 and 3 are relatively poor, and their water levels are lower than those of confined aquifer 1. In terms of zoning characteristics, the water levels in zones 1 and 3 decrease approximately linearly. The confined aquifers in these two zones have a high degree of closure and the leakage is not obvious. Since the permeability coefficient of the weak permeable layer in zone 2 is one order of magnitude higher than that of zones 1 and 3, the leakage effect of the confined aquifer in zone 2 is significantly enhanced, and its water level changes slowly. The above results well reflect the seepage law of the multi-level heterogeneous leakage aquifer system and its sensitivity to the distribution of hydrogeological parameters, which has important theoretical significance for the actual groundwater development and utilization and engineering groundwater control.

[0166] 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.

[0167] 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.

[0168] 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.

[0169] 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 The steps for the functions specified in one or more boxes.

[0170] 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.

[0171] 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. An analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system, characterized in that: The method comprises the following steps: Establish the stable seepage control equation of the multi-level aquifer and adjacent weak permeable layer system and transform it into a set of ordinary differential equations in matrix and vector form; The ordinary differential equations are solved by matrix eigenvalue analysis method to obtain the general solution of the water level in the confined aquifer. Based on the continuity conditions of water level and flux at the interface of heterogeneous parameter partitions, a recursive relation of the constant coefficient vector in the general solution is constructed; The recursive relationship is extended to the lateral boundary of the system, and combined with the lateral boundary conditions, a global analytical solution is obtained; Write a calculation program for the global analytical solution and use it to calculate actual groundwater seepage problems.

2. The analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system according to claim 1 is characterized in that: The process of establishing the stable seepage control equation of the multi-level aquifer and adjacent weak permeable layer system and converting it into a set of ordinary differential equations in matrix and vector form includes the following steps: Construct the seepage control equations of the multi-layer leaky aquifer system, including the horizontal seepage of the confined aquifer and the vertical seepage of the weak permeable layer; Taking the water level of the confined aquifer as the boundary condition for vertical seepage in the aquitard, the analytical expression for vertical seepage is derived. Based on the analytical expression of vertical seepage, the expression of aquitard overflow represented by the water level of the confined aquifer is derived. Taking the water level of the confined aquifer as the dependent variable, a set of ordinary differential equations is established and expressed in the form of matrix and vector operations.

3. The analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system according to claim 2 is characterized in that: The process of constructing the seepage control equations for a multi-layered leaky 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 conductivity of the nth aquifer (n, m) in the horizontal partition m is T (n,m) [L 2 T -1 ], the permeability coefficient of the nth weak permeable layer (n,m) in the horizontal partition m is K′ (n,m) [LT -1 ], thickness is b n [L], L stands for length, T stands for 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 horizontal Darcy flow in aquifer (n, m) and vertical Darcy flow in aquitard (n, m) are expressed as: Among them, h (n,m) =h (n,m) (x)[L] represents the water level in the aquifer (n,m), x m-1 ≤x≤x m , h′ (n,m) =h′ (n,m) (x,z n )[L] represents the position z in the aquitard (n,m) n The water level at and They represent the source and sink terms through the interface between the aquifer and the aquitard respectively:

4. The analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system according to claim 3 is characterized in that: The process of converting the steady seepage control equation into a system of ordinary differential equations in matrix and vector form includes the following steps: In the vertical direction, the hydraulic head at the interface between the aquifer and the aquitard is continuous, and the expressions are: In the horizontal direction, at the interface between the aquifers (n,m) and (n,m+1), the hydraulic head and flux are continuous, and the expressions are: in, and Respectively represent x m the left and right sides; The water level at the top boundary of the aquitard (1, m) is recorded as At the left and right boundaries of the system, considering the constant water level boundary conditions, we have: in, and Respectively represent the constant water level of the left and right boundaries, n = 1, 2, ..., N; The boundary condition for vertical one-dimensional seepage in a weak permeable layer (n, m) is: h' (n,m) (x,0)=h (n,m) (x)(7a); h′ (n,m) (x,b n )=h (n-1,m) (x)(7b); The general solution of equation (2) is: Substituting formula (8) into formulas (3a) and (3b), we obtain: By substituting formula (9) into formula (1), the coupled ordinary differential system of the stable water level of the confined aquifer is obtained, which is expressed in matrix-vector form as follows: in, It contains h (n,m) 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:

5. The analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system according to claim 4 is characterized in that: The process of solving a system of ordinary differential equations by matrix eigenvalue analysis consists of the following steps: make Then formula (10) becomes The particular solution of equation (12) is: The matrix eigenvalue analysis method is used to derive the general solution of the homogeneous part of equation (12); specifically, the following eigenvalue problem is solved to obtain the matrix ξ m The eigenvalue of and the eigenvector Where I is the identity matrix; the eigenvector The eigenvector matrix V is composed as follows m : The overall solution of equation (12) is: in, (2N×1) is the constant coefficient column vector to be determined, U m (x) is a (N×2N) matrix whose non-zero elements and eigenvalues Related, expressed as:

6. The analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system according to claim 4 is characterized in that: The process of constructing the recursive relation for the constant coefficient vector in the general solution includes the following steps: According to the continuity conditions represented by formulas (5a) and (5b), the constant coefficient column vector of adjacent inhomogeneous parameter partitions is obtained: The relationship is: in, Oh m =Φ m+1 (x m ) -1 F m (x m (19a); in, is a (N×2N) matrix, whose non-zero elements are represented as: Using the recursive rule, we can further obtain the following equation (18): in, 7. The analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system according to claim 1 is characterized in that: The process of extending the recursive relationship to the lateral boundary of the system and combining the lateral boundary conditions to obtain the global analytical solution includes the following steps: The recursive relationship is extended to the lateral boundary of the system, and the constant coefficient vector of the lateral boundary partition is obtained in combination with the given boundary conditions; the constant coefficient vector obtained from the lateral boundary conditions is substituted into the recursive relationship to obtain the analytical solution constant coefficient vector of each horizontal partition, and the analytical solution for seepage calculation of multi-level heterogeneous leaking aquifer system is constructed.

8. The analytical method for seepage calculation of a multi-level heterogeneous leaking aquifer system according to claim 7, characterized in that: The process of obtaining a global analytical solution 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 constant water level boundary conditions expressed by formulas (6a) and (6b), we obtain: in, and is a column vector composed of the constant water levels on the left and right boundaries, and is a (N×2N) matrix, whose non-zero elements are represented as: Combining formula (22), formula (23a) and formula (23b), we get the following about the constant coefficient column vector The expression is: What you want Substituting into formula (20) we can get the constant coefficient column vector corresponding to each parameter partition: Then Substitute into formula (16) to obtain the aquifer water level The closed solution of .

Images