#### Abstract

To further investigate the rheological consolidation mechanism of soft soil ground with vertical drains, the fractional-derivative Merchant model (FDMM) is introduced to describe the viscoelastic behavior of saturated clay around the vertical drains, and the flow model with the non-Newtonian index is employed to describe the non-Darcian flow in the process of rheological consolidation. Accordingly, the governing partial differential equation of the ideal sand-drained ground with coupled radial-vertical flow is obtained under the assumption that the vertical strains develop freely. Then, the numerical solution to the consolidation system is conducted using the implicit finite difference method. The validity of this method is verified by comparing the results of Barron’s consolidation theory. Furthermore, the effects of the parameters of non-Darcian flow and FDMM on the rheological consolidation of ground with vertical drains are illustrated and discussed.

#### 1. Introduction

In the southeastern coastal and inland areas of China, soft clay foundations exist widely. Consolidation with vertical drains has been proved to be an effective method for the reinforcement of such soft soil foundations and has always been an important research topic in the field of soil mechanics. In 1948, Barron [1] established the basic differential equation for axisymmetric consolidation and derived the analytical solution under the extreme assumptions of equal strain and free strain, which has been widely applied in practice to predict the consolidation behavior. Yoshikuni and Nakanodo [2] perfected the definition and assumption of free strain, introduced the flow continuity condition of well column, and deduced the theoretical solution of sand-drained ground considering the well resistance under the condition of free strain. Hansbo [3] and Xie and Zeng [4] developed Barron’s theory, respectively, deduced the analytical solutions that satisfy the basic equation of radial consolidation and definite conditions, which further perfected the analytical theory of consolidation with sand drains under the condition of equal strain. Since then, a group of researchers [5–12] have conducted a series of studies on the consolidation problems of sand-drained ground and, respectively, provided the consolidation solutions of sand-drained ground under different conditions, which took into account the influence of variable load, partially penetrated vertical drains, stratification of the ground, well resistance, smearing effect, and other factors. However, in those above studies, the deformation of the soil skeleton was modeled as a linear elastic material, and the rheological properties of soft clay were neglected. Therefore, there sometimes exists notable difference between the theoretical predictions of consolidation and the measured values.

In order to consider the rheological properties of the soil layer, many researchers use various component rheological models, which are the combinations of Hooke’s elastic element (spring) and Newton’s viscous element (dashpot), to simulate the deformation of soil skeleton structure and establish the corresponding rheological consolidation theory considering the viscoelasticity of soft soil, but most of them have been focused on the one-dimensional (1D) consolidation [13–17]. For the consolidation with vertical drains, it can be traced back to the 1960s in China, Qian et al. [18], extended the Barron’s theory to the viscoelastic soil based on the standard Merchant model (SMM). Subsequently, Zhao [19] used the generalized Voigt model to derive the widely applicable analytical solution of sand-drained ground which considered the viscoelasticity of the soil. Liu et al. [20] gave the analytical solutions of viscoelastic theory under the condition of free strain and equal strain based on the SMM and established a relatively perfect viscoelastic consolidation theory considering the influence of well resistance and smearing. Wang and Xie [21], Liu et al. [22], and other researchers also considered the influences of different loading forms, unpenetrated soil layers, and semipermeable boundaries based on the SMM, and further revealed the influence of viscoelastic parameters on the consolidation mechanism. Nevertheless, in the past few decades, many authors [23, 24] pointed out that the fractional-derivative model (FDM) is more suitable for describing the viscoelastic behavior of various real materials. Gemant [25] verified the applicability of the fractional-derivative constitutive model to viscoelastic soil through experiments, which provided the basis for the development of fractional order in geotechnical engineering. Subsequently, some researchers [26–31] began to use the spring-pot element (or Abel’s dashpot) in which the first derivative of strain with respect to time is replaced by the fractional-derivative, the fractional-derivative kelvin model or Merchant model, in which the dashpot element is replaced by the spring-pot element, to fit the rheological experimental data of different soft soils. Liu and Yang [32], Liu et al. [33], and Wang et al. [34, 35] also introduced the fractional-derivative kelvin model or Merchant model to study the 1D consolidation of viscoelastic soil. Zhang et al. [36] put forward the fractional-derivative merchant model (FDMM) to study the interaction between rectangular plates and viscoelastic soil layers. However, to the author’s knowledge, the consolidation of viscoelastic soil with vertical drains has been seldom simulated using a fractional-derivative model.

