Analysis of Storage Effects in the Recirculation Zone Based on the Junction Angle of Channel Conﬂuence

: In this study, numerical simulations using the Environmental Fluid Dynamics Code model were conducted to elucidate the effects of ﬂow structures in the recirculation zone on solute storage based on the junction angle. Numerical simulations were performed at a junction angle of 30 ◦ to 90 ◦ with a momentum ﬂux ratio of 1.62. The simulation results revealed that an increase in the junction angle caused the recirculation zone length and width to increase and strengthened the development of helical motion. The helical motion increased the vertical gradient of the mixing layer and the mixing metric of the dosage curves. The recirculation zone accumulated the solute as a storage zone, which formed a long tail in the concentration curves. The interaction between the helical motion and recirculation zone affected the transverse mixing, such that the transverse dispersion had a positive relationship with the helical motion intensity and a negative relationship with the recirculation zone size. Transverse mixing exhibited an inverse relationship with the mass exchange rate of the recirculation zone. These results indicate that the transverse dispersion is replaced by mixing due to strongly developed storage zones. coefﬁcient and mass exchange rate.


Introduction
The prediction of mixing characteristics of solute inputs in rivers affects decisionmaking management in the operation of hydraulic structures for the safety of local residents. Solute mixing caused by shear flow is commonly analyzed using the advection-dispersion equation based on Fick's law [1]. However, the complex topography of rivers and hydraulic structures such as dams and weirs generates storage zones that induce non-Fickian mixing [2]. Storage zones in rivers mainly consist of surface transient storage in recirculation zones and hyporheic transient storage [3]. Previously, many researchers have used a one-dimensional (1D) storage model to simulate the mixing phenomena in these types of storage zones [4][5][6]. To increase the accuracy of these storage zone models, determination of the model parameters is important. Therefore, many previous studies have conducted parameter calibrations using concentration curves obtained from tracer experiments [7][8][9][10]. However, these calibrated parameters are associated with high uncertainties because the three-dimensional (3D) characteristics of solute mixing in actual tracer tests are translated into results for 1D storage model parameter values.
The inflow of tributaries into main channels is one of the main sources of solute transport. Therefore, analysis of tributary mixing characteristics is important for water quality management. Jung et al. [11] conducted field studies on river confluences to prove that the discharge ratio of the main and tributary channels changes the secondary current intensity, which affects solute transverse mixing. This demonstrates the relevance of the flow characteristics between tributaries and main channels to near-field mixing analysis in confluences. Figure 1 presents a conceptual image of the flow and solute mixing characteristics near a channel confluence. In the channel confluence, a shear layer is formed and a recirculation zone is developed by the decrease in pressure on the inner side of the confluence [12], as shown in Figure 1. A circular motion is formed as a result of discharge from the tributary channel at the outer side of the confluence. The complex characteristics of the flow in the confluence change with the junction angle, discharge ratio between the main and tributary flows, discordance of the channel bed at the entrance to the confluence, geometric symmetry, density contrasts, and potentially, the width-depth ratio [13][14][15][16].
Research on the analysis of flow characteristics in confluences has been conducted based on field measurements, experimental research, and numerical simulations. Schindfessel et al. [17] demonstrated that turbulent kinetic energy increased with an increase in the tributary discharge in the mixing layer developed in the recirculation zone. Yuan et al. [18] reported that helical motion became stronger and elongated when the tributary discharge was greater than the discharge of the main channel. Gualtieri et al. [19] demonstrated that the momentum transfers and the mixing interface width in the confluence hydrodynamic zone depend on the flow conditions. Yuan et al. [20] conducted analysis of the helical motion in confluences and found that flow structures originating from circular cells affect sediment transport. Other researchers have conducted numerical simulations using the Reynolds-Averaged Navier-Stokes (RANS) model to analyze 3D flow characteristics based on field measurements and laboratory experiments [21][22][23][24][25]. Additionally, detailed turbulence structures in the mixing interfaces of confluences have been analyzed through detached eddy simulations using the RANS model [26][27][28][29]. The junction angle influences near-field mixing by creating a recirculation zone and helical cells downstream of the confluence. Near-field mixing at confluences can be divided into rapid mixing, where transverse mixing is completed within a short distance, and slow mixing, where transverse mixing is maintained further downstream [30]. Lane et al. [31] found that the presence of channel-scale circulation at a confluence transforms slow mixing into rapid mixing. Rapid mixing occurs when a tributary flow flows into a main channel with a large junction angle. In such cases, a difference in concentration between the main and tributary channel creates a mixing interface in which the vertical gradient is increased, leading to the stimulation of vertical mixing and rapid completion of transverse mixing. In this context, Ramón et al. [14] stated that the momentum flux ratio and junction angle are important hydraulic factors that influence the location of the mixing interface. Furthermore, weak density contrast at the confluence increases the mixing area by distorting the mixing interface, leading to an increase in transverse mixing [32]. The effect of density contrast on mixing is more dominant farther downstream from the confluence than differences in the velocity of the main and tributary channel [33]. Lewis et al. [16] demonstrated that a tributary flow with strong momentum flux induces helical motion, which moves the mixing interface toward the outer bank, causing transverse mixing to increase. The strong momentum flux affects not only the helical motion but also the flow separation, which causes solute storage. Constantinescu et al. [26] demonstrated that strong Kelvin-Helmholtz instabilities induce vortex pairing in the mixing layer in the recirculation zone, which causes sediment to be actively mixed in the recirculation zone. Despite several studies related to mixing in confluences, it is still difficult to define the effects of storage zones induced by the recirculating flow in confluences on transverse dispersion.
The goal of this study was to analyze the effects of recirculation zones in channel confluences on solute storage. This research focused on the influence of the confluence angle on solute mixing and associated storage effects. The flow structure at a confluence and its effects on solute mixing were simulated using the quasi-three-dimensional model Environmental Fluid Dynamics Code (EFDC). The EFDC model has been used frequently in many riverine, estuary, flow, and solute mixing simulations, and it has been proven to be reasonably accurate [34][35][36][37]. The EFDC was applied to various junction angles to determine their effect on the size of the recirculation zone. An analysis of solute mixing through helical motion and recirculation zones was also conducted. Additionally, the retention time changes of solutes caused by junction angle changes were investigated based on the results of mixing characteristics in confluences. The momentum flux in a tributary channel may cause the flow to detach as it joins the main channel, creating a recirculation zone as a result of low pressure. When two water bodies with the same density are combined into one, the momentum flux ratio is calculated using the following equation:

