#### Abstract

The fractal geometry, anisotropy, discontinuity, and non-Darcy flow of tight reservoirs exert a significant effect on well production performance. In this study, the reservoir fractal geometry is represented by exponential functions on the basis of microseismic data, while the discontinuity of the fractures is presented as a nonequilibrium effect. The impact of the nonequilibrium effect and the low velocity non-Darcy flow on the temporal scale of the wellbore pressure is predicted herein. Results showed that the time scale analysis accurately simulates gas flow in a tight reservoir. The wellbore pressure gradually increases, whereas the pressure in the matrix lags when the nonequilibrium effect is considered. The wellbore pressure is affected in the early period by the nonequilibrium effect. However, at the later stage, the pressure in the matrix is mainly affected by the non-Darcy flow. When the non-Darcy flow is dominant, the pores without gas flowing through are better presented.

#### 1. Introduction

Tight gas reservoirs present an attractive option for energy exploitation due to their abundance, resource richness, and adaptability to the energy situation. Nevertheless, extraction of gas from tight reservoirs is challenging. The most economic and viable approach for tight gas reservoir development is multistage fracturing of horizontal wells. Fracturing produces a multispatial scale reservoir structure [1, 2], which in turn significantly affects fluid flow. Flow theories within conventional reservoir models [3, 4] cannot accurately describe the tight gas flow, since they do not account for the heterogeneous reservoir structure. Hence, rigorous characterization of the reservoir microstructure is essential in order to improve the accuracy of tight gas reservoir production models. Microseismic measurements can provide a reasonable description of the reservoir structure at the microscale. The multiscale structure of tight gas reservoirs can then be scaled up as a fractal geometry based on the microseismic data [5, 6].

Microseismic monitoring technology has been widely used in visualization experiments [7] of fractured horizontal wells using video techniques. The rock damage during hydraulic fracturing acts as a source for microseismic waves. These waves can be captured by microseismic monitoring equipment. Subsequently, the shape, size, distribution, and direction of the fractures can be monitored. During the visualization experiments [7], the video is scaled up as a low velocity non-Darcy flow and dynamic coupling, whereas the fractal geometry is determined based on the sound wave propagation. Since the microstructures of various mediums differ from each other, the flow velocity from the matrix to the natural fracture needs a long time to reach the velocity from the natural fracture to the matrix. However, in the traditional model [6, 8–11], the crossflow coefficients from the matrix to the natural fracture are equal to each other, as they are derived based on the shape factor. Thus, the fractal geometry is not considered in the coupling coefficients. The relationship between the microstructure and the shape factor is not addressed. The crossflow reaches instant equilibrium when a pressure difference is established between mediums.

During fracturing, the microseismic event distributes unevenly through the heterogeneous rocks [12–14]. Nevertheless, the fractures propagate in a self-similar manner. Therefore, the heterogeneity is best captured by a fractal geometry model, which significantly improves the model prediction. A continuous function of a spatial coordinate, such as a linear, power, or exponential [15] function, is employed to represent the fractal geometry [16, 17]. The fractal-based model requires a smaller number of data points, is less computationally elaborate, and is more accurate than traditional models [18–22]. The functions describe the spatial variation of the porosity, the permeability, and the crossflow coefficient. The permeability and the crossflow coefficient in the exponential function vary slower than those in the power function, but faster than the linear function beside the inner boundary. The fracture is upscaled by the function. For instance, the permeability and the crossflow coefficient vary with space slower for the linear function. In this case, the reservoir is well fractured. However, the linear function [23] does not match the microseismic data well due to the viscosity and filtration of the fracture fluid. For the power function, the permeability tends to infinity. Such drawbacks are overcome by the exponential function. Nevertheless, the integration between exponential functions and microseismic events has not been reported previously. Moreover, the heterogeneity has only been accounted for in one direction [5, 6, 15, 24–31]. The gas flow model in consideration of second dimensional heterogeneity has not been reported recently, and the effect of filtration on the production has not been considered. Furthermore, heterogeneity does not conform to field data [32–34], and the production data is not estimated accurately.