In addition, many permeability tests [37, 38] indicate that water flow in soft soil may deviate from Darcy’s law under low hydraulic gradients. Hansbo [37] established a segmented function model to describe non-Darcian flow. Liu and Jiao [39] also introduced this flow model to modify the traditional consolidation theory with sand drains. However, Hansbo’s model is represented by a piecewise function whose cut-off point is not well determined, and its parameters are relatively more. It is believed that the behavior of non-Newtonian fluid of pore water in soft clay is the main factor for the occurrence of non-Darcian phenomenon. Therefore, Swartzendruber [40] proposed a continuous exponential flow equation describing non-Darcy law, and the expression iswhere is the flow velocity; is the hydraulic gradient; and are the slope and intercept of the asymptote of the -*i* curve, respectively, and is similar to the permeability coefficient of Darcy’s law; and is the non-Newtonian index. If becomes zero, non-Darcian flow equation described by equation (1) will degenerate into Darcy’s law.

The mathematical description of this flow model with non-Newtonian index is relatively simple, and it has fewer parameters and easy to determined. Moreover, it has a good fitting with the experimental data of Hansbo [37]. Based on equation (1), Li et al. [41] preliminarily analyzed the influence of its parameters on the 1D consolidation behavior of soft soil layer.

In order to investigate the applicability of the fractional viscoelastic model in the consolidation theory of sand-drained ground, the FDMM based on the fractional derivative defined by Caputo [42] is introduced to describe the rheological properties of saturated clay. And the flow model with the non-Newtonian index, i.e., equation (1), is considered to describe the relationship between flow velocity and hydraulic gradient. Then, under the assumption of free vertical strain, the partial differential governing equation of the ideal sand-drained ground is derived, and its numerical solution is obtained by the implicit finite difference method. At last, the influences of the parameters of non-Darcian flow and FDMM on the rheological consolidation process are investigated.

#### 2. Fractional-Derivative Merchant Model

A spring-pot element (see Figure 1) is defined by Koeller [43], and it describes the stress-strain relationship as follows:where *E* is the elastic modulus; is the viscous time of the viscoelastic body; *F* is the viscosity coefficient; and represents the fractional differentiation of Caputo’s definition, which is defined aswhere is the fractional order, *t* is the elapsed time after loading, and is the Gamma function. If or , the spring-pot element is, respectively, degenerated into a spring element with elastic modulus *E* or a dashpot element with viscosity coefficient *F*. Therefore, if takes a value between 0 and 1, the spring-pot element will represent the mechanical behavior of the viscoelastic materials.

The FDMM as shown in Figure 2 is obtained by using the spring-pot element to replace the dashpot element of the Merchant model, and its creep compliance iswhere is a single parameter Mittag-Leffler function. If and , then equation (4) degenerates into , i.e., the creep compliance expression of the integer-order or classic Merchant model.

Then the stress-strain equation of fractional viscoelastic model with small strain can be expressed as follows:wherein represents the effective stress.

#### 3. Analysis Model

The problem studied herein is shown in Figure 3. A homogeneous layer with thickness *H*, and with a pervious upper boundary and an impervious lower boundary, has been completely consolidated under its self-weight. Vertical drains, set evenly in the soil layer, are usually modeled as an equivalent axisymmetric system as shown in Figure 3(b), in which and are the radii of the vertical drain and the effect zone, respectively. Then, for the ease of mathematical derivation, the assumptions made in this study are parallel to those adopted in the theory of consolidation by vertical drains for free strain condition [1, 2, 20]. And the prominent assumptions are listed as follows:(i)The vertical loads are initially carried by the excess pore water pressure *u* [1, 2, 20].(ii)The deformation of each point in the layer is completely free, and the load does not change its distribution due to differential settlement on the ground [1, 2].(iii)The compressive strains within the soil mass only occur in a vertical direction, ignoring the deformation of the soil in the horizontal direction [1, 2].(iv)The uniformly distributed vertical load is applied instantaneously, and the stress-strain relationship of the soil is simulated by FDMM. The flow of pore water in the consolidation process can be described by equation (1), and the parameters of FDMM and flow model are constant during the consolidation process.

