#### Abstract

The rapidly increasing demand for the consumption of natural gas has attracted the interests to store natural gas in aquifer reservoir. However, natural gas injected into the aquifer reservoir, which could cause ground surface deformation and mechanical integrity destruction of caprock. Taking the aquifer gas storage of S trap as the research object, according to the geological structure and hydrogeological information, a coupling large-scale hydromechanical model is established to evaluate the damage risk of the gas reservoir in S aquifer. The proposed methodology is based on the development of fluid-solid coupling and application of FEM. The different failure mechanisms of S aquifer gas storage caprock can be evaluated on the basis of the tensile failure criterion and Mohr-Coulomb shear failure criterion. To analyze the change of caprock in gas injection and production process more clearly, a reference model is defined as an ideal calculation condition to discuss the mechanical response, pore pressure variation, and surface deformation characteristics of the caprock during injection and production. On this basis, the second scheme of sensitivity analysis is defined. The pressure injection rate, reservoir parameters, in situ stress, and other factors are considered, respectively, and the influence of different input parameters on mechanical stability and surface deformation of caprock is analyzed. Finally, the mechanical stability is analyzed and combined the above two criteria to predict the upper limit injection pressure of S. Simulation results show that the permeability and in situ stress have a significant influence on ground surface deformation and mechanical integrity of caprock, but Young’s modulus and Poisson’s ratio can be ignored; the upper limit pressure coefficient of S is 1.908.

#### 1. Introduction

Since the end of the 20th century, the global natural gas market has fallen short in terms of supply and demand. For example, the marginal profit of the European natural gas market has fallen to negative values in 2020, while China’s natural gas consumption is about 320 billion cubic meters, about a 4% increase than last year [1]. With the increasing demand for natural gas in China, accelerating the construction of gas storage facilities and improving the capacity of natural gas peak sharing and responding to emergency needs has become an important guarantee to promote social and economic development. Such problems can be solved effectively by the use of underground gas storage, which is the most reliable means of natural gas supply and peak regulation for the development and utilization of this resource [2, 3]. However, the process of the development and injection of gas storage is bound to have geomechanical effects, while the influence of rock mechanical and permeability properties on natural gas storage and carbon dioxide sealing is still unknown [4–6].

A large amount of natural gas is injected into reservoirs continuously, which will damage the mechanical integrity and sealing of the caprock [7]. Therefore, it is necessary to quantitatively characterize the degree of deformation and damage of the caprock during the processes of injection and production and to correctly evaluate its mechanical integrity [8]. At present, most research on the stability of aquifer gas storage caprock has focused on the field of geological sequestration of CO_{2} [9, 10], with little research on the analysis of caprock mechanical stability in aquifer gas storage. Because there is a wealth of theoretical models available for use in numerical simulations [11–13], such studies have been widely used in the analysis of the mechanical stability of caprock. Due to the complexity of the lithology and geological environment of caprock, different rock mass constitutive models describe different effects on the deformation and potential failure of caprock [14–16]. In the numerical simulation of the mechanical integrity of caprock, the following three constitutive models are mainly used: (1) the elastic constitutive model, (2) the plastic constitutive model, and (3) the discontinuous medium model. Among them, the elastic model is widely used to simulate caprock mechanics during CO_{2} storage owing to its ease in carrying out multifield coupling calculations with a high degree of coupling [17, 18]. Based on the elastic constitutive model, a large number of hydromechanical coupling geomechanical models have been developed to analyze the mechanical integrity of caprock during the processes of injection and production [19–23]. Current research shows that the mechanical integrity of caprock will be affected by many factors during these processes [24, 25]. For example, the injection rate, the initial in situ stress, and the injection temperature have the greatest effects on the mechanical integrity and stability of caprock [26–28]. However, there are few reports on the mechanical integrity and stability of caprock under the action of multiperiod alternating injection and production. Hence, it is necessary to establish a reasonable geomechanical model to analyze the mechanical integrity and stability of caprock under such conditions.

