A method for simulating time-lag compression interlayer body permeability coefficient based on ground stress change principle

Abstract

The present application relates to a kind of based on the principle of stress change of compression time lag interlayer body variable permeability coefficient simulation method, comprising: the grid space discrete processing of compression time lag interlayer body is carried out, obtains several discrete grid unit;For the discrete grid unit of different position of preestablished, constructs water balance equation;The water balance equation of each discrete grid unit is simultaneously solved to form difference equation group, and the difference equation group is iteratively solved, realizes the compaction / expansion process simulation of each node unit of interlayer body.The present application can reflect the coupling influence of compression time lag interlayer body compaction, expansion and its permeability coefficient, perfect the simulation theory mechanism of compression time lag interlayer body, improve the accuracy and applicability of ground subsidence numerical simulation.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of groundwater system simulation, in particular to a method for simulating variable permeability of compressible time-lag interlayer based on the principle of ground stress change. BACKGROUND

[0002] Land subsidence, also known as ground subsidence or land collapse, is a kind of downward movement caused by the consolidation compression of underground loose strata, which threatens the sustainable development of natural environment and economy and society. When the subsidence phenomenon occurs, the release of water in the aquifer is accompanied by the change of effective stress. The parameters such as water storage property, porosity and permeability coefficient of the clayey soil will change after compression deformation. Due to the high complexity of the compression process, a large number of scholars have conducted extensive research on the mechanism and trend prediction of land subsidence. The numerical simulation method based on physical mechanism is an important means to study and predict the evolution law of land subsidence. This method couples the groundwater seepage model and the soil model to more completely describe the interaction process of groundwater and soil compression.

[0003] In the field of groundwater simulation, MODFLOW is widely used due to its comprehensive functions. The CSUB package in the latest version of the software integrates all the functions of the coarse-grained medium in the aquifer, the non-compressible time-lag interlayer, the compressible time-lag interlayer, and the release of compacted water and the storage of expanded water of the pore water in the aquifer based on the principle of ground stress change to simulate them together. Compared with the previous versions such as IBS1, SUB, SUB-WT, etc., the CSUB package is more perfect in mechanism. During the simulation process, the water storage property of the compressible time-lag interlayer will change with the change of the effective stress due to the compaction, expansion and change of the overburden load, but its permeability is not affected by the change of its porosity, which is inconsistent with the actual situation. In fact, the compaction of the interlayer will cause a significant change in the thickness of the interlayer, resulting in a decrease in the porosity of the compressible medium, which in turn will affect the compaction process and drainage process of the interlayer. Under certain conditions, not considering the change of the permeability coefficient may lead to the porosity of the interlayer being simulated as negative, i.e. the distortion result of "overcompaction". SUMMARY

[0004] The purpose of the present application is to provide a method for simulating variable permeability of compressible time-lag interlayer based on the principle of ground stress change to solve the problems existing in the prior art.

[0005] To achieve the above purpose, the present application provides the following scheme:

[0006] A method for simulating variable permeability of compressible time-lag interlayer based on the principle of ground stress change, comprising:

[0007] The compressed time delay sandwich body is discretely processed in a grid space to obtain a plurality of discrete grid units;

[0008] The water balance equation is constructed for the discrete grid units at different positions;

[0009] The water balance equation of each discrete grid unit is constructed to form a difference equation group, and the difference equation group is iteratively solved to realize the compaction / expansion process simulation of each node unit of the sandwich body.

[0010] Optionally, the discrete grid units at different positions include a starting boundary unit, an end boundary unit and an internal node unit.

[0011] Optionally, the water balance equation of the starting boundary unit is:

[0012]

[0013] wherein, is the current vertical permeability coefficient of the first node unit; is the water head of the groundwater grid unit j at the end of the mth calculation period; is the water head on the first node unit at the end of the mth calculation period; Δz1 is the length of the first node unit, i.e. the starting boundary unit; Δz2 is the length of the second node unit; Δt is the time step; is the selected water storage rate on the first node unit within the period; is the average vertical permeability coefficient between the current first node unit and the second node unit; is the water head on the second node unit at the end of the mth calculation period; z j is the reference elevation corresponding to the calculation of the total ground stress and the effective stress, and the model is taken as the bottom plate of the groundwater grid unit; is the total ground stress at the bottom plate of the grid unit at the end of the mth calculation period; is the pre-consolidation stress on the first node unit at the end of the m-1th calculation period; is the elastic water storage rate on the first node unit within the period; is the effective stress on the first node unit at the end of the m-1th calculation period.