Methods and
where M r is the momentum flux ratio, Q m and Q t are the flow rates of the main and tributary channels, respectively, and U m and U t are the cross-sectional averaged velocities of the main and tributary channels, respectively. As the momentum ratio and discharge ratio increase, the width and length of the recirculation zone also increase, and the recirculation zone range is affected by junction angle changes. The length and width of the recirculation zone can be predicted based on the discharge ratio between the main and tributary channels, as well as the junction angle. Gurram et al. [38] presented an empirical formula to calculate the width and length of a recirculation zone based on experimental results. The recirculation zone characteristics for subcritical flow can be estimated according to the junction angle as where L s and b s are the length and width of the recirculation zone, respectively, as depicted in Figure 1, α is the junction angle, Q r is the discharge ratio defined as Q m /(Q m + Q t ), and Fr d is the Froude number downstream of the confluence defined as U/ gh 0 in which U is the cross-sectional averaged velocity and h 0 is the downstream water depth. The recirculation zone becomes a storage zone for solutes and affects the mixing characteristics downstream of the junction [26]. Therefore, the accurate prediction of the size of the recirculation zone is important.

Transverse Mixing at a Confluence
Under high-momentum-flux-ratio conditions, the shear plane caused by tributary inflows induces helical motion, as shown in Figure 1. This helical motion influences transverse mixing in confluences, which can be divided into slow mixing or rapid mixing according to the hydraulic characteristics. The width of a solute cloud can be defined as six times the standard deviation of a concentration curve under the assumption of the Fickian mixing [39]. When the concentration changes are less than 5% throughout the channel width, transverse mixing is considered to be complete [40]. When the distance between the solute input source and transverse mixing completion point (L m ) is five to 10 times the width of the channel, it is considered to be rapid mixing [30]. In the case of slow mixing, the distance between the source and completion point is 100 times the channel width [31,41]. The transverse concentration changes caused by tributary inflows can be evaluated using a standard deviation defined as where σ r is the dimensionless mixing metric, and σ us and σ ds are the standard deviations of transverse concentration variations in the upstream and downstream cross-sections, respectively. If the value of σ r is one, then the concentration changes throughout the upstream and downstream widths of the channel are the same, and transverse mixing does not occur. In contrast, if the value of σ r is zero, then transverse mixing is considered to be complete [41]. For a transient solute input from a tributary, transverse mixing can be analyzed by integrating time-varying concentration curves over time and converting them into steadystate conditions. The dosage curve is the time-integrated value of the concentration passing through one point, which can be calculated using the following equations [42]: where θ is the time integration of the concentration passing a point, C is the depth-averaged concentration, Θ is the total dosage in a cross-section, η is the dimensionless cumulative discharge, u and v are the depth-averaged velocity in the x and y-direction, respectively, q is the cumulative discharge, and h is the water depth. The transverse dispersion coefficient can be estimated from the dosage curves shown in Equation (5) using the moment method as where D T is the transverse dispersion coefficient, σ 2 θ is the variance of the dosage curves, and W is the channel width.
The recirculating flow in a confluence can form a storage zone that traps the solute inflow from a tributary. Such a storage zone increases the residence time of solutes and generates a long-tailed concentration curve [43]. The mixing characteristics in rivers change as a result of these storage zones, which affects far-field mixing results. Therefore, many studies have adopted the following storage model to analyze far-field mixing considering storage zones [44]: ∂C s ∂t where C f and C s are the concentrations in the flow and storage zones, respectively, K f is the longitudinal dispersion coefficient for the flow zone, ε = A s /A, A s is the cross-section area of the recirculation zone, A is the entire flow area, k is the mass exchange coefficient, P is the boundary length between flow and recirculation zone, and T is the residence time in the storage zone. After a solute has been trapped in the recirculation zone for a sufficient amount of time, the concentration of the flow zone (C f = 0) becomes zero, and the following equation can be derived [6]: where C s0 is the initial concentration in the storage zone. By using Equation (13), the parameters for the storage zone model can be estimated by analyzing changes in the solute concentration-time curves in the recirculation zone.

