Applications of Coupled Explicit–implicit Solution of Swes for Unsteady Flow in Yangtze River

In engineering practice, the unsteady flows generated from the operation of hydropower station in the upstream region could significantly change the navigation system of waterways located in the middle-lower reaches of the river. In order to study the complex propagation, convergence and superposition characteristics of unsteady flows in a long channel with flow confluence, a numerical model based on the coupling of implicit and explicit solution algorithms of Shallow Water Equations (SWEs) has been applied to two large rivers in the reach of Yangtze River, China, which covers the distance from Yibin to Chongqing located upstream side of the Three Gorges Dam. The accuracy of numerical model has been validated by both the steady and unsteady flows using the prototype hydrological data. It is found that the unsteady flows show much more complex water level and discharge behaviors than the steady ones. The studied unsteady flows arising from the water regulation of two upstream hydropower stations could influence the region as far as Zhutuo hydrologic station, which is close to the city of Chongqing. Meanwhile, the computed stage–discharge rating curves at all observation stations demonstrate multi-value loop patterns because of the presence of additional water surface gradient. The present numerical model proves to be robust for simulating complex flows in very long engineering rivers up to 400 km.


Introduction
The Yangtze River ranks first importance in China and it serves as the key navigation route between the east and west part of the country.To build a safe and efficient Yangtze River transportation system is of significant national importance.As a result, several hydropower stations have been built in the upstream region to regulate water flow aiming to improve the navigation status of the middle and lower reaches.After completion of the Three Gorge dam, the navigation conditions between Chongqing and Yichang City have improved considerably, while little flow information has been found upstream of Chongqing.Thus, it would be of practical interest to study the unsteady flow propagation and influence along the river route due to the water regulation of two upstream hydropower stations, which are located at Jinsha River and Min River, respectively.Meanwhile, the sediment transport and alluvial process of Yangtze River are also affected.A site map showing the Yangtze River reach, key gauging stations and two upstream hydropower stations are shown in Figure 1.The unsteady flows demonstrate much more complex phenomena than the steady ones due to the behavior of multi-values in the water surface slope and stage-discharge relationships.Thus, most empirical formula derived from the steady flow assumptions cannot be directly used in the unsteady flow conditions.One breakthrough was made with the establishment of Saint-Venant equations for 1D flow.Later extension of the work to 2D significantly expanded the unsteady flow simulations in a more practical scenario.The advance of computer technology contributed greatly to the development of numerical solution algorithms of Saint-Venant and 2D Shallow Water Equations (SWEs).Meanwhile, various commercial software packages based on the Finite Difference method (FDM), Finite Volume method (FVM) and Finite Element method (FEM) have been widely used for both the research and consultancy purposes, such as HEC-RAS, MIKE 21, FlO-2D, Delft 3D, etc. [1].
There have been consistent documented works on the novel development of SWEs solution algorithms and their practical applications in engineering rivers.For example, one early work in this field could be represented by Ambrosi et al. [2], who performed a numerical simulation on the delta of Po River for a section of 20 km long and 3-day flow cycle.The channel was characterized by complex geometry which consisted of narrow bend and steep gradient in the longitudinal bed profile.Then Kesserwani and Liang [3] presented a second-order Runge-Kutta discontinuous Galerkin scheme for shallow flow simulations involving wetting and drying over complex domain topography.Recently, Liang et al. [4] formulated the 2D SWEs and solved them in an arbitrary curvilinear coordinate system.The model was applied to the study of Malpasset dam-break flood in the Reyran river valley and tidal oscillations in the Severn Estuary and Bristol Channel.Vacondio et al. [5] used a highly robust GPU-parallelized numerical scheme to produce accurate and fast simulations of flooding in large domains with an area of 180 km 2 .Morales-Hernandez et al. [6] proposed a hybrid 1D-2D modeling strategies for fast computation of large rivers flooding for the Tiber River near Rome and they found the computational efficiency was speeded up to 15 times.Similar numerical work was also reported by Chen et al. [7] for a real flood simulation in Jiakouwa flood detention basin, China.Besides, quite a few studies have also been done on the unsteady flow propagation due to the water regulation from upstream hydropower stations.Zhou et al. [8] presented a 2D morphological model for the unsteady flow and multiple grain size sediment transport downstream of the Three Gorges Dam.Pu et al. [9] used a dynamic adaption length approach for sediment transport after dam-break in a prototype regional river.The study was carried out for a 63 km long alluvial river channel in the middle reach of Yangtze River, China.Skublics and Rutschmann [10] used the 2D hydrodynamic model for a 270 km section of Bavarian Danube and compared the river flooding patterns prior to and post the implementation of river regulation measures.More challenging 3D simulations were documented by Chen et al. [11] for a 124 km riverreservoir system from Lewis Smith Dam tailrace to Bankhead Lock & Dam, Alabama, by combining with a thermal temperature model.In recent years, the meshless Smoothed Particle Hydrodynamics (SPH) method has also been developed for simulating the river and flow confluence and it has the The unsteady flows demonstrate much more complex phenomena than the steady ones due to the behavior of multi-values in the water surface slope and stage-discharge relationships.Thus, most empirical formula derived from the steady flow assumptions cannot be directly used in the unsteady flow conditions.One breakthrough was made with the establishment of Saint-Venant equations for 1D flow.Later extension of the work to 2D significantly expanded the unsteady flow simulations in a more practical scenario.The advance of computer technology contributed greatly to the development of numerical solution algorithms of Saint-Venant and 2D Shallow Water Equations (SWEs).Meanwhile, various commercial software packages based on the Finite Difference method (FDM), Finite Volume method (FVM) and Finite Element method (FEM) have been widely used for both the research and consultancy purposes, such as HEC-RAS, MIKE 21, FlO-2D, Delft 3D, etc. [1].
There have been consistent documented works on the novel development of SWEs solution algorithms and their practical applications in engineering rivers.For example, one early work in this field could be represented by Ambrosi et al. [2], who performed a numerical simulation on the delta of Po River for a section of 20 km long and 3-day flow cycle.The channel was characterized by complex geometry which consisted of narrow bend and steep gradient in the longitudinal bed profile.Then Kesserwani and Liang [3] presented a second-order Runge-Kutta discontinuous Galerkin scheme for shallow flow simulations involving wetting and drying over complex domain topography.Recently, Liang et al. [4] formulated the 2D SWEs and solved them in an arbitrary curvilinear coordinate system.The model was applied to the study of Malpasset dam-break flood in the Reyran river valley and tidal oscillations in the Severn Estuary and Bristol Channel.Vacondio et al. [5] used a highly robust GPU-parallelized numerical scheme to produce accurate and fast simulations of flooding in large domains with an area of 180 km 2 .Morales-Hernandez et al. [6] proposed a hybrid 1D-2D modeling strategies for fast computation of large rivers flooding for the Tiber River near Rome and they found the computational efficiency was speeded up to 15 times.Similar numerical work was also reported by Chen et al. [7] for a real flood simulation in Jiakouwa flood detention basin, China.Besides, quite a few studies have also been done on the unsteady flow propagation due to the water regulation from upstream hydropower stations.Zhou et al. [8] presented a 2D morphological model for the unsteady flow and multiple grain size sediment transport downstream of the Three Gorges Dam.Pu et al. [9] used a dynamic adaption length approach for sediment transport after dam-break in a prototype regional river.The study was carried out for a 63 km long alluvial river channel in the middle reach of Yangtze River, China.Skublics and Rutschmann [10] used the 2D hydrodynamic model for a 270 km section of Bavarian Danube and compared the river flooding patterns prior to and post the implementation of river regulation measures.More challenging 3D simulations were documented by Chen et al. [11] for a 124 km river-reservoir system from Lewis Smith Dam tailrace to Bankhead Lock & Dam, Alabama, by combining with a thermal temperature model.In recent years, the meshless Smoothed Particle Hydrodynamics (SPH) method has also been developed for simulating the river and flow confluence and it has the great advantages on computing the mixed super/sub-critical flows.A 1D SPH-SWE approach was proposed by Chang et al. [12] to accurately handle various subcritical and supercritical flow conditions in channel junctions.Kao and Chang [13] adopted a 2D SPH-SWE modeling approach to simulate the dambreak-induced flood and inundation.
Nowadays more and more hydropower stations have been put into operations.Their discharged unsteady flows should generate considerable influence on the downstream navigation, water supply and sediment transport.However, very few studies have been done on very large scale rivers with complex flow confluence, flood wave propagation and superposition.Therefore, the present work aims to use a highly efficient numerical model to investigate the water regulation effect after the completion of Jinsha and Min River hydropower stations.The engineering simulation deals with a practical river reach of 400 km long up to the Chongqing City.The model accuracy will be validated by both the steady and unsteady flow data which were collected on the key gauging stations along the Yangtze River.