[0014] Optionally, the water balance equation of the end boundary unit is:

[0015]

[0016] Optionally, the water balance equation of the internal node unit is:

[0017]

[0018] Optionally, the vertical permeability coefficient is calculated by using the Kozeny-Carman equation.

[0019] Optionally, the average vertical permeability coefficient between two node units is calculated by using the harmonic mean method.

[0020] Optionally, the difference equation group is:

[0021] [A] m [h] m =[r] m

[0022] wherein [A] m is a NN×NN three-diagonal symmetric matrix; [h] m is a one-dimensional unknown head vector with NN elements; and [r] m is a one-dimensional known vector with NN elements.

[0023] Optionally, the iterative solution of the difference equation group comprises:

[0024] In each iterative solution process, two layers of nesting are performed, wherein the outer layer of nesting is to update the length and the vertical permeability coefficient of the node unit according to the current compaction / expansion amount of the node unit, and the inner layer of nesting is to solve the matrix equation constructed according to the current length and the vertical permeability coefficient of the node unit;

[0025] The iterative solution process continues until the length of the node unit and the head value of the node unit meet the preset solution accuracy requirement.

[0026] The beneficial effects of the present application are:

[0027] The present application provides a compression time lag interlayer body variable permeability coefficient simulation method based on the principle of ground stress change. First, the compression time lag interlayer body is subjected to grid spatial discretization processing to obtain a plurality of discrete grid units. Second, the water balance equation is constructed for the discrete grid units at different positions. Finally, the water balance equations of the discrete grid units are simultaneously solved to form a difference equation group, and the difference equation group is iteratively solved to simulate the compaction / expansion process of each node unit of the interlayer body. When simulating the compression time lag interlayer body, the present application can reflect the coupling effect of the compaction, expansion and permeability coefficient of the compression time lag interlayer body, perfect the simulation theory mechanism of the compression time lag interlayer body, and improve the accuracy and applicability of the numerical simulation of land subsidence. BRIEF DESCRIPTION OF DRAWINGS

[0028] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the following will briefly introduce the drawings needed in the embodiments. Obviously, the drawings described below only constitute some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained from these drawings without creative labor.

[0029] Figure 1 A flow chart of a variable permeability simulation method of a compression time lag interlayer body based on the principle of ground stress change for an embodiment of the present application;

[0030] Figure 2 A one-dimensional space discrete graph of a compression time lag interlayer body for an embodiment of the present application;

[0031] Figure 3 A model setting graph for an embodiment of the present application;

[0032] Figure 4 A comparison graph of MODFLOW-CSUB and the simulated settlement amount of the present application for an embodiment of the present application;

[0033] Figure 5 A comparison graph of MODFLOW-CSUB and the vertical permeability coefficient at the center node of the interlayer body of the present application for an embodiment of the present application;

[0034] Figure 6 A comparison graph of MODFLOW-CSUB and the change of the water head value at the center node of the interlayer body of the present application for an embodiment of the present application. DETAILED DESCRIPTION

[0035] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments only constitute some embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0036] In order to make the above-mentioned purposes, features and advantages of the present application more obvious and easy to understand, the present application will be further described in detail below with reference to the drawings and specific embodiments.

[0037] As shown in Figure 1 , the present embodiment proposes a variable permeability simulation method of a compression time lag interlayer body based on the principle of ground stress change, which comprises:

[0038] S1, performing grid space discrete processing on the compression time lag interlayer body to obtain a plurality of discrete grid units;

[0039] S2, constructing a water balance equation for the discrete grid units at different positions.

[0040] S3, the water balance equation of each discrete grid unit is constructed into a difference equation group, and the difference equation group is iteratively solved to realize the compaction / expansion process simulation of each node unit of the interlayer body.

[0041] Further, step S1 is specifically that the compression time delay interlayer body is discretized into one-dimensional NN grid units (such as shown in FIG. 1) by using a rectangular grid. Figure 2

[0042] Further, in step S2, the discrete grid units at different positions include a starting boundary unit, a terminal boundary unit and an internal node unit. The positions are as shown in FIG. 2. The different positions in the embodiment are not specifically classified because the water exchange equations between the interlayer body and the aquifer at the interlayer body and aquifer junction and between the central node units are different. According to the different water balance equations, the equations are divided into three types of equations at the interlayer body and aquifer junction (upper and lower boundaries) and between the internal node units. Figure 2

