Heterogeneous aquifer groundwater tidal response characteristic calculation method based on analytic solution

Abstract

The invention discloses a heterogeneous aquifer groundwater tidal response characteristic calculation method based on an analytic solution, and the method comprises the steps: constructing a multi-aquifer tidal response mathematical model, and comprehensively considering the complex characteristics of a multi-layer structure and horizontal heterogeneity; and deducing a hydrodynamic response analytical solution of the multi-layer heterogeneous aquifer system under the tidal action by using a matrix eigenvalue analysis method. According to the method, synchronous and accurate calculation of groundwater tidal response water level amplitude and phase difference can be achieved under the multi-layer heterogeneous condition, and the method is suitable for a complex system with any number of aquifers and parameter partitions. According to the method, the limitation of a traditional analysis method in multi-layer heterogeneous aquifer tidal response calculation can be efficiently solved, and the calculation precision and efficiency are remarkably improved. The method provides a reliable technical tool and theoretical support for underground water dynamic research, resource evaluation and engineering underground water control in coastal and offshore areas.

Description

Technical Field

[0001] The invention relates to the field of hydrogeology, and in particular to a method for calculating groundwater tidal response characteristics of a heterogeneous aquifer based on an analytical solution. Background Art

[0002] Groundwater systems in coastal or offshore alluvial plains usually have typical multi-layer structural characteristics. The deep cover layers in these areas are often composed of multiple aquifers and adjacent weak permeable layers stacked alternately, forming a complex multi-layer aquifer structure. Due to the complexity and variability of the sedimentary environment, these strata usually show obvious heterogeneity in the horizontal extension direction. This horizontal heterogeneity makes the groundwater flow characteristics more complicated. In addition, affected by the action of ocean tides, aquifers in coastal or offshore areas often show significant tidal response characteristics. The propagation process of tidal head fluctuations in groundwater systems is affected by multiple factors such as multi-layer structure and heterogeneity, and is therefore highly complex. The accurate calculation of this response characteristic is of great significance for understanding the dynamic behavior of groundwater systems, assessing water resource security, and conducting disaster prevention and mitigation research. However, there is currently a lack of effective calculation methods based on rigorous analytical solutions for the tidal response characteristics of multi-layer heterogeneous aquifer systems. Most of the existing analytical methods are limited to simplified models under single-layer heterogeneity or multi-layer homogeneity conditions, and cannot accurately reflect the tidal response dynamics of multi-layer heterogeneous aquifer systems. In response to this technical gap, the present invention proposes a method for calculating the tidal response characteristics of groundwater in heterogeneous aquifers based on analytical solutions. This method takes advantage of the mathematical advantages of analytical solutions, comprehensively considers the complex characteristics of multi-layer structures and heterogeneity, and realizes the rapid and accurate calculation of the tidal response characteristics of multi-layer heterogeneous aquifer systems. Compared with traditional methods, the present invention not only improves the calculation efficiency, but also significantly improves the accuracy of the calculation results, providing new technical means and theoretical support for research and application in related fields. Summary of the invention

[0003] The purpose of the present invention is to provide a method for calculating the tidal response characteristics of groundwater in a heterogeneous aquifer based on an analytical solution. The method systematically considers the complexity of the multi-layer structure and its heterogeneous characteristics by constructing a mathematical model of the tidal response of a multi-layer aquifer. The analytical solution of the hydrodynamic response of a multi-layer heterogeneous aquifer system under tidal action is derived using a matrix eigenvalue analysis method. The present invention can achieve synchronous and accurate calculation of the tidal response water level amplitude and phase difference of each aquifer under multi-layer heterogeneous conditions. Compared with traditional analytical methods, the method of the present invention significantly improves the ability to describe the tidal response process of groundwater under complex heterogeneous conditions. The method is suitable for the study of groundwater dynamics in coastal and offshore areas, and provides solid technical support and theoretical basis for hydrogeological surveys, resource assessments, and engineering groundwater control of aquifer systems.

