Abstract
A governing equation of drawdown was put forward to describe the one-dimensional large-strain consolidation behavior of an aquitard without consideration of the creeping effect. An analytical solution was derived to characterize the drawdown variation in the aquitard subjected to sudden hydraulic head decline in adjacent confined aquifer. The characteristics of the groundwater dynamics and water balance in the aquitard have been analyzed based on the analytical solution. A comparison analysis of results has been made between the large-strain theory and the classical small-strain theory. The type-curve fitting method was used to determine the hydrogeological parameters, on the basis of the observed variations of aquitard deformation with time. The analytical solution was thus validated by a comparison with the observed experimental results. It is found that the water drainage of aquitard is obviously delayed in response to the water head decline in the adjacent aquifer. All delayed water release from the aquitard terminates when the consolidation time reaches the value of . The aquitard deformation predicted by the large-strain theory is less than that given by the small-strain theory, and the prediction discrepancy of these two theories increases with the increasing soil compressibility.
1. Introduction
A multilayered aquifer system usually consists of aquifers of relatively high permeability alternating with aquitards of relatively low permeability in between [1]. Aquitards may extend over large land areas throughout the earth, covering entire basins or alluvial plains, and sometimes they even extend into the sea, such as the Dakota aquifer system in the United States and the Yangtze Delta in China [2–4]. Abundant groundwater stored in the aquitards can be used for water supply when the aquifers are developed [5]. The role of aquitards in the subsurface groundwater system has been increasingly recognized in many respects. The groundwater dynamics in the aquitard is closely related to the groundwater exploitation, groundwater contamination, underground storage utilization, and land subsidence [6, 7].
Aquitards consist of fine sediments such as clay and silt in low-energy depositional environments [8, 9]. Aquitards, whose compressibility is often much larger than aquifers, are capable of releasing substantial quantities of water from the original storage in the process of dissipation of excess pore water pressure [10]. The laboratory measurement and regional groundwater flow modeling have demonstrated that the permeability of aquitard is not zero [8, 11]. The hydraulic conductivity in an aquitard is usually less than 10−8 m/s, and the permeability varies over a large range from 10−15 to 10−23 m2 [12].
The understanding of the groundwater dynamics and water balance is of great importance to the prediction, evaluation, and control of land subsidence, and it is also meaningful to the groundwater resource evaluation and exploitation, as well as the analysis of contaminant and heat transport in the aquifer system [13–15]. Many studies have recently been carried out in terms of the groundwater flow and depletion in the aquitard. A deconvolution method was proposed [16], by which the drawdown in the confining layer can be calculated on the basis of the water head in neighboring aquifers. In particular, this water head may fluctuate in an arbitrary manner, which is not necessarily caused by a prescribed pumping regime [16]. As the water head varies in the confining layer, the water flow through it can be described by Darcy’s law, which gives a direct estimation of water depletion. Konikow and Neuzil [3] presented a simplified method to estimate the groundwater depletion from the confining layers, and the computational errors were generally small. An analytical approach was proposed by Li et al. [5] to evaluate the depletion of groundwater in the confining layers throughout the entire exploitation history and also within limited periods, especially when there are few data on the drawdown history of the pumped aquifer. Alternatively, numerical models can be used to analyze the groundwater depletion in the confining layers by providing abundant information of the aquifer-confining layer system [17, 18]. However, it should be pointed out that the aforementioned studies were mostly done based on Terzaghi’s theory of one-dimensional consolidation for saturated clays.
Since the hydraulic properties of the confining layers of low permeability are usually not well understood, it is difficult to manage effectively the groundwater resources and have an accurate description of contaminant transport in the confined aquifer systems [19]. In general, the hydraulic conductivity and specific storage of an aquitard can be determined by the laboratory and in situ experiments [20–22]. For example, Burbey [23] used a graphical technique to determine the hydrogeological parameters by taking advantages of the time-subsidence data in a pumping test. In terms of this graphic technique, the straight-line portion of semi-log plot of time-compaction data represents the aquitard compaction. A type-curve fitting method was proposed by Zhou et al. [4] to calculate the hydraulic conductivity and specific storativity with the observed flux in response to the sudden decline of water head in the adjacent confined aquifer. Zhuang et al. [19] proposed a type-curve fitting method to estimate the hydraulic conductivity and specific storativity of an aquitard on the basis of the aquitard compaction variation and a piecewise linear function describing the drawdown history of the underlying aquifer. Konikow and Neuzil [3] presented a general relationship between the porosity and the hydraulic conductivity by considering the clay content of the aquitard and a relationship between and specific storativity by considering the consolidation degree. The parameters of aquitard were assumed to be constant during the consolidation process in the aforementioned studies.
Researchers have also investigated the groundwater dynamics and water balance in the small-strain consolidated aquitard [4, 6]. For the large-strain consolidated aquitard (e.g., spatially the newly formed sedimentary soft soil), the variation of hydraulic properties is nonlinear during the consolidation of aquitard [24–26]. The small-strain consolidation theory is not suitable when the deformation of aquitard is larger than 10% of its thickness [27–29]. In this paper, an analytical solution was proposed to account for the drawdown variation in an aquitard subjected to large-strain consolidation, given that the water head in the adjacent confined aquifer suffers from a sudden decline. The characteristics of the groundwater dynamics, water balance, and aquitard deformation were accordingly analyzed. Additionally, a type-curve fitting method was presented for the determination of hydrogeological parameters, and the obtained parameters were validated by the experiments.
2. Theoretical Background
2.1. Mathematical Model and Analytical Solution
A multilayered aquifer-confining layer system was considered in this study. The system consists of an aquitard of infinite lateral extent, which are bounded by two aquifers lying, respectively, above and below the aquitard. The initial thickness of aquitard is . The aquitard is assumed to be homogeneous and saturated and to have a unified thickness (see Figure 1).
In the multilayered aquifer-confining layer system, the scalar equivalent hydraulic conductivity of each aquifer exceeds that of the aquitard by at least two orders of magnitude. The seepage in the aquitard follows the path of the least resistance; that is, water flows downward in the vertical direction and enters the aquifer with a sharp refraction at the interface between the aquitard and the aquifer. Therefore, it is reasonable to consider the water flow in the confining layer to be essentially vertical [16, 30]. The Lagrangian coordinate was used, which was assumed to be positive in vertically downward direction, with the coordinate origin located at the top surface of the aquitard (see Figure 1).
With the soil particles and pore water assumed to be incompressible, the equations for one-dimensional large-strain consolidation of saturated clay in the Lagrangian coordinate system are given as follows [31]: where is the total vertical stress; is the specific gravity of solid particles, and is the unit weight of water; is the void ratio at time ; is the initial void ratio. where is the total pore water pressure; is the hydrostatic pore water pressure; is the long-term departure of pore water pressure from the hydrostatic pore water pressure due to a hydraulic gradient across the soil layer, and is the excess pore water pressure.
The equilibrium of fluid phase requires
And the effective stress principle is given as where is the effective vertical stress.
Darcy’s law is given as where and are the velocities of the pore water and soil particle relative to the datum plane; is the hydraulic conductivity of clay in the vertical direction.
Equation (5) shows that the water flow velocity can be decomposed into two components, including one part arising from a long-term hydraulic gradient across the consolidating soil layer and the other part relating to the excess pore water pressure gradients. The governing equation of continuity of pore water flow is
The governing equation of consolidation is just the combination of (5) and (6).
The nonlinear variations of the compressibility and permeability of soil in large-strain consolidation process are assumed to, respectively, abide by (8) [32] and (9) [26], and soil creeping here is not taken into account. where is the consolidation coefficient of aquitard, is the initial hydraulic conductivity of the aquitard, and and are the coefficients of compressibility of the aquitard at small- and large-strain states, respectively. It should be noted that the small-strain consolidation theory does not consider the variation of hydraulic properties of consolidating aquitard, and it is not suitable when the deformation of aquitard is larger than 10% of its thickness, while the large-strain consolidation theory considers the nonlinear variations of the compressibility and permeability of the aquitard.
The relationship between void ratio and total pore water pressure is given as where is the initial effective vertical stress.
The decrease of excess pore water pressure equals the increase of drawdown when the total stress is constant. Therefore, the relation of void ratio with the drawdown (i.e., the hydraulic head decline) of the aquitard is obtained as where is the specific storativity of the aquitard. Substituting (8), (9), and (11) into (7) gives
The water head is assumed not to vary in the aquifer above the aquitard, and it is assumed to decline by in the confining aquifer underlying the aquitard. Equation (12) is solved at the following initial and boundary conditions.
The analytical solution of drawdown variation with the large-strain consolidation theory is facilitated by the use of the variable transformation and separation (Appendix).
2.2. Groundwater Dynamics in the Aquitard
According to Darcy’s law, the flux per unit horizontal area at the position and time in the aquitard can be expressed as
Substituting the expression of into (17) gives
The cumulative fluxes at the surface and bottom of the aquitard are, respectively, determined to be
2.3. Groundwater Balance of the Aquitard
The bottom of the aquitard is fixed, so the aquitard settlement can be expressed as
Integrating (21) by replacing with (16) yields the settlement at any point in the aquitard
The settlement at the top surface of the aquitard is then given as
Furthermore, the final settlement of the aquitard is identified by letting
The cumulative fluxes at the top and bottom of the aquitard are given in a dimensionless form, respectively, as where denotes the dimensionless time. According to the water balance principle, the water released from the entire thickness range of the aquitard is equal to the difference of the cumulative fluxes per unit horizontal area at the top and bottom of the aquitard.
The water balance results of the aquitard are shown in Figure 2. The curves of and represent the cumulative fluxes at the top and bottom surfaces of the aquitard, respectively. The curve represents the cumulative water release of the aquitard. Three characteristic evolution stages are observed in Figure 2. At the initial stage (i.e., ), the cumulative flux is zero at the cross section . The cumulative flux increases significantly at the cross section , owing to the cumulative water release of the aquitard. At the middle stage (i.e., ), the cumulative water release of the aquitard increases slowly. The cumulative flux increases slightly, and the increase rate tends to be stable at the cross section . Also, the increase rate of cumulative flux tends to be steady at the cross section , indicating that the underlying aquifer is recharged not only by the released water of the aquitard but also by the leakage from the upper aquifer. At the last stage (i.e., ), the cumulative flux at the cross section is larger than the cumulative water release of the aquitard. The main supply source at the cross section is the cumulative flux at the cross section . The curves of and are parallel to each other, indicating that the water release tends to be completed.
2.4. Identification of the Parameters of the Aquitard
To determine the hydrogeological parameters of the aquitard while the water head declines by in the adjacent confined aquifer, we express (23) in a dimensionless form for simplicity
It is assumed that the soil particles and pore water are incompressible. Therefore, the settlement at the top surface of the aquitard equals the cumulative water release of the aquitard.
The logarithmic forms of (23) and (27) are, respectively, expressed as
The standard curve of is presented in Figure 3. The second terms at the right side of (30) and (31) are constant. Therefore, in the logarithmic plot, the test curve of the flux is analogous to the type curve of the dimensionless flux at (see Figure 4 for reference). The curve of is superposed onto the standard curve of , and meanwhile, the axes of the two graphs are kept to be parallel to each other. There is a horizontal difference of and a vertical difference of between these two curves. The matched point is chosen to be the intersection point of the standard curve with the test curve of . The coordinates of matched point , , , and are then substituted into (32), (33), and (34) to determine the hydrogeological parameters.
3. Results and Discussion
The proposed analytical solutions in this paper for calculating the drawdown and depletion of the aquitard are evaluated by solving a large-strain consolidation problem. The drawdown of underlying confined aquifer is 10 m, and the water head of ground surface water (or unconfined aquifer) is assumed to be unchanged. A comparison analysis is also performed between the results obtained by the large-strain and small-strain consolidation theories. The large-strain theory herein refers to the analytical solution developed in this paper, and the small-strain consolidation theory is the Terzaghi’s one-dimensional consolidation theory. The initial thickness of the aquitard is 10 m. The aquitard parameters used in the large-strain consolidation analysis are m/s, , 0.02, 0.03, 0.04 m−1. The corresponding compressibility coefficients are 1, 2, 3, and 4 MPa−1. In the small-strain analysis, the values of specific storativity and hydraulic conductivity are the same as those used in the large-strain consolidation analysis.
The spatial distributions of drawdown in the aquitard with different specific storativity values (, 0.02, 0.03, 0.04 m−1) at years were calculated using the large- and small-strain theories. The results in Figure 5 indicate that the drawdown in the aquitard increases with the decreasing specific storativity. The excess pore water pressure dissipates more quickly at a relatively small specific storativity. It is also observed that the drawdown obtained by the large-strain theory is less than that obtained by the small-strain theory. The dissipation of excess pore water pressure in the large-strain consolidation is slower than that in the small-strain consolidation. The difference of predicted drawdown between the large- and small-strain theories increases with the increasing specific storativity.
The drawdown at different positions (, 3, 5, 7, 9 m) in the aquitard is plotted in Figure 6. The drawdown increases initially and arrives at a stable value eventually. The variation rate of drawdown is large at the beginning and then gradually decreases to zero. It is also interesting to note that the drawdown increases with the increasing value at a given time. The pore water releases more quickly at relatively high value and specific storativity. The drawdown predicated by the large-strain consolidation theory is smaller than that obtained by the small-strain consolidation theory, and the prediction difference between the large- and small-strain consolidation theories is insignificant. The void ratio decreases with the increasing effective stress in the aquitard. Meanwhile, the hydraulic conductivity and specific storativity of the aquitard decrease, which is ascribed to the decrease of void ratio.
The spatial distribution of water flow velocities in the aquitard is determined by (18) at different specific storativity values (, 0.02, 0.03, 0.04 m−1) (see Figure 7). The results suggest that the water flow velocity in the upper part of the aquitard is lower than that in the lower part. The water flow velocity in the upper part is relatively small given a large specific storativity, but the case is opposite for the water flow velocity in the lower part. The water flow velocity predicated by the large-strain theory is smaller than that predicted by the small-strain theory, and the prediction difference increases with the increasing value and specific storativity.
The variation of water flow velocity at different positions in the aquitard is described in Figure 8. The water flow velocity is zero at the beginning in the upper part of the aquitard and then increases to a steady value. For the lower part, the flow velocity increases initially and then decreases to the same steady value. Due to the decreasing hydraulic conductivity and specific storativity of the aquitard considered in the large-strain theory, the pore water flow velocity is lower than that obtained by the small-strain theory, and the predicted water flow velocities are the same when the water flow becomes stable in the aquitard.
The cumulative water release within unit area of the aquitard equals the settlement at the top of the aquitard, which is calculated by (23). As indicated in Figure 9, the depletion of the aquitard increases gradually to a steady value. The cumulative water release in unit area of the aquitard, which is given by the large-strain theory, is, respectively, 0.48, 0.91, 1.30, and 1.65 m3, while is, respectively, 0.01, 0.02, 0.03, and 0.04 m−1. The depletion of the aquitard increases with the compressibility coefficient. In addition, the predicted cumulative water release of the large-strain theory is obviously less than that of the small-strain theory, and the prediction difference of aquitard depletion increases with the increasing specific storativity value.
4. Experimental Verification
4.1. Laboratory Apparatus
A physical model test has been carried out to study the groundwater dynamics, water balance, and deformation characteristics of the aquitard in response to an abrupt water head decline in the adjacent confined aquifer. A schematic description of the testing device is given in Figure 10. It is composed of the main body of consolidation container, deformation monitoring device, flow monitoring device, and water supply device.
This apparatus can measure the deformation and flux of clay layer during testing. The consolidation container is a sealed organic glass cylinder with the height, outer, and inner diameters being 50 cm, 19 cm, and 18 cm, respectively. The soil specimen consists of three parts: the lower sand layer (4.5 cm in thickness), the middle layer filled with soil for testing (17.4 cm in thickness), and the upper sand layer (5.7 cm in thickness). The middle layer representing an aquitard is silty clay, and the basic physical and mechanical properties are given in Table 1.
4.2. Test Procedures
The water head in the clay layer is initially constant. The water head in the lower sand layer declines suddenly by at the beginning and then remains unchanged during testing. We measured the drainage flux at the bottom of the aquifer and the subsidence of clay layer. The testing steps are given as follows. (1)Initial inspection of the testing device. The consolidation container must be ensured to be sealed, and the inlet and outlet pipes must not be obstructed. The testing instrument must be ensured to operate normally(2)Soil sample preparation. The inner wall of consolidation container is cleaned and smeared with a layer of grease or Vaseline; hot water is injected into the consolidation container through the overflow tank and the water table is raised slowly; no bubble is allowed in the filter layer during the water injection process; the consolidation container is successively filled with sand layer, clay layer, and sand layer, and the device to measure the settlement is placed at the top of and below the clay layer after each soil layer is consolidated under its self-weight for a duration of 24 hours(3)Sample saturation. Turn off all the valves except the exhaust valve, and pump out the air in the container with a vacuum pump; increase the negative pressure within the consolidation container gradually at a given pressure sequence (e.g., –0.02, −0.04, −0.06, −0.08, and− 0.1 MPa); the time interval for maintaining each pressure level is 12 hours(4)After the soil layers are saturated, water is added slowly to reach the target elevation; open the valve of inlet pipe, and supply water into the device; mount the water supply tank at the designed height, and fill in it with water; discharge the excess water from the overflow nozzle; adjust the position of dial gauge, fasten it, and record the initial reading; observe the readings of the dial gage until they do not vary. This process generally takes a few days(5)Record the readings of the dial gage during testing. Start the test by opening the valve of outlet pipe; decrease the hydraulic head below the confined soil layer by m at the beginning of test, and then maintain it to be constant; at the same time, record the readings of the dial gage at time sequence of 5 s, 10 s, 30 s, 1 min, 2 min, 4 min, 8 min, and 12 min, and measure the flux from the outlet of overflow tank using an electronic balance (the water discharge is held with a measuring cup on an electronic balance)(6)Observe the readings persistently until they do not vary, and the test is then terminated; the observation duration is 1440 minutes; the flux becomes stable at the end of observation
4.3. Experiment Results and Analysis
The flux per unit area is determined as the flux at the bottom of the clay layer divided by the horizontal cross section area of cylinder (i.e., 254.34 cm2). The values of are listed in Table 2. The deformation of clay layer was also measured by the dial gage. The data of settlement of the clay layer obtained in the test are given in Table 3.
The measured values of are plotted against time in Figure 4. The type curve of is then superposed onto this figure, by keeping the axes of the two graphs parallel to each other (see Figure 4). The matched point is the intersection point of the type curve with the experimental curve of , and its coordinates in the two systems are mm, , min, and . The parameters of the clay layer can be obtained by substituting these coordinates into (32)–(34). The hydraulic conductivity , specific storativity , and hydraulic diffusivity are 4.25 × 10−4 cm/min, 3.23 × 10−3 cm−1, and 0.1315 cm2/min, respectively.
The flux within the observation time at the bottom of the clay layer is obtained by substituting the hydrogeological parameters estimated by the type-curve fitting method into (18). The main advantage of the type-curve fitting method lies in its ability to make full use of all pumping test data and to improve the accuracy of calculation. Figure 11 compares the predicted and measured fluxes. It is seen that predicted by the analytical solution in this study agrees well with the experimental results.
Based on the one-dimensional Terzaghi’s consolidation theory (i.e., the small-strain consolidation theory), the hydrogeological parameters of the aquitard can also be determined using the type-curve fitting method [33]. The hydraulic conductivity, specific storativity, and hydraulic diffusivity are determined to be × 10−4 cm/min, × 10−3 cm−1, and cm2/min. The hydraulic diffusivity estimated by the large-strain theory is the same as that given by the small-strain theory. The hydraulic conductivity and specific storativity determined by the large-strain theory are larger than those determined by the small-strain theory. The values of hydraulic conductivity and specific storativity estimated by the large-strain theory in fact refer to the hydraulic parameters of the clay layer at the beginning and they decrease with the decreasing void ratio. The compressibility coefficient at the bottom of the clay layer in the large-strain theory is a constant (=2.8 × 10−3 cm−1) as estimated by (8). The hydraulic conductivity at the bottom of the clay layer is 3.2 × 10−4 cm/min as estimated by (9), which is close to that obtained by Darcy’s law ( × 10−4 cm/min). The equivalent hydraulic conductivity given by the large-strain theory is 3.68 × 10−4 cm/min [26], which is close to that predicted by the small-strain theory.
5. Conclusion
The governing equation of drawdown was proposed for the description of the one-dimensional large-strain consolidation of saturated and homogeneous clays, with the creeping effect being ignored. An analytical solution accounting for the drawdown variation in a largely deformed aquitard has been derived, which is subjected to abrupt water head decline in the adjacent confined aquifer. Meaningful results have been obtained on the basis of the assumptions in (8) and (9) relating to the compressibility and permeability of soil. The characteristics of the groundwater dynamics and water balance of the aquitard have been examined with analytical solution. A type-curve fitting method is presented for determining the hydrogeological parameters. The main conclusions are given below: (1)When the water head in the confined aquifer declines by , the response of groundwater dynamics in the aquitard is found to be obviously delayed. The dissipation of excess pore water pressure predicted by the large-strain theory is slower than that predicted by the small-strain theory(2)An expression for computing the water flow velocity in the aquitard has been derived. Then, an analytical solution has been proposed to account for the cumulative water release and settlement of the aquitard. A water balance equation has also been established to describe the interrelations among the recharge, discharge, and water release of the aquitard. All delayed drainage of the aquitard is expected to be completed when the consolidation time arrives at . The soil compressibility is one important factor affecting the discrepancy of predicted deformations between the large- and small-strain consolidation theories. This discrepancy is found to increase with the increasing soil compressibility(3)The hydrogeological parameters of the aquitard (i.e., the hydraulic conductivity, specific storativity, and hydraulic diffusivity) are estimated by fitting the curves of analytically predicted and experimentally observed settlement. The hydraulic conductivity and specific storativity determined by the large-strain theory are higher than those given by the small-strain theory. In the large-strain theory, the hydraulic conductivity and specific storativity of the aquitard are considered to decrease with the decreasing void ratio. The hydrogeological parameters estimated by the small-strain theory are constants which are assumed to be on average equivalent to their true values in the consolidation process. The flux at the bottom of the clay layer, which is obtained from the analytical solution, is in good agreement with the experimental observation, implying that the estimated parameters are able to characterize the groundwater flow in the aquitard
The hydrogeological parameters of the aquitard can be estimated by the type-curve fitting method, together with the data of flux per unit horizontal area at the top and bottom of the aquitard. The type-curve fitting method can also be used to measure in situ hydrogeological parameters.
Appendix
Derivation of the Solutions to (12), (13), (14), and (15)
The solution of drawdown variation for the large-strain consolidated aquitard in response to the sudden decline of water head in the adjacent confined aquifer is facilitated by the use of a transform defined as so
Substituting (A.1) into (11), (12), (13), and (14) gives the transformed governing equation and conditions as follows
Because the boundary condition in (A.3d) is nonhomogeneous, it is necessary to simplify it to be
Substituting into in the governing equation, as well as the initial and boundary conditions in ((A.3a)–((A.3d), leads to
The analytical solution of the above mathematical model (i.e., ((A.6a)–((A.6d)) is obtained by the method of separation of variables as
Substituting (A.5) and (A.7) into (A.4), the solution of the mathematical model (i.e., ((A.3a–(A.3d)) is obtained as
Substituting (A.8) into (A.2) gives the form solution of (12), (13), (14), and (15) in the main text as
Notations
| : | Lagrangian coordinate |
| : | Consolidation coefficient of the aquitard at small-strain states |
| : | Consolidation coefficient of the aquitard at large-strain states |
| : | Void ratio |
| : | The initial void ratio |
| : | Specific gravity of solid particles |
| : | Hydraulic conductivity of the clay in the vertical direction |
| : | Initial hydraulic conductivity of the aquitard |
| : | Thickness of the aquitard |
| : | Coefficients of volume compressibility of the aquitard at small-strain states |
| : | Coefficients of volume compressibility of the aquitard at large-strain states |
| : | Excess pore water pressure |
| : | Cumulative flux |
| : | Cumulative fluxes in a dimensionless form |
| : | Flux per unit horizontal area in the aquitard |
| : | Aquitard settlement |
| : | Specific storativity of the aquitard |
| : | Dimensionless settlement |
| : | Drawdown in the aquitard |
| : | Time |
| : | Dimensionless time |
| : | Total pore water pressure |
| : | Long-term departure of pore water pressure from the hydrostatic pore water pressure |
| : | Hydrostatic pore water pressure |
| : | Velocities of soil particle relative to the datum plane |
| : | Velocities of the pore water relative to the datum plane |
| : | Total vertical stress |
| : | Effective vertical stress |
| : | Initial effective vertical stress |
| : | Unit weight of water |
| : | Drawdown in the confined aquifer. |
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 study is financially supported by the National Natural Science Foundation of China (Grant no. 41702253, 41572209) and the National Key Research and Development Program of China (Grant no. 2016YFC0402803).