[0043] In the embodiment, S2, the water balance equation is constructed for the discrete grid units at different positions, and the Kozeny-Carman equation is introduced to calculate the change of the vertical permeability coefficient. Specifically, the following is included:

[0044] For the compression time delay interlayer body, the compaction water release and the swelling water storage of the interlayer body based on the principle of ground stress change can be expressed by the following one-dimensional water head diffusion equation:

[0045]

[0046] Wherein: Kv' is the vertical permeability coefficient (L / T) of the compression time delay interlayer body; h is the water head (L) in the interlayer body; z is the vertical spatial coordinate (L); S s is the water storage rate (-); t is the time (T); γ w is the water bulk density (M / L2 / T2); σ' is the effective stress (M / L / T2).

[0047] For any mth calculation period, the water entering the unit is positive, and the numerical dispersion of the partial differential equation is carried out according to formula (1). First, the water balance difference equation of the equivalent interlayer body node unit i=1:

[0048]

[0049] Wherein: is the vertical permeability coefficient (L) of the first node unit at present; is the water head (L) of the groundwater grid unit j (assuming that the equivalent interlayer body is located on the groundwater grid unit) at the end of the mth calculation period.​​ is the water head on the 1st node element at the end of the mth calculation period; Δz1 is the length (L) of the 1st node element; Δz2 is the length (L) of the 2nd node element; Δt is the time step (T); is the selected water storage rate (1 / L) on the 1st node element during the period, which is the elastic water storage rate or the inelastic water storage rate according to whether the water head of the node element exceeds the pre-consolidation water head; is the average vertical permeability coefficient (L) between the current 1st node element and the 2nd node element; is the water head on the 2nd node element at the end of the mth calculation period; z j is the corresponding reference elevation when calculating the total ground stress and the effective stress, which is taken as the bottom of the groundwater grid element (L); is the total ground stress (in terms of water column height) at the bottom of the grid element at the end of the mth calculation period (L); is the effective stress (based on the bottom of the grid element, in terms of water column height) on the 1st node element at the end of the mth calculation period (L); is the pre-consolidation stress (based on the bottom of the grid element, in terms of water column height) on the 1st node element at the end of the m-1th calculation period (L); is the elastic water storage rate (1 / L) on the 1st node element during the period; is the effective stress (based on the bottom of the grid element, in terms of water column height) on the 1st node element at the end of the m-1th calculation period (L).

[0050] It is obtained that:

[0051]

[0052] For node NN, similarly, we have:

[0053]

[0054] For the 1st node element 1<i<NN-1, the water balance difference equation is established as:

[0055]

[0056] It is obtained that:

[0057]

[0058] The elastic and inelastic water storage rates of each node element in the above formulas are calculated as:

[0059]

[0060] wherein: C c and C ris the non-dimensional compression and recompression index of the delay interlayer (−); e i is the current void ratio of the i-th node element (−); is the effective stress (in terms of water column height) used in the current calculation of the water storage rate parameter of the i-th node element; is the effective stress (in terms of water column height) on the i-th node element at the end of the m-th calculation period (L).

[0061] The Kozeny-Carman equation is used to calculate the vertical permeability coefficient Kv' during the simulation process, which is:

[0062]

[0063] where K is the permeability coefficient of the sediment (L / T), and θ is the porosity of the sediment (−), so:

[0064]

[0065] where: is the vertical permeability coefficient (L / T) of the i-th node element at the end of the m-1-th calculation period; is the porosity (−) of the i-th node element at the end of the m-th calculation period; is the porosity (−) of the i-th node element at the end of the m-1-th calculation period.

[0066] Referring to equation (9), the permeability coefficient is related to the porosity of the interlayer. When the interlayer is compacted or expanded, its porosity changes, which in turn changes the size of the permeability coefficient. The change in the permeability coefficient will then affect the water release process of the interlayer (as you can also find from the water balance equation), thus producing a mutual feedback coupling effect.

[0067] The average vertical permeability coefficient between two node elements is calculated using the harmonic mean method:

[0068]

[0069] Further, in step S3, the iterative solution of the difference equations includes:

[0070] In each iterative solution process, there are two layers of nesting, where the outer layer of nesting is to update the length and vertical permeability coefficient of the node element according to the current compaction / expansion amount of the node element, and the inner layer of nesting is to solve the matrix equation according to the current length and vertical permeability coefficient of the node element;