Descriptions of EFDC
The hydrodynamic module of the EFDC includes quasi-3D flow analysis and solute transport models. The flow analysis model adopts the RANS equation as a governing equation that considers hydrostatic pressure instead of solving the vertical momentum equation [45]. Although this hydrostatic pressure assumption is limited in terms of reproducing vertical flow, this model, with its quasi-three-dimensional characteristics, was determined to be appropriate for this type of research, even though confluences have vertical flow. This flow analysis model adopts the Smagorinsky model [46] and the Mellor-Yamada 2.5 model [47] as turbulence closure methods for estimating the horizontal and vertical eddy viscosities, respectively. The model solves the governing equations using a combination of finite volume and finite difference scheme on an orthogonal curvilinear coordinate system in the horizontal direction. In the vertical direction, the equations are transformed using the sigma coordinate system. The model employs the mode splitting method, which separates the computation into external and internal mode. In the external mode, free surface elevation is obtained from the depth-averaged governing equations. In the internal mode, velocity shear is computed using the results by the external mode. More details about the numerical formulations can be found in Hamrick [45].
The solute transport model is coupled with the flow analysis model, and the governing equation of the solute transport model is given as where c is the time-averaged concentration, m x and m y are the metric coefficients, H = ζ + h, ζ is the bottom elevation, and ε h and ε z are the horizontal and vertical turbulent diffusion coefficient, respectively. The turbulent diffusion coefficients in Equation (14) are considered as the eddy viscosity according to the Reynolds analogy [40], which has been adopted for simulations without density differences between the main and tributary channels [48,49]. However, it has been reported that simulation results applying the Reynolds analogy may lead to overestimation of the transverse mixing [50].

