Abstract

Carbonate reservoirs typically have complex pore structures, so the production wells typically have high production in the early production stage, but they decline rapidly. It is highly challenging to achieve accurate interpretation results. In this paper, a new and practical methodology for production data analysis of fractured and fractured-vuggy carbonate reservoirs is proposed. Firstly, analytical solutions to characterize the different multipore media and simulate transient production behavior of fractured and fractured-vuggy carbonate reservoirs during the transient flow regime are presented. Then, a new function and that related to the dimensionless production rate is introduced, and a series of new decline type curves are drawn to make a clear observation of different flow regimes. In addition, the effects of the storativity ratio, interporosity flow coefficient, skin factor, and dimensionless radial distance of external boundary on production performance are also analyzed. Finally, two example wells from the fractured and fractured-vuggy carbonate reservoirs are used to perform rate decline analysis with both the Blasingame type curves and the new type curves. The validation of the new method is demonstrated in comparison to the results of well test interpretation. The results show that the curves of vs. are -shaped for dual-porosity reservoirs and -shaped for triple porosity reservoirs and also indicate that the interpreted parameters such as permeability, skin factor, storativity ratio, and interporosity flow coefficient using new decline type curves are aligned well test interpretation. In correlation with other traditional well test analysis, this approach effectively reduces the multisolution probability of interpretation.

1. Introduction

Carbonate reservoirs, which approximately account for 60% of the world oil reservoirs, represent a significant amount of oil and gas reserves and have a great potential to be exploited [1]. With further exploration and development, these fractured (see Figure 1(a)) and fractured-vuggy (see Figure 1(b)) carbonate reservoirs have been attracting considerable attention in recent years [25]. The complex structures of the reservoir are characterized by the presence of different porous media and transmissibility [6]. Production wells of carbonate reservoirs typically have high production in the early production stage but they decline rapidly, leading to the difficulties of further exploitation [711]. At present, depletion development is mainly applied to extract oil and gas from carbonate reservoirs, and the pressure drawdown is quite large to have a considerable production rate. Nevertheless, it is difficult to predict the production and analyze production data for fractured and fractured-vuggy carbonate oil reservoirs.

In order to characterize the complex pore structures and flow features of fractured and fractured-vuggy carbonate oil reservoirs, the concept of multiple-porosity is proposed. For the dual-porosity, Barenblatt et al. [12] hypothesized that there are two continuous media, matrix and fracture. The porosity of fracture is much less than that of matrix, but the permeability of fracture is much greater than that of matrix. Warren and Root [13] proposed the idealized dual-porosity model and solved the governing equation of pseudosteady flow analytically. They assumed that the fracture system evenly divides the matrix rock. The matrix system which is in a hexahedral shape is homogeneous and isotropic. Kazemi [14] proposed a new model which is assumed that the reservoir is composed of horizontal fractures and matrix, and the fracture system evenly divides the matrix system. De Swaan [15] developed a model and assumed that the configuration relationship between the matrix system and the fracture system is similar to that in the Warren-Root model, except that the matrix rock blocks are round sphere. For the triple-porosity, Abdassah and Ershaghi [16] firstly proposed the triple-porosity model, which considers that the reservoir is composed of three media: matrix, fracture, and vug. The model is assumed that the fracture communicates directly with the wellbore, and the matrix and vug simultaneously have an unsteady flow to the fracture. Jalali and Ershaghi [17] considered the heterogeneity of the reservoir and analyzed the unsteady pressure behavior of such reservoirs. Al-Ghamdi and Ershaghi [18] developed the triple-porosity/dual-permeability model, which proposed two types of flow units to differentiate between microfractures and macrofractures. Based on the theories and models of these pioneers, in-depth research of the dual-porosity and triple-porosity reservoirs can be found in the recent literature [19, 20].