The discontinuity of fractures promotes different velocities inside natural and hydraulic fractures. However, the flow in the hydraulic fracture is presented by the cubic law and discrete fracture model, which are widely employed in the former literature [35–37]. The crossflow from the natural fracture to the hydraulic fracture is instantly established. However, in fact, the establishment needs a long time and the pressure between the fractures is discontinuous. Instead, some gas may be stored inside the intersections between natural and hydraulic fractures. Consequently, an effect is exerted on the coupling by both the pressure and the stored gas, which is called the nonequilibrium effect. This phenomenon is mainly utilized to describe the viscoelastic fluid flow [38–41]. Different components are equivalent to mediums. The elasticity is described by an additional item in the coupling condition between mediums. The wellbore storage is also the nonequilibrium coupling between the wellbore and the reservoir [11, 22, 42, 43]. The stored gas performs as an elastic fluid when flowing out of the well. The flow velocity from and to the wellbore is not the same instantly. However, few studies have been reported on the stored gas inside the stimulated region.

In addition, an assignable effect is exerted by the low gas flow velocity in the unstimulated region [44, 45] due to non-Darcy flow [11, 46]. Hence, the flow time in the hydraulic fracture lags behind the flow in the natural fracture. Besides, the non-Darcy flow is also attributed to the gas-solid interaction [47]. In most tight gas linear flow models, the non-Darcy flow is presented as the kick-off pressure gradient [48, 49], which was also observed in some experiments [47, 50, 51]. Nevertheless, the propagation of the boundary between Darcy and non-Darcy flows has not been determined yet. Also, the variation of the kick-off pressure gradient with time and space was never taken into account. Lastly, the pressure gradient is sometimes negative in the unstimulated region. Thus, the time needed for the gas to pass the unstimulated region is inaccurately estimated, and the changing pattern in the production is easily missed. This drawback is overcome by the Ikoku-Ramey model [52].

Last but not the least, during production, the velocity and the pressure reduction rate of each medium vary greatly, which is the multiple time scale fluid flow. However, only the production at each time point is described by the traditional flow model [11, 22]. The flow processes, such as fluid flow in different porous mediums and the crossflow, are not described, especially the rapid process. The temporal scale analysis is an important algorithm to describe the flow processes. Moreover, the production reduction is slowed down by the change in the production system. Under unsuitable production system parameters, the fracture may close. Thus, the tight reservoir cannot be further developed. Therefore, the production system is optimized [53] with the aid of the time scale and fed back to the field with the aid of big data technology.

In this study, a detailed description of the reservoir heterogeneity extracted from microseismic measurements is provided. Then, a second dimensional heterogeneous, anisotropic, and dynamic coupling gas flow model is proposed. The crossflow coefficients are different from each other, which is dynamic coupling. They are derived by integrating the Navier-Stokes equation. In this way, the microstructure is upscaled more accurately. Due to the great flow velocity difference between mediums, the crossflow is established after the pressure difference between mediums. Therefore, the coefficients also decrease with space in two directions. The nonequilibrium flow is presented based on the wellbore storage effect. Additionally, the non-Darcy flow is captured by a power function, with an . Since this phenomenon is more evident around vertical wells [54–56], non-Darcy flows in a power function are mainly considered for vertical wells in commercial simulation packages. In contrast, the effect of the non-Darcy flow in horizontal wells has not yet been investigated [57, 58]. Since the flow model with second dimensional heterogeneity, anisotropy, and gas diffusion is difficult to solve semianalytically, a temporal scale analysis approach published earlier by our group is used to solve the model. The result in this work is time continuous. Therefore, the stability and convergence problems caused by the time discretization are avoided, and the trend of pressure with time is described. The wellbore pressure is divided into many parts. The gas flow temporal scale, *t _{i}*s, is the -axis of the time scale distribution diagram, which corresponds to different mediums and flow processes. The time scale is also the gas flow time constant [59]. The time when the dimensionless pressure increases to 63.2% of its final value is presented by the time scale. The trend and propagation of pressure are described by the time scale, while the contribution of temporal scales,

*α*s, is the -axis of the diagram. The contribution of time scales

_{i}*α*s is the changing range of each pressure part and reveals which region is developed. If , the corresponding temporal scale does not exist. The peak in the temporal scale distribution graph means the corresponding factor influences the flow significantly. Consequently, the production system is determined based on the main factor and changed at the start or end of a certain flow period. Lastly, a case study for the validation of the model is presented here. The effect of the heterogeneity of the crossflow coefficient and the non-Darcy flow on the gas flow temporal scales is discussed.