Grid Sensitivity Analysis Results
Our flow analysis results for channel confluence using the EFDC were evaluated through comparisons to the experimental results presented by Shumate [51], which have also been adopted for validating simulation results at confluences in previous studies [21,23,52]. In this experimental study, the tributary approached the main channel with a 90 • approach angle. The flow velocity and water depth were each measured using Acoustic Doppler Velocimetry (ADV) with ±1% accuracy and a point gauge with 1 mm accuracy. We conducted a numerical simulation for the case with a discharge ratio (Q r ) of 0.25, which is a value that develops a prominent recirculation zone, to determine the mixing characteristics in the recirculation zone. For grid sensitivity analysis, the grid size was uniform, and the number of grids for the channel width (W/dx) is listed in Table 1. The suitability of the grid size was evaluated by calculating the mean absolute percentage error (MAPE) in the water depth and recirculation zone size for depth-averaged flow fields (L s and b s ).  Figure 2 presents the MAPE calculation results based on the grid size. In the case of W/dx = 10, a recirculation zone is not developed in the junction based on the relatively large grid size compared to the other cases, so the error is 100%. In the case of W/dx = 15, the length and width of the recirculation zone have MAPE values of 28.5% and 32.3%, respectively. As W/dx changes to 20 and 25, the recirculation zone length errors are 6.6% and 4.9%, and the width errors are 15.8% and 5.3%, respectively. The changes in water depth decrease significantly in the case of W/dx > 15, and in the case of W/dx = 20, the error ranges from 1.8% to 2.6%. Based on these results, the numerical simulations were conducted using a grid size of W/dx > 20.

Model Settings and Validations
Numerical simulations were conducted to analyze changes in solute mixing with flow and mixing changes stemming from junction angle modifications. The channel used for numerical simulation was designed based on the laboratory experimental channel considered by Shumate [51], and the length scale was increased to 10 times the original size. Through simulations using the scaled-up channel, flow and mixing patterns in small-scale streams, which are often found in urban areas covered with concrete (n = 0.013) with a small width-to-depth ratio, were reproduced. According to the Froude similitude, the channel length, width, and discharge values used in our experiments are listed in Table 2.
The confluence channel used for simulation was set to have the same width of W = 9.14 m for the main and tributary channels in all cases. The downstream boundary water depth (h 0 ) for the numerical simulations was 2.96 m, which conforms to the geometric ratio. The junction angle was changed in 15 • intervals with five cases, from 30 • to 90 • . Based on the results of grid sensitivity analysis, the grid size was set to be W/dx > 20, and the number of vertical layers was set to 10 for the flow and mixing simulations by applying quasi-three-dimensional characteristics. The computational grids are presented in Figure 3. Table 2. Simulation conditions for flow and solute mixing analysis in channel confluence.  . Time step (dt) was set to 0.005 s to satisfy the Courant-Friedrichs-Lewy (CFL) condition. The upstream boundary conditions for the main channel and tributary were set to a constant discharge, and the downstream boundary condition was set to a constant water depth, as listed in Table 2. The discharge of the main and tributary channels was set to Fr d = 0.37 for a discharge ratio (Q r ) of 0.25, such that the numerical simulations would have a separation zone downstream of the junction. Based on these hydraulic conditions, the momentum flux ratio (M r ) in this study was 1.62, which is included in the range for previous studies (0.29-2.52) [16,30]. The free-slip condition was employed for the wall boundary. The kinematic boundary condition was applied for the free surface and bottom boundary, and w becomes 0 at the bottom and water surface.
The EFDC model specifies the shear stress at the bottom cell layer for computation of the internal mode as where τ b x and τ b y are the bottom shear stress; u b and v b are the horizontal velocity component at the bottom cell layer; c b = κ 2 /[ln(∆z/2z * 0 )] 2 is the bottom stress coefficient given assuming a logarithmic velocity profile where ∆z is the bottom cell layer thickness, and z * 0 is the dimensional bottom roughness height [45]. In this study, z * 0 was set to 1 mm. For the solute mixing simulations, a scenario with an accidental spill was considered for the analysis of solute storage effects in the recirculation zone. Therefore, a fully mixed neutrally buoyant solute was introduced for 5 s at a concentration of 100 ppm from the tributary channel. The no-flux boundary condition was applied to all boundaries for the solute mixing simulations.
To verify the accuracy of our simulations using grids following the Froude similitude, the experimental results from Shumate [51] were used for performing water depth and vertical velocity profile comparisons to the simulations, as shown in Figures 4 and 5, respectively. A comparison with the experimental data and simulation data for the water depth for the Case C90 at the points of y/W = 0.17, 0.50, and 0.82 revealed an MAPE of 2.7% to 3.3%. The vertical variations of velocity were properly reproduced compared to the measurements except for the result shown in Figure 5a. The discrepancy at x/W = −1.0, y/W = 0.17 was caused by the hydrostatic assumption, which is insufficient to present the vertical pressure gradient in highly distorted velocity gradient due to the flow separation. Despite the discrepancies, the numerical simulation of the channel with scale enlargement following the Froude similitude adequately reproduced the experimental results.