Review of Numerical Model and Solution Method
Two general numerical solution schemes are widely used to solve the Saint-Venant Equations and Shallow Water Equations, i.e., explicit and implicit methods.The former is more sensitive to the solution of previous time step and also constrained by the Courant condition, while the latter uses a linked matrix for all the unknowns in the next time step and is more sensitive to the boundary conditions.Although the implicit scheme is unconditionally stable in theory, adequate ratios should be observed between the temporal and spatial steps to achieve most satisfactory performance.Different discretization algorithms such as FDM, FEM or FVM have been used either in the explicit or implicit solution methods.
The proposed coupled explicit-implicit solution method, which will be used to investigate the unsteady flow propagation and confluence in 400 km long Yangtze River reach, has been detailed in a previous Journal of Applied Fluid Mechanics (JAFM) work by Zheng et al. [14] for a dam-break flow induced alluvial process.Here only a brief review of the model is summarized since our main focus in the present study is the robust engineering application of long rivers and complex flow regimes.In summary, the following Shallow Water Equations (SWEs) are used as the governing principle: where η = water surface (m); t = time (s); h = flow depth (m); u = velocity vector (m/s); n = bed roughness (s/m 1/3 ); g = gravitational acceleration (m/s 2 ); and A H = horizontal eddy viscosity coefficient (m 2 /s).Following Zheng et al. [14], a coupled solution algorithm combining the implicit Eulerian-Lagrangian method (ELM) and explicit Characteristics-based Splitting method (CBS) can be used to accurately and efficiently solve the SWEs for unsteady flows in an engineering scale.The coupled solution scheme first uses the ELM approach from n + 1 time step along the streamline to n time step in a reverse tracking order.To improve the tracking, an explicit CBS method should be used within each computational cycle to couple with the implicit algorithm to refine the water level and velocity field.In a practical computation, at each time step the temporary flow field u* (m/s) is first obtained by using the implicit ELM and then corrected by using the explicit CBS.The principles of this coupled model are based on the fact that the streamline tracking of implicit computations could lead to the deviation of real and modeled traced lines in complex unsteady and transient flows.Hence, an explicit iterative solver such as CBS has to be used within the intermediate time step to interpolate the flow depth and velocity field.It should be noted here that the explicit process is used within local time step of the correction phase to enhance the iteration efficiency.The numerical schemes of ELM and CBS are summarized as follows.
In the ELM computation, for a spatially-discretized triangular mesh, the momentum equation is solved by the following procedure after the water level at time n + 1 has been obtained where Ω = integration domain; N = shape function; is the velocity flux (m 2 /s); ∆t = time step (s); and 0.5 ≤ θ ≤ 1 is the explicit-implicit weighting coefficient.In Equation ( 3), P and Q are respectively defined as For an advection-diffusion type equation, the traditional solution method based on the Characteristics-Galerkin could lead to complex numerical algorithm and low computational efficiency.Here, a more promising CBS method is adopted to overcome these issues in spite of being comprised on the conditional stability.In a 2D domain it is represented as where k = advection-diffusion coefficient; S = source term; and indexes i and j represent the spatial x and y directions, respectively.One distinct advantage of the above coupled method is that by simply increasing the amount of explicit computation, the information lost in the implicit computation, such as the detailed flow field, can be recovered.The advantage of using the coupled model is that the implicit scheme allows a CFL value of 1-10, while the explicit one only allows a value of 0.2-0.4considering both the computational efficiency and accuracy.To maintain the consistency of coupled algorithm in the temporal domain, the explicit iterations should be bounded by the implicit ones in time [14].Besides, for the implementation of numerical boundary conditions, there are three types of the boundary utilized in this model, which include the open boundary at the flow inlet and outlet, land boundary where both the slip and non-penetration conditions should be satisfied, and moving boundary related to the change of the flow transport.More detailed treatment procedures on the boundary conditions can be found in Chen et al. [7].
Since the model is also applied to the upstream region of the Yangtze River, which mostly demonstrates the characteristics of mountainous rivers with steep bed slope and abrupt change of flow regimes, the numerical scheme is expected to have good capacity to capture the shock waves.This is different from the numerical model of Zheng et al. [14] which was only used for the study of benchmark open channel experiment.Therefore, artificial viscosity has to be added to the CBS solutions to eliminate undesirable numerical oscillations induced by the shock-waves.This concept was early proposed by Von Neumann and Richtmyer in the middle of twentieth century.By using the artificial viscosity coefficient µ a (kg/m•s), the flow variable φ in Equation ( 6) should be updated by