_{i}#### 2. Fractal-Based Tight Gas Flow Model

##### 2.1. Physical Model

The physical model in this work is a tight gas reservoir, based on the microseismic data and an earlier work. The model is presented in Figure 1 and is composed of three regions: the hydraulic fracture (region 1), the stimulated area (region 2), and the unstimulated area (region 3) [60]. The zero point is the intersection of the wellbore and the hydraulic fracture.

##### 2.2. Mathematical Description and Solution

###### 2.2.1. Flow in Region 1

The production rate of the well is assumed to be constant. Coupling between regions 1 and 2 conforms to the nonequilibrium effect. The intersection between the hydraulic fracture and the natural fracture is shown in Figure 2. The NS and continuous equations are expressed as follows: where is the flow regime, is the normal direction of , coordinates and are along , is the temperature, is the conventional gas constant, is the compressibility coefficient, is the gas concentration inside the natural fracture, and is the flow velocity inside the natural fracture. From (1), we get

The solution of (2) is where is a constant. Then, where is a function of and , as given by

Integrating (4) along , the flow rate from the natural fracture to the hydraulic fracture, , is expressed as where is the nonequilibrium coefficient in the hydraulic fracture.

Here, is the number of pores in a unit area. Gas flow inside the finite conductivity hydraulic fracture is given by the following equation: where is the storage coefficient of the wellbore.

###### 2.2.2. Flow in Region 2

Region 2 is composed of the inorganic matrix and fracture networks. The natural fracture, the inorganic matrix, and the kerogen pores have a fractal geometry and a multiscale distribution. Hence, the coupling between the flow mechanisms is dynamic. Based on the microseismic data, the heterogeneity in the flow mechanisms is described by exponential functions. The flow equations are as follows: where is the crossflow coefficient from the natural fracture to the matrix and is the crossflow coefficient from the matrix to the natural fracture. They are different from each other due to the dynamic coupling. The NS equation is integrated through pores [61]. In this way, the crossflow coefficient is obtained.

###### 2.2.3. Flow in Region 3

Region 3 is composed of a homogeneous inorganic matrix. A low velocity non-Darcy flow evolves from this region to region 2. The non-Darcy flow velocity, (m/s), is given by where is the coefficient of the non-Darcy flow and is the exponent of the non-Darcy flow. In the case of Darcy flow in the matrix, . The tighter the matrix is, the higher the power, , is. is the permeability of the matrix. The rate of production, , is expressed as where is obtained as follows:

The flow equation in region 3 is as follows:

##### 2.3. The Relation between Microseismic Data and Heterogeneity

The second dimensional heterogeneity in the model is considered and given by continuous functions, , along the direction, and , along the direction. Taking and as examples, we get where and are the natural fracture permeability and porosity at the zero point, respectively. The permeability and porosity are expressed as where is determined based on the cubic or Poiseuille law, is the area of a single pore, and is calculated as follows [5, 6]: where is the total number of microseismic events in the studied region and and are the probability of event pairs closer than and in the and direction, respectively.

From equations (16) to (20), and can be expressed as

##### 2.4. Microseismic Data Collection

The microseismic data is collected from the field [11], as shown in Figure 3(a). Based on the distance between the microseismic events, and are plotted in Figures 3(b) and 3(c) versus and and matched with the exponential pattern, respectively. The correlation dimension is derived based on Figures 3(b) and 3(c) and is considered in the tight gas flow model.

**(a)**

**(b)**

**(c)**

#### 3. Case Study

##### 3.1. Model Validation

Parameters of the tight gas reservoir considered in this study were published previously [11] and are shown in Table 1. The model validation is conducted by comparing the temporal scale analysis results and the semianalytical solution [22] of a previous work. Figure 4(a) shows that the temporal scale analysis result is in good agreement with the semianalytical solution. The error analysis is shown in Figure 4(c). In the early period, the minor differences mainly arise from numerical dispersion. The error decreases with time initially and then increases at the dimensionless time of around 200. In the later period, the error originates from the assumption of the linear flow. Afterwards, the error decreases again.

**(a)**

**(b)**

**(c)**