This paper is aimed to analyze how fluid injection rate, in situ stress coefficient, upper pressure coefficient, and mechanical properties of caprock affect caprock stability and sealing capacity. Based on the geological structure of the S Aquifer storage, a large-scale 2D axisymmetric hydromechanical finite element model is established firstly, which considered the whole site (reservoir, caprock, underburden, and overburden). Then, the developed hydromechanical coupling model is used to model 5 scenarios of hydromechanical impact of the fluid pressure injection on the effective stress field in the reservoir-caprock and discussed the stress sensitivity of caprock. Finally, the mechanical integrity and sealing failure of S caprock are evaluated according to the Mohr-Coulomb criterion and tensile failure criterion.

#### 2. Caprock Failure Tendency Assessment

This section introduced the theory of hydromechanical coupling and the failure criterion of caprock.

##### 2.1. Effective Stress Principle of Rock

Rock is a kind of porous medium, which is composed of solid skeleton, interconnected pores, and fluids (gas and water) stored in the framework pores. The pore pressure is the pressure borne or transmitted by the fluid in the rock, and the effective stress is the stress transmitted by the contact surface between the rock-solid particles. Many laboratory experiments have shown that the pore pressure has different effects on deformation and failure of the fluid fully or partially saturated porous solid [29, 30]. Both theoretical and experimental studies have shown that failure is controlled by the effective stress. In this paper, the stress direction of tensile stress is defined as negative, and the pore pressure of the compressive stress is positive, which can be defined as follows [31]: where is effective stress, is the total stress, is the Kronecker symbol, is pore pressure of water, is pore pressure of gas, and is a coefficient related to saturation and surface tension. Since the test data is difficult to obtain, it is generally assumed that, obviously, , the Biot’s coefficient, which can be defined as where is the volumetric compression modulus of rock and is the compression modulus of solid particles.

The expression of effective stress can be expressed as where is the average pressure of the two fluids.

For the injection of fluid or gas, many literatures [32, 33] have shown that a general reduction of the effective stress might lead to caprock failure, which created the potential flow paths for the gas to migrate towards the surface, i.e., to a tremendous diminution of caprock integrity. Two main mechanical failure mechanisms might occur during fluid injection.

##### 2.2. Tensile Failure Criterion

Tensile failure occurs when the effective stress is significantly reduced and below the tensile strength of rock. According to the theory of the first strength of materials, the tensile failure criterion of caprock is defined as follows [10]: where is the minimum principal effective stress (MPa) and is the tensile strength (MPa). Most sedimentary rocks have relatively low tensile strength. Wu et al. [34] showed that the rock matrix tensile strength is usually assumed to be null, as a conservative assumption, it is assumed that tensile failure will occur if the minimum effective stress is greater than zero. Therefore, it is relatively easy to predict tensile failure.

##### 2.3. Shear Slip Failure Criterion

The shear stress is the stress component, which acts along the fracture plane (as shown in Figure 1). It is written using the principal stress components as follows: where is shear stress of shear plane, is the effective normal stress, is the cohesion of rock, and is the angle of rock.

**(a)**

**(b)**

As shown in Figure 2, the strength criterion determined by equation (5) can be represented by the failure line, whose slope is and its intercept on the -axis is .

During the gas injection stage, the pore pressure of reservoir-caprock’s rock could increase, which leads to the decrease of effective stress. The Mohr circle of stress will move to the left and gradually approach the strength envelope (as shown in Figure 2). When the stress Mohr circle is tangent to the strength envelope, the shear failure will occur; if the stress state is lower than the failure line, the caprock is stable.

#### 3. Geomechanical Model Setup

This section describes the geomechanical model of gas storage, including geometry, fluid injection schedule, material properties, initial and boundary conditions, and constitutive laws. The geomechanical modeling was performed with a finite element program.

##### 3.1. Geometry

S is one of the traps located in the Bohai Bay Basin where there is rich in water-bearing trap resources and has carried out oil and gas exploration for more than 50 years. It has the location and technical advantages of constructing aquifer gas storage around the sea of Bohai. Since 2002, a number of aquifer trap targets have been preliminarily screened. S is an anticline structure and the closure is about 45 m. There is no fault in the trap (as shown in Figure 3). The direct caprock is gray and purplish-red mudstone, the thickness of single layer is 3-14.5 m, and the cumulative thickness is 106.5 m. The reservoir is mainly composed of sandstone, which buried depth is about 1125-1320 m, the thickness of single layer is 5-28 m, and the cumulative thickness is 71 m. The average porosity and permeability are 34% and 896 mD, respectively. Considering the safety of injection and production, the site selection evaluation and the geomechanical problems could present in the process of gas storage injection and production (as shown in Figure 4), it is necessary to carry out the research on sealing and mechanical integrity evaluation of S aquifer caprock.