[0071] The iterative solution process continues until the length of the node element and the water head value of the node element meet the preset solution accuracy requirements.

[0072] Specifically, in the present embodiment, S3, the equations of each node element are combined to form a difference equation set, and the solution is implemented to realize the simulation of the compaction / expansion process of each node element of the thick interlayer body.

[0073] The equations (3), (4) and (6) of each node element are combined to obtain a symmetric matrix equation:

[0074] [A] m [h] m =[r] m (11)

[0075] Wherein: [A] m is a NN×NN three-diagonal symmetric matrix; [h] m is a one-dimensional unknown head vector of NN elements; [r] m is a one-dimensional known vector of NN elements. The matrix equation can be solved by Gaussian elimination method.

[0076] Equations (3), (4) and (6) represent water balance equations at different positions (boundary or internal node elements). When there are multiple node elements, the equations of each element are combined to obtain the equation set, wherein the A matrix is the coefficient matrix of the node element head h. The specific meaning in the symmetric matrix equation can be seen from the following formula.

[0077] Each element in the matrix equation is:

[0078]

[0079] Since the length Δz i , the vertical permeability coefficient Kv' i , the water storage parameter S sk,i , etc. will change in the simulation, the matrix equation (11) must be solved iteratively. Each iteration process needs two layers of nesting, the outer layer of nesting is to update Δz i and Kv' i according to the current compaction / expansion amount of the node element, and the inner layer of nesting is to construct the matrix equation according to the current Δz i and Kv' i to solve. The iterative solution process will continue until Δz i and the node element head value meet the solution accuracy requirement.

[0080] The present embodiment simulates the process of thick interlayer body drainage caused by the gradual decrease of aquifer head, which is divided into 1 stress period, and the total simulation time is 1000d. Using a time step multiplier equal to 1.05, it is divided into 100 calculation periods of different lengths. The model setting is as follows: Figure 3As shown, the model mesh consists of 1 layer, 1 row, and 3 columns. The top elevation of the model is 0 meters, the bottom elevation is -1000 meters, and the row and column widths are both 1 meter. A delayed interlayer is located in the middle mesh, and the head values ​​of the constant head mesh cells on both sides are 0 meters. The permeability coefficient of the middle mesh is set to 1.0E+6 to ensure that the head value is also 0 meters. Detailed model parameters are shown in Table 1. For comparison, the constant permeability coefficient method of the MODFLOW6-CSUB package was also used to simulate the aquifer in this embodiment.

[0081] Table 1 Model Parameters

[0082] Parameter Value Grid cell initial water head (m) 0 Wet soil specific gravity (dimensionless) 1.7 Saturated soil specific gravity (dimensionless) 2 Coarse material porosity (dimensionless) 0.2 Interlayer elastic water storage (1 / m) 1.0e-5 Interlayer inelastic water storage (1 / m) 1.0e-2 Interlayer porosity (dimensionless) 0.45 Interlayer initial water head (m) 10 Interlayer preconsolidation water head (m) 10 Interlayer thickness (m) 10 Number of discrete elements 19 Initial permeability (m / d) 0.001

[0083] Achieved application effects:

[0084] The purpose of this comparative test case is to examine the difference in simulation results between the method in this embodiment and the constant permeability coefficient method under CSUB conditions for a thick interlayer under significant compaction and near-hydraulic head equilibrium. Figure 4 The figure shows a comparison of the settlement calculation results of the present invention and MODFLOW-CSUB. By the end of the simulation period, the settlement calculated by the two methods was similar, but the result calculated by MODFLOW-CSUB showed that the compaction was faster. This indicates that the variable permeability coefficient method has a significant impact on the compaction process without changing the final compaction amount of the interlayer. Changes in the vertical permeability coefficient will significantly delay the compaction process of the interlayer. Figure 5 , Figure 6 The figures show the changes in vertical permeability coefficient and hydraulic head at the central node of the interlayer. It can be seen that under the variable permeability coefficient method, the vertical permeability coefficient decreases rapidly at the beginning of the simulation due to the higher compaction rate of the interlayer. However, as the compaction process ends, the vertical permeability coefficient remains essentially unchanged. Meanwhile, under the constant permeability coefficient method, the dissipation rate of the hydraulic head at the central node of the interlayer is significantly faster, exhibiting a rapid decrease followed by a relatively stable state. Although the hydraulic head values ​​are quite similar at the end of the simulation (0.0019 m and 0.0469 m for the constant and variable permeability coefficient methods, respectively), the dissipation time is longer under the variable permeability coefficient method. These results indicate that the variable permeability coefficient method can fully reflect the slow compression process caused by changes in permeability coefficient during the simulation without changing the final compaction amount of the interlayer. Therefore, the method is reasonable and consistent with reality.