In recent years, empirical, numerical, and analytical methods have been widely used for production data analysis of fractured and fractured-vuggy carbonate reservoirs [2024]. Because of its relative simplicity, empirical methods are commonly applied in the actual oil field to conduct history matching and make forecasts for production wells [21, 25, 26]. Arps [25] decline curve analysis is a classic empirical method used mainly to predict recoverable reserves and future production rates. However, this method is based on empirical observations of production decline. Furthermore, the decline rate greatly dominates the decline curve analysis. Conventional methods have been the use of numerical simulation for performing production data analysis of fractured and fractured-vuggy carbonate oil wells [22, 27, 30, 31]. However, it is time-consuming in gridding and simulating, which makes numerical models less attractive. Analytical models not only capture the reservoir architecture of carbonate reservoirs but also have much higher computational efficiency than numerical models, which have recently attracted more and more attention [2531]. Fetkovich [32] expanded Arps’s approach to cover the transient flow condition using analytical flow equations. However, this method is only suitable for a well producing at constant bottom-hole pressure (BHP) condition while the Pressure-Volume-Temperature (PVT) properties are not considered. Blasingame et al. [33] proposed a rate decline analysis method that accounts for variable BHP and production rate conditions. This allowed the analysis of well production rate with variable BHP and variable production rate, by reducing them to an equivalent constant-rate data. In addition, the PVT properties are also incorporated in the models. Therefore, Blasingame type curves have been widely used for field data analysis. Agarwal et al. [34] expanded Fetkovich’s and Blasingame’s work to develop a new type decline curve based on well test definitions, as opposed to the dimensionless variables in the Blasingame method. By such, an approach to the rate decline analysis is to reduce the multisolution to certain content. In practice, the boundary-dominated regime should be observed on the type curve to avoid multisolution. Based on these theories and methods, extensive researches and field applications have been carried in the literature [3539]. However, when the wells exhibit a long transient flow regime and the boundary-dominated regime could not be observed, especially the flow behavior characteristics of fractured and fractured-vuggy carbonate reservoirs, these methods are generally led to multisolution. Therefore, proposing an effective method to capture the production performance of wells during the transient flow regime is of significant importance [40].

Firstly, a novel methodology is presented to analyze production data of wells from fractured and fractured-vuggy carbonate reservoirs during the transient flow regime. Then, a traditional dual-porosity and triple-porosity model is applied to model the fractured and fractured-vuggy carbonate reservoirs, respectively. With Laplace transform, the analytical solution of the mathematical model is obtained. In order to develop the new type curves (see Table 1), functions related to the dimensionless production rate are introduced. Furthermore, a series of new type curves are drawn to make a clear observation of different flow regimes. The effects of some critical parameters on production performance are also analyzed with the proposed model. In addition, the new method is benchmarked against the well test interpretation. Finally, two field cases are used to demonstrate the practicability of the method. The main contribution of this paper is the provision of a novel approach for production data analysis of fractured and fractured-vuggy carbonate oil wells during the transient flow regime. Compared to other methods, the accuracy of the proposed method for production data analysis is higher. This new method provides another choice for field engineers and related research and further reduces the multisolution probability by comparing among the different interpretation methods. Although the traditional multiple media model has been studied in this paper, this new method is also applicable to other complex fractured reservoirs. The corresponding research work will be carried out in the next step.

2. Methodology

2.1. Physical Model and Assumptions

It can be seen from Figure 2(a) that the fractured carbonate reservoir is modeled with a dual-porosity model, which is structured by fracture and matrix systems. The single-phase oil flows from the matrix system to the fracture system, and then to the wellbore (see Figure 2(b)). The fractured-vuggy carbonate reservoir is modeled with a triple-porosity model, which is structured by fracture, vug, and matrix systems (see Figure 2(c)). The oil firstly flows from the vug system to the fracture system, then from the matrix system to the fracture system, eventually to the wellbore (see Figure 2(d)). The matrix-fracture interporosity flow and the vug-fracture interporosity flow in the reservoirs are described by the pseudosteady state model [9, 16, 19]. The radial cylindrical dual-porosity and triple-porosity reservoirs are considered in which a single production well is located at the center, completely penetrating the formation. Some simplifying physical model assumptions for the derivation of the governing equation are as follows: (1)Total compressibility (rock and single-phase fluid) is a constant with a low value(2)Isothermal and Darcy flow is applicable for the models(3)Capillary and Gravity forces are neglected(4)The initial pressure is uniformly distributed in the formation(5)Changes of fluid density and viscosity are neglected(6)The matrix, fracture, and vug systems are homogeneous, respectively

2.2. Dimensionless Mathematical Model

