Abstract
Multiseam coalbed methane (CBM) exploitation can not only reduce single-well investment but also increase the length of service of the well and significantly enhance the CBM economic recovery of the entire basin. To compare with and further to guide the actual project of CBM production, this study proposed a conceptual gas-water two-phase separate flow model for single coal seam considering the solubility of gas. This mathematical model was solved analytically by separation of variables and verified through history matching of the production data from the No. 3 seam of 1# test well of Jincheng and then applied to investigate the effect of gas solubility on the gas pressure. Furthermore, based on the coupled two-phase separate flow model of single seam, another two-phase separate flow model for the development of multicoal seam development was established. Similarly, the analytical solution of this model for multicoal seam layers was matched with the in situ data of TS-1 well of Liupanshui coal mine. It is found that the height difference and pressure difference between the two seams play key roles in the multiseam CBM development comprehensively.
1. Introduction
To improve vertical well coalbed methane (CBM) production, a great amount of well drainage control methods have been researched, such as vertical drainage of single seam layer, multiseam fracturing drainage of multicoal seam layers, and single-seam fracturing drainage of multicoal seam layers [1]. Through investigation of the multiseam CBM exploitation in Powder River Basin, the U.S. Department of Energy [2] has found that multiseam CBM exploitation can not only reduce single-well investment but also increase the length of service of the well and significantly enhance the CBM economic recovery of the entire basin [3]. Moreover, it is found that some wells adopting multiseam CBM exploitation have been working for more than 10 years in the earlier well group in the Panzhuang block of the southern Qinshui Basin, with the gas production remaining around 1000 m3/d [4]. Thus, the effective development of multicoal seam layers will be helpful to obtain high yield in vertical CBM well production [4, 5].
Based on multicoal seam development, a variety of factors which can affect mining such as burial depth, coal thickness, gas content, porosity, permeability, structural setting, and damage mechanism of upper exposed producing layers have been studied [6–10]. Besides, in the technical perspective, several kinds of technology such as well test analysis, numerical simulation, and physical simulation have been applied to the development of multicoal seam [1, 11]. For instance, after the comparison between the theory and the field data from coal reservoir in southern Qinshui Basin, China, Li and Li [1] put forward that the main controlling factors influencing the multicoal seam development consist of the height of the critical production gas, the grades of reservoir pressure, the difference of ability of supply for liquid, and the permeability; in the aspect of numerical simulation, the laws of the geological mining factors to the upper and lower pillar’s stability in multicoal seam strip mining were studied by FLAC3D and the function relation between the stress increasing coefficient of upper and lower pillars with the mining depth, mining widths, interlayer spacing, mining thickness, property of interstratified rock, and the spatial relationship was obtained [11].
In order to closely simulate and further to guide the practical production of CBM wells, the modeling of multicoal seam development has become the research emphasis naturally. To select favorable regions for multiseam CBM exploitation, Zhang et al. [3] have conducted a multiobjective fuzzy matter-element model to evaluate the potential of multiseam CBM exploitation in the study area based on the correlation analysis. On the other hand, in general, pore fractures of low-rank coal reservoirs are filled with water, especially in a shallow burial depth of 1000 m [12]; thus, the study of CBM exploitation process is actually the process of two-phase flow, on which many scholars have already conducted large amount of researches [13–15]. Ibrahim and Nasr-El-Din [14] have developed an analytical model for saturated CBM reservoirs where two-phase flow occurs from the start of production, and the model can be used as a prediction or a history matching tool for CBM performance. Clarkson and McGovern [16] developed a model for simulating single-phase flow of water above desorption pressure and two-phase flow of gas and water below desorption pressure, but effects of transient flow were ignored. As a supplement, using a revised DDA model, the model of transition from single-phase to multiphase flow during transient flow for CBM reservoirs has been conducted [13]. Due to the difficulty of obtaining the analytical solution, numerical simulation is always adopted to solve the model founded in the study of two-phase flow; for example, Ziarani et al. [17] have investigated the effect of sorption time on CBM recovery through numerical simulation. However, there is almost no mathematical model focused on the relationship of the fluid pressure, production rate and mining time, mining distance on the development of multiseam layers at two-phase flow stage, not to mention the analytical solution of the fluid pressure and velocity, and so on.
It is generally accepted that CBM mainly exists in coal seam in adsorbed state, but more and more studies reveal that dissolved gas should not be ignored [18, 19]. Based on this conclusion, lots of experiments for methane solubility are conducted [12, 20]. Samples taken from Hailar Basin of Xinjiang, China, were physically simulated of dissolved methane content [12]. The experimental results indicate that the solubility of methane is influenced by temperature, pressure, coal reservoir water salinity, etc. But few coupling seepage models considering the solubility factor are established and discussed in theoretical research actually.
This paper proposes a two-phase separate flow model for the mining of multicoal seam layers. Firstly, considering the solubility of gas, a two-phase separate flow model of single seam is proposed for the multiphysical process. Secondly, this coupled mathematical model is analytically solved by separation of variables and applied to compare with the field data of the No. 3 seam of 1# test well of Jincheng, Shanxi, for the verification of model accuracy. Thirdly, the impacts of gas solubility on the gas pressure are analyzed. Furthermore, based on the coupled two-phase separate flow model of single seam, another coupled mathematical model of two-phase separate flow for the development of multicoal seam layers is established. Similarly, the analytical solutions of this model for multicoal seam layers is matched with the field data of TS-1 well of Liupanshui coal mine, China [10]. Finally, the impact of interlayer difference on the development of multicoal seam layers is demonstrated.
2. Gas-Water Two-Phase Separate Flow Model for Single-Seam CBM Development
Before the establishment of the two-phase flow model, assume that the gas and water two-phase flows are two completely separate fluids with different velocities and properties for the convenience of obtaining the analytical solutions. However, the gas solubility in water is considered in the following models.
2.1. Gas-Water Flow Continuity Equations
In the following fluid continuity equations at the two-phase flow stage, subscripts g and w are used to represent the gas phase and water phase, respectively.
The continuity equation for gas is
The continuity equation for water is where is the coal porosity and is the gas solubility in water. is the density of water and is the density of gas under formation conditions. is the velocity of gas and is the velocity of water. and denote the fluid strength of gas and water, respectively. is the Hamiltonian operator. and are the saturations of gas and water, respectively, and the pores are fully saturated by gas and water:
2.2. Gas-Water Flow Motion Equations
The gas-water coupled seepage equations of motion are [21] where and are the absolute velocities of water and gas, respectively. is the velocity of the coal solid. is the absolute permeability of porous media. and are the relative permeabilities of water and gas, respectively. and are the viscosities of water and gas, respectively. and are the pore pressures of water and gas within porous media, respectively. g is the gravity acceleration, and H refers to the vertical elevation.
According to Equations (1)–(4), ignoring the second-order phase, the governing equations of the gas-water two-phase flow are obtained.
For the water phase,
For the gas phase, where
is the volumetric strain of the coal solid.
2.3. Deformation Field Equations of Coalbed
According to the theory of elasticity, the coalbed maintains the force balance during CBM mining, and the differential equation of equilibrium is where is the total stress tensor and is the body force tensor.
Through the principle of effective stress, the relationship between the total stress tensor and the effective stress tensor can be written as where is the effective stress tensor. denotes the pore pressure of coalbed. is the Kronecker function, and it is given by
As described above, the coalbed experiences elastic deformation in the CBM mining process; thus, the generalized Hooke law can be used to characterize the stress-strain relationship of the coal. where is the strain tensor, is the Biot coefficient. is the shear modulus of elasticity, and is the lame constant, and they can be written as where is the modulus of elasticity and is the Poisson ratio.
Suppose the coal solid produced only small deformation in the CBM mining process, the following equation can be obtained through the continuity in the deformation process as where and are the displacement components of coal solid.
From Equation (14), we can obtain where is the displacement vector.
Combined with Equations (10)–(19), the governing equations are simplified into
2.4. Fluid Pressure Model at Two-Phase Flow Stage
Combined with Equations (7), (6), and (16), the fluid pressure model in gas-water two-phase flow becomes
This model can be further simplified into
2.5. Analytical Solutions
Before the analytical process, the relationship between the gas density and gas viscosity is derived.
In the theoretical analysis of coal mining, the volume coefficient is usually introduced due to the fact that the density of ground standard state is easier to measure, (,, is the ground standard state density).
Firstly, the compressibility coefficient is introduced. Under isothermal conditions, the compressibility coefficient of the fluid is expressed as
where is the volume of fluid and is the pressure. Since the dimension of pressure is Pa or MPa, the unit of compressibility is or .
Integrating both sides of Equation (21) yields where the density corresponds to the reference pressure under standard atmospheric pressure.
The expression between viscosity and pressure is as follows: where and are the viscosities under pressure and atmospheric pressure , respectively. is the pressure coefficient, which is proportional to void volume and inversely proportional to absolute temperature.
In this paper, for the convenience of distinction, the viscosities of and are corresponding to the pressures of and , respectively, with .
Combining Equations (20) and (21) yields
Taylor series expansion is carried out to at , which yields and where , which is the Lagrange residual term of , with.
2.5.1. Solving Process
Figure 1 shows the schematic diagram of a single-coalbed in CBM mining process. In the single-coalbed CBM development case, both the values of and from model (II) are 0. Suppose the coal seam is infinite cylinder, the solution of model can be converted into the following definite problems (III) and (IV). where is the initial pore pressure in the coal seam, r is the distance to the center of the well, and t is the time. where .
To simplify the calculation of definite problems (III) and (IV), a general form of definite problem (V) is constructed as
Method of variable separation is used in solving problem (V) as follows.
Firstly, assuming that
Then, separation of variables is made after this expression generated into problem (V); after that, order
Equation (33) represents the zero-order Bessel equation. Then, eigenvalue problems will be solved as
It should be noted that is not an intrinsic value. Multiply Equation (34) with , and we have namely,
Integrate to both ends of Equation (36) from 0 to , and we have
If , then that is to say,
Hence where is constant.
Nevertheless, as there must be
is not an intrinsic value.
Making
Considering as , let us discuss the condition of in Equation (34), by
We can obtain
The solution of Equation (45) is where represents zero-order Bessel function and represents zero-order Neumann functions, which have the approximate expression, respectively. where is Euler’s constant, with a value of 0.577216. According to the natural boundary conditions, the value of is 0. And the boundary condition is
If , then
Yet, the number of zero point of zero-order Bessel function can be several. Set the zero point to ; thus,
Intrinsic function is
Substituting the value of into Equation (32), we have
To overlay , the constant is determined as follows: where
Owing to the orthogonal system with weight constituted by ,
is the generalized Fourier coefficient of , namely,
Substituting into Equation (53), the solution of the definite problem (V) is where
Substitute Equations (51) and (62) into Equation (57); hence,
Given , by the definition, the above equation becomes
Substituting Equation (60) into Equation (57), the analytical solutions of definite problem (V) is
2.5.2. Final Solution
Combining Equation (61) with definite problems (III) and (IV), the solution of gas pressure at the two-phase flow stage is where
The solution of water pressure at the two-phase flow stage is where .
2.6. Model Verification
2.6.1. Field Data
The engineering data is from Jincheng 1# test well of Shanxi province, China. The No. 3 seam experiences single coal seam development, whose mineral depth is from 521.6 m to 527.4 m, the minimum horizontal stress on the main seam tested by hydraulic fracturing is 7.9 MPa~9.4 MPa, and the maximum horizontal stress determined from long spacing sonic and coal reservoir fracture pressure is 13.7 MPa~15.7 MPa. Meanwhile, the actual measured gas content is 19.6 m3/t, and the free gas is in dynamic equilibrium with the adsorbed gas, with gas saturation 73.0%. The single-well drainage area is within km2 .
Table 1 describes the daily actual measured gas production and water production of different time and exploitation stages, among of which I-V represents water flow stage, free gas desorption stage, fluctuation desorption stage, stability desorption stage, and attenuation desorption stage, respectively.
Table 2 shows physical parameters of coal bed; various physical parameters will be validated with the theoretical model.
2.6.2. Verification of Theoretical Model on Fluid Pressure
As the analytical solution of the gas and water pressure has been given by Equations (62) and (64), substituting each parameter value into the theoretical model yields the gas pressure near the wellhead:
As shown in Figure 2, the gas pressure curve near the wellbore calculated by the theoretical model was compared with the field pressure data from No. 3 seam of 1# test well. The figure shows that water will be drained from the beginning in the mining process, while gas production occurred after the first 120 days. Firstly, both the pressure of gas and water are decreasing; secondly, the speed of the fluid decrease eventually tends to 0, and the pressure change approaches flattening. This is because as gas desorption occurs with time increasing, correspondingly, the increasing gas migration velocity results in the gas pressure drop. With the effect of desorption, the amount of desorbed gas gradually reduces, thereby leading to the drop of changing rate of gas pressure.
It is given that the original reservoir pressure is 3.99 MPa, gas pressure is 2.69 MPa, critical desorption pressure is 2.41 MPa, and the depletion pressure 0.3 MPa. The pressure obtained from this two-phase separate flow model can well fit the measured field data.
2.6.3. Verification of Theoretical Model on Gas Flow Rate
The shale gas production rate is defined as [22] where .
Based on Darcy’s law, the relationship between shale gas production and water production is obtained by Xue et al. [23].
Equation (67) can be transformed into the final form of water production rate as
Substituting each parameter of Table 2 into Equations (66) and (68) can yield the results of the water and gas production rate. Figures 3 and 4 represent the performance of the water production rate and gas production rate with time, respectively. The figure shows that in the interval of 0~120 days, the coal seam is at the single-phase fluid flow stage from the in situ data, and no gas is output at this stage; in the interval of 120~295 days, the fluid flow is at a nonsaturated stage, and since the speed at this stage does not seek a mathematical expression, the curve is expressed by the boundary value of the single-phase fluid flow stage and two-phase flow stage; furthermore, the seam enters into a stable desorption stage after 295 days. In general, the overall change of the theoretical curve and measured value is of the same trend, such as the gas production and average water production reduce with time going on, and the rate of change slow down. It is necessary to note that the theoretical curves show the trend of the ideal state, while some factors in engineering will influence the trend of gas production rate, such as the amount of water pumped in hydraulic fracturing will affect the stability of daily gas production.
2.7. The Influence of Gas Solubility on Changes of Gas Pressure
For a more comprehensive consideration, the solubility factor is intentionally added in the model for two-phase flow. Figure 5 describes the effect of the solubility on gas pressure. As is shown in Figure 5, as the value of gas solubility in water becomes larger, the decreasing magnitude of gas pressure increases, and the scope of mining becomes wider. This is due to the fact that the greater the solubility of gas in water in depressurization process, the easier the desorption of gas, and the sooner the permeation rate of gas, thus the greater the pressure drop.
3. Gas-Water Two-Phase Separate Flow Model for Multiseam CBM Development
The case of multicoal seam development is discussed in this section. According to the pressure balance principle, due to the pressure gradient in multicoal seam development, if the pressure of the wellbore is higher than that of the coal seam, intrusion of the fluid would happen; namely, the fluid will flow into the coal seam from the wellbore; otherwise, the coal seam produces gas.
Reflected into the model proposed in this paper, as shown in Figure 6, the two-phase flow for multicoal seam development was studied. Assume that the fluid pressure of coal seam 1 is greater than that of coal seam 2, namely, , in that way, the fluid output from coal seam 1 will flow into coal seam 2. For the modeling of fluid pressure from coal seam 1, it meets the solution of definite problems (III) and (IV); but for coal seam 2, it is equivalent for fluid flowing into the coal seam, with the flow rate where and are the strength difference of water and gas between the two coal seams, respectively. is the permeability of coal seam 1. is the diameter of the wellbore. and are the viscosities of gas and water in the coal seam 1, respectively. and are the pressures of water and gas in the coal seam 1, respectively. and are the pressures of water and gas in the coal seam 2, respectively.
3.1. Solving Process
Definite problems (VI) and (VII) of the two-phase flow in coal seam 1 are constructed as follows: where . where and .
To simplify the calculation, the similar definite problem (VIII) with the above problems is built.
Firstly, aiming to construct a function both meeting equation and boundary conditions, the boundary conditions will experience partial homogeneous treatment.
Substituting Equation (61) into the original problem, the conversion of function is made as
Thus, is converted into a mixed problem as
Separation of variables is used to solve the above problem (IX), whose detailed solution procedure is similar to problem (VI). Eventually, the solution is where is the solution of the triangulation equation of , and eigenvalues can be obtained as
Then, the solution for the definite problem (IX) is
Corresponding to the solution of Equation (79), the expressions of gas pressure and water pressure can be obtained, respectively.
For gas flow where
For water flow where .
3.2. Model Verification
3.2.1. Field Data
As Table 3 shows, the field data is from TS-1 well of Liupanshui coal mine within Tianpanjiang production area in Western Guizhou, China [10].
Table 4 shows physical parameters of coal bed, which is from the field data of No. 9 coal seam; various physical parameters will be validated with the theoretical model.
3.2.2. The Analysis of Noninterfering Coal Seam Mining
The case of noninterfering coal seam mining is discussed in this section. Firstly, each parameter value of both coal seams is given in Table 4. Take , that is, to research the area near wellhead. Besides, it is given that the gas desorption begins from the 80th day. Without considering the interaction between the two coal seams, combined with Equations (62) and (64), the pressure of the two-phase flow can be simplified into where the subscript 1g represents gas phase of No. 52 seam, the subscript 1w represents water phase of No. 52 seam, the subscript 2g represents gas phase of No. 9 seam, and the subscript 2w represents water phase of No. 9 seam.
Figure 7 shows gas and water pressure curves of both the No. 52 seam and No. 9 seam in single coal seam mining. The four curves are given by the analytical solution above. As is shown in the figure, firstly, whatever the coal seam is, the gas pressure is always decreasing from the start of gas production, and the gas pressure of No. 52 seam is always greater than that of No. 9 seam; secondly, the solution of water pressure of both coal seams in the fluctuation desorption stage is not studied. For No. 52 seam, we can find that the water pressure is not always making the same trend during different water flow stages. And finally, in the 500th day, the water pressure of No. 52 seam drops to 4.75 MPa. In addition, for No. 9 seam, the water pressure maintains the downward trend, and the water pressure of No. 9 seam drops to 4.39 MPa in the 500th day.
3.2.3. The Effects of Factors on the Development of Multicoal Seam
The effects of factors on the development of multicoal seam are discussed in this section. The impact of altitude difference between the two coal seams is considered, and the remaining factors were not considered temporarily.
Due to the two-phase flow phenomenon in the CBM process, it should be emphasized that the two coal seams interact each other and intrusion of the fluid occurs between the two coal seams in multicoal seam development in the case that (1)Ifthe fluid of No. 9 coal seam flows to the wellhead; meanwhile, a portion of the fluid flows backward into No. 52 coal seam. (2)If
a portion of the fluid of No. 52 coal seam will not flow to the wellhead but flows vertically downward and backward into the No. 9 coal seam.
Figure 8 demonstrates the relationship of pressure difference and height difference of the two seams, both of which are the key factors to influence the development of multicoal seam. As is shown, for water flow, the red solid line represents that there are no influence of water flow between the two coal seams; namely, both the No. 9 seam and No. 52 seam produce water separately and have no effect on the other coal seam. In region ①, water flows from No. 9 seam into No. 52 seam; however, in region ②+③, gas and water from No. 52 seam would flow into No. 9 seam in the mining process. For gas flow, the blue dotted line means that there is no influence of gas flow between the two coal seams. In region ①+②, gas flows from No. 9 seam into No. 52 seam; however, gas would flow from No. 52 seam into No. 9 seam in region ③.
As Figures 7 and 8 show, the fluid pressure of No. 52 seam is always greater than that of No. 9 seam. By calculation of Equations (78)–(81), it can be obtained that the relationship of the pressure satisfies Equation (85), which means a part of the gas and water would flow backward into the No. 9 coal seam in the mining process.
3.2.4. The Analysis of Multicoal Seam Mining
As Figure 6 shows, the field data is from TS-1 well of Liupanshui coal mine; a portion of the fluid (gas and water) of No. 52 coal seam will not flow to the wellhead but flows vertically downward and backward into the No. 9 coal seam.
Substituting the parameter value of Tables 3and 4 into Equation (80), the expressions of the gas flow can be given as
Plotting according to expression from Equation (86), the contrast curves between the development of single coal seam and multicoal seam can be obtained as follows.
As Figure 9 shows, the gas pressure decreases with time regardless of whether there is interaction between the two coal seams. In addition, the gas pressure for the development of multicoal seam decreases faster than that of single coal seam. What is more, the gas pressure for the development of multicoal seam keeps higher than that of single coal seam until the 200th day, and it is lower after the 200th day. It is because when the intrusion happens, the pressure gradient increases, which can accelerate the gas production, followed by the gas pressure dropping faster.
For water phase, the classification discussion will be used as follows: when d, the pressure of No. 52 is greater than that of No. 9 coal seam, and the water intrusion occurs from No. 52 to No. 9 seam; and when d, water intrusion will happen from the No. 9 to No. 52 seam, which does not match the time scope of in this paper.
Thus, substituting the parameter values of Tables 3 and 4 into Equation (82), the expressions of the water flow can be given as follows:
For water flow of No. 52 seam
For water flow of No. 9 seam
According to Equations (87) and (88), the contrast curves of water pressure between the development of single coal seam and multicoal seam can be obtained as follows.
As Figure 10 shows, the water pressure decreases with time regardless of whether there is interaction between the two coal seams. In addition, at the two-phase flow stage, the water pressure for the development of multicoal seam decreases faster. What is more, the gas pressure for the development of multicoal seam keeps lower than that of single coal seam in the whole mining of 500 days. It is because a large amount of drainage occurs in single water flow stage, and plenty of water intrusion will affect the drainage of No. 9 coal seam while considering the impact of multicoal seam; however, the amount of drainage is relatively small when the intrusion happens, and the water pressure gradient will increase in No. 9 coal seam, which can accelerate the drainage, followed by the water pressure dropping faster.
4. Conclusions
Two models of two-phase flow for the development of single coal seam layer and multicoal seam layers were developed, respectively, to simulate the fluid flow in coal seam gas reservoirs in the mining process. These models were then analytically solved by separation of variables and applied to analyze the field data from different CBM formations (the 1# test well of Jincheng and the TS-1 well of Liupanshui coal mine, China). The effects of gas solubility on the gas pressure were investigated. Besides, the impact of interlayer difference on the development of multicoal seam layers was explored. Based on these studies, the following conclusions can be made.
Firstly, the two-phase models for the development of single coal seam and multicoal seam, respectively, are proved to be correct after matching the field data well, respectively. Moreover, the gas pressure is influenced by the gas solubility; if the value of gas solubility increases, the gas pressure will drop increasingly.
Secondly, the key factors incorporating the height difference and the pressure difference of the two seams are proved to influence the multicoal seam development comprehensively. For the field from the TS-1 well of Liupanshui coal mine, during the seal layer mining process, the pressure of the upper No. 52 seam is greater than that of the lower No. 9 seam; thus, the gas and water will intrude into the No. 9 coal seam and stimulate the production and overall shorten the mining time.
Nomenclature
| : | Gas solubility in water |
| : | Saturations of water |
| : | Density of water |
| : | Saturations of gas |
| : | Density of gas |
| : | Density of gas under standard atmospheric pressure |
| : | Velocity of water |
| : | The velocity of gas |
| : | Fluid strength of water |
| : | Fluid strength of gas |
| : | Hamiltonian operator |
| : | Velocity of the coal solid |
| : | Absolute permeability of porous media |
| : | Relative permeability of water |
| : | Relative permeability of gas |
| : | Viscosity of water |
| : | Viscosity of gas |
| : | Viscosity of gas under standard atmospheric pressure |
| : | Pore pressure of water |
| : | Pore pressure of gas |
| g: | Gravity acceleration |
| H: | Vertical elevation |
| : | Volumetric strain of the coal solid |
| : | Total stress tensor |
| : | Body force tensor |
| : | Effective stress tensor |
| : | Pore pressure of coalbed |
| : | Kronecker function |
| : | Strain tensor |
| : | Biot coefficient |
| : | Shear modulus of elasticity |
| : | Lame constant |
| : | Modulus of elasticity |
| : | Poisson’s ratio |
| : | Displacement component in -direction |
| : | Displacement component in -direction |
| : | Displacement vector |
| : | Initial pore pressure in the coal seam |
| r: | Distance to the center of the well |
| t: | Time |
| : | Zero-order Bessel function |
| : | Zero-order Neumann functions |
| c: | Euler’s constant |
| : | Generalized Fourier coefficient. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
The authors are supported by the Independent Innovation Project of “Double First Rate” of China University of Mining and Technology (Grant No. 2018ZZCX04) and the National Natural Science Foundation of China (Grant Nos. 42030810, 51934007).