Model Applications I: Steady Flow Simulations
In this section, the numerical model is applied and validated by the steady flow data from long-term hydrological series.The computational domain covers the upstream Xiangjia Dam of Jinsha River and Gaochang Hydrological Station of Min River, up to downstream Changshou region of Chongqing City, which spans along the main reaches of the Yangtze River with a total length of about 400 km.The topographic elevations of the study area as shown in Figure 2a are obtained from the field observations.Due to the complexity of the river channel with abrupt enlargement and narrowing in the cross-sections, triangular meshes are used to discretize the simulation domain.The size of the mesh is around 25 m-50 m with a node number of 291,969 and grid number 562,502.The typical mesh distributions at Zhutuo gauging station, which is located immediately upstream of the Chongqing, is shown in Figure 2b.
The Manning's bed roughness n of the computational area was calibrated by using the principle of minimum navigation depth.The basic roughness value starts from a reference level of 0.02-0.05,with an adjustment of −0.1%-0.1% following the change in the flow discharge.This calibration was carried out based on the long-term hydrological observation data of the Yangtze River.As a result, the Manning's bed roughness n distributions in the computational domain from Yibin to Zhutuo gauging stations are obtained and shown in Figure 3. n is also related to the surface roughness k s and both reflect the roughness features of river bed.However, they have quite different dimensions.k s has the unit of m while the unit of n is s/m 1/3 .
There are three types of boundary setups, i.e., flow discharge, water level and land boundaries, respectively.For the upstream flow discharge, two large rivers such as Jinsha River and Min River use the measured data while other tributary rivers use the annual mean flows.The downstream outlet adopts the water level as boundary conditions, for example, the normal and dry-season water levels at Three Gorges Dam area are 175 m and 155 m, respectively.The land boundary is described when the computed water depth is below a threshold minimum value.Then both the slip and non-penetration boundary conditions are imposed on these computational cells.For simplicity, the influence of seepage to and from the groundwater is not considered.The validity of this could be verified by the fact that the total flow discharge between two adjacent gauging stations keeps nearly constant within a certain time period.With regard to the inflow and outflow of nearby tributary rivers, most of them can be reasonably ignored since the flow discharges of these small rivers are several orders below the main Yangtze flow.However, in the model simulations, several larger tributary rivers such as Heng River, Tuo River, He River and Jialing River, have been included.
According to the inflow data provided by the main hydrological stations along the Yangtze River, which are measured in Yibin, Luzhou and Zhutuo at 1400 m 3 /s, 2060 m 3 /s and 2160 m 3 /s, respectively, the steady flow simulations have been carried out for a sufficient long time until the simulated flow conditions become stable.In the computations, the steady flow status was evolved from the initial unsteady flow and flow discharges along the channel were used as the criterion to decide whether the flow has reached the steady status or not.The model was run for several times of the flood wave travelling time from the upstream to the downstream end.Thus, the computed flow discharges along the route remained more or less the same with variations being below 1%.
The computed and measured base water levels are shown in Figure 4 and some representative flow velocity fields are provided at key observation sites in Figure 5, where the highlight area illustrates the location of nearby gauging stations.It is shown from Figure 4 that the maximum difference between the measured and modeled water levels is below 0.12 m, with most of them being below 0.05 m and the relative errors to the local measured flow depth are well within 5%.The good agreement is an indication that the numerical model is accurate and reliable for the long river channel routing thus could be used with confidence for more complex unsteady flow simulations in the next section.(a)   Water 2017, 9, 91 9 of 18