[0004] To achieve the above functions, the present invention designs a method for calculating the tidal response characteristics of groundwater in a heterogeneous aquifer based on an analytical solution. For a multi-layer heterogeneous aquifer, the following steps S1 to S5 are performed to complete the calculation of the tidal response characteristics of groundwater in a multi-layer heterogeneous aquifer:

[0005] Step S1: dividing the multi-layer heterogeneous aquifer into multiple aquifers and weak permeable layers, and dividing the multi-layer heterogeneous aquifer into multiple heterogeneous parameter partitions in the horizontal extension direction, constructing the control equation of the unstable movement of groundwater in the multi-layer heterogeneous aquifer, and converting it into a group of ordinary differential equations in matrix and vector forms;

[0006] Step S2: using matrix eigenvalue analysis method to obtain analytical expressions of ordinary differential equations;

[0007] Step S3: Based on the continuity conditions of the water level and flux at the interface of the heterogeneous parameter partition, a recursive relation of the constant coefficient vector in the analytical expression is constructed;

[0008] Step S4: Extend the recursive relationship to the lateral tidal water level boundary of the multi-layer heterogeneous aquifer, combine the lateral boundary conditions, first obtain the constant coefficient vector of the heterogeneous parameter partition of the lateral boundary, then derive the constant coefficient vector of any heterogeneous parameter partition, and finally construct a global analytical solution;

[0009] Step S5: Write a calculation program for the global analytical solution and use it to calculate the groundwater tidal response problem in an actual multi-layer heterogeneous aquifer.

[0010] Beneficial effects: Compared with the prior art, the advantages of the present invention include:

[0011] 1. Compared with the existing analytical methods, the present invention can flexibly deal with groundwater flow systems composed of any number of aquifers and supports the setting of any parameter partitions in the horizontal extension direction. Whether it is complex heterogeneous conditions of the formation or multi-layer interactive overflow scenarios, this method can meet the analytical calculation requirements and significantly expand the scope of application of traditional methods.

[0012] 2. The present invention is highly adaptable and can flexibly adapt to a variety of hydrogeological conditions through parameter settings, including specific situations such as homogeneous aquifers, completely impermeable conditions of weak permeable layers, and single-layer leaking aquifers. No matter how complex the heterogeneous conditions are, it can provide accurate tidal response characteristic calculation results. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 is a flow chart of a method for calculating groundwater tidal response characteristics of a heterogeneous aquifer based on an analytical solution provided in an embodiment of the present invention;

[0014] Figure 2is a conceptual diagram of groundwater tidal response in a multi-layer heterogeneous leaking aquifer provided according to an embodiment of the present invention;

[0015] Figure 3 is a spatial distribution diagram of water level amplitude in response to tide in a confined aquifer provided in an embodiment of the present invention;

[0016] Figure 4 It is a spatial distribution diagram of water level phase difference in tidal response of a confined aquifer provided according to an embodiment of the present invention. DETAILED DESCRIPTION

[0017] The present invention will be further described below in conjunction with the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and cannot be used to limit the protection scope of the present invention.

[0018] The method for calculating the tidal response characteristics of groundwater in a heterogeneous aquifer based on an analytical solution provided in an embodiment of the present invention is as follows: Figure 1 For a multi-layer heterogeneous aquifer, the following steps S1 to S5 are performed to complete the calculation of the corresponding characteristics of groundwater tides in the multi-layer heterogeneous aquifer:

[0019] Step S1: dividing the multi-layer heterogeneous aquifer into multiple aquifers and weak permeable layers, and dividing the multi-layer heterogeneous aquifer into multiple heterogeneous parameter partitions in the horizontal extension direction, constructing the control equation of the unstable movement of groundwater in the multi-layer heterogeneous aquifer, and converting it into a group of ordinary differential equations in matrix and vector forms;

