Numerical Study of Discharge Adjustment Effects on Reservoir Morphodynamics and Flushing Efﬁciency: An Outlook for the Unazuki Reservoir, Japan

: The Unazuki Reservoir is located on the Kurobe River, which is inﬂuenced by a catchment with one of the highest sediment yields in Japan. Due to a sufﬁciently available discharge during ﬂood events, annual sediment ﬂushing with full water-level drawdown (i.e., free-ﬂow sediment ﬂushing) is conducted to preserve the effective storage capacity of the reservoir. Nevertheless, the upstream half of the reservoir (i.e., study segment) suffers from the excessive deposition of coarser sediments. Remobilization of these coarser materials and their transportation further downstream of the reservoir is a priority of reservoir owners for sustainable reservoir functions, such as ﬂood-risk management and hydroelectric energy generation. In this paper, an already conducted sediment-ﬂushing operation in the Unazuki Reservoir is simulated, and its effects on sediment scouring from the study segment of the reservoir together with changes in bed morphodynamics are presented. A fully 3D numerical model using the ﬁnite volume approach in combination with a wetting/drying algorithm was utilized to reproduce the hydrodynamics and bed changes using the available onsite data. Afterwards, the effects of discharge adjustment on the morphological bed changes and ﬂushing efﬁciency were analysed in the study segment using an additional single-discharge pulse supplied from upstream reservoirs. Simulation results showed that an approximately 75% increase in the average discharge during the free-ﬂow stage changed the dominant morphological process from deposition into an erosive mode in the study segment. If the increase in discharge reaches up to 100%, the ﬂushed volume of sediments from the target segment can increase 2.9 times compared with the initiation of the erosive mode.