Effects of Junction Angle on Flow Characteristics
In this study, the reproducibility of the flow structures associated with the junction angle of the designed confluence was evaluated based on numerical simulation results. Figure 6 presents the depth-averaged velocity and streamlines in the recirculation zone developed by the tributary inflow. The inflow of the tributary causes flow separation and the maximum velocity is generated by flow contraction at the right bank. The stagnation zone appears in the upstream junction corner and the velocity at this point decreases. The size of the recirculation zone increases as the junction angle increases, which has been reported in previous studies [12,23,38,53].    (Figure 8b). The results also reveal that helical motion develops strongly as the junction angle increases. Additionally, the strength of the helical motion increases near the left bank (0 < y < 2 m), where the recirculation zone is developed as the junction angle decreases. This indicates that the transverse transport of the solute in the recirculation zone may be weakened as the junction angle increases.
In Figure 9, the length (L s /W) and width (b s /W) of the recirculation zone, based on the depth-averaged velocity fields shown in Figure 6, are compared to the results obtained using Equations (2) and (3), as suggested by Gurram et al. [38]. The values of L s /W and b s /W are significantly reduced in the range of 45 • -75 • in the simulation results. These L s /W changes are similar to those obtained using Equation (2), but in the C30 case, the recirculation zone length is overestimated in the simulation results. The results of b s /W exhibit some underestimation compared to the calculations using Equation (3). These discrepancies in the simulation results could have occurred as a result of the hydrostatic assumption in the EFDC model. However, despite these limitations of the model, the overall results of the numerical simulations only exhibit minor differences compared to the empirical equations based on experiments, and the trend of the recirculation zone area changes with changes in the junction angle was reproduced reasonably accurately.

Solute Transport Results
The effect of junction angle on solute mixing was evaluated. Figure 10 presents the depth-averaged concentration distributions of the simulation results for t = 2.4 min and 4.1 min. When the solute cloud crosses the recirculation zone, a relatively small amount of solute inflows into the recirculation zone compared to the flow zone. The boundary layer created between the recirculation and flow zones disturbs solute trapping by the recirculation zone. Once the solute is trapped in the recirculation zone, it persists there for a relatively long time. Comparisons of the concentration in the recirculation zone revealed that the concentration in the C30 case was the highest after 2.4 min as a result of a relatively weakened boundary layer effect as the junction angle decreases. In contrast, after 4.1 min, the solute in the C90 case is retained the longest, resulting in a concentration that is higher than in the other simulation cases. The C30 case exhibits different mixing characteristics than the C60 and C90 cases because the solute input from the tributary bounces off the right bank, resulting in higher concentrations. In Figure 11, the concentration distributions and streamlines of the v−w velocity, in the same cross-section shown in Figure 7, are depicted for t = 2.4 min. All cases reveal results in which the concentration at the bottom of the channel is relatively higher than that at the upper water surface. These results were obtained because the u velocity is lower at the bottom of the channel, resulting in less solute dispersion compared to the water surface, as shown in Figure 7. Additionally, the clockwise helical motion, which creates a strongly distorted mixing interface, causes rapid mixing with high discharge from the tributary. The distortion of the mixing interface results in an increase in the effects of u w on turbulent diffusion [18]. As the junction angle increases, the vertical gradient of the mixing layer also increases as a result of the strongly developed helical cell. The solute cloud transported near the water surface by helical motion spreads to the outer bank. In contrast, near the channel bottom, the solute cloud moves to the inner bank where the recirculation zone is located. Then, the cloud is stored to accumulate the solute. The transverse mixing of solute inflows from the tributary was analyzed using dosage curves, which were calculated using Equation (5). Figure 12 presents the transverse change in the dimensionless dosage (θ/Θ) plotted versus the dimensionless cumulative discharge (η). The transverse distributions of θ/Θ gradually become flatter toward the downstream cross-section, except for the C30 case. In the C30 case in particular, θ/Θ increases at the right bank as a result of wall reflection. In the C45 case, wall reflection at the right bank causes an increase in θ/Θ from x = −3W. In the C75 and C90 cases, the effects of wall reflection rarely appear, even as far as x = −5W; in these cases, θ/Θ locally increases near the left bank at x = −W and x = −2W as a result of the recirculation zones. In the C90 case, the local increase in θ/Θ is maintained up to x = −3W.
The mixing metrics defined in Equation (4) can be utilized to evaluate transverse mixing. Figure 13 presents the changes in the dimensionless mixing metric (σ r ), which was normalized using the standard deviation of the dosage curve at x = 0 (σ us ), as shown in Figure 12a. Upstream of x = −2W, the value of σ r was not affected by the junction angle, except for the C30 case. The C30 case shows that the deviation of the transverse dosage is minimized as a result of active transverse mixing, causing σ r to decrease quickly. However, as a result of wall reflection, σ r increases downstream of x = −2W, where transverse mixing is completed. These results indicate that as the junction angle decreases, transverse mixing becomes more active, causing a decrease in the transverse concentration deviation. In cases where the junction angle is below 45 • , the tributary flow with high momentum causes the solute cloud to be transported to the far bank. Then, the solute is positioned near the outer bank as a result of wall reflection, resulting in an increase in the transverse concentration deviation.