##### 3.2. Numerical Model

The method of hydromechanical coupling with single-phase fluid flow in porous media is employed. Note that this paper mainly focuses on the mechanical effect caused by the change of total pore pressure, rather than the distribution of water and gas phase saturation [36]. According to the basic geological conditions of S, a 3D conceptual model of aquifer gas storage is established (as shown in Figure 5), and the vertical injection well is simulated by the axisymmetric method (as shown in Figure 6). The use of the axisymmetric model enables us to consider a very large lateral extension, allowing an accurate representation of hydraulic boundary conditions. This model extends vertically from the ground surface to a depth of 2690 m, and the horizontally is around 12 km to simulate laterally infinite conditions. It includes four typical horizon layers, i.e., overburden, caprock, aquifer reservoir, and underburden. The storage reservoir has a thickness of 30 m, and the thickness of the caprock is 36 m. The top of the overburden is placed at the depth of 1124 m. The injection wellbore at the symmetry axis of the reservoir (as shown in Figure 5). Eight-node hexahedral displacement and pore pressure elements are used to mesh the model. A structured technique mesh is used, which near the injection well with a minimum mesh cell of 10 cm in the reservoir and caprock, and becomes coarser further away (as shown in Figure 7). At the bottom and lateral boundaries, no flow and the normal displacements are fixed.

##### 3.3. Rock Mechanical Properties

The rock mechanical properties are obtained according to the available well logging data of the S as well as laboratory experiments. Table 1 lists the rock mechanical properties and fluid-flow properties. This paper mainly discusses the mechanical integrity and ground surface deformation of caprock. Therefore, the reservoirs, overburden, and underburden are regarded as linear elastic materials, which ignored the cohesion and friction.

##### 3.4. Initial Conditions

As shown in Figure 6, the model is cut into 8 equal parts, which is used to simulate the influence of different reservoir and caprock permeability areas on the seal property of caprock in the process of injection and production. The model is a semiopen boundary, with a closure of 45 m and a reservoir depth of about 1200 m. The parameters of the injection area are consistent with the reservoir. The vertical well is injected along the reservoir depth range line, and the initial formation pressure of the reservoir is 11.30 MPa. The initial gas injection pressure is about 21.26 MPa. The coefficient of lateral pressure is 1.0.

##### 3.5. Boundary Conditions

As the lateral boundaries of the geomechanical domain were set far from the injection point, the normal displacements applied to the lateral boundaries were fixed to zero. The vertical displacement at the bottom of the geomechanical model was also fixed to zero (as shown in Figure 6). No flow is considered at the top and bottom boundaries. The right boundary is the hydrostatic boundary. The excess pore pressure is used to calculate the in situ stress in order to reach a mechanical equilibrium between the initial stress fields and the applied boundary conditions of the model.

#### 4. Application to S Aquifer Case

This section analyzed two numerical simulation scenarios: the first is defined as the ideal calculation condition, which is taken as the reference scenarios to discuss the mechanical response, pore pressure change, and ground surface deformation characteristics of the caprock during gas injection and production; the second is the sensitivity analysis scheme, which discusses the influence of different input parameters on the stability of caprock and ground surface deformation.

##### 4.1. Modeling Scenarios

According to the geological information and the feasibility research scheme of S, this paper mainly discussed the following 5 scenarios (As shown from Tables 2–7).

##### 4.2. Reference Scenario No.1

###### 4.2.1. Calculation Results of A0

Before gas injection and production, the distribution of the initial in situ stress field should be simulated to provide initial conditions. Based on the equilibrium results of the initial stress field, the vertical stress, horizontal stress, vertical displacement, pore pressure nephogram, and half-year injection nephogram are taken for comparative analysis, which is shown in the following nephograms.

