Abstract

The first part of this paper presents the major drawbacks of the traditional methods for generating joints in Particle Flow Code 2D (PFC2D). Violent oscillations in the postpeak shear stress and shear-induced dilation in the normal direction occur in specimens generated by directly removing bonds in joints and using the discrete fracture network (DFN) method. The specimens generated by the additional wall method can be used to simulate realistic shear mechanical properties in the direct shear test, but it is difficult to achieve a uniform initial stress distribution within the specimen due to the constraint of particle motion. The second part of this paper explores an improved method to generate realistic joints based on the particle grouping technique and the smooth joint model, and the validity of this method is verified by a set of numerical direct shear tests. The numerical results show that the proposed joint generation method can effectively eliminate the oscillation of the postpeak shear stress and shear-induced dilation in the normal direction. In addition, the mechanical behaviours of the rough jointed rock mass correspond well with the theoretical results obtained from Patton’s and Barton’s models. The proposed model can also simulate the asperity degradation of rough jointed rock masses.

1. Introduction

In engineering construction, the material properties of joints are important factors in jointed rock masses [14]. Empirical models based on the roughness of the joint surface have become a focus of research on jointed rock masses [59]. Two of the most widely used empirical formulas based on the roughness of the joint surface are the bilinear strength criterion proposed by Patton [5] and the shear strength empirical expression proposed by Barton and Choubey [7].

With the development of numerical modelling techniques and computational software, Particle Flow Code (PFC) has been successfully applied in research on jointed rock masses. Instructive results have been obtained from numerical modelling [1015]. Especially, Ivars et al. [16] describe synthetic rock mass (SRM) modelling and use the bonded particle model for rock to represent intact material and the smooth joint contact model (SJM) to represent the in situ joint network. Their method improves the accuracy of rock mass simulation and greatly expands the application of PFC in jointed rock mass. By using this software, the interactions of particles are transmitted through virtual contacts and bonds. Particle displacements are not restricted, allowing large deformations of rock to be simulated. The interparticle bond can break readily when the external loading stress exceeds the bonding strength, reflecting the evolution of rock damage [1721].

The rationality of joint generation in PFC can be verified by simulations of direct shear tests. Two methods are commonly used for joint generation in PFC2D: the first approach is to remove the bonds between particles in the joints, and the second method is to employ the smooth joint model on particles located within the joints based on the discrete fracture network (DFN). However, numerical results of direct shear tests on rock masses with a horizontally persistent joint have indicated that neither of the two methods can reproduce realistic shear properties of joint rock. When using the first method, the generated joint surface is often irregular, and the shear force direction is determined by the contact direction because particles are generated randomly. Thus, the postpeak shear stress exhibits unrealistic oscillations with increasing shear displacement. In the second method, the mechanical properties can be properly simulated only when the shear displacement is smaller than a certain value, beyond which the postpeak shear stress oscillates unrealistically. This oscillation is mainly caused by the unbalanced contact forces of different sections of rock block resulting from the nonuniform grain distribution along the joint surface. When the shear displacement is larger than a certain value, particles may move across the joint surface under the vertical loads and touch other particles on the opposite side to form new contacts. Consequently, new contact bonds will be installed, and a stress concentration will form, leading to oscillations in the postpeak shear stress.

To solve the above problems, Chiu et al. [12] and Chiu et al. [22] modelled rock joints by artificially specifying a joint roughness coefficient. Although their method eliminated the variations in the dilatancy angle, the shear strength still oscillates, and the evolution of shear damage cannot be modelled properly. Bahaaddini et al. [13] proposed a method to generate joints with additional walls. In their method, an additional wall is first created along the joint, and then particles are generated on both sides of the wall. After reaching the initial equilibrium, this artificial wall is removed from the model, and a smooth joint model is used on the contact along the joint. However, it is difficult to generate a set of joints within the rock block when using this method of joint generation because the use of additional walls within the granular packing constrains particle motion, leading to a nonuniform initial stress distribution.

