Abstract
River networks and estuaries are very common in coastal areas. Runoff from the upper stream interacts with tidal current from open sea in these two systems, leading to a complex hydrodynamics process. Therefore, it is necessary to consider the two systems as a whole to study the flow and suspended sediment transport. Firstly, a 1D model is established in the Pearl River network and a 3D model is applied in its estuary. As sufficient mass exchanges between the river network and its estuary, a strict mathematical relationship of water level at the interfaces can be adopted to couple the 1D model with the 3D model. By doing so, the coupled model does not need to have common nested grids. The river network exchanges the suspended sediment with its estuary by adding the continuity conditions at the interfaces. The coupled model is, respectively, calibrated in the dry season and the wet season. The results demonstrate that the coupled model works excellently in simulating water level and discharge. Although there are more errors in simulating suspended sediment concentration due to some reasons, the coupled model is still good enough to evaluate the suspended sediment transport in river network and estuary systems.
1. Introduction
As a link between marine environments and rivers, estuaries are characterized by a variety of complex and complicated processes [1]. Therefore, there has been a growing interest in the field of estuarine and coastal sediment dynamics, including sediment transport, geomorphological process, and the transport of heavy metal by sediment particles [2, 3]. In the researches on hydrodynamics and suspended sediment transport in this region, numerical model is no doubt a highly efficient and low-cost method to describe estuarine flow and sediment transport processes in recent decades. Normally, the 1D model is more suitable than the 3D coastal model to simulate the flow and suspended sediment transport in river network due to the complex river network geometries and complex structures like weirs, barrages, and dams. The 3D model is more useful in evaluating the coastal flow and sediments movement due to the vertical circulations. However, if both the upstream river network geometries and the downstream tidal currents in the estuary are complicated, application of individual 1D or 3D model for the river and the estuarine area will result in either an inadequate representation of these processes (in case of 1D) or a computationally inefficient method (in case of 3D) [4]. Therefore, it is necessary to take upstream river network and downstream estuarine system as a whole in one modelling system, to resolve the relevant physical processes in a consistent way.
A lot of efforts have been made recently to develop coupled modelling systems, which were designed for different estuaries [5–7]. Some of the coupled models were used to predict flood effects [8], some to investigate mass fluxes and transformations of nutrients [9], and some to prevent further water pollution [10] and so on. Some researchers developed the “domain decomposition” technique to increase computational efficiency [6]. This technique creates efficient model grids by “stitching” together different partial grids, which enable the high resolution and low resolution areas to be combined and modelled. However, this technique may not be stable due to the existing nested grids. Zhou et al. [10] provide a method that the 1D and 3D schematisations do not need to overlap. The upstream water level of the 3D model is used for the 1D outflow open boundary and the downstream outflow of the 1D model is provided to the upstream inflow of the 1D model [10]. The 1D model and the 3D model run simultaneously and the interface information will be saved for offline computation. Offline coupling requires huge storage space. In order to resolve this, an online coupling method is developed. This method resolves the real-time exchange and interaction between the coastal waters and river network, allowing for accurate and mass conserving modelling of complex coastal water and river network system. It is easy to implement as it only replaced file I/O (read) by communication operations (receive) [11].
The Pearl River Estuary (PRE) and its upstream river network form the largest river system in the southern China. Historically a rich and fertile area, it has now become one of the most populated and most developed regions in China. However, as the economy grows in the region, serious problems also occur at the Pearl River network and its estuary, such as salinity intrusion [12], riverbed erosion, and water pollution. Therefore, coupled 1D–3D hydrodynamic modelling, such as horizontal fronts, vertical stratification, and salt intrusion, was developed to simulate the hydrodynamics and salinity in this area [4]. Suspended sediment transport, however, was not involved in this coupled model. Recently, a coupled model has been used to investigate the transport of pollutants and water quality parameters, to provide a scientific basis for the improvement of the water quality of the Pearl River and to prevent further water pollution [10]. In addition, the models were also used for the construction of the nutrient budgets for the entire area and the characterization of the key biogeochemical processes in the transfers and transformations of nutrients [9]. However, previous works on the coupled models fail to provide mathematics derivation on the relations of water level between 1D and 3D models. Moreover, most coupled models focus on the hydrodynamics or water quality in the river network and its estuary, but little attention has been paid to the suspended sediment transport. In fact, the transport of most organic carbon from land to oceans is closely related to river sediment transport. Hence, understanding the estuarine dynamics of suspended sediment is also crucial to resolve estuarine problems [11, 13]. Therefore, this paper provides a method to establish a coupled model of hydrodynamics and suspended sediment transport by strict mathematical derivation, rather than physical interpretation in the river network and its estuary.
2. Methods
The model coupled the 1D river network model with the 3D coastal model. It takes into consideration the complex topography of the river network and the vertical gradients of the coastal hydrodynamics. The model is applied in the PRD, which is divided into two parts. The 1D model is set up in the Pearl River network and the 3D model is applied in the PRE to investigate the transport of water flow and suspended sediment in the delta.
2.1. One-Dimensional Hydrodynamics and Suspended Sediment Model
Based on the Saint-Venant equations and the nonequilibrium transport equation for the suspended sediment, a numerical model of junction-channel for the river network is established by using the “junction-control” method for water level and suspended sediment concentration. Firstly, the river network is schematized as river channels, junctions, and cross-sections. Secondly, the variables at each junction are calculated by solving the governing equations. Finally, the water levels, discharges, and suspended sediment concentrations at each cross-section in each river channel are calculated. More details about the hydrodynamics difference scheme and transforms of the finite-difference equations can be found in Zhang et al. [12]. The algorithm for suspended sediment concentration in the river network is similar to that for unsteady flow. The details are introduced in the following parts.
The governing equation for the nonequilibrium transport of suspended load [14] is expressed as
The riverbed transfiguration equation is where is the discharge, is the cross-section-averaged suspended sediment concentration, is the cross-section-averaged suspended sediment transport capacity, is the cross-sectional area, is the water surface width, is the mean elevation of the cross-sectional water surface, is the mean elevation of the riverbed, is the cross-section-averaged velocity, is the settling velocity of suspended sediment, is a recovery coefficient relating to the suspended sediment of the saturation, is the bed load sediment transport rate, is the cross-sectional mean grain diameter of bed load, and is the dry density of suspended sediments.
For suspended sediment transport equation, the upwind finite difference scheme and Preissmann four-point implicit difference scheme are employed here. According to the direction of flow, the difference equations are written as where , , , and are known coefficients of the difference equations.
In a tidal river, flow direction changes with ebb and flood current. The flow at a single river channel can be classified by four types (Figure 1): (a) downflow, , ; (b) upflow, , ; (c) faced-flow, , ; (d) depart-flow, , ; and are the discharge at the start and end cross-section of a single channel.
For the downflow, it is single direction flow. If the sediment concentration at the up-boundary node () is known, the sediment concentration at each node can be obtained by the following formula: where and are the coefficients; they can be gotten by (3) and written as
For the upflow, it is also single direction flow. If the sediment concentration at the down-boundary node () is known, the sediment concentration at each node can be obtained by the following formula: where and are the coefficients; they can be gotten by (4) and written as
For the faced-flow, it is double direction flow. The single river can be divided into two single direction flow river segments. If the sediment concentrations at the up- and down-boundary nodes (, ) are known, the sediment concentration at each node can be obtained by formulas (5) and (7).
For the depart-flow, it is another double direction flow. Firstly, the position of the stagnant point () must be determined, the sediment concentration at the stagnant point () may be obtained by the transport equation, and then the single river can be divided into two single direction flow river segments in which interboundary is the stagnant point; finally the sediment concentration at each river junction can be obtained by the formulas (5) and (7). The position of the stagnant point is calculated by the following formula:
Finally, the sediment concentration at the stagnant point is
2.2. Three-Dimensional Hydrodynamic and Suspended Sediment Model
The 3D hydrodynamic and suspended sediment model is established based on ECOM model, a free-surface, fully nonlinear, primitive-equation estuarine and coastal ocean model. The equations under the orthogonal curve coordinate include water continuity equation and the momentum equations: where is the Coriolis coefficient, is the referenced density, is the local density, is the gravitational acceleration, and are the turbulent viscosity coefficients of turbulent momentum in vertical and horizontal, respectively, and , , and are the velocity in , , and direction in the Cartesian coordinate system.
The diffusion process of salinity is also necessary to consider for the model. The conservation equations are described as where and are the diffusion coefficients of the vertical and horizontal, respectively, and is potential temperature.
The transport equation of the suspended sediment is written as where is the sediment concentration and is the settling velocity of sediments.
The details of ECOM model can refer to Zhu et al. [15].
2.3. Coupled 1D and 3D Hydrodynamic and Suspended Sediment Model
In this paper, the relationship expression of water level between 1D model and 3D model is deduced by strict mathematics method. The 1D and 3D models meet at outlets of river network and estuary, and the hydroinformation is transferred through interfaces, by which the additional hydrodynamic equation of the 1D and 3D coupled numerical model can be derived. The 1D and 3D boundary conditions at interfaces are satisfied through iteration.
The water level and discharge at the interfaces of the 1D model and 3D model can be reasonable considered equal: where and are water level and discharge at the th interface of the 1D model, respectively, and and stand for water level and discharge at the th interface of the 3D model, is the discharge of the grids, is the number of the grids of the interface, and is the number of the layers of water depth.
Due to the difference of water levels at grids of 3D model, the down-boundary of the 1D model is the averaged water levels of the grids in the 3D model:
Based on the method of the 1D hydrodynamics model, the governing equation can be deduced to matrix of water level at junctions as below: where is the coefficient matrix, is the water levels at junctions, and is constants.
However, in a coupled model, the interfaced water levels are not known, which lead to containing unknown term and constant term. Therefore, the unknown part should be separated from the matrix : where is the water levels at the interfaces, is the corresponding coefficient matrix, and is the constant term.
Gauss elimination method is used to simplify the coefficient matrix to a diagonal matrix :
Then, the water levels at junctions connected to the interfaces can be expressed by the water levels at interfaces of 1D model and 3D model: where and are coefficient and constant after unitization, respectively.
Utilizing the relations between the water level and discharge at junctions [12], the discharge at the beginning and the end of a single river channel can be deduced as below:
Equation (19) is substituted to (20) and (21), respectively; then the discharge at the interfaces of 1D and 3D models can be totally represented by the water level at the interfaces, and are the coefficients of the equations. Consider the following:
After simplifying (22), the relations between water level and discharge can be established in 1D model. Consider the following: and is the coefficient of equations and is constant.
Finally, taking into account (14), an iterative equation for water level at interface is established as below:
At the beginning of the calculation, the initial discharges are used as the upper boundary condition of the 3D hydrodynamics model and can get the water levels . Then, substituting these water levels into (24) is to get the approximate values , which can be regarded as the downstream boundaries of 1D hydrodynamics if water levels satisfy the convergence condition ( is a given error). Otherwise, have to be substituted into (23) to get discharges for 3D model calculation again. This iterative process stops until the convergence condition is satisfied.
As for the suspended sediment model, two additional conditions are needed to be considered. Firstly, the sediment fluxes at the interfaces of the 1D model and 3D model should be equal:
Secondly, the spatial averaged sediment concentration of 3D model should be closed to that of 1D model:
2.4. Study Area and Data Processing
The Pearl River Delta (PRD) is located in the northern of the South China Sea. It is a dynamically complex estuarine system, which encompasses Pearl River network and PRE. The Pearl River, a subtropical river, is the second largest Chinese river in terms of annual water discharge (over 3 × 1011 m3). The channel network is formed by a large number of intricately interlaced deltaic branches, including three principal rivers: Xijiang, Beijiang, and Dongjiang, and some small rivers draining the PRD. It is a catchment area of 26,820 km2 and a total length of 1,600 km [16]. Meanwhile, as the distance from one junction to another is 0.68–1.07 km/km2, the Pearl River network is considered as one of the world’s most complicated river networks, which is formed by a large number of intricately interlaced branches [12]. In the river network, the mean suspended sediment concentration is about 0.284 kg/m3 [17], and annual mean sediment load in 1957 to 2006 is 75.3 × 109 kg/yr, with an annual average sediment load of 66.83 Mt/yr in Xijiang alone [18].
These fluxes pass through the river network and flow into the PRE through the eight outlets Hutiaomen and Yamen discharging into the Huangmao Bay, Modaomen and Jitimen discharging into open sea directly, and Humen, Jiaomen, Hongqimen, and Hengmen discharging into the Lingding Bay [4]. The PRE has a broad bell shape, in which the Lingding Bay is like an inverted funnel with the narrow neck in the north and wide mouth opening to the south [19]. As the gravitational circulation develops in the PRE, the residual currents are landwards in deep channels and seawards in the western part of the PRE [20]. Normally, the circulation and water properties in the PRE are caused by various forcing mechanisms, including river discharges, tides, monsoon winds, coastal currents, and buoyancy forcing associated with the mixing of freshwater and saline water. The tides in the PRE mainly come from the Pacific oceanic tidal propagation [21], with a mean tidal range between 1.0 and 1.7 m. Thus PRE can be classified as a microtidal estuary. Propagating from offshore towards the estuary, the tidal range increases gradually and reaches a maximum near the river outlets [19]. Salinity in the PRE seldom exceeds 32 psu even during the dry season and is highly stratified in the vertical with a very low salinity surface layer. The distributions of salinity have a significant seasonal variability [7]. Further, the carrying capacity of sediments in the PRE is basically controlled by tidal flow, wave, fluid, and particle features. In the PRE, the very fine-grained silt and clay are contained in the suspended sediments, within which most particles are less than 100 μm in diameter and the median particle diameter is about 8 μm [22].
2.5. The Setup of the Coupled 1D and 3D Model
The model is applied to the Pearl River network and its estuary. The river network in 1D model is discretized into 347 channels, 13 boundary rivers, and 1,850 cross-sections with the interval varying from 0.2 to 3.0 km and 220 junctions. The length of the river network is about 1,600 km. The upstream boundaries include Shizui at the Tanjiang River, Gaoyao at Xijiang, Shijiao at Beijiang, Laoyagang at the Liuxi River, and Boluo at Dongjiang. The downstream boundaries include Dahu at Humen, Nansha at Jiaomen, Wanqingsha at Hongqili, Hengmen at Hengmen, Denglongshan at Maodaomen, Huangjin at Jitimen, Hutiaomen at Hutiaomen, and Huangchong at Yamen. These eight outlets are also the interfaces of the 1D and 3D model. The 3D estuarine model covers the area from the eight outlets to the −30 m isobaths (Figure 2). More than 57,000 orthogonal curvilinear grids are used in the 3D model with spatial resolution from 0.05 km to 1 km and 10 equidistant sigma levels in vertical.
The hourly observed runoff at the upstream stations is the upper boundary conditions of the coupled model. The salinity at those stations is definitely zero. At the open sea part of the coupled model, the tidal levels are evaluated by the harmonic analysis from the most important tidal constituents of the PRE, with the uniform values of temperature and salinity from the observations. Consider the following: where is the mean water level, is the amplitude of partial tide, and is the phase of partial tide.
The simulation is initialized with zero water levels and velocities in the bathymetry of 1998-1999 and runs about 30 days. The results of the model are used as the new initial conditions and the simulation hot-started with the results. The processes are repeated several times until the error between the last two simulations satisfies a given condition. These results are then used as the initial conditions of the coupled model. Two different periods, including July 1998 (wet season) and December 1998 (dry season), are used for model calibration. In both the 1D and 3D model domains, the nonuniform parameters, such as bed roughness parameter and particle diameter, have been adjusted in order to optimize the water level, discharge, and suspended sediment concentration over the upstream river network and its coastal waters. The specific key parameters of the coupled model are shown in Table 1.
3. Results
The coupled model uses the field measurements in July 1998 (wet season) and December 1998 (dry season) at eighteen stations (Figure 2) as calibration. The model skill score () is employed to evaluate quantitatively the model performance of water level, discharge, and suspended sediment concentration against the measured time series:
This method is widely utilized by numerical modeling system [12, 23] and the performance levels are categorized as follows: is excellent, 0.5–0.65 is very good, and 0.2–0.5 is good; if , it means a poor fit [24].
The model skill scores for water levels are fairly high (Figures 3(a) and 4(a)), indicating that the coupled model can accurately simulate water level in river network. Among the 18 stations, most of them performed excellently with only 4 stations scoring between 0.5 and 0.65 in the wet season (Figure 3(a)). In the dry season, 14 stations score more than 0.65 and the rest score between 0.5 and 0.65 (Figure 4(a)). In terms of the discharge, although the skill scores are not as excellent as that of water level, the model skill scores still show a very good skill of simulating discharge. In the wet season, the skill scores for most stations are greater than 0.5, implying a very good skill, although 4 stations score between 0.2 and 0.5 (Figure 3(b)). The coupled model also performed well in simulating discharge in dry season. Among the 18 stations observed, the skill scores of 13 stations are greater than 0.5 (Figure 4(b)), also showing a very good skill. Based on the agreement between the simulated data and the measured data throughout the river network, it can be concluded that the coupled model can simulate hydrodynamics process well in the river network and its estuary. Although the coupled model results do not perform quite as well as water level and the discharge, the simulated suspended sediment concentration in the Pearl River network also seems to be in a good agreement with the observed data. Figures 3(c) and 4(c) show that most of stations get their skill scores between 0.2 and 0.5 for both the wet and dry season. In the dry season two stations score below 0.2 and in the wet season three. There are further two stations in the wet season and one in the dry season that score greater than 0.5. Generally speaking, the simulated suspended sediment concentrations do not match the observed excellently but still satisfy the demand.
(a)
(b)
(c)
(a)
(b)
(c)
Figure 5 shows the results of the skill scores of the velocity at different layers in the PRE. In general, the skill scores of the velocity range from 0.4 to 0.6, showing a good or a very good performance. What is more, the skill scores of the velocity at the surface layer are a little bit better than those at the bottom. As for the suspended sediments, although the skill scores of the suspended sediment concentration (Figure 6) are relatively smaller compared with the skill scores velocities (Figure 5), the skill scores of all the stations are above 0.2, indicating that the model does a good job to simulate the suspended sediment concentration in the PRE. The skill scores of suspended sediment for some stations are even greater than 0.5, showing a very good performance. The PRE is a sedimentary system characterized by intricate depositional structures in space and time [11]. Due to the intensive anthropogenic activity, the bathymetry of the entire Pearl River network has changed dramatically. In particular, because of the large-scale sand excavation, the impacts of the anthropogenic activity cannot be neglected in the investigation of the erosion and deposition in the river network and its estuary. Figure 7 quantifies the impacts of the sand excavation in the bathymetry evolution. It compares the sediment erosion between the sand excavation and the natural conditions from Sixianjiao to Chaolianzhoutou in 1998. Figure 7 provides an insight into the amount of sediment in reaches under the sand excavation and the natural erosion and deposition. It is clear that the amount of the sand excavation is close to that of sediment erosion in Reach 1 and Reach 2 (Figure 8(a)) and larger than the amount of sediment erosion in Reach 3 and Reach 4. In Reach 4, the amount of sand excavation is almost 2.5 times that of the natural erosion. Therefore, the amount of sand excavation has to be taken into account in the comparison between the observed and the modelled sediment erosion and deposition in the river network. Figure 8 shows the thickness of simulated sediment erosion and deposition in the aforementioned reaches in 1998. The results show that erosion generally occurs in these channels in a year, although the river channel sedimentation plays a dominant role in the wet season. There is some discrepancy between the simulated and the measured river bed (Figure 8(b)), and the results still show that the model can be used to evaluate sediment transport in this region.
(a)
(b)
4. Discussion
The results show that the model skill score for sediment concentration is lower than that for water level and discharge. There are some potential reasons for the relative larger discrepancy between the simulated and observed value. First of all, more parameters in sediment simulation need to be selected than that of water level and discharge simulation. Normally, the roughness coefficient is the most important parameter in hydrodynamics simulation. However, large domain simulation in suspended sediment transport needs more data and parameters, such as the particle diameter and sediment settling rate. These parameters have large influences on the suspended sediment concentration simulation. However, due to the lack of data, some of these parameters are treated as relatively uniform distribution over the Pearl River network. This leads to error when simulating sediment transport. Secondly, the bathymetry of the Pearl River network changes dramatically due to the intensive sand excavation during the last three decades. Therefore, the asynchrony of hydrographic and bathymetry measurement may directly lead to the inaccuracy of the simulation results, especially in suspended sediment transport simulation. In fact, in our model, the field measurements were taken in 1998, but the bathymetry survey was done in 1999. The asynchrony of hydrographic and bathymetry measurement should be another reason for the simulation errors. Finally, some complex impacts, such as wind and discharge let-out by nearby reservoir, are not considered in this model, which probably affect the results of the coupled model.
This model was then used to simulate the flow division in the PRD, because flow division at some important junctions is important for the flow and sediment transport in distributary channel networks. Makou and Sanshui Stations are located at the apex of the delta. Their water discharge enters the Xijiang and the Beijiang, respectively. Table 2 shows the diversion ratio of flow at the two main junctions in the Pearl River network: (a) Makou and Sanshui and (b) Tianhe and Nanhua. The diversion ratios of flow between Makou and Sanshui Stations are 75.7% and 24.3%, respectively, showing the disequilibrium of the flow distribution in the waterway network. At the same time, the diversion ratios of flows at Tianhe and Nanhua are relative equal with the 54.2% in Tianhe and 45.8% in Nanhua. Although the flow flux displays an obvious seasonal variation, the diversion ration of flow is relative stabile. The flow flux at Makou Station is 1185 × 108 m3 in the wet season and 429 × 108 m3 in the dry season. The diversion ratio of flow is 74.5% in the wet season and 79.3% in the dry season. The diversion ratios of flows at the Tianhe and Nanhua Stations are also relatively stable in the wet season and the dry season. It can be observed that diversion ratio of the suspended sediment transport is similar to the flow ratio (Table 3). More sediment enters the Xijiang River as the flow distributes at the Makou and Sanshui station. Although the amount of sediment is much lower in dry season than that in the wet season, the majority of sediment still enters Xijiang River. Table 4 shows the annual flow flux of the eight outlets of the PRD. It can be noticed that Modaomen, with the annual water flux of 623 × 108 m3, is the largest contributor of the flow flux to the sea, accounting for 24.8% of the total flow flux from the delta to the sea. Humen comes second, accounting for the 21.4% of the total flow flux to the sea. Flow flux from Jiaomen outlet is 515 × 108 m3, accounting for 20.5%. Flow fluxes from the rest of the outlets are relatively small, which are 205 × 108 m3, 301 × 108 m3, 64 × 108 m3, 79 × 108 m3, and 187 × 108 m3, respectively. The annual sediment fluxes from outlets to the sea are in general consistent with the individual flow ratio. Modaomen is still the most important outlet to transport sediment from delta to the sea, accounting for nearly 30% of the total sediment flux to the sea. Following Modaomen, there are Humen and Jiaomen, both of which transport more than 40% sediment from PRD to the sea. The changes in sediment transport in the Pearl River delta and its estuary have great impacts on the water quality in this region. Therefore further study should be conducted to improve the water quality of the river. Nowadays, the application of nanotechnology for the remediation of contaminants may give promising results in future [25, 26]. An insight of the flow and sediment transportation mechanism in the river and estuary system has been obtained through the paper, and a reliable basis to control water pollution of the Pearl River Delta can be proposed combined with the present studies.
5. Conclusions
In deltaic estuaries, the exchanges of the mass flux between river network and estuary require the researchers to consider these two systems as a whole. Therefore, this paper provides a mathematical method to connect the 1D river networks with the 3D estuarine model to study the flow and suspended sediment transport in this region. Based on the classic Saint-Venant equations and the nonequilibrium transport equation, 1D flow and suspended sediment model is established. Based on the orthogonal curve coordinate and the ECOM model, 3D estuarine model is established to describe the complex horizontal and vertical variability of the flow and suspended sediment transport. The relationship expression of water level at interfaces of 1D model and 3D model is strictly deduced by utilizing the relations between the water level and discharge at junctions. The 1D and 3D hydrodynamics models are then coupled by iterative computations. Finally, the 1D and 3D suspended sediment models are also coupled by adding two additional conditions. The advantage of this coupled method is that the schematisations do not overlap at the interface of the two systems, which can improve the stability of modelling.
The coupled model has been applied successful in the Pearl River network and its estuary, which is a complicated system. The simulated results are compared with the field measurements to evaluate the accuracy of the coupled model. The results show that the coupled model is capable of simulating water levels, discharges, and the suspended sediment concentration in the river network, as well as the velocity and suspended sediment at the shallow estuary. Generally speaking, the coupled model is elaborate enough to simulate water levels, discharges, and velocity in the study area. Due to the limit of the data and the simplification of the model, although the simulated suspended sediment is not as good as water level and discharge, it is still acceptable to investigate the suspended sediment transport in river network and estuary system.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research was financially supported by the “Natural Science Foundation of China” (NSFC, Project nos. 41006046 and 41376094), the “Joint Research Projects NSFC-NWO” (Project no. 51061130545), and the “Commonweal Program of Chinese Ministry of Water Resources” (Project no. 201301072).