At the beginning, there is a fast gas flow out of the hydraulic fracture, which increases the wellbore pressure. Then, the gas starts to flow out of the natural fracture. The natural fracture permeability beside the inner boundary is higher. The rate of wellbore pressure increase slows down. At a later stage, the gas struggles to flow from the matrix. However, more gas flows out of the natural fracture and the hydraulic fracture. Consequently, the wellbore dimensionless gas pressure increases faster.

The temporal scale diagram is presented in Figure 4(b). Since there is high permeability in the fracture, the gas flows in the fracture faster at a dimensionless time equal to 10^{-5}. The left peaks in the figure correspond to the fracture. The flow velocity in the hydraulic fracture is higher. Therefore, the peak is quite low but slightly higher than the peak around (Figure 3(b)). The peak around corresponds to the natural fracture. Since the natural fracture is heterogeneous, the area with very low permeability behaves as a “nonsealing fault.” The gas struggles to flow out of the “nonsealing fault.” More peaks and dips appear around the temporal scale to . The permeability and porosity are distributed continuously, and thus, the peaks connect to each other. Since the outer boundary is tight, the right side of each peak gets higher, and the height of the peaks increases with time. The peaks around the temporal scale to correspond to the matrix beside the outer boundary. More gas remains in the matrix and takes a much longer time to flow to the wellbore. Consequently, the peaks at the longer temporal scale, , are much higher and less frequent than the ones occurring at a shorter temporal scale.

##### 3.2. Field Study

The comparison between temporal scale analysis results and field results is shown in Figure 5. The comparison is reasonable. First, the gas flows rapidly out of the hydraulic fracture. The production rate is high in the early stages. Then, the gas flows from the natural fracture and the matrix at a slower rate, leading to a rapid production decline. At a later stage, as the pressure of the fracture declines, gas crossflow occurs. Lastly, more gas flows to the hydraulic fracture and the well production rate declines more slowly.

##### 3.3. The Nonequilibrium Effect

Figure 6(a) shows the wellbore dimensionless gas pressure at different values of the nonequilibrium effect coefficient, . The pressure derivative is shown in Figure 6(c). The gas flow velocity to the wellbore is presented by the pressure derivative. During the gas flow inside the hydraulic fracture, the gas pressure inside the hydraulic fracture decreases. As a consequence, the pressure derivative decreases. At the dimensionless time of around 0.3, the gas flows out of the natural fracture. Due to the high permeability area beside the inner boundary, the pressure derivative increases sharply. As a result of the tight outer boundary, less gas flows from the outer boundary, and the derivative reaches a plateau. At the dimensionless time 10, since the gas flows out of the matrix, the derivative increases slightly faster. Afterwards, the gas flows from region 3 and the derivative rises sharply again. Eventually, the gas struggles to flow from the outer boundary and the derivative increases slowly. Only the pressure at the early period is affected by . At higher , more gas is stored at the intersection between the natural fracture and the hydraulic fracture. The gas flows into the wellbore in the early stage. Consequently, in the early period, the wellbore dimensionless pressure is lower. Afterwards, more gas at the intersection flows out. As a result, the pressure derivative is about 0.3, for the case . The pressure derivative is about 0.7, for the case , while it is about 0.5, for the case . The gas in the matrix flows into the intersection. The rate at which the pressure increases slows down, and the derivative decreases faster at higher . The pressure at the intersection recovers. The dimensionless gas pressure curve approaches the case . Minimal effect is exerted on the approaching time by . The pressure derivative increases to about 1 in the late stage.

**(a)**

**(b)**

**(c)**

The temporal scale diagram is shown in Figure 6(b). The flow time is 10^{-4} in the hydraulic fracture and 10^{-3} to 0.01 in the natural fracture. At a higher , the gas at the intersection of the natural fracture and the hydraulic fracture flows to the wellbore, instead of the gas in regions 2 and 3. The gas pressure propagates to regions 2 and 3 for a longer time. Therefore, in the left side of Figure 6(b), the peaks shift to the right. In addition, less gas inside the natural fracture contributes to production. Consequently, the peaks before the time scale 10^{-2} become lower. At a later stage, the pressure propagation accelerates. As a result, the time scale of the peak between the time scales 0.01 and 0.05 is less significant. The coupling time between regions 1 and 2 is about 0.1, and the dimensionless coupling time between natural and hydraulic fractures is about 0.1. At higher , more gas flows from the intersection between regions 1 and 2, instead of region 3. Therefore, the corresponding peak is higher. When gas flows from the intersection, less gas flows out of the matrix beside the inner boundary. A more significant effect is exerted on the inner boundary than on the outer boundary. The gas dimensionless flow time in the matrix is 1 to 50. The peaks at the right side of the time scale 1 are lower, especially the left part of the peak. Due to the existence of the nonequilibrium effect, the development time of the matrix lags.