Model Applications II: Unsteady Flow Simulations
To investigate more complex confluence and propagation of the unsteady flows over extremely long river channels such as 400 km long due to the upstream water regulation operations, the numerical simulation would provide a more robust tool as compared with the laboratory experiment, which has both facility unavailability and scaling effect.In this section, the proposed coupled model is used to simulate this prototype process, based on four different inflow discharge scenarios from two upstream hydropower stations, i.e., Xiangjia Dam at Jinsha River and Gaochang at Min River (as shown in Figure 1 or Figure 2a).The first three cases correspond to the low inflow discharges while the last one is of medium-large inflow.
The four case studies were based on the prototype observations of the water levels measured on the upstream location and other hydrological gauging stations along the middle and lower reaches of the Yangtze River.(1): The first case, during 11-15 February 2013, corresponds to an inflow discharge of 1420 m 3 /s-1440 m 3 /s at Jinsha River, where the water level variation is relatively stable, and 705 m 3 /s-1230 m 3 /s at Min River, where the maximum water level variation is around 0.63 m (2): The second case, during 22-26 February, corresponds to an inflow discharge of 1380 m 3 /s-1540 m 3 /s at Jinsha River, where the maximum water level variation is 0.15 m, and 800 m 3 /s-1640 m 3 /s at Min River, where the variation is 0.99 m (3): The third case, during 26 February-3 March, corresponds to an inflow discharge of 1400 m 3 /s-1930 m 3 /s at Jinsha River, where the maximum water level variation is 0.84 m, and 765 m 3 /s-1590 m 3 /s at Min River, where the variation is 0.91 m (4): The last case, during 20-27 October 2014, corresponds to an inflow discharge of 2900 m 3 /s-5860 m 3 /s at Jinsha River, where the maximum unsteady flow level variation is 2.25 m, and 1700 m 3 /s-3180 m 3 /s at Min River, where the maximum variation is 1.05 m.For all four cases, the upstream water regulations by the hydropower station generated significant effect in the far reaching downstream areas.The detailed hydrographs recorded on both upstream control sections and several downstream gauging stations along the Yangtze River channel are shown in Figure 6a-d, respectively.It should be noted that the water levels of six stations in Figure 6 are measured from different reference datum planes as indicated in the figure caption.