[0085] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Various modifications and improvements made to the technical solutions of the present invention by those skilled in the art without departing from the spirit of the present invention should fall within the protection scope defined by the claims of the present invention.

Claims

1. A method for simulating the permeability coefficient of a body with a compressive time lag interlayer based on the principle of ground stress change, characterized in that, The method comprises the following steps: carrying out grid spatial discretization processing on the compression time lag sandwich body to obtain a plurality of discrete grid units; constructing water balance equations for the discrete grid units at different positions; solving the difference equation set by iteration to realize simulation of the compaction / expansion process of each node unit of the sandwich body; the discrete grid units at different positions include a starting boundary unit, an end boundary unit and an internal node unit; the water balance equation of the starting boundary unit is: in, This represents the current vertical permeability coefficient of the first node element; For the first Groundwater grid unit at the end of each calculation period water head; For the first At the end of the calculation period, the first Water head at each node unit; The length of the first node element, i.e., the starting boundary element; This is the length of the second node element; For time step; For the first time period The water storage rate to be selected on each node unit; This is the average vertical permeability coefficient between the current first node unit and the second node unit; For the first At the end of the calculation period, the first Water head at each node unit; This serves as the reference elevation for calculating total ground stress and effective stress. For the first The total geostress at the bottom plate of the grid cell at the end of each calculation period; For the first At the end of the calculation period, the first Preconsolidation stress on each node element; For the first time period Elastic water storage rate on each node unit; For the first At the end of the calculation period, the first Effective stress on each node element; the water balance equation of the end boundary unit is: in, Let be the current vertical permeability coefficient of the Nth node element; For the first The water head at the Nth node unit at the end of each calculation period; This is the length of the Nth node unit, i.e., the starting boundary unit; The length of the (NN-1)th node unit; The water storage rate to be selected at the NNth node unit within the time period; This represents the average vertical permeability coefficient between the current (NN-1)th node unit and the NNth node unit; For the first The water head at the (NN-1)th node unit at the end of each calculation period; For the first Preconsolidation stress on the NNth node element at the end of each calculation period; The elastic water storage rate of the NNth node unit within the time period; For the first The effective stress on the NNth node element at the end of each calculation period; the water balance equation of the internal node unit is: wherein, Kz is the average vertical hydraulic conductivity between the current ith node element and the i-1th node element; hi is the water head at the end of the i-1th node element; L is the length of the ith node element, i.e., the starting boundary element; L is the length of the ith node element, i.e., the starting boundary element; L is the length of the i-1th node element; R is the tentative water storage rate in the ith node element during the time period; Kz is the average vertical hydraulic conductivity between the current ith node element and the i+1th node element; hi+1 is the water head at the end of the ith node element; hi+1 is the water head at the end of the ith node element; σp is the pre-consolidation stress at the end of the ith node element; σp is the pre-consolidation stress at the end of the ith node element; R is the elastic water storage rate in the ith node element during the time period; σe is the effective stress at the end of the ith node element; and σe is the effective stress at the end of the ith node element.

2. The method according to claim 1, characterized in that, the vertical permeability coefficient is calculated by using the Kozeny-Carman equation.

3. The method according to claim 2, wherein the method is characterized by, The average vertical permeability coefficient between two node units is calculated by using the harmonic mean method.

4. The method according to claim 1, wherein the method is characterized by, the difference equation set is: in, for A tridiagonal symmetric matrix; for A one-dimensional head vector with elements to be determined; for A one-dimensional known vector with n elements.

5. The method according to claim 1, wherein the method is characterized by, solving the difference equation set by iteration comprises: in each iteration solving process, two layers of nesting are carried out, wherein the outer layer of nesting is to update the length and vertical permeability coefficient of the node unit according to the current compaction / expansion amount of the node unit, and the inner layer of nesting is to solve the matrix equation constructed according to the current length and vertical permeability coefficient of the node unit; the iteration solving process continues until the length of the node unit and the water head value of the node unit meet the preset solving accuracy requirement.

Images