In this paper, a novel approach is proposed to perform production data analysis of wells from fractured and fractured-vuggy carbonate reservoirs during the transient flow regime. The proposed rate decline analysis method is very similar to the well test analysis, and the dimensionless variables are based on the conventional well test definitions, as opposed to the Blasingame dimensionless variables. By such, an approach to the rate decline analysis is to reduce the multisolution to certain content. After some reasonable assumption and simplicity, the conventional mathematical modeling of a well in a dual or triple porosity reservoir has been established in a radial cylindrical system [13, 16, 20], as shown in Appendix A. To solve the modeling, we use the dimensionless effective radius to establish the differential equations. The dimensionless definitions of all kinds of parameters are shown in Appendix B. The dimensionless mathematical modeling is as follows:

For a dual-porosity reservoir [13, 20],

For a triple-porosity reservoir [19],

Definite conditions are shown by taking the triple-porosity reservoir as an example:

Initial condition:

The well is assumed to produce at constant rate production, and the inner boundary can be written as

The outer boundary is assumed to be closed and given by where , , and are dimensionless pressure of matrix, fracture, and vug systems, respectively; , , and are fluid capacitance coefficient of matrix, fracture, and vug systems, respectively; is the interporosity flow factor of matrix system to fracture system; is the interporosity flow factor of vug system to fracture system; is the dimensionless radial distance; is the dimensionless radial distance of side external boundary; is the dimensionless production time; is the skin factor.

2.3. Solution to the Mathematical Model

By the use of Laplace transform and inverse Laplace transform methods, the model can be solved. According to the process of solving the model (see Appendix C), the dimensionless pressure in Laplace space for constant rate production can be expressed by [13, 16, 20, 4042]

For a homogenous reservoir,

For a dual-porosity reservoir,

For a triple-porosity reservoir,

Based on the Duhamel principle [43], the solution in Laplace space for modeling well productivity at constant wellbore pressure production can be obtained by where and are dimension pressure and production in Laplace space, respectively; is the Laplace transform variable; is the modified Bessel function of the first kind, zero order; is the modified Bessel function of the first kind, first order; is the modified Bessel function of the second kind, zero order; K1() is the modified Bessel function of the second kind, first order.

3. New Type Curves

3.1. Description of the New Type Curve Method

In real space, the dimensionless production rate () can be obtained using Stehfest numerical inversion [44] to convert back to . Therefore, we can obtain the standard log-log type curves of , as shown in Figure 3. It can be seen from Figure 3 that the type curves of different dimensionless radial distance of external boundary () overlap during the transient flow regime. In addition, the type curves separated gradually with the increase of . This new type curve matching method not only is accurate and practicable but also reduces the multisolution to a certain extent. Particularly, the production well exhibits a long transient flow regime and a boundary-dominated regime could not be observed.

In order to develop the new type curves, a function related to the dimensionless production rate is introduced [45].

Equation (14) can be rewritten as

According to Equation (15), the standard log-log type curves of can be obtained. Figure 4 shows the type curves of vs. for a homogenous reservoir.

Furthermore, we need to introduce the derivative of the function to form the rate decline analysis type curves.

The differential form of Equation (14) is written as

The derivative of the function can be written as

Taking Equation (16) into Equation (17), one can obtain where is the dimensionless production rate and is the inverse Laplace transformation.

According to Equation (18), the standard log-log type curves of can be obtained. Figure 5 shows the type curves of vs. for a homogenous reservoir.

Combining Equation (13) with Equations (15) and (18), we can obtain the standard log-log type curves of rate decline analysis of , , and vs. . Figure 6 shows the type curves of rate decline for a homogenous reservoir. When the theoretical curve is fitted with the measured curve, the permeability, skin coefficient, and other parameters can be interpreted.

3.2. Type Curves and Flow Regime Recognition for Multiple Media

Based on the solutions of the proposed models, the corresponding type curves for multiple media can be obtained. In this section, type curves for production dynamics are plotted and the flow regimes are analyzed as well.

Figure 7 shows the standard type curves of rate decline analysis for a well producing located at a dual-porosity reservoir. An entire transient flow process is clearly shown, and the following four main flow periods can be recognized in Figure 7.

Period I is the skin effect and early transition stage. The curves of , , and vs. exhibit a downward line.

Period II is the interporosity flow stage of the matrix system to the fracture system. The curve of vs. assumes a -shaped hump shape, which represents the interporosity flow of the matrix system to the fracture system. This period is controlled by the storativity ratio of fracture () and the interporosity flow coefficient of the matrix system to the fracture system ().