Storage Effects in the Recirculation Zone
To analyze the mixing properties in the recirculation zone, the concentration-time curves and skewness coefficient in the center of the recirculation zone are plotted in Figure 14. In Figure 14a, the rising limb of the concentration curve increases steeply while the falling limb decreases gradually, which is a typical property of mixing in the storage zones [55]. The strongly developed recirculation zone, which depends on an increase in the junction angle, leads to an increase in the maximum concentration in the storage zone. However, in the C30 case, the maximum concentration is higher than those in the C45 and C60 cases because the recirculation zone area is smaller than in those cases, allowing the passage of a relatively large amount of solute, as shown in Figure 10e. Figure 14b presents the results for the skewness coefficient of the concentration-time curves in Figure 14a. One can see that the skewness coefficient increases significantly as the junction angle increases, especially between 30 • and 60 • . These results demonstrate that the increase in the tail caused by solute trapping is stabilized at junction angles above 60 • .  Figure 15 presents the transverse changes in concentration in cross-sections containing the center of the recirculation zone at t = 2.4 min and t = 4.1 min. At t = 2.4 min, the results reveal that a concentration dip occurs in the recirculation zone, and the lowest concentration tends to increase as the junction angle decreases. In contrast, at t = 4.1 min, the C90 case yields the highest concentration. Additionally, the maximum concentration decreases as the junction angle decreases. These results indicate that the concentration accumulated in the storage zone increases with the junction angle and that the retention time of the solute also increases.

Solute Mixing Downstream of Junction
In this study, the effects of the junction angle on the mixing characteristics at the channel confluence and retention of the solute cloud in the storage zone were analyzed. First, the reproducibility of the flow pattern at the confluence was evaluated by changing the junction angle. Following the experimental research by Gurram et al. [38], the simulation results in this study also demonstrated that as the junction angle increases, the length and width of the recirculation zone increase. Additionally, the simulation results reproduced a clockwise helical cell, which was developed at the outer bank downstream of the confluence. Field measurements by Rhoads and Kenworthy [54] revealed similar results, where one secondary cell developed downstream of the junction under the condition of M r > 1. As the junction angle decreased, the center of the helical cell moved slowly toward the bottom of the channel ( Figure 8) and showed similar results to Shakibainia et al. [23], where numerical simulations revealed a decrease in the u velocity. These previous studies verify the applicability and validity of the simulation results derived in this study.
Previous studies have reported that transverse mixing is related to the momentum flux ratio (M r ). Lewis et al. [16] and Lewis and Rhoads [41] reported that transverse mixing increases with M r . Additionally, the results presented by Pouchoulin et al. [30] support the positive correlation between transverse mixing and M r based on field measurements, where D T /hu * was 1.34, M r = 1.15 and 2.21, and M r = 1.54. However, analysis of the transverse dispersion associated with flow structures by junction angle has been limited in previous studies. In this research, simulation results demonstrated that transverse mixing depends on two main features: helical motion and recirculating flow. Therefore, the relationship between the two main flow features caused by the junction angle and transverse dispersion coefficients is presented in Figure 16. In Figure 16, D T /hu * is calculated using Equation (8), resulting in a range of 0.32-0.82. These results show that D T /hu * increases with the relative intensity of helical motion versus the u velocity 1 W W 0 I dy/u, where u is the depthaveraged u velocity (Figure 16a). In contrast, D T /hu * decreases with an increase in the length (Figure 16b) and width of the recirculation zone (Figure 16c). Since the recirculation zone size changes over flow depth, these results are valid for the recirculation zone size based on the depth-averaged flow field. As the junction angle increases, the helical motion intensity decreases because the increase in the u velocity was more dominant than that in the v velocity as a result of flow contraction caused by the occurrence of a recirculation zone. Therefore, one can see that the interrelationship between the helical motion and recirculation zone caused by a change in the junction angle affects transverse mixing.