In this study, we explore an effective numerical method for joint generation in PFC2D. The remainder of this paper is organized as follows. In Section 2, the numerical model configuration of the direct shear test is presented, and the input microscopic parameters are calibrated. In Section 3, a systematic study of the commonly used joint generation methods is presented, and the advantages and disadvantages of each method are discussed. To overcome the disadvantages of these methods, a new method based on the particle grouping technique and a smooth joint model is discussed in Section 4. In Section 5, simulations of direct shear tests are performed on two types of specimens with standard shaped joints, which are generated according to the standard JRC (joint roughness coefficient) joint curve or triangular sawtooth joint curve. The obtained numerical results are compared with the results obtained using Patton’s and Barton’s models. Section 6 presents the discussion, and Section 7 presents the conclusions of this study.

2. PFC2D Model Configuration of Direct Shear Tests

2.1. Calibration of the Model Parameters

In PFC2D, the rock structure is discretized as a collection of disc-shaped particles with random distributions of particle positions and radii. The particle contacts can be described by either the contact bond (CB) model or the parallel bond (PB) model. In the PB model, the particle motion and interaction are governed by Newton’s second law. The interparticle bond will break when the normal force or shear force exceeds the bond strength, leading to damage evolution in the joints. In addition, the PB model also permits the transfer of moment between particles. Therefore, the PB model is employed for the tests presented in this study to simulate the mechanical behaviour of rock.

The choice of input microscopic parameters can considerably influence the accuracy of numerical modelling. In this paper, the input parameters of the numerical model are calibrated against the experimental results on Hawkesbury sandstone [23] (the experiments were conducted at the University of Wales, Australia). In the simulation, the rock specimen has a width of 84 mm, a height of 42 mm, a minimum particle radius , and a maximum particle radius of . This rock sample consists of approximately 12,817 randomly distributed particles. The elastic modulus of the numerical specimen is proportional to the elastic modulus of the contact and the elastic modulus of the PB. The shear strength of rock increases as the normal and tangential strengths of the PB increase. Poisson’s ratio increases as the ratio of the normal to the tangential stiffness of the contact increases. The input microscopic parameters used in the current study are listed in Table 1. Uniaxial compression tests are conducted on five randomly generated specimens, and the obtained mechanical parameters are compared with the experimental data in Table 2. The numerical tests accurately reproduce the mechanical characteristics of the experiments.

2.2. Numerical Modelling of Direct Shear Tests

The numerical simulation of a direct shear test on a persistent jointed rock specimen is one of the most effective ways to evaluate the accuracy of a joint generation method. The model configuration of the direct shear test used in the current study is shown in Figure 1. The specimen consists of 13,389 particles and has a height of 40 mm and a width of 100 mm. The input microscopic parameters are listed in Table 1. The loading frame is composed of frictionless walls and is divided into upper and lower parts. The lower part is fixed in space, whereas the upper part can move with movements in walls 1, 2, and 3. The movement of wall 3 is used to apply normal stress, whereas walls 1 and 2 are used to apply a constant shear velocity on the specimen. To maintain a quasistatic state during shearing, the shear velocity is set as 0.002 m/s, and the numerical time step is set to 8.158 × 10−8 s. Thus, the lateral wall moves 1.63 × 10−10 m at each step. During the simulation, the shear and normal displacements can be obtained by monitoring the average horizontal displacement of walls 1 and 2 and the vertical displacement of wall 3. The average shear stress along the shear plane can be obtained by monitoring the horizontal forces acting on walls 1 and 2.

3. Evaluation of Existing Generation Methods

3.1. Removing Bond Method

Cundall [24] proposed a method to generate a stripe of unbonded particles within the joint (Figure 2). This method has been adopted and discussed by many scholars [11, 25, 26].

For a horizontal persistent joint, when the friction coefficient of the joint approaches zero, the shear strength also decreases to zero. To determine the accuracy of this method, the friction coefficient of the joint (unbonded particles) is set to 0.03, 0.2, and 0.3, and the normal confining stress is set to 1 MPa. The relationships of shear displacement with shear stress and normal displacement are shown in Figures 3 and 4, respectively.

Figure 3 illustrates that the shear stress decreases suddenly after reaching the peak value. Then, it oscillates unrealistically, which is not consistent with the experimental observations. In addition, the normal dilatancy effect of the jointed rock specimen can be clearly observed in Figure 4, which is not the mechanical response of a horizontal persistent jointed rock. Although the particle size (maximum radius of 0.40 mm) is small compared to the rock specimen (100 mm in width), the joint surface formed by the unbonded particles has a degree of initial roughness due to the random generation of particles, as shown in Figure 2. During the shear process, the shear stress gradually concentrates on certain pairs of particle contacts along the joint plane because the joint surface has a jagged shape and the deformation resistance is uneven along this plane. Then, these pairs of particles rotate around each other and move away from the original joint plane. These complex particle movements finally result in the normal dilatancy effect and the oscillations in the shear stress. In addition, the increase in the joint width will only reduce the peak shear stress and cannot effectively eliminate the stress oscillations and normal dilatancy effect [24].

