Abstract
Salt precipitation is generated near the injection well when dry supercritical carbon dioxide (scCO2) is injected into saline aquifers, and it can seriously impair the CO2 injectivity of the well. We used solid saturation () to map CO2 injectivity. was used as the response variable for the sensitivity analysis, and the input variables included the CO2 injection rate (), salinity of the aquifer (), empirical parameter , air entry pressure (), maximum capillary pressure (), and liquid residual saturation ( and ). Global sensitivity analysis methods, namely, the Morris method and Sobol method, were used. A significant increase in was observed near the injection well, and the results of the two methods were similar: had the greatest effect on ; the effect of and on was negligible. On the other hand, with these two methods, had various effects on : had a large effect on in the Morris method, but it had little effect on in the Sobol method. We also found that a low had a profound effect on but that a high had almost no effect on the value.
1. Introduction
CO2 emissions from human activities have exacerbated global warming. In recent years, scientists have tried to use the carbon capture and storage (CCS) method to mitigate the greenhouse effect. In comparison to the other potential storage types (including deep saline aquifers, depleted oil reservoirs, and unminable coal seams), deep saline aquifers are widely distributed and have a high storage potential; thus, they are considered to be effective storage sites for CO2 geological storage [1]. A schematic diagram of CCS is presented in Figure 1. Extra pressure buildup near the injection well, which indicates the production of salt precipitation, is monitored in CO2 geological storage site projects [2, 3]. Through laboratory tests, Bacci et al. found that salt precipitation changed a rock's porosity (), that even small variations in will cause considerable permeability () fluctuations, and that a decrease in permeability can result in an increase in injection pressure and a decrease in injectivity.
Due to the lack of field test data, the analysis of parameters affecting is always obtained from laboratory experiments and numerical simulation. Ott et al. conducted a core scale drainage test to compare the salt precipitation of the core with two different pore structures, they found that the porous rock is more prone to form salt precipitation than single-pore rock, and they believe that the impairment of injectivity depends on the mobility of the brine phase, which was based on several core flood experiments [4]. By comparing different CO2 injection rates (), Peysson et al. and Ott et al. found that capillary backflow is almost negligible at higher CO2 injection rates and that the appearance of salt precipitation is very limited [5–7]. Tang et al. used brine with different salinities () to carry out a CO2 flood experiment, and they found that has a significant impact on injectivity loss [8]. Pruess and Müller used numerical simulation to perform sensitivity studies on precipitation, they found that capillary pressure () causes an increase in , and they suggested that the effect of the can be alleviated by injecting pure water in advance [9]. Guyant et al. stress that a high permeability reservoir under a low CO2 injection rate has the most salt precipitation [10]. Wang and Liu used TOUGH2/ECO2N to analyze the influence of , and they found that will increase when decreases, increases, or increases [11, 12]. The results of the above studies are focused on local sensitivity analysis of a single factor and study with a single factor with specific changes; however, the global sensitivity evaluation of each factor has not yet been studied.
A global sensitivity analysis method can test the interaction between parameters, and it can test the effect of multiple parameters, which change concurrently on the response variable. The conventional approach to performing global sensitivity analysis is the Morris sensitivity test method, the Fourier amplitude sensitivity test (FAST), and the Sobol sensitivity test method [13–15]. Global sensitivity analysis has been widely used in hydrological models and design models [16, 17]. Global sensitivity analysis requires a large number of model calculations, and the salt precipitation model for geological storage is very complex; thus, it is very difficult to calculate the global sensitivity of salt precipitation. Zheng et al. performed a sensitivity analysis of the CO2 transport process in deep saline aquifers and compared the results of Morris, Sobol, and other global sensitivity analysis methods, and they found that the sensitivity orders of the input variables are different for different response variables and that the computational burden of the Sobol method is too large [18, 19]. Wainwright et al. performed a sensitivity analysis using iTOUGH2 for /brine migration and pressure buildup. Global sensitivity analysis of CO2 sequestration always focuses on the distribution of CO2 and site pressure changes, and these data are used for project monitoring; however, few studies have performed global sensitivity analysis of the phenomenon of salt precipitation [20].
The traditional Sobol method extracts the parameters first, and then, it applies the parameters into a formula or model to calculate the response value. If a model has a clear mathematical formula, it is feasible to extract many samples, such as a polynomial model for the Sobol calculation based on the experimental results [21]. However, physical models are often very complicated, and it is difficult to obtain a clear mathematical formula.
Zheng et al. selected parameters to calculate response variables based on the TOUGH2/ECO2N model, and they found when the number of samples is small, the sensitivity coefficient may be negatively affected by the calculation; thus the number of samples they used was 8129 [19]. The traditional operation of the TOUGH2/ECO2N module consists of four steps: modifying the input parameters, entering the run instructions, calculating, and viewing the output results. Each set of parameters is run once, and the computational complexity of the Sobol sensitivity is too large. A surrogate model is an engineering method used when an outcome of interest cannot be directly measured; thus a model of the outcome is used as an alternative. Zhang and Sahinidis used the surrogate model from the polynomial chaos expansion (PCE) to conduct uncertainty analysis, but they mainly considered the uncertainty between the residual water saturation and the injection rate [22]. Wu et al. proposed the use of Kriging surrogate models for uncertainty and sensitivity analysis in nuclear engineering [23]. Palar and Shimoyama used the Kriging method to solve expensive multiobjective design optimization problem [24]. Zhao et al. use Kriging surrogate model as an optimization approach for identifying the release history of groundwater sources [25]. Hou et al. used this surrogate model to conduct the sensitivity and uncertainty analysis [26], it is more precise and allows more sampling times to be calculated so that the results of the Sobol global sensitivity analysis are more accurate.
In this paper, a simple radial flow model was established, and the numerical simulation of the salt precipitate model was performed by TOUGH2/ECO2N [12]. To analyze the global sensitivity of the salt precipitation, we chose as the response variable of the global sensitivity analysis method, which was one of the results of the numerical simulation. Input variables included the CO2 injection rate (), salinity of the aquifer (), empirical parameter , air entry pressure (), maximum capillary pressure (), liquid residual saturation in the relative permeability function (), and liquid residual saturation in the capillary pressure function (). The global sensitivity analysis was carried out by the Morris method and the Sobol method, which are qualitative and quantitative methods, respectively [13, 15]. The impacts of each parameter and parameter combination on were investigated. The Morris method could quickly select the insensitive parameters, and the number of samples used in the calculation was small. The Sobol sensitivity method was more complicated to calculate, but it calculated the quantitative sensitivity and sensitivity orders of all the variables.
Salt precipitation is always generated near the injection well when dry scCO2 is injected into deep saline aquifers, and it can seriously impair CO2 injectivity. The effect of salt precipitation on the injectivity of a well is represented by solid saturation () near the injection well. It is necessary to analyze the factors that affect the value of near the injection well and then determine the main parameters and the negligible parameters that affect the value. We used two global sensitivity analysis methods, and the main method used in this paper was the Sobol global sensitivity analysis method, which was used to determine the sensitivity orders of all the input variables [15]. The Kriging surrogate model was established to calculate the response variable , and this model replaced the operation of TOUGH2/ECO2N and reduced the simulation time to increase the efficiency of the Sobol sensitivity analysis method [27]. The Monte Carlo algorithm was used to calculate the Sobol global sensitivity values for all input parameters, and the Monte Carlo algorithm reduced the computation load and improved the calculation accuracy [28, 29].
2. Numerical Simulation of Salt Precipitation
The depth of a formation layer suitable for CO2 geological storage is generally more than 800 m; the temperature and the pressure of the layer must exceed the critical value of CO2 (31.1 centigrade and 7.38 MPa, resp.). During the process of supercritical CO2 (scCO2) being injected into deep saline aquifers, the scCO2 is considered a nonwetting phase and the brine is considered a wetting phase. If the nonwetting phase is injected into the wetting phase, the water in the brine evaporates into the dry CO2 stream, and salt precipitation is formed when the brine reaches its solubility limit. Accumulation of salt precipitation near the injection well blocks the CO2 flow path, resulting in a decrease in near the injection well and reduced well injectivity.
Verma and Pruess provide a simple conceptual model to characterize the variation of with near the injection well [30]. These authors use a tubes-in-series model to simulate the natural condition of porous media rock (Figure 2), and they obtain the change in near the injection well with according to the tubes-in-series model:where is the initial permeability and is the fraction of original porosity at which permeability is reduced to zero. Figure 3 shows an example relationship from (1), for parameters of and .
When the nature of rock is determined, the change in rock is mainly controlled by . Therefore, studying the change in near the injection well is equivalent to studying the change in . The following sections describe the process to establish a simple two-dimensional radial flow model to study the global sensitivity of .
2.1. Establishment of a Salt Precipitate Model
Assuming that CO2 was injected into a homogeneous, anisotropic deep saline aquifer at a constant injection rate of 1 kg/s, the thickness of the storage layer was approximately 100 m, the injection well radius was 0.3 m, and the radial distance was 2500 m. The specific radial mesh is given in Table 1. In the longitudinal direction, the model was split into 22 layers. The injection well was 11.5 m from the bottom of the model. The thickness of the top and bottom layers was set to 2 m, the section between the injection well and the upper boundary was evenly divided every 5 m, and the area between the injection well and the lower boundary was also set for every 5 m. Using the MESHMAKER module to establish a radial model, the model is shown in Figure 4, and the meshing is shown in Figure 5 [31].
Horizontal permeability () and vertical permeability () were considered for the horizontal and vertical anisotropy. The liquid capillary pressure () equation used the van Genuchten model [32]:where is the liquid saturation and is the liquid saturation when liquid is saturated and is the residual water saturation. The liquid relative permeability () equation used the van Genuchten-Mualem model [32, 33], and the gas relative permeability () equation used Corey model [34]:where is the residual gas saturation. The parameters of the model are presented in Table 2. The plots of the capillary pressures and relative permeabilities are given in Figures 6 and 7.
2.2. Results of the Numerical Simulation
The ECO2N module was used to simulate the migration of the CO2 in deep saline aquifers and obtain the value of , , and . The time of the numerical simulation was set to 100 days, half a year, and one year. The results of the numerical simulation are given in Figures 8, 9, and 10.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
Figure 8 shows the contour map of near the injection well, and since the vapor content of the gas is small, the gas saturation can be approximated by the saturation of CO2. With the increase in time, CO2 migrates gradually and tends to accumulate upward due to buoyancy. The value of the injection well is approximately 0.9, and it is higher than in other areas. Figure 9 shows the contour map of near the injection well. The value of near the injection well is the largest, which indicates that the salt precipitates accumulate mainly near the well area. With the increase in the injection time, the range of salt precipitation gradually increases, and the direction of the salt precipitate is consistent with the CO2 flow. The values of and near the injection well almost reach the maximum value with the increase in time, and near the injection well almost approaches zero, indicating that the injection well becomes a dry zone. The increase in near the injection well means that the permeable pores are being blocked, which affects CO2 injectivity, causes pressure buildup in the injection well, as shown in Figure 10, and may increase the risk of leakage during the CO2 injection process.
To assess the effect of salt precipitate model parameters on injectivity, of the injection well was selected as the response variable, and the global sensitivity of the model parameters was analyzed.
3. Global Sensitivity Analysis
3.1. Morris Qualitative Global Sensitivity Analysis
Morris proposed a data screening method that can select parameters that have a low impact on the results and reduce the number of analysis variables [13]. In this paper, the Morris method was used to analyze the sensitivity of salt precipitate model parameters, the range of parameters selected for the Morris sensitivity analysis as shown in Table 3, and the selection of additional parameters reference TOUGH2/ECO2N [12]. The basic steps of the Morris method were as follows:(1)Each parameter was assumed to be scaled to take on values in the interval .(2)Each parameter took on values from , where was the number of sampling points in the interval . In this paper, we took 5 sampling points in the interval , and the parameters took on values from .(3)The seven parameters in this paper had values within the respective sampling interval and were constructed in an 8 × 7 order matrix () according to
In matrix , 0 represents the original value of the parameter, 1 represents the value after the parameter was changed, each column represents a parameter, and only one parameter changes between each of the adjacent two rows.(4)The values of the eight were calculated by using each row of parameters in the salt precipitate numerical model.(5)The sensitivity of each parameter was calculated using the following equation:
where is the sensitivity of each parameter, is the difference between two adjacent rows of the matrix , and is the only parameter different between two rows of the matrix .
Because a deviation in one of the matrices appeared in the Morris method, we selected multiple samples and calculated the average value. The number of samples () in this paper was four. The Morris global sensitivity method is expressed as the mean and standard deviation . means the value of the global sensitivity, and means the degree of interaction between parameters.
3.1.1. Results of the Morris Global Sensitivity Analysis Method
The results of the Morris sensitivity analysis are shown in Figure 8. The descending ranking of the sensitivity model parameters is , , , , , , and . The effect of on the is the highest, and the higher the , the more the salt that can be precipitated in the brine. The also had a significant impact on , and the smaller the was, the more the brine stayed near the injection well, and the more the water vapor diffused into the dry CO2 flow, causing more salt precipitation. , , and had a similar sensitivity value, which was approximately half of the value of the and , where the impact of was very small, and the influence of was negligible.
The influence of the interaction between parameters on the response variable can be taken from the value of the ordinate in Figure 11; it can be divided into four levels: the greatest influence (), the second level (, ), the third level (, ), and the lowest influence ().
3.2. Sobol Quantitative Global Sensitivity Analysis Based on the Kriging Surrogate Model
The Sobol method is a quantitative global sensitivity analysis method based on a variance analysis proposed by Sobol in 1990 [15]. The sensitivity of a single input parameter can be evaluated by calculating the contribution of that parameter to the output variance, and the cross sensitivity of multiple input parameters can also be evaluated by calculating the contribution of multiple input parameters to the output variance. The basic steps of the Sobol method were as follows:(1)The model was decomposed by (8) according to the single parameter and the combination parameter.(2)Then, the total variance of the model was decomposed by (9) according to the variance of the single parameter and the variance of the combination parameter.(3)The total variance of the model was calculated by using (10), and the total variance characterized the effect of all parameters on the model output.(4)The variance of the parameter was calculated using (11), and the variance characterized the effect of a single parameter or a combination of parameters on the model output.(5)The variance ratio calculated by (12) was the global sensitivity coefficient and the total sensitivity coefficient of the parameter was given by (13) in the following:where is the first-order sensitivity coefficient, and this coefficient described the contribution of the independent effect of the parameter on the sensitivity. is the second-order sensitivity coefficient, which describes the contribution of the parameter interaction to the sensitivity. is the total sensitivity coefficient, and this coefficient is the sum of every step of sensitivity coefficients of the parameter. In (13), characterizes the effect on variance with all parameters except and their interaction.
If the model has a clear mathematical expression , Sobol sensitivity can be calculated from (12) and (13), but if the model does not have a clear mathematical expression, the sensitivity is calculated by the Monte Carlo method [29, 35]. In this paper, we used the Monte Carlo method to calculate the Sobol sensitivity because it was difficult to establish a definite mathematical expression for . If the number of samples is too low, then the sensitivity coefficient would be negative. The Sobol sensitivity calculation usually has two drawbacks: the input parameters and the running time of the model. If the number of samples is too small, then the sensitivity coefficient would be negative. To reduce the computational load while sampling a large number of samples, a surrogate model was proposed instead of the salt precipitate model in TOUGH2/ECO2N to calculate . The surrogate model was a mathematical model that, by fitting the discrete data, could help to predict the output of the physical model.
3.2.1. Kriging Surrogate Model and the Monte Carlo Method
The Kriging surrogate model is a semiparametric model proposed by Krige in 1951, and it consists of a parametric model and a nonparametric model. The parametric model is a regression analysis model and the nonparametric model is a stochastic distribution model [36]. The main reason to use the Kriging surrogate model for the Sobol sensitivity calculation is to construct a surrogate model based on the existing sample data and then use the surrogate model to calculate the model output. In this paper, the Kriging model was built using the Matlab DACE toolbox [27]. Of the 120 sets of data selected from the numerical simulation results described in Section 2, 100 sets of data were used to establish the surrogate model and 20 sets of data were used to verify the accuracy of the model. The test result is provided in Table 4.
The mean square error (MSE) equaled 0.000109, and it was observed that the solid saturation values fitted with the Kriging model, meeting the accuracy requirements. The coefficient of determination () equaled 0.999861, which was very close to 1, indicating that the Kriging model had a high reference value.
The variable range and the probability distribution of each parameter should be given before using the Monte Carlo method to calculate the Sobol sensitivity. It was considered that the probability distributions of all parameters were a uniform distribution. To reduce the correlation between the samples, sampling parameters were rearranged using the Latin hypercube sampling method and using Matlab to edit the program to reduce the statistical correlation coefficient (Spearman correlation coefficient) [28]. The number of samples () was 8000, and the number of calculations for the model was , where was the number of parameters in the model. The modified Monte Carlo method is presented below [37]:where , superscripts (1) and (2) represent two dimensional Latin hypercube sampling arrays of , is the sampling point of space, and denotes that the th column data of the matrix is taken from array (2) and the other columns data are taken from array (1). Finally, the first-order sensitivity coefficient, the second-order sensitivity coefficient, and the total sensitivity coefficient can be calculated according to (12) and (13).
3.2.2. Results of the Sobol Global Sensitivity Analysis Method
The results of the Sobol sensitivity analysis are shown in Figures 12 and 13. Figure 12 shows the first-order and total sensitivity coefficients of the salt precipitate model parameters. has the greatest influence on , which is much higher than the other parameters. Figure 12(b) shows the sensitivity coefficients for removal, and the sensitivity coefficients and are close to zero, indicating that these two parameters had little effect on the solid saturation. The sensitivity coefficients of each parameter are arranged in ascending order according to Figure 12: , , , , , , and , which is similar to the Morris method except for . Figure 13 shows the contribution of the parameter interaction to the sensitivity, and it can be seen from the figure that the combination of and has the greatest influence on the . In addition to the combination of and , the impact of other combination parameters is small (Figure 13(b)).
(a)
(b)
(a)
(b)
4. Discussion
4.1. Sensitivity Analysis of Parameters in the Salt Precipitate Model
Wang and Liu used a similar radial model to analyze the sensitivity of salt precipitate model parameters; the response variable for the analysis is , and the variable parameter is the capillary pressure equation parameter (, , , and ), where each parameter takes five different values for the model calculation, and a significant change in was observed [11]. This analysis was a representative single-factor local sensitivity analysis, where only the impact of each variable on the response variable was analyzed, but the sensitivity of each factor could not be determined qualitatively or quantitatively. The effects of parameter interaction on the response variable could not be determined (Figure 14).
(a)
(b)
(c)
(d)
The Morris global sensitivity analysis method quickly found the minimum sensitivity parameter (), but the sensitivity values of the other three parameters (, , and ) were approximately the same (Figure 11). The increase in the number of calculations led to a different sensitivity of the arrangement. The Sobol quantitative global sensitivity analysis clearly gave the sensitivity coefficient of these four parameters, from large to small , , , and , where the sensitivity coefficient of and was close to zero (Figure 12). For the interaction between the parameters, both methods showed that the interaction between these four parameters and other parameters had little effect on the response variable. For large-scale and longtime salt precipitate numerical simulations, it was feasible to set and as empirical constants.
4.2. Effect of the CO2 Injection Rate on Solid Saturation
Many researchers have studied the effect of on by numerical simulation; they consider that an increase in decreases the , and their results are shown in Figure 15 [9, 38–40]. The Morris method and Sobol method give a different sensitivity arrangement for . In the Morris method, is one of the main factors causing the change in , which is consistent with the conclusion of Pruess et al. [9, 38–40]. However, in the Sobol method, the influence of the on is much smaller than that of , , and .
To study the effect of on near the well, we chose more values for the numerical simulation (ranging from a low injection rate of 1.5 kg/s to a high injection rate of 100 kg/s). We modified the injection model by enlarging the radius of the model and refined the mesh near the injection well to set up a wide range of . This paper referred to a TOUGH2/ECO2N case to divide the grid and set the model parameters (Tables 5 and 6) [12]. We selected 11 different values for the simulation, and the results are shown in Figure 16.
Figure 16 shows the relationship between different values of and near the injection well. The range of salt precipitation increases with the increase in . The values of at the injection well decrease with the increase in , and at a low , the brine in the unit near the injection well evaporates more completely. In the low range (1.5 kg/s~15 kg/s), it is clear that near the injection well decreases with increasing , which is consistent with the conclusions of Pruess et al. This finding is also consistent with the results of the Morris method ( has considerable influence on ) [9, 38–40]. However, in the high range (15 kg/s~100 kg/s), the value of near the injection well generally does not change. The variation in under high values of was not investigated in the previous sensitivity analysis of salt precipitation [9, 38–40].
It was important to distinguish between at the injection well and near the injection well; the injection well was a boundary condition, and its value could be inaccurate. The response variable for the sensitivity analysis was near the injection well. If the value of was small, then the effect of on near the well was significant. If the value of was large, then the influence of on near the well was negligible. This conclusion also explained the difference between the Morris and Sobol results. In this paper, the range of was 5 kg/s~25 kg/s, the sampling number and the sampling frequency of Morris method were small (4 in this paper), most of the values (75% in this paper) were in the low range, and the effect of on near the well was significant (Figure 11). However, the Sobol method required a large sampling number (8192 in this paper) to calculate the sensitivity coefficient, the selected values were uniform, and the effect of on the was not significant (Figure 12(b)).
5. Conclusions
The radial model was used to simulate the salt precipitate near the injection well during CO2 injection into deep saline aquifers. Numerical simulation results showed that evaporation led to salt precipitate near the injection well. The increase in reduced the injectivity and led to extra pressure buildup near the injection well. In this paper, of the unit near the injection well was selected as the response variable. The Morris global sensitivity analysis and the Sobol global sensitivity analysis based on the Kriging surrogate model were carried out and the main conclusions obtained from this study were as follows:(i)The effect of on was the greatest, and a decrease in could decrease the salt precipitation, such as injection of pure water into the saline aquifer. The influence of , , and on was secondary: was related to the pore distribution of the formation, and were related to the pore size of the formation, and these three parameters could be characterized as the heterogeneity of the formation. and had similar results in sensitivity analysis, and the influence of these two parameters on were negligible. These two parameters did not need to be changed in the numerical simulation of CO2 geological storage.(ii)In the results of the Morris method, and the interaction with other parameters had considerable influence on . However, the results of the Sobol method showed that had little effect on and that the interactions of and , , and had considerable influence on , while the interaction of and other parameters had a low effect on .(iii)The high second-order sensitivity coefficient may have caused the difference between the Morris method and the Sobol method, and the results of the Sobol method, with more sampling numbers, were more reliable than those of Morris method. In this paper, the Morris method and Sobol method had different results related to the sensitivity of , but the second-order sensitivity coefficient of was not large, which indicated that the effect of could not simply be characterized by sensitivity.(iv)There was a critical value () during CO2 injection (15 kg/s in this paper); when was less than , had considerable influence on near the injection well, and when exceeded , its influence could be negligible.
The results of global sensitivity analysis to assess salt precipitation for CO2 geological storage in deep saline aquifers provide possible help for future field text and numerical simulations. Salt precipitation mechanism should be considered in the future work, which could affect the salt precipitation model.
Nomenclature
Salt Precipitate Model Parameters| : | Permeability |
| : | Initial permeability |
| : | Gas relative permeability |
| : | Liquid relative permeability |
| : | Formation pore distribution empirical parameter |
| : | Pressure |
| : | Air entry pressure |
| : | Capillary pressure |
| : | Liquid capillary pressure |
| : | Input absolute value of maximum capillary pressure |
| : | Critical injection rate of |
| : | Injection rate of |
| : | Water saturation |
| : | Solid saturation |
| : | Liquid saturation when liquid is saturated |
| : | Residual water saturation |
| : | Residual gas saturation |
| : | Residual water saturation in relative permeability function |
| : | Residual water saturation in capillary function |
| : | Gas saturation |
| : | Temperature |
| : | Salinity of saline aquifers |
| : | Porosity |
| : | The fraction of original porosity at which permeability is reduced to zero. |
| : | Variance |
| : | The sensitivity of each parameter |
| : | The number of parameters in the model |
| : | Number of samples |
| : | The number of sampling points in the interval |
| : | The number of matrix selections |
| : | Sobol sensitivity coefficient |
| : | Parameter |
| : | The value of the Morris global sensitivity |
| : | The degree of interaction between parameters. |
| : | Coordinate axis direction |
| : | Parameter type. |
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The research was supported by the National Natural Science Foundation of China (no. 51179060), the Education Ministry Foundation of China (no. 20110094130002), the Fundamental Research Funds for the Central Universities (no. 2015B26514), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (no. KYCX17_0470), and the Fundamental Research Funds for the Central Universities (no. 2017B648X14).