Period III is the entire radial flow stage of the matrix and fracture systems. The curves of and vs. converge to a level straight line.

Period IV is the external boundary response stage. The curves of and vs. ultimately converge to a straight line with a slope of “-1,” which is the reflection of the pseudosteady state flow.

Figure 8 shows the standard type curves of rate decline analysis for a well producing in a triple-porosity reservoir. The following six main flow periods can be recognized:

Period I is the skin effect and early radial flow stage.

Period II is the interporosity flow stage of the vug system to the fracture system. The curve of vs. assumes a -shaped hump shape, which represents the interporosity flow of the vug system to the fracture system. The shape and location of -shaped “hump” are controlled by the storativity ratio of vug () and the interporosity flow coefficient of the vug system to the fracture system ().

Period III is the entire radial flow region of the matrix and fracture systems. The curves of and vs. converge to a level straight line.

Period IV is the interporosity flow stage of the matrix system to the fracture system. As stated above, this period is controlled by the storativity ratio of fracture () and the interporosity flow coefficient of the matrix system to the fracture system ().

Period V is the entire radial flow region of the vug, matrix, and fracture systems. The curves of and vs. converge to a level straight line.

Period VI is the external boundary response stage. The curves of and vs. ultimately converge to a straight line with a slope of “-1.”

When comparing the type curves of a dual-porosity reservoir with those of a triple-porosity reservoir, the difference is that the curves of vs. is -shaped for a triple-porosity reservoir, with the left -shaped being the reflection of the interporosity flow of the vug system to the fracture system and the right -shaped being the reflection of the interporosity flow of the matrix system to the fracture system.

3.3. Characteristics of Transient Production Dynamics

Based on the proposed model, the effects of some critical parameters on production performance are analyzed.

Figures 9 and 10 illustrate the effects of the storativity ratio () and interporosity flow coefficient of the vug system to the fracture system () on the production dynamics of a well in fractured-vuggy carbonate reservoirs, respectively. The represents the relative capacity of fluid stored in the vug system. The bigger is the response of relatively abundant reserves in the vug system. The mainly affects the shape of the left -shaped “hump,” as shown in Figure 8; a bigger leads to a higher -shaped hump. The mainly affects the location of the left -shaped “hump” (Figure 9). A smaller leads to the later the time of interporosity as depicts the starting time of the interporosity flow of the vug system to the fracture system.

Figures 11 and 12 show the effects of the storativity ratio () and interporosity flow coefficient of the matrix system to the fracture system () on type curves, respectively. The represents the relative capacity of fluid stored in the matrix system. As shown in Figure 10, the mainly affects the shape of the right -shaped “hump.” A bigger leads to a higher -shaped hump, which means comparatively abundant reserves in the matrix system. The mainly affects the location of the right -shaped “hump” (Figure 11). Because represents the starting time of the interporosity flow of the matrix system to the fracture system, therefore, the smaller the is, the later the time of interporosity is.

The influence of skin factor on type curves is also analyzed, as shown in Figure 13. It can be seen that a smaller skin factor leads to a higher location of the dimensionless rate curve in the early period. In addition, the smaller the skin factor, the slower the rate declines. This is because the incremental value of skin factor results in the increasing additional filtration resistance and the skin effect transition period will last longer. Moreover, the larger the skin factor, the slower the pressure wave propagates to the external boundary and the faster the rate declines.

Figure 14 shows the influence of dimensionless radial distance of external boundary on type curves. One can find that mainly influences the late flow periods. The smaller the , the faster the rate declined after the pressure wave propagated to the closed boundary. This is because the external boundary response would be delayed when the increased. Therefore, the production drop should decrease to maintain the constant BHP of the production well.

4. Model Verification and Field Application

In this section, we validated the proposed models using two field cases from fractured carbonate oil reservoir and fractured-vuggy carbonate oil reservoir.

4.1. Validation Case 1: A Fractured Carbonate Oil Well

This field example is a vertical well completed in a fractured carbonate oil reservoir in the Tarim Basin of China. Table 2 shows the basic formation and well parameters that are crucial for well test interpretation and rate decline analysis.

The well was put into production on October 15, 2011. During the initial production stage, the well had a production rate of 183.3 m3/d. From October 2011 to March 2015, the total production rate declined from 183.3 m3/d to 83.7 m3/d. The approximate accumulated oil production rate was  m3. The formation pressure was 56.72 MPa.