The above analyses demonstrate that the realistic mechanical behaviour of horizontal persistent jointed rock masses cannot be reproduced by the removing bond method due to the random distribution of particles in PFC2D.

3.2. Smooth Joint Model

Pierce et al. [27] proposed the smooth joint model to simulate the mechanical properties of a joint. This model can eliminate the negative influence of particle contact orientation in the direction of the reacting force on the joint. The smooth joint model is employed on jointed contacts whose particles are on both sides of the joint on the joint surface. The particles belonging to these jointed contacts slide along the joint surface, and the original contact direction is neglected.

3.2.1. Input Parameters of the DFN Smooth Joint Model

The stiffness parameters of a smooth joint are typically calibrated via numerical and experimental results. Park and Min [28] found that the stiffness of the smooth joint model in PFC can be uniquely determined by the joint stiffness obtained experimentally; specifically, it can be directly calculated from the experimental joint stiffness and particle distribution of a numerical specimen. The realistic joint surface is continuous, as illustrated in Figure 5(a). If particles are regularly generated in PFC without overlaps and the resulting length of the smooth joint is the same as that of the realistic joint, as shown in Figure 5(b), the stiffness of the smooth joint model is the same as that obtained experimentally. If the particles are randomly distributed in PFC (Figure 5(c)), the resulting smooth joint length differs from the real joint length, leading to a difference in the numerical and experimental joint stiffnesses. In this case, the numerical stiffness of the smooth joint model can be obtained by multiplying the experimental stiffness value with a correction coefficient, which can be calculated by the following equation:where is the length of the smooth joint section in PFC, is the number of contacts in the smooth joint model, and is the contact radius of the smooth jointed contact. , where is the radius multiplier and is typically set to 1.0 and and are the radii of two particles in the smooth jointed contact. is the length of the jointed surface in the equivalent continuum model, is the total number of continuous jointed surfaces, and is the geometric length of the continuous jointed surface.

The numerical stiffness parameters of the smooth joint model obtained from experiments on Hawkesbury sandstone [23] and the experimental analyses of Bahaaddini et al. [13] are shown in Table 3. The correction coefficient can be calculated by (1) after the generation of a numerical specimen.

3.2.2. Mechanical Properties of the DFN Smooth Joint Model

To verify the rationality of the DFN smooth joint model in simulating the mechanical behaviour of a joint, this model is used to generate a rock specimen with a horizontal persistent joint for numerical simulations of a direct shear test. The microscopic parameters used in the simulation are provided in Table 1, and the stiffness of the joint is provided in Table 3. The friction coefficient of the joint is set as 0.2, 0.3, or 0.57, and the normal confining stress is set as 3 MPa. Various parameters, such as the shear and normal displacement, shear stress, particle motion on both sides of the joint surface, and number of particles moving across the jointed plane, are monitored during the simulation. The results are presented in Figure 6. The shear stress and displacement curves correspond well with the experimental results in the initial stage, and no particle moves across the jointed plane. After reaching a certain shear displacement, particles begin to move across the jointed plane, with the total number increasing from 0 to 10, 11, and 14 in Figure 6(c). The shear stress oscillates, and an unrealistic dilatancy effect is observed from the normal displacement.

By default, the smooth joint model will be applied to any newly formed particle contact if particles are located on both sides of the jointed plane. However, particles initially belonging to the smooth joint may move from one side to the other side of the jointed plane when the shear displacement is sufficiently large. Such movement leads to a stress concentration on the newly formed contact, with macroscopic responses of shear stress oscillation and normal displacement dilatancy. These phenomena have been analysed in detail by Bahaaddini et al. [13].

3.3. Addition Wall Method

Bahaaddini et al. [13] proposed a method to generate joints by additional walls, as shown in Figure 7. In this method, the model is first divided into an upper part (A) and a bottom part (B) by a joint surface, which is replaced by additional walls. After the generation of particles in regions A and B reaches an initial equilibrium, the additional walls are removed, and the contacts along the jointed surface are set to the smooth joint model.