**(a)**

**(b)**

At the point with radius *r* and depth *z* and at the time *t*, the pore water pressure is *u*, the volume strain is , and the vertical and radial components of the flow velocity are and , respectively. The vertical and radial components of the hydraulic gradient are and , respectively, and , , .

The continuous condition of volume change can be expressed as

Substituting equations (1) and (5) into equation (6) yields

The initial and boundary conditions of the problem are as follows:

For the convenience, the dimensionless parameters are introduced as follows:

Then, equations (7)–(10) can be turned intowhere , , and .

#### 4. Numerical Solution of Consolidation Governing Equation

The implicit finite difference method is introduced to obtain the numerical solution of equation (11) under the conditions of equations (12) and (13). That is, the soil layer is uniformly divided into *N* layers by step in the radial direction within the range of and dispersed into *M* layers by step in the vertical direction within the range of . Meanwhile, an increment along the time axis is taken as . Therefore, the equation (11) can be discreted as

Thus,wherewherein the value of is and the value of *j* is .

The initial condition of equation (12) and the boundary condition of equation (13) can be discretized as

For the impervious boundary at the edge of the influential zone of the sand drain, set a symmetric virtual node *N* + 1 after node *N*, and let to get

For the impervious boundary at the bottom of the sand-drained ground, set a symmetric virtual node *M* + 1 after node *M*, and let to get

For the corner boundary node (*N*, *M*), let , , and the following equation can be obtained:

Thus, equations (15)–(21) compose a closed system of equations, which can be obtained by the iterative method. After calculating the dimensionless excess pore water pressure *U*, the overall average degree of consolidation *U*_{p} defined by the dimensionless excess pore water pressure *U*, which represents the dissipation rate of the overall pore water pressure in soil layer, can be obtained:

Assume that *U* is linear within the interval along *R*-direction and within the interval along *Z*-direction, thus

By substituting equation (23) into equation (24), the overall average degree of consolidation at the time can be expressed as follows:where .

#### 5. Solution Validation

The FDMM degenerates into a linear elastic model if and , and if , the non-Darcian flow described by equation (1) does into the Darcy. Thus, if , , and , then the present problem will be reduced to the concerns of Barron’s axisymmetric consolidation theory with ideal sand wells. Taking , , , the iterative error of equation (15) is 10^{−8}, then the numerical solution of this equation is obtained by MATLAB programming for the special cases of and , respectively. The overall average degree of consolidation *U*_{p} obtained is shown in Figure 4. Meanwhile, for the sake of comparison, the axisymmetric consolidation solution of the ideal sand wells by Barron under the assumptions of free strain and equal strain are also shown in the figure. It can be seen that, at the initial stage of consolidation, the numerical solution in this study is consistent with the analytical solution by Barron’s theory under the assumption of free strain and is slightly larger than the consolidation solution that is under the assumption of equal strain. However, with the development of time or the increase of drain spacing ratio *n*, the results under the two extreme assumptions tend to be consistent. This shows that the proposed numerical solution in this study is effective.

#### 6. Influence of Non-Newtonian Index *I*_{0}

Firstly, the influence of *I*_{0} on the dissipation of dimensionless pore water pressure *U* in sand-drained ground is investigated. Figure 5 shows the influence of non-Newtonian index *I*_{0} on the distribution of pore water pressure *U* along the radial and vertical direction, where , , , , , and *I*_{0} varies from 0 to 5. Similar to the theory of 1D consolidation [41], at a given time, the pore water pressure *U* calculated based on the non-Darcian flow described by non-Newtonian index is greater than that on Darcian flow (), and the difference becomes more obvious with the increase of *I*_{0}, i.e., non-Darcian flow described by non-Newtonian index delays the dissipation of pore water pressure compared with Darcian flow. Further analysis of Figure 5(a) shows that, in the effect zone of the sand drain, the retarding effect of non-Darcian flow on the dissipation of pore water pressure *U* is related to the distance *R* from the center of the sand drain, that is, the phenomenon becomes more obvious with the increase of the distance from the sand well. It can be observed in Figure 5(b), at the same time, the distribution of the pore water pressure along the vertical direction also increases with the increase of distance from the drainage surface. In contrast, unlike the distribution of pore water pressure *U* along the radial direction, when a “certain distance” is exceeded, the distribution of pore water pressure in the vertical direction remains basically unchanged. And with the increase of *I*_{0}, this “certain distance” decreases relatively. The reason is that the vertical drainage path is long in the sand-drained ground, and the influence of dissipation of the pore water pressure *U* is only significant in the areas which are close to the drainage surface. That is to say, the vertical flow of pore water in the consolidation process can be ignored when the sand well is longer, and the more the flow of water deviates from Darcy’s law, the more significant the effect of radial drainage.