Introduction
Sediment deposition occurs in reservoirs as a consequence of the construction of dams on rivers [1]. Sediment deposition and the subsequent loss of effective storage capacity reduce the life span of dams and have negative effects on reservoir services such as flood control, hydropower generation, irrigation and water supply, resulting in significant economic losses [1][2][3]. Various measures have been used to control excessive sedimentation in reservoirs by means of deposition control and removal of already deposited sediments aimed at increasing their lifespan [1,4]. Amongst practical measures, drawdown flushing (i.e., free-flow sediment flushing) plays an important role in reservoir storage capacity recovery and conservation because it represents an efficient hydraulic technique for sediment removal [5]. Free-flow sediment flushing involves a complete decrease in the water level by opening low-level outlets to temporarily develop accelerated riverine flow (i.e., free flow) within the reservoir. Depending on the geometry of the reservoir, this accelerated flow scours a flushing channel in the deposited sediment ( [2,3,5]. Subsequently, fine and coarse sediments that are usually deposited at the head of the reservoir are remobilized due to the increased bed shear stresses. In Japan, the sediment load coming from catchments is generally high because of the dominant geological and hydrological conditions. Therefore, a significant amount of sediment is transported in Japanese rivers and traps, especially in mountainous areas where a notable number of reservoirs have been constructed. Consequently, the incoming sediment loads into reservoirs fed by these rivers are high. The Kurobe River, located in Toyama Prefecture, is one of the most important rivers in Japan owing to a cascade of five reservoirs along this river and the significant amount of hydropower generated by water stored in these reservoirs. Thus, the reservoir owners are interested in employing coordinated full-drawdown flushing to enhance sediment scouring in the Unazuki Reservoir, the most downstream reservoir on the Kurobe River, and the Dashidaira Reservoir located just upstream of it. A special focus is placed on specific areas that encounter excessive depositional issues to ensure sustainable reservoir functionality (e.g., hydropower generation and flood control) [6]. Supplying sediments, especially finer sediments downstream of the Unazuki Reservoir, can secure the fertility of the existing alluvial fan, which is one of the largest alluvial fans in Japan and effectively contributes to protecting the river course from degradation and associated embankment failure. In addition, supplied fine sediments play a major role in preserving the coastline around the river mouth (i.e., Shimo-Nikawa), which is historically subject to severe wave-induced erosion. However, the amount of coarser sediments transported to downstream reaches of the Unazuki dam is lower than the amount of finer sediments.
Although 1D and 2D numerical models are used extensively to simulate sedimenttransport issues, they may not be able to present the required details in some cases. 1D models are not appropriate for simulating the detailed flow field, and the bed changes notably when the flow condition varies between shallow and deep in different parts of the reservoir and when sandbars appear and move, which usually occur during a flushing operation. Simulation of a complex 3D flow field that includes secondary currents is beyond the ability of 2D and quasi-three-dimensional (i.e., 3D) numerical models utilizing the hydrostatic pressure assumption. However, it is obvious that secondary currents contribute significantly to the natural sediment-transport processes in meandering reservoirs [7]. Therefore, the velocity variation over the flow depth, i.e., helical flow, has a substantial effect on sediment transport (e.g., in channel bends) and makes the use of 3D numerical models inevitable as a step forward in reservoir sediment management studies [6]. Threedimensional (3D) numerical models have been used by various researchers to simulate the flow field and sediment-transportation issues in dam reservoirs [8][9][10][11][12]. Recently, 3D numerical models were used to analyse the sediment dynamics in reservoirs during a flushing operation at both the physical model scale [13][14][15] and prototype scale [6,[16][17][18]. Assessing the performance of 3D numerical models for simulating various sediment management strategies in reservoirs under different boundary conditions will provide valuable outcomes for further development and wider application of these models for integrated sediment management frameworks.
The computational fluid dynamics (CFD) software package SSIIM (Simulation of Sediment Movements in Water Intakes with Multiblock Option, Trondheim, Norway), which was developed at the Norwegian University of Science and Technology (NTNU), was used in this study. SSIIM is a fully 3D numerical model employing Reynolds-averaged Navier-Stokes equations (RANS), which calculates the flow field based on the mass and momentum conservation equations utilizing different turbulence closure approaches. This CFD program was utilized with enhanced grid generation features to simulate the sedimentation and flushing channel evolution in reservoirs, including large reservoirs, which presented satisfying outcomes [19][20][21][22]. Nonetheless, it should be mentioned that simu-lating the sediment-flushing operation in a steep full-scale reservoir such as the Unazuki Reservoir is complex due to the temporal variation in hydraulic boundary conditions and riverbed sediment distribution, even over short distances; therefore, associated studies in the literature are scarce. Recently, Esmaeili et al. (2018 and) conducted a study focusing on potential measures for increasing the sediment-flushing efficiency in the Dashidaira Reservoir located upstream of the Unazuki Reservoir [6,17]. They showed that one of the feasible measures that noticeably affects the reservoir-bed morphology and can potentially increase sediment scouring during the flushing operation is adjusting the operational discharge by using artificially additional discharge during the free-flow condition. However, their study case was focused on a much narrower reservoir with an already fixed location of the main channel, while a number of scouring channels can be observed in the study segment in the Unazuki Reservoir. Because of the wider geometry in this segment of the Unazuki Reservoir compared with that of the Dashidaira Reservoir, the water depth is much shallower, and flushing flow shows a different behaviour that results in different erosional patterns. In addition, the incoming sediment load introduced from the upstream boundary into the Unazuki Reservoir represents a more complex operational condition.
In this study, the SSIIM model was utilized to simulate the 2012 free-flow flushing operation of the Unazuki Reservoir, which suffers from excessive deposition of coarse sediments in the study segment. Sediment deposition continues during coordinated sediment flushing when the flushed sediments from the Dashidaira Reservoir enter the Unazuki Reservoir. Therefore, dam owners were interested in assessing potential approaches (e.g., adjusting the flushing discharge) for mitigating this problem. To achieve this goal, the bed changes index (BCI) was introduced to represent the dominant morphological processes (i.e., erosion and deposition) quantitatively in different subareas in the study segment. Afterwards, the influence of operational discharge adjustment scenarios using additional discharge pulses during the free-flow phase was numerically modelled in the study segment to investigate how it affects bed morphology and whether it can increase the total volume of flushed sediment (TVFS) and flushing efficiency (FE), which is defined as the ratio of the flushed sediment volume to the water-used volume. Thus, the relationship between the flushing efficiency and the amount of water used for flushing was established. It should be noted that the additional discharge can be supplied from upstream reservoirs, and its onsite test was performed tentatively.

Site Description
The origin of the Kurobe River is located at Washiba Mountain with an elevation of 3000 m, located in the northern Japanese Alps, and flows into Toyama Bay in the Japan Sea ( Figure 1a) [23]. The length of the river is 85 km, and its catchment area is approximately 682 km 2 . The bed slope varies between 1% and 20%, which characterizes the river as a steep river. In the catchment area of the Kurobe River, the average annual rainfall and total sediment yield are 4000 mm and 1.4 × 10 6 m 3 /year, respectively, which are both among the highest values in Japan [24]. The completion of the Unazuki gravity dam with a height of 97 m and its multipurpose reservoir occurred in 2000. The main purpose of constructing the Unazuki dam was flood control, whereas water supply and hydroelectricity generation are other functions of the dam. A discharge of 700 m 3 /s can be regulated by the dam when a flood occurs with a peak discharge of 6900 m 3 /s. The gross and effective storage capacities of the Unazuki Reservoir are 27.7 and 12.7 million m 3 , respectively [24]. The mean annual sediment load and the ratio of the total storage to the mean annual runoff in the Unazuki Reservoir are 0.96 × 10 6 m 3 and 0.014, respectively [25]. The Unazuki dam is one of the first dams in Japan constructed with sediment-flushing facilities. This dam is controlled directly by the Ministry of Land, Infrastructure, Transport and Tourism of Japan (i.e., MLIT). To mitigate negative environmental impacts on downstream areas, coordinated sedimentflushing operations have been performed every year since 2001 in the Unazuki Reservoir and another upstream reservoir (i.e., Dashidaira Reservoir, approximately 7 km farther) during the first major flood event in the rainy season (i.e., between June and August). The coordinated sediment-flushing operation during flood events was implemented by authorities since aquatic animals have already adapted to the perturbation caused by natural floods [4]. It should be taken into account that a flood not only provides sufficient energy to transport flood-borne sediments through the reservoir with a lowered water level but also scours previously deposited sediments from the reservoir [6]. The deposited sediments in the study segment of the reservoir could not be eroded effectively during previously conducted sediment-flushing operations, which highlights the necessity for modifying the flushing operation with current arrangements (e.g., adjusting the flushing discharge). Consequently, excessive deposition of coarse-bed materials occurs in this segment and subsequent storage loss causes an increase in flood inundation risks during torrential floods and threatens the sustainable goal of hydroelectricity generation. If coarse-bed materials could be scoured and displaced further downstream towards the dam, a notable chance would be given for a final flushing through the reservoir within the next flushing operations, especially because the flushing facilities of the dam are able to sluice coarse sediments. Therefore, stakeholders are interested in scouring these coarse sediments from the study segment and transporting them into the downstream river reach, which would also be ecologically beneficial through maintaining a more sustainable habitat structure for aquatic animals. Figure 1b illustrates the measured longitudinal profile of the Kurobe River along the Unazuki Reservoir, which clearly reveals continuous bed aggradation while annual sediment-flushing operations are ongoing.

Field Observations and Data Arrangement
The measured bed levels before the flushing operation in June 2012 and the locations of cross-sections A-A to H-H, which were used for further assessment of morphological bed changes caused by the flushing, are shown in Figure 1c. The study segment in the Unazuki Reservoir has been divided into two subareas, namely, Areas I and II. Sampled bed material is coarse in both subareas, but distinctive coarser material can be found in Area I since it is located further upstream in the the reservoir. In addition, the inflow discharge and water-level fluctuations during the 2012 flushing operation are illustrated in Figure 1d. The preliminary drawdown occurred from 7 to 27 h after the start of the operation, and then, the free-flow state occurred between 27 and 34 h. The water-level drawdown magnitude from the beginning of the operation to the free-flow state was approximately 21 m at the dam site. Figure 1e demonstrates that the water level decreased and sandbars appeared within the study segment during the free-flow stage. Moreover, based on the available onsite samples, six sediment sizes.
Ranging between 118.3 mm and 0.37 mm were considered to be representative grain sizes in the bed of the reservoir. Figure 1f illustrates an aerial photograph taken at an early stage after the flushing operation in the study segment in the Unazuki Reservoir.
In addition, Table 1 shows the average sediment-size distributions in different cross-sections prior to the flushing operation [17].

Numerical Model
The fully 3D numerical model SSIIM solves the continuity equation together with RANS equations in three dimensions and within a non-orthogonal grid to compute the water motion for turbulent flows [26]. The Reynolds stress term is modelled using the standard k-ε turbulence model with constant empirical values [27]. The finite-volume approach over a control volume is applied as a discretization method to transform the partial differential equations into algebraic equations. The convection term in the RANS equations is solved using the second-order upwind scheme (i.e., SOU) [26]. Although other schemes can be used in the model, the conducted studies in the literature showed that the application of the SOU scheme can present higher quality results when a complex flow field, such as secondary currents and velocity variation over the flow depth within channel bends, plays a major role [6,21].
The unknown pressure field in the RANS equations is calculated by employing the semi-implicit method for pressure-linked equations (i.e., SIMPLE method) [28]. Consequently, the calculated pressure distribution is non-hydrostatic. An implicit free-water surface algorithm was used in the computations and was based on the pressure gradient between a cell and the neighbouring cells. Then, the free-water surface is computed using the local difference between the water surface elevation in the cell, z p , and the ith neighbour cell, z i [21,29].
where ρ is the water density, g is the acceleration due to gravity, p p is the pressure in the cell and p i is the pressure in the ith neighbour cell computed by the SIMPLE method. Thus, application of Equation (1) represents the elevation difference between two cells as a function of the pressure difference. The equation can be used for each cell and is solved implicitly. Since each cell has eight neighbouring cells and the water level in all cells is unknown in the implicit scheme, an iterative method is used, taking a weighted average of the values of neighbouring cells into account. The weighted average coefficient for each neighbouring cell (a i ) is a function of the Froude number, flow direction and location of neighbouring cells [21]: where Fr is the Froude number, w is the dot product of → r and → u , → r is the direction vector pointing from the centre of a cell to the centre of the neighbour cell aimed at taking into account the upstream/downstream effect and → u is the velocity vector of the cell. This coefficient is then used to discretize the following equation: This implicit and iterative approach is a robust and stable method that can also be used in connection with estimates of sediment transport and morphological bed changes under unsteady flow conditions. The use of an implicit discretization scheme allows large time-step sizes to be employed in the model [11,21].
In the SSIIM model, the grid is non-orthogonal, unstructured and adaptive, and it moves vertically with changes in the bed and free-water surface elevation. During the computations, only the water body is modelled, resulting in a reduction in simulation time. The water surface is recomputed after each time step, and the employed wetting/drying algorithm enables the model to have a varying number of grid cells (e.g., in the vertical and lateral directions) with respect to the water depth in the computational domain [30]. The wetting/drying algorithm causes the cells to dry up and disappear from the computational domain if the water depth is smaller than a user-defined lower boundary. If the water depth becomes greater than the lower boundary in the dry area, this algorithm regenerates a cell. This approach makes it possible to have a dynamic grid that can move in the lateral direction, allowing changes in the bed and water level to be accurately modelled.
The Dirichlet boundary condition (logarithmic velocity distribution) was used for the water inflow, whereas the zero-gradient boundary condition was specified for the water outflow and sediment concentration calculation. For the boundary condition at the bed and walls, where there is no water flux, the empirical wall law introduced by Schlichting (1979) was utilized [31]. Additionally, bed roughness in the form of dunes and ripples is taken into account using an empirical formula introduced by Van Rijn, which employs the characteristics of sediment-size distribution and bed form height within the computational domain as follows [32]: where k s is the bed roughness, d 90 is the characteristic sediment size, y is the water depth, and ∆ is the bed form height, which is calculated as follows: where d 50 is the median sediment size, τ is the shear stress, and τ c,i is the critical shear stress for the i-th sediment fraction. The sediment-transport computation for simulating the morphological changes is divided into suspended sediment and bed-load transport. The suspended sediment transport is calculated by solving the transient convection-diffusion equation.
where U j is the water velocity, w is the fall velocity of the sediments, c is the sediment concentration over time t and spatial geometries (i.e., x and z), and Γ T is the turbulent diffusivity coefficient that is calculated by the following equation: where Sc is the Schmidt number representing the ratio of turbulent eddy viscosity (i.e., ν T ) to the diffusion coefficient and is set to 1.0 as the default value in the SSIIM model [26].
To compute the equilibrium suspended sediment concentration, used as the boundary condition in the cells close to the bed, an empirical formula developed by Van Rijn is used [33]: where a is a reference level set equal to the roughness height, d i is the diameter of the i-th fraction, ν is the kinematic viscosity, τ is the shear stress, τ c,i is the critical shear stress for d i , which was calculated from Shield's curve, and ρ w and ρ s are the densities of the water and sediment, respectively. The obtained sediment concentration can be converted to an entrainment rate and used for modelling the resuspension of sediments from the bed. When sediments are prescribed in a bed cell using this formula, the sediment continuity is not satisfied for the cell. Then, the continuity defect can be used to change the bed elevation. Additionally, time-dependent computation can model how the bed geometry changes as a function of sediment erosion and deposition [19,29]. In the SSIIM model, the bed load can be simulated using the Van Rijn formula or alternatively by the Meyer-Peter-Müller (MPM) formula [34,35]. The MPM formula is more appropriate for steep rivers, which mainly transport coarse sediments close to the bed, which is the case for the Kurobe River.
As in the SSIIM program, the k-ε turbulence model is used, and the bed shear stress is calculated as follows [36]: where τ bed is the bed shear stress, c µ is a constant in the k-ε turbulence model that equals 0.09, ρ is the water density, and k is the turbulent kinetic energy. For computing the critical bed shear stress, the particle Reynolds number is calculated using each sediment fraction. The particle Reynolds number is then used to extract the critical bed shear stress (i.e., τ c ) through the Shields curve. When the existing bed shear stress exceeds the critical bed shear stress, the erosional process is initiated. Subsequently, vertical bed-level changes are calculated with respect to the continuity defect in the cell close to the bed: where δz 0 is the vertical bed-level movement, A is the bed-cell area in a horizontal plane, c f is a conversion factor between the sediment flux and bed-level movement that is a function of the solid fraction of the bed material, S i is the sediment inflow, and S o is the sediment outflow quantities [36]. The model accounts for side slope effects by utilizing a reduction function of the critical shear stress for incipient motion by introducing the formula of Brooks together with a sand slide algorithm [12,37]. The sand slide algorithm corrects the bed slope if it exceeds a defined critical angle of repose of the sediments during excessive erosion and thereby accounts for side bank erosion. In fact, the implemented sand slide algorithm acts as a limiter when erosion continues and the bed slope increases [12,38].
In the model, different types of non-uniform grain-size distributions can be implemented over the computational grid, which enables the model to use the measured bed sample data and nest the values following the computation of a sediment-transportation process [26].

Model Setup and Calibration
The bathymetric surveys performed before the flushing, shown in Figure 1c, were used to make the computational grid of the Unazuki Reservoir. The mesh cell sizes in the streamwise and transversal directions were 10-20 m and 5-10 m, respectively. There were 162 cells in the streamwise direction and 55 cells in the transverse direction. Refining the mesh and decreasing the time step (i.e., 20 s) had a negligible effect on the computed flow velocity field and bed changes within the study segment. The bed material density was supposed to be 2650 kg/m 3 . To introduce the hydrodynamic boundary condition for simulation, the water levels and inflow discharge fluctuations shown in Figure 1d were used. Moreover, a non-uniform bed material size distribution with spatially varying fractions was introduced to the model using the six representative sediment sizes shown in Table 1. Indeed, the computational domain was divided into a number of small segments, and each segment had its own non-uniform sediment-size distribution.
Owing to the good performance of the MPM bed-load, sediment-transport formula for flushing simulations of Alpine reservoirs and the Dashidaira Reservoir with conditions almost similar to those of the Unazuki Reservoir (e.g., steep slopes and a wide variety of sediment-size distributions), this sediment-transport formula was selected for further simulations in the Unazuki Reservoir [6,22]. Because the Unazuki Reservoir is located downstream of the Dashidaira Reservoir and the sediment-flushing operation in the Unazuki Reservoir is connected to the Dashidaira Reservoir, a reference case has been built to conduct the sensitivity analysis and calibration process for the Unazuki Reservoir using the adopted values of the main empirical parameters in the Dashidaira Reservoir in the first stage. Both qualitative and quantitative assessments were utilized during the calibration process. A qualitative assessment was conducted in the first stage to check whether noticeable erosion downstream of Area II could be captured by the 3D model, mainly by using Figure 1f for comparison. A simulation case was considered for second-stage assessment if the mentioned erosion could be observed in the simulated bed topography. In the second stage, a quantitative assessment was conducted to compare the dominant morphological process (i.e., deposition or erosion) in each area with the field observations and to measure the deviation of the simulated bed levels from the measured onsite data. Consequently, the BCI parameter that reveals the dominant morphological process in each target area compared with a reference case is defined as follows [6]: where z i_ms is the measured or simulated bed level after the flushing operation at each grid node and z i_reference is the measured or simulated reference bed level at the corresponding node, which is used for comparison purposes to provide information about erosion or deposition over a specific zone of the computational grid. Furthermore, n is the number of grid nodes considered for comparison purposes. Positive and negative BCI values represent bed-level aggradation and degradation, respectively. In addition, mean absolute error (MAE) of the simulated bed levels from the measured onsite data was calculated.
A case was considered for further assessment if it represented deposition in Areas I (e.g., BCI > 0.05) and erosion in Area II (e.g., BCI < −0.035) as dominant morphological processes and if the calculated MAE for each area was smaller than 1 metre. In the third stage, the updated parameter of each successful case of the second stage was selected and implemented. The estimated incoming sediment load into the Unazuki Reservoir, including the flood-borne sediments and the flushed sediments from the Dashidaira Reservoir, was extracted from the outcomes of an adopted and widely used 1D numerical model for the Kurobe River and employed as the upstream boundary condition [5]. In this way, the sediment supply upstream of the Unazuki Reservoir was taken into account for numerical simulations. The incoming sediment load was introduced as a volumetric sediment concentration for each sediment-size fraction during the entire simulation process; therefore, the sediment load changes proportionally to the water discharge when various scenarios are implemented later. This approach can represent a more realistic condition in terms of sediment load coming into the reservoir since the transported sediment load increases when flood discharge increases. The model calibration process involved sensitivity analysis for roughness, active layer thickness, water content of the bed material, and critical angle of repose. The chosen values were set to 0.5 m, 1 m, 40% and 32 degrees, respectively, for the study segment in the Unazuki Reservoir. In addition, the existing constant in Equation (9) was set to 0.035 instead of 0.015, and the lower boundary for water depth in a cell should be considered dry and removed from the computational domain adjusted to 0.3 m. These values could satisfy the aforementioned qualitative and quantitative criteria for the calibration process. Figure 2 shows the computational grid adjustment at the middle of drawdown (i.e., t = 18 h; Figure 2(a1)), during the free-flow conditions (i.e., t = 32 h; Figure 2(a2)), and the corresponding surface water velocity field (i.e., Figure 2(b1,b2)). Figure 2(b1) reveals the gradual shifting of the flushing flow direction from the left bank towards the centre and right bank of the study zone. The onset of flow bifurcation and water flow concentration in two distinctive narrow flushing channels during the drawdown stage (i.e., t = 18 h) can also be seen in this figure. The flow bifurcation in addition to the water-level drawdown contributes to the further development of these two narrow channels. During the free-flow condition (i.e., t = 32 h), velocities increase to 2.5 m/s, and supercritical flows emerge in several zones of the flushing channels that formed (Figure 2(b2)).

Evaluation of Morphological Bed Changes in the Reservoir Caused by the Flushing Operation
When measured bed levels before flushing are used as the reference case, employing measured bed levels after flushing resulted in BCI parameters equal to 0.11 m and −0.07 m in Areas I and II, respectively. BCI values were 0.11 m and -0.1 m, respectively, when simulated bed levels after flushing were employed with the same reference case, which shows a better agreement between the numerical outcomes and measurements in Area I than in Area II in terms of predicting the dominant morphological process [17].
In reality, both sediment deposition and erosion occur in the study segment in the Unazuki Reservoir. The deposition tendency is a result of the significant amount of flushed sediment from the upstream Dashidaira Reservoir, which flowed into the Unazuki Reservoir. While erosion is the dominant process in areas closer to the dam, deposition occurs at the head of the reservoir, where the river enters the reservoir and bed shear stresses and flow velocities decrease as the result of the cross-sectional widening. The comparison of the bathymetric measurements before and after the flushing operation revealed a sediment deposition volume of 19.89 × 10 3 m 3 in the study segment in the Unazuki Reservoir (i.e., areas I and II), while the calculation in the calibrated model showed a sediment deposition volume of 7.22 × 10 3 m 3 in this segment. The accuracy of the simulated average bed changes was 0.62 m in Area I and 0.52 m in Area II. Owing to the existence of wide and shallow sections in the study segment, it was found that different simulated bank erosional patterns in some sections made a major contribution to the mentioned discrepancy between the numerical outcomes and measurements. The inherent uncertainties associated with the estimated incoming sediment load into the Unazuki Reservoir, especially the bed load, in addition to the bathymetric measurements that were not performed immediately after the flushing operation, may exacerbate the discrepancy. As shown in Figure 1f, there are several patternless longitudinal scouring channels in the form of channel networks. These scouring channels appear as the result of channel transitions from the left side towards the centre and later towards the right side of the reservoir over the transition area accompanied by flow bifurcation during the sediment-flushing operation (see Figure 2).
In Figure 3   While almost same operational instructions are employed, the measured bed changes after each flushing event in the study segment is far different than others as the incoming sediment is variant. As can be seen from Figure 3, whereas deposition is dominant in a large portion of the study segment after 2010 and 2012 flushing operations, erosive mode is dominant in the case of 2013 operation. Additionally, this figure illustrates distributed scouring channels without fixed geometrical characteristics (e.g., location, size, depth, etc.), which appear after each flushing operation and highlights the complexity of simulation of bed changes pattern when distinctive flushing channels do not exist. Esmaeili et al. (2017) concluded that the empirical sediment-transport formula (i.e., MPM formula) has a better performance in the areas covered with the coarse materials, which is the case for Area I. Additionally, because of emerging a very shallow flow condition in the study segment, especially during the free-flow stage, flushing flow can be deflected easily due to the small disturbances, as there is no distinctive flushing channel to attract and stabilize the flushing flow [39,40]. As a consequence, a small difference between numerical model outputs at the early stages of simulation and the real operational characteristics can finally lead to noticeably different bed morphology and transported sediment volume compared with the measurements. As the reservoir is steep and characterized by high sediment yield and dynamic bed changes, small floods can easily disturb the morphology of the riverbed in the study segment, which has similarities with braided streams. Therefore, it is difficult to simulate the bed morphology with high accuracy that agrees well with the actual measurements, which are not conducted immediately after the flushing operation [17].
However, extracted TVFS value in each area and BCI parameter for each cross-section in Table 2 show that the model has captured the dominant sediment deposition mode in cross-sections located in Area I and upstream half of Area II and also it could simulate locally a correct tendency for erosion or deposition except in cross-sections G-G. Since the main target of this study is providing an outlook for tracking the governing behaviour of the study segment and quantifying the probable outcomes under different discharge scenarios through a comparative study, the simulated bed topography of 2012 flushing event, adopted as the calibrated case, was used as a reference case for comparison purposes.  Figure 4 demonstrates the measured bed levels before and after the flushing operation and compares them with simulated bed levels of the reference case in selected cross-sections from Area I (i.e., C-C and D-D) and Area II (i.e., F-F and H-H). It also depicts bed variation in the same cross-sections at different stages of the flushing operation (i.e., middle of draw down and also beginning and end of the free-flow stage) to reveal the governing morphodynamic processes in these sections as representative locations for Areas I and II. As can be seen, simulated bed levels at cross-sections C-C and D-D located in Area I and cross-section F-F located in the upstream half of Area II have higher level of agreement with measurements (e.g., Figure 4(a1-c1)) excluding the side bank areas. This trend is observed with a lower agreement for cross-section H-H located in the downstream half of Area II. On the other hand, simulated bed variations in different stages show how bed levels dynamically change within the study segment where local erosive mode can turn into deposition mode and vice versa (e.g., Figure 4(a2,b2)). Figure 5a illustrates the variation in average sediment size (i.e., D 50 ) in different subareas in the study segment before and after the flushing operation. D 50 in different subareas in the study segment before and after the flushing operation was obtained based on the measured field data shown in Table 1 and the simulation results in the calibrated case (i.e., reference case), respectively. The results show that in Area I (i.e., up to crosssection D-D) D 50 decreases, while in Area II, it even increases slightly. This indicates that eroded coarser bed materials from Area I and the incoming sediment load into the Unazuki Reservoir are transported further downstream and deposited in Area II in the study segment because of its wider geometry and subsequent velocity reduction, especially when the water level rises after closing the bottom outlets. Since deposited sediments in Area II still have the chance to be eroded during the free-flow stage where the flow velocity and subsequent driving force are high, the increase in the D 50 of Area II is not considerable. Figure 5b shows the flushed sediments from bottom outlets during an annual flushing operation.

Discharge Adjustment Effect on Morphological Bed Changes and Flushing Efficiency
One of the feasible scenarios in the Unazuki Reservoir is the introduction of an additional discharge pulse during the free-flow stage of the flushing operation (i.e., the ADF scenario). An additional single-discharge pulse could be supplied from reservoirs located upstream of the Unazuki Reservoir. This additional discharge could enhance sediment entrainment by increasing the average flow velocity and induce bed shear stress in the reservoir. The ADF scenario is currently applicable, and preliminary tests for this scenario are being implemented in the Dashidaira Reservoir located upstream of the Unazuki Reservoir [6]. Figure 6a shows the single-discharge pulses that could be supplied in different ADF scenarios. According to the literature on the narrower Dashidaira Reservoir, another scenario was also tested by utilizing pulses of discharge during the free-flow condition (i.e., PDF scenario). In this scenario, the additional water used was introduced in the form of two discharge pulses, with the higher discharge pulse located in the second half of the free-flow stage, as illustrated in Figure 6b. This scenario was also tested to clarify how this method affects the bed morphology and whether it can outperform the ADF scenario for developing larger scouring areas in the study segment, such as that in the Dashidaira Reservoir. Figure 6c shows the discharge rates and original water level (i.e., Q and h, respectively) together with modified discharge rates under different ADF scenarios in the Unazuki Reservoir. Under the different ADF scenarios, the original boundary conditions of the water level (i.e., recorded during the 2012 flushing operation) have been retained, but the original boundary conditions of the discharge during the free-flow stage of the 2012 flushing operation change depending on which ADF scenario is used. For instance, ADF 150 indicates that 150 m 3 /s of additional single-discharge pulse has been added to the original discharges during the free-flow stage of the 2012 flushing operation. The effects of introducing additional single-discharge pulses under various ADF and selected PDF scenarios on the FE and the TVFS are illustrated in Figure 6d. The horizontal axis shows the ratio of average discharge during the free-flow stage using different ADF scenarios (i.e., Q 2 ) to the average discharge during the free-flow stage when no additional discharge is introduced in the reference case (i.e., Q 1 ). FE 2 and FE 1 are the flushing efficiencies when an ADF scenario is employed and when no additional discharge is employed in the reference case, respectively. As shown in Figure 6d, the change from the depositional mode to erosional mode in the study area of the Unazuki reservoir occurs when the average discharge during the free-flow stage increases by approximately 75% (i.e., Q 2 / Q 1 = 1.75). This indicates that the total water volume used for the flushing operation should be increased by approximately 11% in the Unazuki Reservoir. Afterwards, both FE and TVFS increased sharply, according to the increase in the used average discharge during the free-flow stage. Moreover, when Q 2 /Q 1 is approximately 2 under the ADF 230 scenario and the water-used volume increases by approximately 15%, TVFS will increase by approximately 290% compared with the initiation of the erosional mode. Regarding the diagram shown in Figure 6d, when the average discharge during the free-flow stage increased by approximately 140% (i.e., ADF 300 scenario), the FE variation had the highest value. Under these conditions, the total water-used volume for the flushing operation increased by approximately 19%. When the average discharge and the associated driving force of sediments are artificially increased during the free-flow stage using the provided single discharge pulse, the volume of incoming sediment load into the Unazuki Reservoir increases. If the enhanced driving force is not sufficient, it can lead to more sediment deposition and a lower TVFS and FE compared with the case using a smaller discharge pulse (e.g., ADF 270 versus ADF 250 scenarios).
In Figure 7, the calculated bed changes under the ADF 150, ADF 250, and ADF 300 scenarios are depicted. As can be observed from Figure 7a-c, introducing the additional discharge contributed to the mobilization and redeposition of bed materials and incoming sediment load within upstream zones in Area I and the formation of scouring channels in Area II. In the study segment of the reference case, 7.22 × 10 3 m 3 sediment deposition was computed. The sediment deposition volume was approximately 8.2 × 10 3 m 3 under the ADF 150 scenario. Under this scenario, two narrow and deep scouring channels appear along the left and right banks in the study segment, while there is no serious scouring in other places. When the magnitude of the introduced discharge pulse further increases (e.g., under the ADF 250 and ADF 300 scenarios), the geometry of scouring channels along the right and left banks develops, and new scouring channels are established in Area I and over the downstream half of Area II. The magnitude of erosion is generally higher in Area II because this area is initially covered with finer sediments, which can be eroded more easily than coarser sediments under a specified discharge magnitude. Therefore, under the ADF 250 and ADF 300 scenarios (i.e., Figure 7b,c), the TVFS was approximately 45.8 × 10 3 and 55.8 × 10 3 m 3 , respectively. In Figure 8, the simulated bed levels in the reference case and under the ADF 150 and ADF 300 scenarios are plotted in cross-sections A-A, C-C, E-E, and G-G, respectively. The location of these cross-sections can be found in Figure 1c. The figure reveals how the ADF 300 scenario contributes to greater scouring magnitudes by utilizing larger discharge magnitudes compared with the ADF 150 scenario.  To evaluate the morphodynamics of the study segment compared with the reference case here (i.e., simulated bed topography after the 2012 flushing operation) more quantitatively, the BCI parameter was used for subareas between different cross-sections, and the results are presented in Table 3. Table 3 also shows how D 50 varies under different scenarios within different subareas. As can be seen in Table 3, when the discharge pulse during the free-flow stage is not sufficiently high (e.g., under the ADF 150 scenario), sediment erosion only occurs locally at the entrance zone (i.e., before cross-section A-A where BCI is −0.09). Then, the eroded sediments along with transported sediments into the reservoir would be deposited just further downstream of Area I, e.g., in the inner part of the channel bend, where the flow velocity decreases because of the reservoir geometry (i.e., gradual widening of the channel) as shown in Figure 1e. Therefore, local excessive deposition and a decrease in the bed slope occur in Area I. Additionally, the transported sediments are deposited in the upper half of Area II due to the aforementioned velocity reduction. The aforementioned processes can reduce the induced bed shear stress on the deposited sediments, which can further magnify the excessive deposition in this area, as shown in Figure 7a. Positive BCI values in all subareas between cross-sections A-A to F-F are consistent with this figure.
Negative BCI values in two downstream subareas (i.e., between cross-sections F-F and H-H) indicate that flushing operation under the ADF 150 scenario results in local erosion in the mentioned subareas because finer sediments exist there. In addition, average sediment size (i.e., D 50 ) decreases in all subareas compared with the reference case already shown in Figure 5, except in the aforementioned downstream subareas where a small increase in D 50 occurs. This outcome can be attributed to the transportation of coarser bed materials to these two subareas during the water-level drawdown and free-flow stages, which settle there when bottom outlets are closed and the water level rises in the reservoir. When the magnitude of the discharge pulse increases during the free-flow stage (i.e., under the ADF 250 and ADF 300 scenarios), according to the BCI values shown in Table 3, erosion propagates further downstream of Area I (i.e., up to cross-section D-D) but also erosion in the upstream zone of Area I (i.e., segment from reservoir entrance to cross-section B-B) would be much higher than that occurring under the ADF 150 scenario. In the same way, bed-level degradation in Area II would be much more significant compared with that occurring under the ADF 150 scenario. However, because of the higher volume of sediments coming into the reservoir in addition to the higher erosion of coarser bed materials from Area I and their transportation to the downstream subareas, D 50 slightly increases in almost all subareas compared with the reference case. Table 3. Average bed-level changes assessed by the BCI parameter compared with the reference case (i.e., calibrated model for the 2012 flushing operation) and median sediment-size variation in the target subareas in the Unazuki Reservoir under different ADF scenarios and a selected PDF scenario. Negative and positive signs of TVFS indicate deposition and erosion, respectively.

Scenario ADF 150
Segment  While the dominant morphological process within the study segment changes from bed aggradation under the ADF 150 scenario to degradation under the ADF 250 and ADF 300 scenarios, TVFS will increase accordingly compared with the ADF 150 scenario, which can be favourable owing to the larger sediment removal from the study segment of the reservoir. It should also be considered that inflowing discharge during the free-flow stage is concentrated in the already formed patternless and concise scouring channels during the drawdown stage, which mainly contributes to deepening and slow widening of these channels [20]. This deepening and widening procedure is more effective in areas covered by small sediment sizes (e.g., Area II), which is revealed in Table 3.
In terms of the PDF scenarios used, the first discharge pulse (i.e., P1) was set with discharge and duration equal to 300 m 3 /s for 3 h, while the second pulse (i.e., P2) had two variants, e.g., 600 m 3 /s for 2 h and 400 m 3 /s for 3 h (i.e., PDF P1 300 3-P2 600 2 and PDF P1 300 3-P2 400 3). These PDF scenarios were selected among different compositions of discharge pulses since they presented the highest TVFS values, and the total volume of water used was equal to that consumed under the ADF 300 scenario. The numerical outputs in Table 3 show that employing the PDF P1 300 3-P2 600 2 scenario results in a TVFS equal to 55.4 × 10 3 m 3 , which is similar to the TVFS obtained under the ADF 300 scenario, while under the PDF P1 300 3-P2 400 3 scenario, the TVFS was equal to 43.6 × 10 3 m 3 . Thus, the PDF P1 300 3-P2 600 2 scenario was considered for further assessments as an alternative to the ADF 300 scenario. In this case, the extracted BCI parameter reveals that in the entrance area, before cross-section A-A, significant erosion occurs (e.g., BCI is −0.55), while in other subareas in Area I, sediment deposition is the dominant morphological process. Furthermore, considerable erosion occurs in Area II, similar to the ADF 300 scenario. In other words, the main advantage of implementing this scenario is providing sufficient driving force to slightly move forward the eroded coarser materials from the entrance zone of Area I before redeposition in the lower parts of this area, which can be a favourable outcome to some extent. In Figure 9, simulated bed levels at the entrance of Area I and cross-section B-B under the introduced PFD scenario and under the ADF 300 scenario are depicted. This figure quantitatively shows the performance of employing discharge pulses for effectively eroding the coarser bed materials from the entrance area behind cross-section A-A, which is consistent with the results shown in Table 3. Considering Figure 9 in addition to the extracted results in Table 3, one can generally conclude that the application of the PDF scenario did not show significant advantages compared with the ADF 300 scenario in the study segment in the Unazuki Reservoir. Thus, the owner of the Unazuki Reservoir and authorities of the Kurobe River can consider these outcomes for updating sediment management strategies (e.g., starting the operation from the Unazuki Reservoir and increasing the drawdown speed) as a practical step forward to mitigate the sedimentation challenges in this reservoir.

Conclusions
The following outcomes are obtained from the present work:

•
In the study segment in the Unazuki Reservoir, a portion of the flushed sediment load released from the upstream Dashidaira Reservoir in addition to the flood-borne sediments will settle. Even a full water-level drawdown is not able to resuspend and scour these particles effectively, since sufficient driving force cannot be supplied by the ongoing regulations of operational flushing discharges. The study segment in the Unazuki Reservoir is wide, which contributes to developing shallow flow conditions, especially during the free-flow stage. Since there is no distinctive flushing channel to attract and stabilize the flow in this segment, flow deviation across the reservoir occurs during different stages of the flushing operation, which results in the development of a set of distributed scouring channels. The simulated surface flow velocity, bed changes of the reservoir and arrangement of the previously mentioned scouring channels were almost consistent with onsite observations and the aerial photograph taken soon after the flushing operation, which emphasizes the importance of picking onsite data before they are disturbed by the dynamic conditions of the reservoir. However, simulated bed changes showed some differences from bathymetric measurements conducted long after flushing. Thus, measurements of bed changes in the shortest possible time after completing the flushing operation can positively contribute to developing a more accurate model calibration. In this context, detecting the flushing flow directions by recording the entire operation in addition to measuring the real surface velocities can also be useful for a more accurate model calibration.
• If the average discharge during the free-flow condition increases by approximately 75% compared with the reference case, it can contribute to changing the dominant morphological process from deposition to erosion in the study segment of the reservoir, which is a valuable step forward to conducting the flushing operation in a more efficient manner. Subsequently, the total water volume used for flushing purposes increases by approximately 11% in the Unazuki Reservoir. If the single discharge pulse becomes larger under more intensive ADF scenarios (e.g., ADF 250 and ADF 300 scenarios), both TVFS and FE increase sharply due to excessive erosion, especially in the lower half of the study segment covered with finer sediments. In addition, the application of ADF scenarios with a high magnitude discharge pulse (e.g., ADF 300 scenario) can increase the erosion of coarser sediments from the upper half of the study segment, including the reservoir entrance area, and transport them further downstream, which is especially beneficial for remedying the excessive sedimentation within the study segment. Furthermore, this approach effectively contributes to fine sediment removal from the lower half of the study segment. In this context, the effect of providing a longer duration for the free-flow stage with higher magnitudes of discharge pulse can be considered to supply more driving forces aimed at eroding and flushing the coarse sediment sizes for further development of the current study. Outcomes show that discharge adjustment, in the form of two discharge pulses with a higher discharge pulse set in the second half of the free-flow stage, can potentially improve the erosion of coarser materials from the entrance area in the Unazuki Reservoir. Another notable effect would be on erosion of finer sediments located in the lower half of the study segment, similar to the ADF 300 scenario.