This method has the following three disadvantages:(1)For an irregular joint, many additional walls must be generated. During the cycles to reach the initial equilibrium, nearly every wall must move to bring the particles to equilibrium and achieve a uniform initial stress. Thus, if the joint is composed of many walls, it is difficult to maintain the uniformity of stress on each wall and form a uniform stress field inside the specimen.(2)For multiple groups of persistent joints (Figure 8(a)), a group of parallel walls should be placed, which means that the movements of nearly all the particles are restricted by these parallel walls during specimen generation. Thus, it is difficult to achieve an ideal equilibrium state because of the uniform stress field, as shown in Figure 8(c). The stress values range from 0 Pa to −3.8e6 Pa. This disadvantage is more discernible when the joint density is higher.(3)For multiple groups of nonpersistent joints, as shown in Figure 8(b), if additional walls are fixed, the same problem of achieving an ideal equilibrium arises because of the uniform stress field, as shown in Figure 8(d). The stress values range from −2.0e5 Pa to −3.4e6 Pa. If additional walls are allowed to move, certain walls would drift away from their designed location under the collision of kinetic particles during specimen generation.

After the additional walls are removed and the initial confining stress is applied, the smooth joint model is only employed on partial segments of the joint surface because particles on both sides are generated randomly without fully touching along the joint under the confining stress. Thus, the confining stress is concentrated within the contact area, as shown in Figure 8(c).

4. New Joint Generation Method

Because of the disadvantages of existing methods, particularly for generating multiple groups of persistent or nonpersistent joints, a novel joint generation method is proposed based on the particle group technique in PFC. At the beginning of specimen generation, the particles can move freely. After the specimen is formed, the particles are divided into several groups according to the regions partitioned by the joints, and the smooth joint model is employed on the contacts located on the boundaries of two groups of particles through an entire loop for contacts. During the cycles, a newly formed contact is identified by contact searching, and the smooth joint model is implemented if it is a contact on the boundaries of two groups of particles. The process of generating a jointed specimen for a direct shear test can be divided into the following eight steps:(1)Generate the specimen boundary: Walls are generated according to the size of the specimen. Because the specimen is used in a direct shear test, the vertical walls on both sides of the specimen are divided into upper and lower parts, as shown in Figure 9(a).(2)Generate the initial granular assembly: Particles are generated randomly inside the boundaries according to the initial porosity, as shown in Figure 9(b). The initial void ratio is set to 3.0%, the minimum radius of the particles is set to 0.25 mm, and the maximum radius is set to 0.4 mm. The stiffnesses of the particles and walls are set accordingly; the wall stiffness is 10% higher than the particle stiffness. Particles can move freely under the loading of the linear contact force, and thus, the particle velocities are set to zero every 10 cycles. The initial specimen is obtained after the initial force balance is reached and the externally overflowing particles are removed, as shown in Figure 9(c).(3)Apply the specified isotropic stress: A small isotropic stress (e.g., 0.01 MPa with a relative tolerance error of 0.003) is applied to the specimen. This isotropic stress is applied by the movement of the walls on which the stresses are monitored. When the relative tolerance errors of the monitored stresses are less than 0.003 compared to the target stress of 0.01 MPa, the movements of walls are stopped, and the initial isotropic stress is applied. This stress is used to eliminate the particle interlocking forces.(4)Elimination of floating particles: A floating particle is defined as a particle in contact with less than three surrounding particles. A new particle is generated at the location of the removed particles, and its position and radius are adjusted to ensure that there are at least three contacts around one particle in the specimen.(5)Particle grouping: Particles are divided into upper and lower parts by the joint, and different groups of particles are marked by different colours, as shown in Figure 9(d).(6)Employ the PB model and smooth joint model: The contacts whose particles have a surface distance of less than 5 × 10−5 m are identified. If the particles are located at the interface of two groups, the smooth joint model is employed for the contacts; otherwise, the PB model is employed. The sliding direction of the smooth joint model is assigned based on the relative positions of the contact and the joint plane.(7)Set the initial distance of the contacts: The initial distances of all the contacts are set to the current surface distance of particles for both the PB model and smooth joint model, and the contact force and displacement are set to zero.(8)Update the contacts during cycles: A newly formed contact is searched and identified; the smooth joint model is implemented if two particles in contact belong to different groups of particles and the location of this contact is within the effective joint range. The proposed method of joint generation was presented via the flow chart to corresponding steps, as shown in Figure 10.