After gas injection into aquifer a half year, both horizontal and vertical stresses are almost unchanged (as shown in Figures 8 and 9), which indicated that injection rate of pressure has little effect on caprock stress. With the injection pressure accumulating, the pore pressure of the reservoir increases obviously (as shown in Figure 10), which resulted in the decrease of effective stress of the reservoir gradually. The increase of formation pressure caused the expansion effect of the reservoir. Furthermore, it leads to the tendency that caprock may occur ground surface uplifting deformation of caprock. After half a year of injection, the maximum uplift displacement is about 6 cm (as shown in Figure 11). The maximum horizontal displacement is about 4 cm (as shown in Figure 12), which occurs in the caprock of the injection area.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

###### 4.2.2. Results of a Half-Year Injection

After a half-year injection, three horizontal and one vertical survey lines were selected. The horizontal survey lines are taken from the reservoir caprock interface, the middle part of the reservoir, and the upper part of the caprock. The vertical survey lines are taken from the whole model (), respectively. The excess pressure is defined as the difference between final and initial pore pressure. Figure 13 gives the excess pore pressure distribution curve along the radial direction. It can be seen that the top of the pore pressure of caprock is almost unaffected by the injection pressure, while the pore pressure distribution along the survey line of the reservoir caprock interface is completely consistent with the measuring line in the middle of the reservoir. The maximum excess pore pressure is mainly distributed in the injection near zone and, then, gradually decreases and tends to be stable. When it is close to the outer boundary, the excess pore pressure decreases rapidly, and the effect of overpore pressure is obvious, and the response range is about 3 km, which is slightly larger than the closed boundary.

Figure 14 gives the excess pore pressure distribution curve along the vertical direction. It can be seen that the pore pressure of reservoirs increases during the injection, the pore pressure within the range of 50 m of the bottom of caprock has a certain effect, while the effect at the top of the caprock is not significant. The pore pressure in the middle of the caprock tends to decrease due to the low permeability of the caprock and the deformation of the caprock. The pore pressure also increased to a certain extent (100 m) at the bottom top.

Gas injection and production have a significant influence on the formation of pore pressure, and the excess pore pressure will direct influence on the formation effective stress state and the caprock integrity. Figures 15 and 16 give the distribution curve of the formation effective stress (vertical and horizontal) along the radial and vertical direction. The gas injection made the formation pressure increased, which resulted in the effective stress influence of formation is obviously reduced. The maximum stress reduction value occurs in the reservoir injection near the zone (as shown in Figure 15). This working condition is a special case with a lateral pressure coefficient of 1. Before gas injection, the initial horizontal stress and vertical stress are the same. However, the magnitude of the vertical stress reduction is greater than that of the horizontal stress after a half year, because there is no constraint effect on the ground surface and the lateral deformation receives constraint effect. After the injection of half a year, the pore pressure distribution curve is basically consistent with the distribution trend of formation effective stress, but the formation stress is reduced and the pore pressure is increased. Due to the influence of gas injection, in the caprock-reservoir-bottom system, the effective stress change amplitude of caprock and bottom layer is obviously smaller than that of the reservoir (as shown in Figure 16). The caprock mainly appears at the bottom, while the bottom layer appears at the top. The maximum stress change amplitude appears on the formation interface, where in the reservoir-caprock interface and the middle reservoir interface.

**(a)**

**(b)**

**(a)**

**(b)**

Figure 17 gives the ground surface deformation along the injection zone. The maximum ground surface deformation is slightly smaller than the maximum caprock deformation, which is about 4 cm, while the maximum horizontal deformation is about 2.5 cm. The maximum vertical deformation occurs at the injection zone, while the maximum horizontal deformation occurs at the water boundary. After 1 year (i.e., after the gas completes an injection-production cycle), the ground surface deformation cannot completely recover to the initial state. The maximum vertical ground surface deformation after one injection-production cycle is about 1 cm (as shown in Figure 17(a)), and the horizontal deformation is also about 1 cm (as shown in Figure 17(b)). According to the theoretical analysis of pore elasticity, it is mainly because of the existence of a certain fluid excess pressure in the formation, which caused the complete dissipation time of excess pressure that is relatively large compared with the production time. In other words, the low permeability caprock and bottom layer take a relatively long time to dissipate, which indicate that the key factor of the ground surface deformation is the excess pressure in the formation.

**(a)**

**(b)**

##### 4.3. Critical Scenario No. 2