Model Applications II: Unsteady Flow Simulations
To investigate more complex confluence and propagation of the unsteady flows over extremely long river channels such as 400 km long due to the upstream water regulation operations, the numerical simulation would provide a more robust tool as compared with the laboratory experiment, which has both facility unavailability and scaling effect.In this section, the proposed coupled model is used to simulate this prototype process, based on four different inflow discharge scenarios from two upstream hydropower stations, i.e., Xiangjia Dam at Jinsha River and Gaochang at Min River (as shown in Figure 1 or Figure 2a).The first three cases correspond to the low inflow discharges while the last one is of medium-large inflow.
The four case studies were based on the prototype observations of the water levels measured on the upstream location and other hydrological gauging stations along the middle and lower reaches of the Yangtze River.From Figure 6, a systematic analysis on the unsteady flood waves has been carried out to illustrate the variation of maximum wave height propagation and deformation, as shown in Tables 1-3, respectively, for the different case studies.Meanwhile, the computed wave heights for all these key control sections are also shown for a comparison.The good agreement indicated that the coupled numerical model is capable of capturing the unsteady flood waves and their deformations in an accurate manner.The results in Tables 1-3 suggested that the unsteady flood waves generated from the water regulation of upstream hydropower stations are a propagating wave travelling over a resistant bed.The wave height decreases during its propagation due to the energy loss, while the wave length increased and its multi-peak pattern gradually changed into a single one.Meanwhile, it could be concluded that the flow discharge and water depth also change along the wave propagation.Due to the river storage effect, unsteadiness of the flow tends to be weakened, which is manifested by the fact that wave height decreases with an increase in the propagation distance.Some local wave heights have been increased, but due to the joining of nearby tributary rivers.In general, the flood wave heights demonstrate a decreasing pattern.It can also be observed in Figure 6a,b that when the regulated water levels on uptream Xiangjia Dam Hydropower Station is relatively stable, the unsteady flood waves found at Yibin confluence behave almost the same as the waves found at Gaochang Hydrological Station.Thus, it may be concluded that the upstream water regulation generates small effect on the unsteady flow patterns of further downstream of Zhutuo gauging station.As long as the flood wave height generated on the Xiangjia Dam is well below a threshold value, the river channel downstream of Luzhou should feel much less of its influence.However, in Figure 6c, it can be seen that when the flood wave height reaches 0.84 m at Xiangjia Dam located on Jinsha River and 0.91 m at Gaochang Hydropower Station located on Min River, an obvious unsteady flow pattern of 0.45 m in wave height can still be observed at Zhutuo station.During the dry season, the regulated water level variations at Xiangjia Dam could significantly increase and this would lead to an increase of the unsteady flood waves along the downstream channels, whose effect could propagate further away.In general, the regulations of water on two upstream hydropower stations can greatly affect the downstream flow patterns in the Yangtze River as far as 400 km away from the source location.
For the further validation of numerical model and demonstration of its robustness to simulate long-river unsteady flows, the computed time history of the water level variations at selected key gauging stations (as marked on Figure 1 or Figure 2a) along the Yangtze River is shown in Figure 7a-f.Since there exist significant water level differences between cases 1-3 and case 4, for each station, the comparison has been shown separately for clarity, such as Figure 7(a1,a2).Meanwhile, the field observation data at these gauging stations are also shown for a comparison for the validation purpose.It is shown that the numerical predictions agree quite satisfactorily with the field data in such a challenging application case.However, relatively larger discrepancies are found in the downstream locations such as Hejiang and Zhutuo gauging stations.This could be attributed to the delay in flood wave data measurement in these areas.The maximum error between the measured and computed water levels in the simulations was found to be around 0.25 m, which happens in the region of flow confluence where the water depth reaches 15-20 m due to the converging of two rivers.Here the flows demonstrate strong 3D characteristics and thus the present 2D model could not accurately address this issue.The present study is for the mountainous rivers and there is almost no flat area.Some hydrological gauging stations were built on the slightly flat topographies, where the computational errors are around 0.1 m under the relatively shallow flow depth of 5-8 m.By considering the difficulties in the prototype observation and long-river numerical modeling, this kind of disagreement should be deemed as satisfactory.
Water 2017, 9, 91 12 of 18 much less of its influence.However, in Figure 6c, it can be seen that when the flood wave height reaches 0.84 m at Xiangjia Dam located on Jinsha River and 0.91 m at Gaochang Hydropower Station located on Min River, an obvious unsteady flow pattern of 0.45 m in wave height can still be observed at Zhutuo station.During the dry season, the regulated water level variations at Xiangjia Dam could significantly increase and this would lead to an increase of the unsteady flood waves along the downstream channels, whose effect could propagate further away.In general, the regulations of water on two upstream hydropower stations can greatly affect the downstream flow patterns in the Yangtze River as far as 400 km away from the source location.
For the further validation of numerical model and demonstration of its robustness to simulate long-river unsteady flows, the computed time history of the water level variations at selected key gauging stations (as marked on Figure 1 or Figure 2a) along the Yangtze River is shown in Figure 7a-f.Since there exist significant water level differences between cases 1-3 and case 4, for each station, the comparison has been shown separately for clarity, such as Figure 7(a1),(a2).Meanwhile, the field observation data at these gauging stations are also shown for a comparison for the validation purpose.It is shown that the numerical predictions agree quite satisfactorily with the field data in such a challenging application case.However, relatively larger discrepancies are found in the downstream locations such as Hejiang and Zhutuo gauging stations.This could be attributed to the delay in flood wave data measurement in these areas.The maximum error between the measured and computed water levels in the simulations was found to be around 0.25 m, which happens in the region of flow confluence where the water depth reaches 15-20 m due to the converging of two rivers.
Here the flows demonstrate strong 3D characteristics and thus the present 2D model could not accurately address this issue.The present study is for the mountainous rivers and there is almost no flat area.Some hydrological gauging stations were built on the slightly flat topographies, where the computational errors are around 0.1 m under the relatively shallow flow depth of 5-8 m.By considering the difficulties in the prototype observation and long-river numerical modeling, this kind of disagreement should be deemed as satisfactory.
(a1) (a2) (b1) (b2) To evaluate the efficiency of model in a more quantitative manner, the Nash-Sutcliffe coefficient [15] has been calculated on the unsteady water levels in Figure 7, as shown in Table 4. Besides, a similar analysis was also carried out for the steady flow simulations as shown in Figure 4, and the Nash-Sutcliffe coefficient was calculated to be 1.0.According to Nash and Sutcliffe [15], an efficiency coefficient of 1.0 corresponds to a perfect match between the modeled and observed data.The present agreement in both steady and unsteady flows shows that the coupled Explicit-Implicit solution model can provide satisfactory predication in practical and complex engineering scenarios.The relationships between the water level and flow discharge are determined by quite a few complex hydraulic factors.A simple hydrograph should demonstrate a single correlation of two variables in a prototype observation.However, for the complex unsteady water level and flow discharge relationships, due to the influence of flood waves, they more or less demonstrate a loop pattern.In the following Figures 8 and 9, we will show the hydrograph at Yibin and Zhutuo gauging stations, respectively, for the four different flood events as shown in Figure 6.The relationships between the water level and flow discharge are determined by quite a few complex hydraulic factors.A simple hydrograph should demonstrate a single correlation of two variables in a prototype observation.However, for the complex unsteady water level and flow discharge relationships, due to the influence of flood waves, they more or less demonstrate a loop pattern.In the following Figure 8 and 9, we will show the hydrograph at Yibin and Zhutuo gauging stations, respectively, for the four different flood events as shown in Figure 6.
It is shown that the water level and flow discharge demonstrate an anti-clock loop pattern, encompassing a relatively large band region.This is due to the additional water surface slope generated by the unsteady flows from the discharge of upstream hydropower stations.That is to say, as the upstream flood wave arrives, the water level increases with the flow discharge.At this moment, the downstream water still maintains its original condition.Thus, a positive water surface slope is generated, which is steeper than the original slope under the same flow rate.On the other hand, when the water retreats from the upstream, a reversed negative water surface slope is generated to reduce the original slope under the same flow discharge.Due to the existence of this additional surface slope, the up-wave tends to be compressed to become steeper while the down-wave tends to be elongated to appear milder.At Yibin gauging station as shown in Figure 8 due to the generation of larger additional surface slope, more flow discharge is found during the up-wave stage than the downwave stage under the same water level.Thus, the loop size of hydrograph is considerably larger than that of Zhutuo station, which is shown in Figure 9.The latter can be approximately regarded as a single relationship of the hydro-rating curve.In Figure 8 and Figure 9, the correlation lines are also plotted passing through the center of the loops as an indication of the averaged water surface slope.It has been found out that the unsteady flow discharge is indeed closely related to the water surface slope, i.e., the change rate of the flood tide and ebb tide.A steeper water surface slope indicates a larger gap between the steady and unsteady flow discharges under the same water level.It is shown that the water level and flow discharge demonstrate an anti-clock loop pattern, encompassing a relatively large band region.This is due to the additional water surface slope generated by the unsteady flows from the discharge of upstream hydropower stations.That is to say, as the upstream flood wave arrives, the water level increases with the flow discharge.At this moment, the downstream water still maintains its original condition.Thus, a positive water surface slope is generated, which is steeper than the original slope under the same flow rate.On the other hand, when the water retreats from the upstream, a reversed negative water surface slope is generated to reduce the original slope under the same flow discharge.Due to the existence of this additional surface slope, the up-wave tends to be compressed to become steeper while the down-wave tends to be elongated to appear milder.At Yibin gauging station as shown in Figure 8 due to the generation of larger additional surface slope, more flow discharge is found during the up-wave stage than the down-wave stage under the same water level.Thus, the loop size of hydrograph is considerably larger than that of Zhutuo station, which is shown in Figure 9.The latter can be approximately regarded as a single relationship of the hydro-rating curve.In Figures 8 and 9, the correlation lines are also plotted passing through the center of the loops as an indication of the averaged water surface slope.It has been found out that the unsteady flow discharge is indeed closely related to the water surface slope, i.e., the change rate of the flood tide and ebb tide.A steeper water surface slope indicates a larger gap between the steady and unsteady flow discharges under the same water level.
Finally, for the last case study, as shown in Figure 6, the flow velocity fields are provided in Figure 10a-d on several control stations of interest to have a better understanding on the complex flow regimes in the main channel and the side bank.It can be realized that the operations of upstream hydropower stations and their water regulation processes could generate significant influence on the downstream river flow patterns.
Water 2017, 9, 91 16 of 18 Finally, for the last case study, as shown in Figure 6, the flow velocity fields are provided in Figure 10a-d on several control stations of interest to have a better understanding on the complex flow regimes in the main channel and the side bank.It can be realized that the operations of upstream hydropower stations and their water regulation processes could generate significant influence on the downstream river flow patterns.