We first plotted pressure and the pressure derivative curves using the original buildup test data. Then, we chose the well-known commercial software of well test analysis Saphir (Ecrin v4.12) developed by the French company Kappa, to conduct well test interpretation. Based on the real field data, we plotted the curves of , , and vs. . We also chose the well-known commercial software of rate decline analysis Topaze (Ecrin v4.12) to perform rate decline interpretation with Blasingame decline type curves. The matching curves of well testing interpretation and rate decline interpretation using a dual-porosity model in Saphir and Topaze are shown in Figures 15(a) and 15(b), respectively. It can be seen from Figures 15(a) and 15(b) that the matching effects are all satisfactory. The main interpretation parameters are shown in Table 2.

At the same time, we plotted the curves of , , and vs. and the proposed method was used to perform rate decline interpretation of the same well. The matching curves of rate decline interpretation by using a dual-porosity model with a closed boundary for side are shown in Figure 14(c). As shown, a good match between the new type curves and field data is obtained. The main interpreted parameters using the new method are also presented in Table 3.

According to the matching results shown in Table 2, there is little difference in interpretation results between well test interpretation and the proposed method. However, there is a big difference in the matching results between well test interpretation and the Blasingame method, especially the characteristic parameters and in the early transient flow period. This is because the production well shows a long transient flow regime and a boundary-dominated regime could not be observed. Compared with the Blasingame method, better history matching results are observed and more reasonably interpreted parameters are obtained, which reveals the new method’s reliability in analyzing production data during transient flow regimes.

4.2. Validation Case 2: A Fractured-Vuggy Carbonate Oil Well

To test the practical applicability of the proposed method, field data obtained from a vertical well completed in a fractured-vuggy carbonate oil reservoir in the Tarim Basin of China is also analyzed. Table 4 summarizes the basic formation and fluid parameters. The well was put into production on January 14, 2010. The initial production rate was 85.84 m3/d. The total production rate declined from 85.84 m3/d to 51.23 m3/d from January 2011 to January 2016. The approximate accumulated oil production rate was  m3, and the formation pressure was 56.91 MPa.

Based on the real field data, we also chose the well-known commercial software of well test analysis Saphir and rate decline analysis Topaze (Ecrin v4.12) to perform the well test interpretation and rate decline interpretation, respectively. The matching curves of well testing interpretation and rate decline interpretation using a triple-porosity model in Saphir and Topaze are shown in Figures 16(a) and 16(b), respectively. It can be seen from Figures 16(a) and 16(b) that both of the matching effects are satisfactory. Table 4 shows the main interpretation parameters. Then, the proposed method is used to perform rate decline interpretation of the same well. The matching curves of rate decline interpretation using the new method are shown in Figure 16(c). A satisfactory matching result is obtained, and the main interpreted parameters are also shown in Table 5.

As can be seen from Table 5, the final matching results between well test interpretation and the new method are less than 3%. However, there is a big difference in the interpretation results between well test interpretation and the Blasingame method, especially the characteristic parameters and in the early transient flow period. This is because the production well is still in transient flow and the reservoir boundaries have not been reached; Blasingame decline curves should not be expected to be applicable. Therefore, the proposed method can accurately undertake production data analysis of wells from fractured and fractured-vuggy carbonate reservoirs during the transient flow regime.

5. Summary and Conclusions