##### 3.4. The Effect of Heterogeneity Conditions

Figure 7(a) shows the wellbore dimensionless gas pressure at different heterogeneity conditions, and the pressure derivative is presented in Figure 7(c). In the first dimensional heterogeneity case, the gas pressure propagates along the direction. The distance to the wellbore has little effect on the gas flow velocity. Therefore, the pressure derivative curve reaches a plateau in the early period. Afterwards, the gas flows out of the matrix and the pressure derivative rises. In the second dimensional heterogeneity case, due to the tighter outer boundary, the pressure in the hydraulic fracture beside the outer boundary does not recover quickly. The initial pressure increases faster, but only for a short time. The increase rate is enhanced at an earlier stage than in the first dimensional case. The pressure at a later stage is still higher than that of the first dimensional case. As the gas pressure struggles to propagate in the tight outer boundary, the pressure curve eventually corresponds to that of the first dimensional case. The fracture permeability decreases along the direction. With the propagation of gas pressure, the flow velocity into the hydraulic fracture decreases. Therefore, the pressure derivative increases with time in the early period. Afterwards, the derivative curve is similar to the first dimensional case. Since less gas flows out of the fracture, the curve is higher than the first dimensional case. Only the pressure between the dimensionless time 1 and 1000 is affected by the heterogeneity crossflow. As the gas pressure propagates to the outer boundary, the gas struggles to flow from the matrix to the fracture. Consequently, the increase rate of the gas pressure is slightly higher only at the beginning of the gas flow into the matrix. The slope of the gas pressure derivative curve gradually increases. Afterwards, because of the lower gas pressure in region 2, more gas flows out of region 3. Thus, the rate declines, and the slope of Figure 6(c) stops further increasing. The gas pressure curve coincides with the one for the case with homogeneous crossflow.

**(a)**

**(b)**

**(c)**

In the second dimensional heterogeneous case, the right side of the peaks around the time scale 10^{-3} to 10^{-2} is higher than that in the first dimensional case. The left part of the peak around 0.1 is steeper. More gas flows to the wellbore. Therefore, a small peak appears around the time scale 0.03. In addition, the gas inside the matrix struggles to flow out at different temporal scales, especially beside the outer boundary. As a consequence, the peak around the time scales 1 to 10 becomes higher and wider, especially to the left of the peak. A small peak appears around the time scale 10. The flowing time in the matrix is 1 to 10. The area with large heterogeneity behaves as a “nonsealing fault.” As the gas pressure propagates through the “nonsealing fault,” the gas in the matrix flows into the “fault” more easily than the gas in the fracture. The dip between the two peaks corresponds to the “nonsealing fault,” while the gas beside the outer boundary also struggles to flow to the wellbore. The dimensionless crossflow time in the matrix beside the outer boundary is 20. Therefore, the peak around 20 becomes lower. Only the peaks to the right of time scale 1 are affected by the heterogeneous crossflow. As shown by the blue line, since the crossflow coefficient beside the inner boundary is higher than that beside the outer boundary, the crossflow rate decreases. The peak around 1 to 10 is consequently lower on its left side. The crossflow mainly occurs beside the inner boundary. A small peak appears around the time scale 1, while the peaks between 10 and 10^{3} become higher, especially to the right of the peaks between 10 and 10^{3}. The flow time in region 3 increases to 10^{3}. At the “nonsealing fault,” the gas in the matrix also has great difficulty in flowing out. Therefore, the dip between 10 and 20 disappears, while the left side of the peak between 20 and 100 becomes flat.

##### 3.5. The Effect of Non-Darcy Flow

The wellbore gas pressure in the case with low velocity non-Darcy flow is presented in Figure 8(a). Only the pressure at a later period is affected by the low velocity non-Darcy flow. At a higher non-Darcy flow exponent, , a lower number of small pores allow gas to flow through. The gas mainly passes through large pores and the driving force becomes lower. The wellbore pressure becomes lower than the case without a non-Darcy flow and becomes steady. At a lower non-Darcy flow coefficient, , the gas struggles to flow from region 3. Therefore, the gas pressure is higher. Without considering the non-Darcy flow, the later gas pressure is overpredicted, and the pores, where gas does not flow through, are not properly presented.