[0020] The number of aquifers and aquitards in the multi-layer heterogeneous aquifer is unlimited, the number of heterogeneous parameter partitions in the horizontal extension direction is also unlimited, and the matrix and vector form of ordinary differential equations takes the complex form of the aquifer tidal response water level amplitude as the dependent variable.

[0021] The specific method of step S1 is as follows:

[0022] A multi-layer heterogeneous aquifer consisting of an arbitrary number of aquifers and an arbitrary number of weak permeable layers is considered, and the multi-layer heterogeneous aquifer can be divided into multiple parameter partitions on a limited horizontal extension length, and then the control equations of the unstable movement of groundwater in the multi-layer heterogeneous aquifer are constructed, including the horizontal movement in the confined aquifer and the vertical movement in the weak permeable layer.

[0023] The multi-layer heterogeneous aquifer is divided into N aquifers and N weak permeable layers. The aquifers are numbered from top to bottom, with the top aquifer numbered n=1 and the bottom aquifer numbered n=N; above each aquifer is a weak permeable layer with the same number; refer to Figure 2, in the horizontal direction, the multi-layer heterogeneous aquifer is divided into M heterogeneous parameter partitions, and the horizontal heterogeneous parameter partitions are numbered from left to right, from m = 1 to m = M; in each horizontal heterogeneous parameter partition, the hydrogeological parameters of the aquifer and the weak permeability layer are considered to be homogeneous and isotropic, and have uniform thickness; the nth aquifer located in the heterogeneous parameter partition m is recorded as aquifer (n, m); the nth weak permeability layer located in the heterogeneous parameter partition m is recorded as weak permeability layer (n, m); the hydraulic conductivity of the aquifer (n, m) is recorded as T (n,m) [L 2 T -1 ], the water storage coefficient is S (n,m) [-], the tidal load efficiency factor is γ (n,m) [-]; the permeability coefficient of the aquitard (n, m) is K′ (n,m) [LT -1 ], the water storage rate is S s ' (n,m) [L -1 ], thickness is b n [L], the tidal load efficiency factor is γ′ (n,m) [-].

[0024] Establish the XOZ coordinate system, refer to Figure 2 ,by Figure 2 The right direction is the positive direction of the x-axis, and the x-coordinate of the interface between the heterogeneous parameter partitions m and m+1 is marked as x m , for each aquitard, define a local vertical coordinate z n (0≤z n ≤ b n ), its positive direction is upward; coordinate z 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 the horizontal Darcy flow in aquifer (n, m) and the vertical Darcy flow in aquitard (n, m) are expressed as:

[0025]

[0026] Among them, h (n,m) =h (n,m) (x,t)[L],(x m-1 ≤x≤x m ) represents the water level at the horizontal position x in the aquifer (n, m) at time t, 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 term through the interface between the aquifer and the aquitard, that is, the overflow per unit area, where 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 The expression is as follows:

[0027]

[0028] Taking the dynamic water level of the confined aquifer as the boundary condition of the water level distribution of the weak permeable layer, the analytical expression of the water level of the weak permeable layer is derived. On this basis, the expression of the overflow of the weak permeable layer represented by the water level of the confined aquifer is further derived;

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

[0030]

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

[0032]

[0033]

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

[0035] The periodic water level at the top boundary of the aquitard (1, m) is denoted as h 0,m (t), and express it in plural form as:

[0036]

[0037] in, Including the amplitude and phase information of the water level, is a unit imaginary number, w is the frequency; at the left and right boundaries of the multi-layer heterogeneous aquifer, the fluctuating water level boundary conditions are considered and a hybrid form is used to comprehensively express it:

[0038]

[0039] in, represents the amplitude of the periodic water level at the left and right boundaries;

[0040] Based on the separation of variables method, the water level control equations of each confined aquifer are established with the water level amplitude in complex form as the dependent variable, and converted into a group of ordinary differential equations expressed in the form of matrix and vector operations to simplify complex calculations.