The procedures for generating multiple groups of persistent or nonpersistent joints are similar to the above generation procedure and employ the group technique for different regions of particles. Consider the specimens with persistent and nonpersistent joints in Figures 8(a) and 8(b) as examples. Particles inside the specimens are grouped according to the distribution of joints, and different groups are assigned different colours in Figures 11(a) and 11(b) such that the locations of joints and groups of particles can be easily identified. During the installation of the smooth joint model in step (6), the joint model is installed by identifying the group to which two contacting particles belong. The other steps are the same as those listed above. Figures 11(c) and 11(d) show the stress distributions that occur after the initial confining stress (1 MPa) is applied to the two specimens generated by the proposed method. The stress values range from −2e5 Pa to −3.2e6 Pa for persistent joints and from −4e5 Pa to −2.6e6 Pa for nonpersistent joints. A comparison of Figures 8(c) and 8(d) illustrates that the stress field is more uniform for the initial confining stress. Thus, the proposed method is more feasible than the existing methods, particularly for partially developed joints.

5. Validation of the Proposed Method

To validate the rationality of the proposed rock joint generation method, specimens with a triangular sawtooth shape [5] and joints following the JRC standard [7] are generated for a numerical direct shear test. The simulation results are compared to the theoretical values. The generated rock specimen has a height of 40 mm and a width of 100 mm. The values of the microscopic parameters and the stiffness of the joint are provided in Tables 1 and 3, respectively, and the friction coefficient of the joint is set to 0.2, 0.3, and 0.58. The shear loading rate is set to 0.002 m/s, and the normal confining stress is set to 1, 2, 3, and 4 MPa.

5.1. Simulation of a Horizontal Persistent Joint

Direct shear tests are performed on joint rock specimens with a horizontal persistent joint. The shear stress and displacement curves are shown in Figure 12. The oscillation of the postpeak shear stress is eliminated, indicating that the stress concentration along the jointed plane has disappeared. Thus, the proposed rock joint generation method accurately simulates the mechanical properties of the joint.

5.2. Simulation of a Triangular Sawtooth Joint

Direct shear tests are conducted on a specimen with a triangular sawtooth joint using the proposed method. Four types of sawtooth joints are generated with the inclination angles of 15°, 20°, 25°, and 30°. The shear strengths of the joints with different inclination angles are compared to the theoretical results obtained using Patton’s model [5] in Figure 13. The simulation result for the sawtooth joint with an inclination angle of 30° and friction coefficient of the joint of 0.58 is shown in Figure 14. Figure 15 shows the damage evolution along the sawtooth joint with an inclination angle of 30° under different normal stresses when the shear displacement reaches 3 mm.

Figure 13 illustrates that the numerical results for the triangular sawtooth jointed specimen correspond well with the theoretical results. Figure 14 shows that as the normal stress increases, the peak shear stress increases, and the normal displacement and dilatancy decrease, which is consistent with the experimental observations. When the shear displacement is 3 mm for the specimen with 30° inclined sawtooth joints, 163, 181, 244, and 258 fractures form under normal stresses of 1, 2, 3, and 4 MPa, as shown in Figure 15. Therefore, the number of fractures increases with increasing normal stress for a given shear displacement, which is consistent with the experimental results.

5.3. Simulation of the Joint with Standard JRC Roughness

Barton and Choubey [7] proposed standard roughness curves; five standard roughness curves with JRCs of 0–2, 6–8, 10–12, 14–16, and 18–20 and a friction coefficient of the joint of 0.58 are selected for this study, and direct shear tests are conducted on the numerical joint specimen. The standard JRC curve is well approximated as a polyline sampled at intervals of 0.5–1.0 mm along the horizontal direction. The particles are divided into several groups by the polyline, and then, the smooth joint model is employed on the contacts at the interface of two particle groups. The shear strengths of different JRC joint surfaces are compared to the theoretical results obtained from Barton’s model [7] in Figure 16. The numerical results for the joint with a JRC of 10–12 are presented in Figures 17 and 18.