Focusing on the ground surface deformation and stress failure of caprock, we selected three measuring points and ground surface survey line to analyze the sensitivity of each influencing factor. The input parameters included caprock elastic modulus, Poisson’s ratio, permeability coefficient, permeability coefficient of caprock, and in situ stress coefficient (in Table 8).

###### 4.3.1. Ground Surface Deformation

A “one factor at a time” sensitivity analysis is conducted to evaluate the weight of each parameter on the ground surface deformation and stress failure of caprock. It consists of varying each input parameter separately and to measure its effect on the output. Each input parameter is alternatively assigned its lower and upper values is shown in Table 8, whereas the other input parameters are fixed to their values as defined according to the basic scenario A0. It can be seen from Figure 18 that the maximum horizontal displacement and maximum vertical displacement of the ground surface slightly increase with the increase of the peak time of the injected fluid pressure difference or the upper limit pressure coefficient. But the increasing extent is different, which indicates that the influence degree of upper pressure coefficient is lower than that of the pressure injection rate.

**(a)**

**(b)**

Figure 19(a) shows the maximum ground surface displacement decreases with the increase of permeability coefficient of caprock, which proved that the permeability of caprock is the key factor about sealing capacity. However, the influence of other caprock parameters is slight (as shown in Figures 19(b) and 19(c)).

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

Figure 19(d) gives the maximum ground surface deformation with different boundary water openness. However, the reservoir is an important position of gas storage, and the influence of the opening degree of boundary permeability changes on the surface is relatively small. Figure 19(e) shows that the influence of lateral pressure coefficient on surface deformation under fluid injection can be negligible.

###### 4.3.2. Failure Sensitivity Analysis of Caprock

Through the analysis of Section 3.2, the bottom of the caprock and the reservoir interface are the areas with great stress changes. The sudden changes of the lithology of reservoir and caprock is easy to be damaged. Three measuring points at the bottom of caprock are selected as the research objects to study the influence of injection rate, upper limit pressure coefficient, caprock permeability, caprock elastic modulus, Poisson’s ratio, reservoir boundary, and in situ stress coefficient on the mechanical integrity of caprock. The calculation results are shown in the following figures. Where in these figures, is the critical internal friction angle of caprock, and is the maximum effective principal stress.

The decrease of internal friction angle will lead to the failure envelope gradually, which indicates that the internal friction angle is the key factor to judge whether the seal failure occurs. Figure 20 shows the relationship between fluid injection rate and the maximum effective principal stress and internal friction angle. With the increase of fluid pressure injection rate, the maximum effective principal stress and friction angle in caprock changes obviously. But it can be seen from Figures 21 and 22 that the maximum effective principal stress and internal friction angle changes greatly when upper limit pressure coefficient and in situ stress coefficients changed, which showed that the influence degree of upper limit pressure coefficient and in situ stress coefficients are more influential than fluid pressure injection rate.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

As the key factor of gas storage safety construction, the influence of in situ stress and upper limit pressure coefficient should not only be considered. Therefore, it is necessary to consider the mechanical response of the mechanical properties of the caprock to the injection and production process. Figures 23–25 give the relationship between caprock mechanical properties and the maximum effective principal stress and friction angle, which proved the mechanical properties of the caprock have a little influence on the mechanical integrity of the caprock, even can be ignored. Figure 26 shows that the influence of reservoir boundary on mechanical integrity of caprock also can be ignored.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

#### 5. Discussion

This section discussed the caprock mechanical stability and upper limit injection pressure of S through basic scenario A0.

##### 5.1. Mechanical Stability of S Caprock

Due to there are too many scenarios in this paper, this section only selected A0 for caprock stability analysis. The selected measuring points 1, 2, and 3 at the bottom of the caprock are analyzed. Figure 27 shows the stress damage during a period of injection and production, which is analyzed in a Mohr diagram. The dotted line is the outermost tangent line of the Mohr circle at each measuring point when the cohesion is 0, and the straight line is the critical internal friction angle to judge whether the caprock is destroyed. For the increase of injection rate, the formation appears excess pore pressure, which leads to the decrease of effective stress of caprock. It can be seen from Figure 27 that with the change of injection rate, the effective principal stress gradually decreases, resulting in the Mohr circle of stress at the bottom of caprock to move to the left. When the injection pressure reaches the peak times, the critical internal friction angle of the caprock gradually decreases and far away from the failure envelope, which indicates that the caprock does not undergo shear failure. In addition, there is no tensile stress produced in the caprock of this paper, which indicates that there is no tensile failure. After a cycle of injection and production, the pore pressure of the reservoir decreases gradually, and the minimum principal pressure stress of the caprock increases, indicating that the formation is gradually restored.