Storage Zone Effects on Transverse Mixing
The recirculation zone at a confluence traps the solute cloud and extends the retention time of the solute, as shown in Figure 10. However, the parameters that represent the storage zone characteristics are difficult to obtain from field measurements. In this study, the parameters for the storage zone model were estimated based on simulation results, and the relationships with transverse mixing in confluences were analyzed.
The main parameters for the storage zone model are the residence time (T), area of the developed recirculation zone (A s ), and relative percentage of the zone area among the total area (ε). Table 3 lists the estimated storage zone model parameters, which were obtained from the recirculation zone in the depth-averaged flow field. T was calculated using Equation (13) based on the averaged concentration-time curves in the recirculation zone. A s was determined as the area of the recirculation zone, and ε was calculated from the area of the horizontal plane containing the recirculation zone. These results indicate that all storage zone parameters decrease as the junction angle decreases as a result of a weakened recirculation zone, as shown in Figure 6. When the junction angle decreases from 90 • to 45 • , T gradually decreases, but the decrease is substantial at 30 • . These results indicate that the solute storage effect of the recirculation zone is greatly reduced when the relative area of the storage zone (ε) reaches approximately 11%. In natural streams containing storage zones, 1D longitudinal mixing is achieved by shear flow dispersion and mixing in the storage zone, as shown in Equation (10). Regarding the relationship between shear flow dispersion and storage zone effects, Choi et al. [56] demonstrated a negative correlation between the mass exchange rate (ε/T) and longitudinal dispersion coefficient in the flow zone (K f ). This inverse relationship can be said to reduce the shear flow dispersion, as mixing by the storage zone is strengthened. Additionally, this study revealed that the existence of a recirculation zone affects transverse mixing. Therefore, to discuss the storage effects on transverse mixing due to a recirculation zone, the relationship between ε/T and D T /hu * is plotted in Figure 17. Similar to the relationship between ε/T and K f , D T /hu * also has a negative correlation with ε/T. As in Fischer et al. [40], the longitudinal shear dispersion has an inverse relationship with the transverse dispersion coefficient; in the present study, the decrease in D T /hu * leads to an increase in longitudinal shear dispersion. Therefore, ε/T and D T /hu * simultaneously reduce the longitudinal shear dispersion, and transverse mixing is replaced by mixing in the storage zone when the recirculation zone is strongly developed.

Conclusions
The paper presents the results of a numerical study about the effect of junction angle on flow structure and transverse mixing at a confluence. Flow and mixing analyses in confluences were conducted under conditions of M r = 1.62 and α = 30 • − 90 • to analyze the storage effects caused by recirculation zones. In this study, the following results were obtained:

•
As the junction angle increases, the strength of helical motion and vorticity magnitude increase, leading to an increase in the vertical gradient of the mixing layer and mixing metric of the dosage curves.

•
The recirculation zone creates a boundary layer that interrupts the inflow of solute and simultaneously accumulates solute. These storage effects of the recirculation zone generate long-tailed concentration curves, where the skewness significantly increases at junction angles of 30 • to 60 • .

•
The interaction between the helical motion and recirculation zone affects transverse mixing, where the transverse dispersion coefficient exhibits a positive correlation with the helical motion intensity and negative correlation with the recirculation zone size.

•
The transverse dispersion coefficient has an inverse relationship with the mass exchange rate of the storage zone. Based on these results, it can be inferred that transverse dispersion is replaced by mixing due to strongly developed storage zones.
The aforementioned analysis results are limited to the case of M r = 1.62 and to the time-averaged flow field. Despite the properly reproduced flow simulation results, the solute mixing simulation results have limitations because the results were not completely validated using experimental results. Therefore, additional detailed research on solute and solids mixing due to the storage zone would require a mass exchange rate based on turbulence structures near recirculation zone boundaries.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study will be shared by the authors if requested.