[0041] The analytical solution of the hydrodynamic response of the aquifer (n, m) and the aquitard (n, m) is expressed in the complex domain as follows:

[0042]

[0043] Substituting equations (6) and (8b) into equation (2), and rearranging, we obtain the following equation:

[0044]

[0045] in,

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

[0047]

[0048] Therefore, the analytical expression of equation (9) is:

[0049]

[0050] Substituting equation (11) into equations (3a) and (3b), the unit area overflow discharge from the aquifer (n, m) to the aquitard (n, m) and (n+1, m) is:

[0051]

[0052] in:

[0053] f (n,m) =K′ (n,m) σ (n,m) csch(σ (n,m) b n )(12b)

[0054] g (n,m) =K′ (n,m) σ (n,m) coth(σ (n,m) b n )(12c)

[0055] By substituting equations (12a) and (8b) into equation (1), we obtain a coupled ordinary differential equation (ODE) system for the dynamic components of the water level fluctuation in the confined aquifer, which is expressed in matrix-vector form as follows:

[0056]

[0057] Among them, T m is a (N×N) diagonal matrix whose diagonal elements are T (n,m) , is included A column vector of elements, All elements are Column vector, S m It is a (N×N) diagonal matrix, and the diagonal elements are S (n,m) , G m is a (N×N) diagonal matrix with the following diagonal elements:

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

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

[0060] G m (N,N)=(g (N,m) -f (N,m) )γ′ (N,m) (14c)

[0061] γ ′ (n,m) represents the tidal load efficiency coefficient of the aquitard (n,m), F m is a tridiagonal matrix with the following diagonal elements:

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

[0063] F m (1,1:2)=[g (1,m) +g (2,m) ,-f (2,m) ](15b)

[0064] F m (N,N-1:N)=[-f (N,m) ,g (N,m)](15c)

[0065] Among them, F m (N,N-1:N) represents the matrix F m The elements in row N, column N-1:N.

[0066] Step S2: using matrix eigenvalue analysis method to obtain analytical expressions of ordinary differential equations;

[0067] The specific solution process is as follows:

[0068] make Equation (13) can be further expressed as:

[0069]

[0070] Among them, γ m is a (N×N) diagonal matrix with diagonal elements γ (n,m) , γ(n,m) represents the tidal load efficiency coefficient of the aquifer (n,m). It can be found that is a constant, so the particular solution of equation (16) is:

[0071]

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

[0073]

[0074] in, and Respectively represent matrices The eigenvalues ​​and eigenvectors of , I is the unit matrix. Let N eigenvalues ​​and their corresponding eigenvectors be The feature vector The eigenvector matrix is ​​formed as follows

[0075]

[0076] According to this equation, the overall solution of equation (16) is:

[0077]

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

[0079]

[0080] in, Representation Matrix The elements in row n and columns 2n-1:2n.

[0081] Step S3: Based on the continuity conditions of the water level and flux at the interface of the heterogeneous parameter partition, a recursive relation of the constant coefficient vector in the analytical expression is constructed to ensure the calculation consistency between different partitions;

[0082] The specific method of step S3 is as follows:

[0083] According to the water level and flow continuity conditions represented by equations (5a) and (5b), the constant coefficient column vector of adjacent heterogeneous parameter partitions is obtained: The relationship is:

[0084]

[0085] in:

[0086]

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

[0088]

[0089] Using the recursive rule, we can further obtain from equation (22):

[0090]

[0091] in:

[0092]

[0093] Step S4: Extend the recursive relationship to the lateral tidal water level boundary of the multi-layer heterogeneous aquifer, combine the lateral boundary conditions, first obtain the constant coefficient vector of the heterogeneous parameter partition of the lateral boundary, then derive the constant coefficient vector of any heterogeneous parameter partition, and finally construct a global analytical solution;

[0094] The specific method of step S4 is as follows:

[0095] The recursive relationship is extended to the left and right lateral boundaries of the heterogeneous aquifer, and the relationship between the heterogeneous parameter partition 1 and the constant coefficient column vector corresponding to M is obtained as follows:

[0096]

[0097] According to the lateral boundary conditions expressed by equations (7a) and (7b), the following equation is obtained:

[0098]

[0099]

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

[0101]

[0102] Combining equations (26) with (27a) and (27b), we obtain the following equation about the constant coefficient column vector The expression is:

[0103]

[0104] What you want Substituting into equation (24) we can obtain the constant coefficient column vector corresponding to each inhomogeneous parameter partition: Then Substituting into equation (20) we can obtain the aquifer water level The closed-form solution of ; then, the amplitude and phase difference of the fluctuating water level are calculated by the following formula:

[0105]

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

[0107] Step S5: Write a calculation program for the global analytical solution and use it to calculate the actual groundwater tidal response problem of multi-layer heterogeneous aquifers, realize the synchronous and accurate calculation of the tidal response water level amplitude and phase difference of any multi-layer heterogeneous aquifer, and provide an efficient and reliable technical tool for practical applications.

[0108] The following is an application embodiment of the present invention:

[0109] To ensure generality, consider a multi-layer heterogeneous aquifer (N = 3) consisting of three confined aquifers and their adjacent weak permeable layers, where the horizontal extension length of the multi-layer heterogeneous aquifer is 150 m and the thickness of the weak permeable layer is b. n =2m, left boundary coordinate x 0 =50m, right boundary coordinate xM =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). The overall characteristics of heterogeneity are set as follows: the parameters in the three partitions of the second confined aquifer are consistent; in the second partition, the parameters of different weak permeable layers and confined aquifers are also consistent, specifically set as follows: The hydraulic conductivity of the confined aquifer T (2,m) =150m 2 / d, T (n,2) =150m 2 / d, T (1,1) =T (1,3) =T (3,1) =T (3,3) =50m 2 / d, water storage coefficient S (2,m) =10 -3 , S (n,2) =10 -3 , S (1,1) =S (1,3) =S (3,1) =S (3,3) =10 -2 , vertical permeability coefficient of weak permeable layer K ( ′ 2,m) =10 -1 m / d, K ( ′ n,2) =10 -1 m / d, K ( ′ 1,1) =K ( ′ 3,1) =K ( ′ 1,3) =K ( ′ 3,3) =10 -2 m / d (n=1,2,3), water storage rate S s ' (2,m) =5×10 -4 m -1 , S s ' (n,2) =5×10 -4 m -1 , S s ' (1,1) =S s ' (3,1) =S s '(1,3) =S s ' (3,3) =10 -4 m -1 In terms of boundary condition setting, it is assumed that both sides of the multi-layer heterogeneous aquifer are periodic water level boundaries with a period of 0.5d and an amplitude of 1m, and the top water level change is not considered. According to the derived analytical solution, the spatial distribution of the tidal response water level amplitude and phase difference of the three confined aquifers is calculated as follows: Figure 3 , Figure 4 shown.

[0110] The embodiments of the present invention are described in detail above with reference to the accompanying drawings, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge scope of ordinary technicians in this field without departing from the purpose of the present invention.

Claims