**(a)**

**(b)**

The above conclusion can be seen more clearly from the contour diagram of the pore pressure of the same time () in the sand-drained ground with Darcian flow () and non-Darcian flow (), in Figure 6, respectively. The contours of the pore water pressure far from the upper drainage boundary are almost linearly parallel to the *Z*-axis, which indicates that the vertical component of the hydraulic gradient is small in this area, so the vertical component of the corresponding flow velocity in this area is also small. The transverse comparison of the two figures also shows that the straight-line segments of the contours under the condition of the non-Darcian flow are longer than that of Darcian flow. That is to say, when the flow model described by non-Newtonian index is considered, the influential distance of the pore water pressure along the vertical direction will be relatively shorter. It indicates that for the layer with the non-Darcian flow, the water flow in the radial direction should be paid more attention to. It is also consistent with Liu and Jiao’s conclusion [39] considering the influence of Hansbo’s flow on the consolidation of sand-drained ground.

**(a)**

**(b)**

In order to further investigate the influence of non-Newtonian index *I*_{0} on the degree of consolidation *U*_{p}, Figure 7 shows the curves of the degree of consolidation *U*_{p} versus the dimensionless time *T* under different values of *I*_{0}. In the Figure, all of the average consolidation curves of non-Darician flow are below the corresponding curve of Darcian flow (). And with the increase of *I*_{0}, the overall average degree of consolidation *U*_{p} at same time becomes less, or it needs longer time to obtain the same value of the overall average degree of consolidation. For instance, corresponding to the overall average degree of consolidation equal to 90%, the dimensionless time *T* calculated based on Darcian flow is 0.44, while the time based on the non-Darcian flow with *I*_{0} = 0.5, 1.0, 3.0, and 5.0 is about 1.25, 2.15, 5.76, and 9.35 respectively, which is 2.84, 4.89, 13.10, and 21.25 times as much as that under the condition of Darcian flow. Therefore, it can be seen that the influence of *I*_{0} on consolidation behavior of sand-drained ground is remarkable, which also illustrates the importance of considering the influence of non-Darcian flow in practical engineering applications. Moreover, in order to compare the influence of radial seepage, Figure 7 also shows the influence curves of the degree of radial consolidation *U*_{R} with *I*_{0} under the same parameter, which are represented by a dashed line in the diagram. It can be seen that at the same time, the overall average degree of consolidation *U*_{p} with coupled radial-vertical flow is slightly larger than *U*_{R} considering only radial flow, and the absolute difference (*U*_{p} − *U*_{R}) between the two is small at the early and final stages of consolidation but larger in the middle stage. The relative difference (*U*_{p} − *U*_{R})/*U*_{p} is larger in the initial stage of consolidation, increases slightly with the increase of *I*_{0}, and decreases gradually with time. That is, the relative difference between the two is significantly obvious at the early stage of consolidation. For example, when *T* = 0.05, the distribution range of relative difference for *I*_{0} = 0∼5 is about 17.6%∼24.6%, while for *T* = 1.0 and 3.0, the distribution range is 2%∼12% and below 3.9%, respectively. Therefore, for the thicker sand-drained ground, the vertical flow can be neglected in the middle and late phase of consolidation, and the influence of radial flow should be mainly considered.

#### 7. Influence of the Parameters of FDMM

##### 7.1. For the Fraction Order *α*

Figure 8 shows the curves of pore water pressure *U* at radial impervious surface and overall average consolidation degree *U*_{p} with the dimensionless time *T*, where , , , , , and *α* varies from 0.1 to 1. Obviously, the four curves all intersect each other, which indicates that the influence of fractional-order *α* on different stages of sand-drained ground is different. As a whole, at the initial stage (*T* is less than 0.10), the larger the fractional order *α*, the smaller the pore pressure *U* at the radial impervious boundary, and the greater the overall average degree of consolidation *U*_{p}, which means the overall dissipation of pore water pressure in the sand-drained ground becomes faster with the increase of *α*. However, as time goes on, the situation in the middle stage of consolidation is just the opposite. This is similar to Liu’s 1D consolidation analysis based on the fractional-derivative kelvin or Merchant model [32, 33]. When *T* is greater than 5, the four curves intersect each other again and almost overlap. At this stage, the degree of consolidation has exceeded 95%, which indicates that the value of *α* has little influence on the later consolidation process of sand-drained ground. Under this scenario, its influence is relatively negligible from the perspective of engineering significance.

**(a)**

**(b)**

##### 7.2. For the Dimensionless Coefficient *V*

As defined in equation (10), , *V* can be regarded as the dimensionless coefficient of viscosity. In order to investigate the influence of viscosity, Figure 9 shows the variation curves of pore water pressure *U* at the radial impervious surface () and overall average degree of consolidation *U*_{p} with time for the cases that , , , , and *V* is 0.001, 0.01, 0.1, and 1 respectively. By comparing Figures 9(a) and 9(b), it can be seen that, at the initial stage of rheological consolidation, with the increase of *V*, the pore water pressure *U* at the impervious surface becomes less, and the overall average degree of consolidation *U*_{p} becomes greater, which means the viscosity of soil speeds up the dissipation of pore water pressure at this stage. As the development of time, approximately at , the four curves intersect each other. Subsequently, the overall dissipation rate of the pore pressure decreases with the increase of *V*. This indicates that at the later stage of consolidation, the viscosity of soil delays the overall dissipation of pore pressure in sand-drained ground. But at this stage, the four curves are almost coincident, i.e., the influence of *V* is relatively negligible. This is also consistent with Liu’s [32, 33] 1D consolidation analyses based on the fractional-derivative model.

**(a)**

**(b)**

##### 7.3. For the Elastic Modulus Ratio *S*

is the modulus ratio of the independent spring and the series spring of the fractional Kelvin body. In order to investigate the influence of *S*, the curves of pore pressure *U* at the radial impervious surface () and *U*_{p} versus time factor *T* for the cases that , , , , and *S* varies from 0 to 10 are shown in Figure 10(b). In Figure 10, is equivalent to disregarding the rheological effect of the soil, because the FDMM is reduced to a linear elastic model if . It can be seen from Figure 10(a) that, with the increase of *S*, the time when the pore water pressure *U* at the point begins to dissipate is delayed, and the dissipation of pore pressure becomes also slower. For example, the dimensionless time *T* when the pore water pressure *U* at the point begins to dissipate for (linear elastic model) is approximately equal to 0.03. For *S* = 0.5, 1, 5, and 10, the corresponding time *T* is 0.05, 0.06, 0.13, and 0.23, respectively. In Figure 10, corresponding to *S* = 0∼10, the time *T* corresponding to the same overall average degree of consolidation *U*_{p} equal to 30% are 0.04, 0.06, 0.07, 0.16 and 0.30, respectively. Obviously, with the increase of *S*, the time to obtain the same degree of consolidation *U*_{p} is also significantly prolonged. This is also similar to the 1D consolidation conclusion based on the traditional integer-derivative model by Li et al. [16] or Liu et al. [17]. That is, the consolidation process of the soft soil ground will slow down with the decrease of the elastic modulus *E*_{1} of the spring in the series Kelvin body.

**(a)**

**(b)**

Compared with Section 6, it can be seen that the influences of non-Newton index *I*_{0} and modulus ratio *S* on pore pressure *U* and overall consolidation degree *U*_{p} shows a similar rule. However, due to their different physical meanings, their influence mechanism is different. The non-Newton index *I*_{0} describes the strength of the non-Darcian characteristics of flow in the foundation. Although the modulus ratio *S* reflects the influence of elastic modulus *E*_{1} of fractional kelvin body but indirectly reflects the viscosity of strength of FDMM. Because the larger *S* is, the smaller *E*_{1} is, and the viscosity of fractional kelvin body is the stronger when the other conditions are the same. The non-Darcian characteristics of flow and the viscosity of soil layer are the main reasons for the dissipation of pore pressure and the slowing of consolidation rate of the foundation.

#### 8. Conclusions

The non-Darcian flow model with the non-Newtonian index is introduced to describe the flow of water in the consolidation process of sand-drained ground, and the FDMM is done to describe the viscoelastic behavior of the soil around the vertical drains. Then, the axisymmetric rheological consolidation equation of the ideal sand-drained ground is established based on the assumption of free strain, and the numerical solution of the equation is obtained by the implicit finite difference method. Through the analysis of the relevant parameters, the following conclusions are drawn:(1)As the non-Newtonian index *I*_{0} increases, the pore pressure dissipation rate in saturated clay layer with sand drains will decrease, i.e., the consolidation of sand-drained ground will slow down when considering the non-Darcy characteristics of flow.(2)Whether the water flow in the pore obeys Darcy’s law or not, the calculated consolidation rate considering the coupled radial-vertical flow is relatively greater than that considering the radial flow only. However, for the thicker soil layer, the vertical flow can be ignored relatively in the middle and late stages of consolidation, and the consolidation of the soil layer depends mainly on the radial flow.(3)The influence of the parameters of FDMM on the consolidation process is different. Among them, the influence of modulus ratio *S* on the consolidation process is more significant. With the increase of *S*, the overall dissipation rate of pore pressure in soil becomes less, and its effect is mainly reflected in the middle stage of the consolidation process. However, the effects of fractional order *α* and the dimensionless viscosity coefficient *V* on the consolidation rate and the dissipation of pore water pressure in sand-drained ground are more complicated, i.e., with the increase of fractional order *α*, the dissipation of the excess pore water pressure will accelerate in the early stage of consolidation, but will slow down in the middle and later stages of the consolidation process. The influence of viscosity coefficient *V* is greatly similar to that of *α*, but with the increase of *V*, the dissipation of the excess pore water pressure will accelerate for a longer time.(4)In general, the solution in this paper is sensitive to the parameters of non-Newtons index *I*_{0} and modulus ratio *S* and is relatively insensitive to the parameters of fractional order *α* and dimensionless viscosity coefficient *V*. In this way, in practical application, the insensitive parameters can be selected first, and only the values of non-Newtonian index *I*_{0} and modulus ratio *S* can be paid attention to.

#### Nomenclature

C_{v}: | Coefficient of consolidation, |

E: | Elastic modulus |

E_{0}: | Modulus of independent spring |

E_{1}: | Modulus of a spring of the Kelvin body |

E_{α}(x): | Mittag-Leffler function |

F: | Viscosity coefficient |

H: | Total thickness of clayey layer |

i: | Hydraulic gradient |

i_{0}: | Non-Newtonian index |

i_{r}, i_{z}: | Radial and vertical components of hydraulic gradient, respectively |

i′, j, k: | Node numbers in the radial direction, vertical direction, and time axis, respectively |

I_{0}, R, Z, V, β: | Dimensionless factors of i_{0}, r, z, F and r_{e}, respectively, as shown in equation (10) |

J: | Creep compliance |

k_{u}: | Permeability coefficient |

M, N: | Number of soil layers in the vertical direction and radial direction, respectively |

m_{v}: | Volume compressibility of the soil layer |

n: | Drain spacing ratio = r_{e}/r_{w} |

P: | Distributed load |

r, z: | Radial and vertical coordinates, respectively, as shown in Figure 3 |

r_{w}, r_{e}: | Radii of the vertical drain and effect zone, respectively, as shown in Figure 3 |

S: | Ratio of modulus |

t: | Time |

T: | Dimensionless time corresponding to time t, , as shown in equation (10) |

u: | Excess pore water pressure |

U: | Dimensionless excess pore water pressure corresponding to excess pore water pressure u, as shown in equation (10) |

U_{p}: | Overall average degree of consolidation with coupled radial-vertical flow |

U_{R}: | Radial average degree of consolidation considering only radial flow |

: | Flow velocity |

, : | Radial and vertical components of flow velocity, respectively |

α: | Fractional order |

η: | Viscous time |

γ_{w}: | Unit weight of water |

ε, ε_{v}: | Vertical strain and volume strain, respectively |

σ′, σ: | Vertical effective stress and total vertical stress, respectively. |

#### Data Availability

All data generated during the study are available from the corresponding author by request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 51578511); their generous support is gratefully acknowledged.