Conclusions
A highly efficient numerical model based on the coupled explicit-implicit Eulerian-Lagrangian solution scheme [14] has been applied to simulate the unsteady flood flow propagation and confluence in the middle and lower reaches of the Yangtze River.To improve the shock wave capturing capacity of mountainous rivers, an artificial viscosity has been added to the previous numerical scheme to remove the numerical oscillations and obtain more accurate solutions.The model is applied to a 400 km long prototype river for the examination of the effect of water regulations from the upstream hydropower stations, which has a profound industry background.The numerical results have been validated by both the steady and unsteady flow data observed on the hydrological gauging stations.We have demonstrated through the simulations that in spite of the continuous damping of the flood wave height, its influence can still reach as far as Zhutuo Station, which is 200 km downstream, during the water regulation operations by the two upstream hydropower plants.With the complex process of unsteady flow propagation and confluence, and the additional water surface slope generated by the rising and retreating waves, the relationships between the water level and flow discharge demonstrate strong loop patterns.This phenomenon is more outstanding when the additional surface slope becomes steeper.The present study could provide useful information on the downstream navigation and water supply practice, when the upstream hydropower plants have been built for the water regulation purpose.
All the computations were carried out on a desktop equipped with Intel ® Xeon ® CPU E5-2620 2.00 GHz and it cost 24 h of CPU time for 168 h of flow simulation.