1. A method for calculating the tidal response characteristics of groundwater in heterogeneous aquifers based on analytical solutions, characterized in that: For a multi-layer heterogeneous aquifer, the following steps S1 to S5 are performed to complete the calculation of the corresponding characteristics of groundwater tides in the multi-layer heterogeneous aquifer: Step S1: dividing the multi-layer heterogeneous aquifer into multiple aquifers and weak permeable layers, and dividing the multi-layer heterogeneous aquifer into multiple heterogeneous parameter partitions in the horizontal extension direction, constructing the control equation of the unstable movement of groundwater in the multi-layer heterogeneous aquifer, and converting it into a group of ordinary differential equations in matrix and vector forms; Step S2: using matrix eigenvalue analysis method to obtain analytical expressions of ordinary differential equations; Step S3: Based on the continuity conditions of the water level and flux at the interface of the heterogeneous parameter partition, a recursive relation of the constant coefficient vector in the analytical expression is constructed; Step S4: Extend the recursive relationship to the lateral tidal water level boundary of the multi-layer heterogeneous aquifer, combine the lateral boundary conditions, first obtain the constant coefficient vector of the heterogeneous parameter partition of the lateral boundary, then derive the constant coefficient vector of any heterogeneous parameter partition, and finally construct a global analytical solution; Step S5: Write a calculation program for the global analytical solution and use it to calculate the groundwater tidal response problem in an actual multi-layer heterogeneous aquifer.

2. The method for calculating the tidal response characteristics of groundwater in a heterogeneous aquifer based on analytical solution according to claim 1 is characterized in that: In step S1, the number of aquifers and aquitards in the multi-layer heterogeneous aquifer is unlimited, the number of heterogeneous parameter partitions in the horizontal extension direction is also unlimited, and the matrix and vector form ordinary differential equation system takes the complex form of the aquifer tidal response water level amplitude as the dependent variable.

3. The method for calculating the tidal response characteristics of groundwater in a heterogeneous aquifer based on analytical solution according to claim 1 is characterized in that: The specific method of step S1 is as follows: The multi-layer heterogeneous aquifer is divided into N aquifers and N aquitards. The aquifers are numbered from top to bottom, with the top aquifer numbered n=1 and the bottom aquifer numbered n=N. Above each aquifer is an aquitard with the same number. In the horizontal direction, the multi-layer heterogeneous aquifer is divided into M heterogeneous parameter partitions. The horizontal heterogeneous parameter partitions are numbered from left to right, from m=1 to m=M. In each horizontal heterogeneous parameter partition, the hydrogeological parameters of the aquifer and the aquitard are considered to be homogeneous and isotropic, and have uniform thickness. The nth aquifer located in the heterogeneous parameter partition m is recorded as aquifer (n, m); the nth aquitard located in the heterogeneous parameter partition m is recorded as aquitard (n, m); the hydraulic conductivity of the aquifer (n, m) is recorded as T (n,m) [L 2 T -1 ], the water storage coefficient is S (n,m) [-], the tidal load efficiency factor is γ (n,m) [-]; the permeability coefficient of the aquitard (n, m) is K′ (n,m) [LT -1 ], the water storage rate is S s ' (n,m) [L -1 ], thickness is b n [L], the tidal load efficiency factor is γ′ (n,m) [-]; Establish the XOZ coordinate system, and mark the x coordinate of the interface between the heterogeneous parameter partitions m and m+1 as x m , for each aquitard, define a local vertical coordinate z n (0≤z n ≤ b n ), its positive direction is upward; Coordinate z 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 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],(x m-1 ≤x≤x m ) represents the water level at the horizontal position x in the aquifer (n, m) at time t, 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 term through the interface between the aquifer and the aquitard, where 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 The expression is as follows: In the vertical direction, the water level at the interface between the aquifer and the aquitard is continuous, and its expression is: 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: 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 express it in plural form as: in, Including the amplitude and phase information of the water level, is a unit imaginary number, w is the frequency; at the left and right boundaries of the multi-layer heterogeneous aquifer, the fluctuating water level boundary conditions are considered and a hybrid form is used to comprehensively express it: in, represents the amplitude of the periodic water level at the left and right boundaries; Based on the separation of variables method, the analytical solution of the hydrodynamic response of the aquifer (n, m) and the aquitard (n, m) in the complex domain is expressed as: Substituting equations (6) and (8b) into equation (2), and rearranging, we obtain the following equation: in, Considering the water level continuity at the interface between the aquifer and the aquitard, the boundary condition for the vertical one-dimensional flow in the aquitard (n, m) is: Therefore, the analytical expression of equation (9) is: Substituting equation (11) into equations (3a) and (3b), the unit area overflow discharge from the aquifer (n, m) to the aquitard (n, m) and (n+1, m) is: in: f (n,m) =K′ (n,m) s (n,m) csch(s (n,m) b n (12b) g (n,m) =K′ (n,m) σ (n,m) food(σ (n,m) b n (12c) By substituting equations (12a) and (8b) into equation (1), we obtain the coupled ordinary differential system for the dynamic component of the water level fluctuation in the confined aquifer, which is expressed in matrix-vector form as follows: Among them, T m is a (N×N) diagonal matrix whose diagonal elements are T (n,m) , is included A column vector of elements, All elements are Column vector, S m It is a (N×N) diagonal matrix, and the diagonal elements are S (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) (14a) G m (n,n)=(g (n,m) -f (n,m) )γ′ (n,m) +(g (n+1,m) -f (n+1,m) )γ′ (n+1,m) (14b) G m (N,N)=(g (N,m) -f (N,m) )γ′ (N,m) (14c) In the formula, γ ′ (n,m) represents the tidal load efficiency coefficient of the aquitard (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) ](15a) F m (1,1:2)=[g (1,m) +g (2,m) ,-f (2,m) ](15b) F m (N,N-1:N)=[-f (N,m) ,g (N,m) ](15c) Among them, F m (N,N-1:N) represents the matrix F m The elements in row N, column N-1:N.