The highlight of this paper is a novel approach for production data analysis of wells from fractured and fractured-vuggy carbonate reservoirs during the transient flow regime. The main conclusions are as follows: (i)All flow regimes, including transient flow and late boundary-dominated flow regimes, can be recognized clearly from the new decline type curves, which reveal the efficiency of the new approach in analyzing production data during transient flow regimes(ii)The type curves for the dual-porosity and triple-porosity reservoirs are quite unique. The curves of vs. are -shaped for dual-porosity reservoirs and -shaped for triple porosity reservoirs(iii)The storativity ratio of the vug system (), interporosity flow coefficient of the vug system to the fracture system (), the storativity ratio of matrix (, interporosity flow of the matrix system to the fracture system (), skin factor (), and dimensionless radial distance of external boundary () plays the most significant role in production performance, which should be all considered in the models(iv)Two field cases from the Tarim Basin are used to show the accuracy and practicability of the new method. Compared with the Blasingame method, better history matching results are observed and more reasonable interpreted parameters are obtained, which reveal the reliability of the new method in analyzing production data during transient flow regimes

Appendix

A.1. Mathematical Modeling

In a cylindrical coordinate system, the governing differential equations for well productivity modeling are written as

For a dual-porosity reservoir [13, 20],

For a triple-porosity reservoir [16],

Definite conditions are shown by taking the triple-porosity reservoir as an example:

Initial condition:

Inner boundary condition:

Exterior boundary condition:

A.2. Dimensionless Definitions

The dimensionless definitions are as follows [20]: (1)For a dual-porosity reservoir

Dimensionless radial distance .

Dimensionless radial distance of external boundary .

Dimensionless production time .

Dimensionless pressure for constant rate production .

Fluid capacitance coefficient of fracture .

Interporosity flow factor from the matrix system to the fracture system . (2)For a triple-porosity reservoir

Dimensionless radial distance .

Dimensionless radial distance of external boundary .

Dimensionless production time .

Dimensionless pressure for constant rate production .

Fluid capacitance coefficient of vug .

Fluid capacitance coefficient of matrix .

Interporosity flow factor from the vug system to the fracture system .

Interporosity flow factor from the matrix system to the fracture system.

A.3. Solution of the Mathematical Model

The triple-porosity reservoir model is taken as an example to solve the mathematical model [13, 16, 20, 40, 41].

Laplace transform is used to solve the established model and defined as

Taking the Laplace transform into Equation (3) to Equation (8), one can obtain

Initial condition:

Inner boundary condition:

Exterior boundary condition:

Through Equations (A.11) and (A.12), one can obtain

Substituting Equations (A.16) and (A.17) into Equation (A.10), one can obtain

in which

Equation (A.18) is the zero-order approximate form of the dimensionless mathematical model of a well in a triple-porosity reservoir, and Equation (9) is the solution of the mathematical model.

Nomenclature

:Oil volume factor, dimensionless
:Total compressibility of fracture, MPa-1
:Total compressibility of matrix, MPa-1
:Total compressibility of vug, MPa-1
:Total compressibility, MPa-1
:Formation thickness, m
:Modified Bessel function of the first kind, zero order, dimensionless
:Modified Bessel function of the first kind, first order, dimensionless
:Modified Bessel function of the second kind, zero order, dimensionless
:Modified Bessel function of the second kind, first order, dimensionless
:Permeability of fracture, mD
:Permeability of matrix, mD
:Permeability of vug, mD
:Pressure, MPa
:Pressure of fracture, MPa
:Dimensionless pressure of fracture system
:Initial formation pressure, MPa
:Pressure of matrix, MPa
:Dimensionless pressure of matrix system
:Pressure of vug, MPa
:Dimensionless pressure of vug system
:Flow rate, m3/d
:Radial distance in reservoir formation, m
:Dimensionless radial distance
:Radial distance of side external boundary, m
:Dimensionless radial distance of side external boundary
:Wellbore radial, m
:Effective wellbore radial, m
:Skin factor, dimensionless
:Production time, hours
:Dimensionless production time
:Laplace transform variable, dimensionless.
Greeks Symbols
:Geometric shape factor of the matrix block, m-2
:Geometric shape factor of the vug block, m-2
:Interporosity flow factor of matrix system to fracture system, dimensionless
:Interporosity flow factor of vug system to fracture system, dimensionless
:Fluid viscosity, mPa·s
:Porosity, fraction
:Porosity of fracture system, fraction
:Porosity of matrix system, fraction
:Porosity of vug system, fraction
:Fluid capacitance coefficient of fracture system, dimensionless
:Fluid capacitance coefficient of matrix system, dimensionless
:Fluid capacitance coefficient of vug system, dimensionless.
Superscripts
:Laplace transform.
Subscripts
D:Dimensionless
m:Matrix
e:External
t:Total
f:Fracture
v:vug
i:Initial
w:Wellbore.

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 they have no conflicts of interest.

Acknowledgments

The authors acknowledge that this study was partially funded by the National Major Science and Technology Projects of China (Grant No. 2017ZX05030002003). We also thank the National Natural Science Foundation of China (Grant Nos. U1762210, 51574258, and 41672132) for financial support.