Figure 16 illustrates that the numerical results for different standard JRC joints correspond well with the theoretical results obtained using Barton’s model [7]. Figure 17 shows that the peak shear stress increases and the normal displacement decreases with increasing normal stress. When the shear displacement is 3 mm, 96, 117, 138, and 184 fractures form under normal stresses of 1, 2, 3, and 4 MPa, respectively, as shown in Figure 18. Therefore, the number of fractures increases with increasing normal stress, consistent with the experimental results.

5.4. Damage Evolution of Joints

The asperity degradation of a joint surface can be observed in the simulations after the shear strength reaches a certain value; the number of cracks increases and the shear stress decreases with increasing shear displacement. The numerical results for the joint specimen with a JRC of 18–20 subjected to a normal stress of 3 MPa are chosen for further analyses, as shown in Figure 19.

The shear process has four stages, as shown in Figures 19(a) and 19(b): (1) The elastic stage: when the shear displacement is less than 0.4 mm, the shear stress increases linearly with increasing shear displacement. No crack appears during this stage. (2) The elastic-plastic stage: when the shear displacement is between 0.4 and 1.37 mm, the shear stress increases with increasing shear displacement, and the number of cracks starts to increase slowly. (3) The softening stage: after the stress reaches its peak, as the shear displacement increases from 1.37 to 1.55 mm, the shear stress decreases rapidly, and the number of cracks increases abruptly. (4) The residual stage: when the shear displacement increases from 1.55 to 3.0 mm, the residual shear stress maintains a nearly constant value. The number of fractures increases gradually from 92 to 200, and the fracture growth rate decreases. As shown in Figure 19(c), almost no fractures occur when the shear displacement is small (e.g., 0.4 mm), and the number of fractures increases gradually with the shear displacement.

6. Discussion

Based on a single joint numerical specimen, a new joint generation method was proposed, and its accuracy was verified. As shown in Figure 20, according to the dip directions of joints in the actual rock mass, the joints can be divided into several sets. Taking one set of joints as an example as shown in Figure 21, the new method for joint generation is explained in this section, and the other sets of joints can be generated in the same way. For each joint, two groups of particles are identified according to their location relative to the joint, that is, the right or left side, within a certain range, whose size is suggested as 2 times the particle diameter, as shown in Figure 22(a). Because the number of joints in one set is large, the identification numbers of particle groups are different not only between two sides of a joint but also among different joints. As shown in Figure 22(b), different colours represent different particle groups. These group numbers are employed to identify the contact model of newly formed contacts in the calculation process. If two newly contacted particles belong to the corresponding groups of a joint, this contact is set to the smooth joint model. The resulting model for a single set of joints is shown in Figure 23.

7. Conclusions

Three existing joint generation methods were analysed by conducting numerical direct shear tests on jointed rock mass specimens in PFC2D. The disadvantages of these methods were demonstrated via analyses of the mechanical behaviours of the generated joints under shear loading. For a horizontal persistent joint generated by the removing bond method and the DFN smooth joint model, the postpeak shear stress oscillates considerably, and the dilatancy of the normal displacement with increasing shear displacement does not correspond well with the experimental results. Although realistic mechanical behaviours of a joint can be simulated if the joint is generated by the method of additional walls, the application of the method is limited in modelling a persistent joint due to the complexity of the joint generation process. When the method of additional walls is used to model multiple groups of joints that are either persistent or nonpersistent, the movements of particles are restricted by the additional walls, leading to difficulties in achieving a uniform initial stress field and generating an ideal jointed specimen.

This study proposed a novel method to generate joints based on the group technique. The main concepts of this method are the division of particles into different groups based on the distribution of joints and the application of the smooth joint model to the contacts located at the interfaces of particle groups. The newly formed contacts during the simulation are forced to follow the same rules such that real-time history and identification are obtained for the formation of contacts during the simulation. Jointed specimens with triangular sawtooth joints and a standard JRC joint were generated by the proposed method and used to conduct a direct shear test. The numerical results were compared with the theoretical values obtained using Patton’s and Barton’s models, and good agreement was obtained. These results verify the validity and rationality of the proposed joint generation method. The proposed model can also simulate the asperity degradation of rough jointed rock masses.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was funded by the Water Resource Science and Technology Innovation Program of Guangdong Province (no. 2014-09) and The Youth Fund of China Natural Science Foundation (no. 5170090818).