**(a)**

**(b)**

The temporal scale diagram is presented in Figure 8(b). Due to the non-Darcy flow, less gas flows out of region 3. The matrix beside the outer boundary contributes less to the production. Therefore, more gas flows out of regions 1 and 2 to the wellbore. The peak around the time scale 1 to 10 is higher to its right, while the peak around the time scale 2 falls greatly. Since the flow velocity in smaller pores is slower, the peak around 100 is shifted to the right and increases in height. Moreover, from time scale 20 to 100, a wider peak appears. The dimensionless flow time in larger pores is 20 to 100. In this case, the gas flows easily through pores with large sizes. Thus, the corresponding peak is lower. The difference between the pore size is better presented by the non-Darcy flow. At lower , the gas flow velocity in region 3 further declines. The peak around 10 is much higher than the case with lower , while the right peaks shift further to the right. Moreover, fewer pores allow gas to flow through. The gas flow time through region 3 is longer, and the distance between the peaks becomes larger. A small peak appears between the major peaks. More gas flows from regions 1 and 2. The peak around 100 becomes higher, while the one at the time scale 1000 becomes lower.

#### 4. Conclusions

This study shows a gas reservoir model of a horizontal well considering a nonequilibrium flow and a non-Darcy flow can be used to describe the reservoir microstructure. A gas flow model considering heterogeneous, anisotropic, and nonequilibrium effects as well as dynamic coupling was established and solved using a temporal scale analysis approach. A case study was analyzed to validate the model, and the effect of heterogeneity of the mass transfer coefficient on the gas temporal scales was discussed. The microseismic data are in good agreement with the field data using an exponential function. The increase in wellbore pressure first slows down, then gradually becomes rapid. The discontinuity is described accurately by the nonequilibrium effect. Pressure and its derivative are affected in the early period by the nonequilibrium effect, whereas they are affected by the non-Darcy flow at a later period. Time of pressure propagation through the matrix lags in the case with the nonequilibrium effect, and the gas flows out of the matrix and natural fracture later. In the case with the heterogeneous crossflow coefficient, the matrix is less developed, especially that beside the outer boundary. Moreover, the heterogeneous flow is better presented in the temporal scale distribution diagram. The non-Darcy flow is well described by a power function for the first time. Pores through which gas does not flow are better handled by the non-Darcy flow exponent. However, at lower non-Darcy flow values, the reservoir is tighter, and the flow velocity further decreases.

#### Symbols

: | Flow regime |

: | Concentration |

: | Gas compressibility coefficient |

: | Nonequilibrium effect coefficient |

: | Total compressibility (Pa^{-1}) |

: | Probability of event pairs closer than in the direction |

: | Probability of event pairs closer than in the direction |

: | Non-Darcy flow coefficient |

: | Anisotropy coefficient |

: | Reservoir height (m) |

: | Permeability (m^{-12}) |

: | Non-Darcy flow exponent |

: | Number of hydraulic fractures |

: | Gas pressure (Pa) |

: | Gas production (m^{3}/s) |

: | Conventional gas constant |

: | Temperature |

: | Time (s) |

: | Temporal scale |

: | Hydraulic fracture width (m) |

: | Length along the hydraulic fracture (m) |

: | Length along the wellbore (m) |

: | Compression factor. |

*Greek Letters*

: | Contribution of each temporal scale |

: | Viscosity (Pa·s) |

: | Transfer coefficient (s^{-1}·Pa^{-1}) |

: | Pseudopressure function (Pa/s) |

: | Porosity. |

*Subscripts*

0: | Initial |

D: | Dimensionless |

f: | Natural fracture |

fm: | From the natural fracture to the matrix |

j: | The index of temporal scales |

lf: | Hydraulic fracture |

m: | Matrix |

mf: | From the matrix to the natural fracture. |

#### Data Availability

No data were used to support this study.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was supported by the Postdoctoral Application Research Project in Qingdao (qdyy20200066), the Postdoctoral Self-Innovative Project (20CX06094A), and the China Scholarship Council.