Figure 1 .
Figure 1.Site map of Yangtze River reach, key gauging stations and upstream hydropower stations.

Figure 1 .
Figure 1.Site map of Yangtze River reach, key gauging stations and upstream hydropower stations.

Figure 2 .
Figure 2. (a) Topographic map of computational area; and (b) mesh distributions of Zhutuo gauging station.

Figure 4 .Figure 3 .Figure 3 .
Figure 4. Measured and computed base water levels along the studied river reach, with blue lines indicating the deepest trench topography.

Figure 4 .
Figure 4. Measured and computed base water levels along the studied river reach, with blue lines indicating the deepest trench topography.

Figure 4 .Figure 3 .
Figure 4. Measured and computed base water levels along the studied river reach, with blue lines indicating the deepest trench topography.

Figure 4 .Figure 5 Figure 5 .
Figure 4. Measured and computed base water levels along the studied river reach, with blue lines indicating the deepest trench topography.

Figure 5 .
Figure 5. Computed flow velocity field on key site along the studied river reach, with red squares indicating the nearby gauging stations: (a) Yibin; (b) Luzhou; (c) Hejiang; and (d) Zhutuo.Figure 5. Computed flow velocity field on key site along the studied river reach, with red squares indicating the nearby gauging stations: (a) Yibin; (b) Luzhou; (c) Hejiang; and (d) Zhutuo.
(1): The first case, during 11-15 February 2013, corresponds to an inflow discharge of 1420 m 3 /s-1440 m 3 /s at Jinsha River, where the water level variation is relatively stable, and 705 m 3 /s-1230 m 3 /s at Min River, where the maximum water level variation is around 0.63 m (2):The second case, during 22-26 February, corresponds to an inflow discharge of 1380 m 3 /s-1540 m 3 /s at Jinsha River, where the maximum water level variation is 0.15 m, and 800 m 3 /s-1640 m 3 /s at Min River, where the variation is 0.99 m (3): The third case, during 26 February-3 March, corresponds to an inflow discharge of 1400 m 3 /s-1930 m 3 /s at Jinsha River, where the maximum water level variation is 0.84 m, and 765 m 3 /s-1590 m 3 /s at Min River, where the variation is 0.91 m (4): The last case, during 20-27 October 2014, corresponds to an inflow discharge of 2900 m 3 /s-5860 m 3 /s at Jinsha River, where the maximum unsteady flow level variation is 2.25 m, and 1700 m 3 /s-3180 m 3 /s at Min River, where the maximum variation is 1.05 m.For all four cases, the upstream water regulations by the hydropower station generated significant effect in the far reaching downstream areas.The detailed hydrographs recorded on both upstream control sections and several downstream gauging stations along the Yangtze River channel are shown in Figure6a-d, respectively.It should be noted that the water levels of six stations in Figure6are measured from different reference datum planes as indicated in the figure caption.