4. The method for calculating the tidal response characteristics of groundwater in a heterogeneous aquifer based on analytical solution according to claim 3 is characterized in that: The specific solution process in step S2 is as follows: make Equation (13) can be further expressed as: Among them, γ m is a (N×N) diagonal matrix with diagonal elements γ (n,m) , is a constant, so the particular solution of equation (16) is: The matrix eigenvalue analysis method is used to derive the general solution of the homogeneous part of equation (16). First, solve the following eigenvalue problem: in, and Respectively represent matrices The eigenvalues ​​and eigenvectors of , I is the unit matrix; let N eigenvalues ​​and their corresponding eigenvectors be The feature vector The eigenvector matrix is ​​formed as follows According to this equation, the overall solution of equation (16) 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, expressed as: in, Representation Matrix The elements in row n and columns 2n-1:2n.

5. The method for calculating the tidal response characteristics of groundwater in a heterogeneous aquifer based on analytical solution according to claim 4 is characterized in that: The specific method of step S3 is as follows: According to the water level and flow continuity conditions represented by equations (5a) and (5b), the constant coefficient column vector of adjacent heterogeneous parameter partitions is obtained: The relationship is: in: in, is a (N×2N) matrix, whose non-zero elements are represented as: Using the recursive rule, we can further obtain from equation (22): in:

6. The method for calculating the tidal response characteristics of groundwater in a heterogeneous aquifer based on analytical solution according to claim 5 is characterized in that: The specific method of step S4 is as follows: The recursive relationship is extended to the left and right lateral boundaries of the heterogeneous aquifer, and the relationship between the heterogeneous parameter partition 1 and the constant coefficient column vector corresponding to M is obtained as follows: According to the lateral boundary conditions expressed by equations (7a) and (7b), the following equation is obtained: in, and is the column vector composed of the amplitude of the left and right boundary water levels, and is a (N×2N) matrix, whose non-zero elements are represented as: Combining equations (26) with (27a) and (27b), we obtain the following equation about the constant coefficient column vector The expression is: What you want Substituting into equation (24) we can obtain the constant coefficient column vector corresponding to each inhomogeneous parameter partition: Then Substituting into equation (20) we can obtain the aquifer water level The closed-form solution of ; then, the amplitude and phase difference of the fluctuating water level are calculated by the following formula: in, express The model, Respectively The imaginary and real parts of .

Images