**(a)**

**(b)**

**(c)**

##### 5.2. Prediction of Injection Upper Limit Pressure

Based on the A0 condition, through the change of the injection pressure and permeability to predict the upper limit injection pressure of S gas storage caprock, the upper limit injection pressure can be judged by the following formula: where is the pressure injected into reservoir, is the water pressure above the caprock, and is the pressure of reservoir pressurization value.

###### 5.2.1. Critical Internal Friction Angle

According to the Mohr-Coulomb shear failure criterion. When the injection pressure is MPa, the Mohr circle change trend of caprock stress is shown in Figure 28. At this time, the maximum critical internal friction angle is about 28.5°, which indicated that shear failure occurred locally in the caprock. The minimum principal compressive stress of the corresponding caprock is shown in Figure 29. The caprock has always been in the state of compressive stress without tensile failure.

###### 5.2.2. Tensile Failure

According to the tensile failure criterion, the maximum injection pressure is simulated, and the corresponding minimum principal stress is calculated. When maximum injection pressure is MPa, the minimum principal compressive stress is close to 0 (see Figures 30 and 31). The tensile failure occurs in the caprock when the compressive stress approaches 0.

###### 5.2.3. Breakthrough the Pressure

Due to the lack of measured breakthrough pressure data of caprock. According to the former Soviet scholar pilip [37], who studied the model of the relationship between permeability and breakthrough pressure through the test where is the penetration pressure difference, unit: MPa; is the permeability of caprock, unit: mD.

According to the related data of the adjacent area, the caprock permeability is 0.001 mD, and the corresponding breakthrough pressure is 10.26 MPa. The critical injected gas pressure of reservoir is about 21.56, which can be calculated by Eq. (6).

As shown in Table 9, under the constraints of the three evaluation criteria, the caprock first happened gas breakthrough, the second is a shear failure, and the last is a tensile failure. In order to ensure the sealing of caprock in the process of gas injection and production, the safest value should be selected as the critical injection pressure. Therefore, 21.56 MPa is selected as the upper limit injection pressure, and the corresponding upper limit pressure coefficient is about 1.908.

#### 6. Conclusions

In the present work, a coupled large-scale hydromechanical model was developed for modeling the mechanical change of caprock. There are five factors that were considered to describe the influence degree on caprock integrity, which included permeability, Poisson’s ratio, Elastic Modulus, in situ stress, and upper limit pressure coefficient. Then, the maximum injection pressure was predicted by the Mohr-Coulomb, tensile failure criterion, and capillary pressure.

Based on the geological structure, logging data, and conventional triaxial tests, a large-scale FEM model has been established for the engineering application. The results indicate that the minimum horizontal principal stress of caprock is in the range of the compressive stress when the maximum injection pressure is less than the breakthrough pressure of caprock, which indicates that there is no tensile failure occurred. The Mohr circle is below the strength envelope line during the injection production process, and for the caprock with undeveloped fractures, there will be no shear failure bad.

The simulated result shows that the in situ stress coefficient is more sensitive to the mechanical integrity of the caprock, followed by the upper limit pressure coefficient and injection rate coefficient, while the elastic modulus and Poisson’s ratio of the caprock have little influence on the mechanical integrity of the caprock, which can be ignored.

Further studies are necessary to be done for a better understanding of the caprock integrity and the sealing capacity. Especially when there are faults in the target trap, it will be included in the future research work.

#### Conflicts of Interest

There is no conflict of interests exists in this paper.

#### Acknowledgments

The authors gratefully acknowledge the support of the Open Research Fund of State Key Laboratory of Deep Rock and Soil Mechanics and Underground Engineering (Grant No. SKLGDUEK2001), the Natural Science Foundation of Heilongjiang Province (Grant No. LH2020D004), and the National Natural Science Foundation of China (Grant No. 42072166).