Flat vs. curved rigid-lid LES computations of an open-channel confluence

Confluences are important junctions in networks of natural or human-made open channels, having a significant impact on the upstream water levels and on the degree of mixing of the incoming flows and fluxes of suspended or dissolved matter (Rice et al. 2008). The complex hydrodynamics of confluences has already been investigated intensively, in the field as well as in laboratory conditions. Laboratory tests often take place in schematized confluences of sharp angled flumes with rectangular cross-sections (see e.g. Best 1987; Weber et al. 2001). However, even the high-resolution data collected in experiments did not allow a quantitative characterization of some important, highly three-dimensional flow structures forming in open-channel confluences (Bradbrook et al. 2000; Huang et al. 2002). Supported by the relatively recent increase of available computational power, numerical tools have been important in producing additional data on the study of the confluence hydrodynamics zone (Bradbrook et al. 2000; Huang et al. 2002; Constantinescu et al. 2011). Within the field of Computational Fluid Dynamics (CFD), numerical simulations of these experiments make use of the available experimental data for validation purposes. In turn, validated numerical models are used to further the knowledge on confluence hydrodynamics (Bradbrook et al. 2000; Lane et al. 2000; Constantinescu et al. 2011; Penna et al. 2018). To simulate a turbulent flow in a confluence, different types of numerical modelling strategies are available. Methods based on the Reynolds-Averaged Navier–Stokes (RANS) equations model the effect of all the turbulent scales on the mean motion by using some turbulence model (Rodi et al. 2013), whereas eddy-resolving methods such as Large Eddy Simulations (LES) resolve the larger turbulent scales in space and time on the computational mesh (Rodi et al. 2013). Since LES is more computationally expensive than a RANS approach, wall boundary layers are often not fully resolved. Either a wall function approach is then adopted (Schindfessel et al. 2015) or a RANS model is applied in the near-wall region (such as the Detached Eddy Simulation (DES) approach described by Rodi et al. 2013). Previous studies have pointed out that the complex flow dynamics of a confluence, the mixing layer instability and upwelling events interpreted from these flows, the formation of small Kelvin–Helmholtz vortices in the shear layer between the two flows, and several coherent vortical and upwelling structures, motivate the use of LES (Bradbrook et al. 2000) or unsteady-RANS (URANS) (Luo et al. 2018). In URANS simulations, the turbulence model is adapted to serve LES-like performance in unsteady regions and use RANS capabilities in the relatively stable regions.

Open-channel flow under certain conditions has been found to be characterized by large-scale turbulent structures acting as a secondary component to the mean flow. As Constantinescu et al. (2011) explain, head losses in the confluence of open-channel flows yield an increase of pressure and a super-elevation of the free-surface upstream of the confluence. The super-elevation causes vertical motion of fluid towards the bed. As flow advances to the bed, it moves laterally away from the mixing interface toward each bank and then upward toward the surface. This motion results in a pair of counter-rotating vortices. However, this topic still remains a source of controversy as noted by the present authors in Pouchoulin et al. (2018). The controversy may be partially attributed to the fact that the number of cells and their orientation are influenced by a wide range of factors (flow and momentum ratio, angle of the confluence, bed discordance, cross-sectional ratio, etc.) (Bradbrook et al. 2000; Biron et al. 2002). Some studies (e.g. Rhoads & Kenworthy (1995) and Bradbrook et al. (2000)) show that two counter-rotating cells can easily evolve into a single cell, depending on the discharge ratio. For the flow-ratio of the present study (q ≈ 0.25), as with that reported in Shakibaeinia et al. (2010), and Pouchoulin et al. (2018), the flow pattern in the downstream main channel is more complex than previous studies have shown and it adds further sources for debate.

Many studies have paid great attention to these secondary currents in confluences. Some of them were carried out in the field (see literature review in Rice et al. (2008)) and others in the laboratory such as Weber et al. (2001), whereas Bradbrook et al. (2000) and Shakibaeinia et al. (2010) conducted their investigation on the secondary currents in confluences by means of numerical models.

Schindfessel et al. (2015) have shown the intermittent character of the upwelling events of open-channel confluences. Therefore, despite the computational cost, eddy-resolving methods (such as LES models (Rodi et al. 2013) or URANS models (Luo et al. 2018)) are highly advantageous for the study of such features and to understand the flow patterns and their interaction with the boundaries of the flow (bed, walls, air, etc.). For the treatment of the free-surface, different numerical modelling strategies are also available (McSherry et al. 2017). Interface-capturing methods use a computational domain that does not only comprise the water phase, but also extends into the air phase above the free-surface. The governing flow equations for each of the fluids (water and air) and an additional equation to capture the free-surface are then solved over a fixed mesh. In the case of the Volume of Fluid (VOF) method, the free-surface is determined by means of a transport equation for the volume concentration of the liquid phase. A more recent interface-capturing method is the Level Set Method (LSM), which was originally developed to simulate the motion, in two or three dimensions, of the interface between two phases of a multi-phase flow. VOF-based modelling of confluence flows has been adopted, among others, by Yang et al. (2013) and Luo et al. (2018). Similarly, LSM applications to open-channel flows can be found, for example, in McSherry et al. (2017).

Contrarily to the fixed mesh interface-capturing methods, interface-tracking methods make use of a moving mesh. The computational domain only comprises the water phase, hence only the governing equations of the water flow need to be solved. The computational mesh is moved after every time step, such that the upper boundary coincides with the free-surface geometry. Note that this requires an appropriate algorithm to readjust the mesh in the entire computational domain, each time the surface has moved. An example of such a method in the context of LES is Kara et al. (2015). It is obvious that the abovementioned approaches to treat the free surface in an open-channel flow lead to a significant increase in computational cost, in comparison to the solution of the water flow in a closed-conduit of a similar domain size. Consequently, a popular approach is to introduce a rigid-lid as the upper boundary of the computational domain. The lid is a fixed and impermeable boundary, at some estimated position of the actual free-surface. To prevent the development of spurious boundary layers (in the absence of wind shear), the lid is treated as a frictionless wall. The pressure field onto the rigid-lid, computed either with a RANS-based or a LES-based model, can then be translated into a virtual field of free-surface elevations, by means of the hydrostatic pressure law assumption (McSherry et al. 2017). To simulate open-channel flows which are characterized by small surface deviations compared to a characteristic mean water depth in the domain (approximately 10% of the depth is the suggested threshold value; see, for example, Rodriguez et al. (2004), Constantinescu et al. (2011) and McSherry et al. (2017)), the rigid-lid is most often a flat surface. Since this rigid-lid approximation suppresses the actual surface deformations, a certain error is introduced in the continuity equation. Typically, a constant discharge is imposed, leading to overestimations or underestimations of the bulk velocity in the zones where the rigid-lid is, respectively, lower or higher than the real water level. Yet, for open-channel flows with a sufficiently low (subcritical) Froude number, this approach is often recommended. A discussion on this issue is carried out, for example, in Constantinescu et al. (2011). However, such an approximation not only brings uncertainty in the assessment of the water levels, but also compromises the linkage between the turbulent structures and the free-surface dynamics. This linkage has been studied in the past for open-channel flows (e.g. Kara et al. 2015). Lane et al. (2000) devoted efforts to the study of the aforementioned relation in the specificities of an open-channel confluence.

The main objective of this paper is to assess the possibilities and limitations of the rigid-lid approximation for LES computations of the flow hydrodynamics in a 90° angled confluence of open channels with concordant, horizontal beds and rectangular cross-sections of equal width. LES simulations (with wall functions) will be made of a flow case for which experimental results have been published in the literature (Weber et al. 2001) and for which the surface level deviations and the Froude number meet the abovementioned recommendations in the literature. The first three simulations will be carried out with flat rigid-lids, at different levels above the bed. The fourth simulation will adopt a curved rigid-lid surface which is a better approximation of the actual free-surface. Both the time-averaged flow characteristics and the turbulence characteristics of the simulations will be analyzed. This work complements the simulations with rigid-lid, VOF and dynamic meshes by others, for the same geometry but with a higher resolution and a large eddy resolving model. Due to advances in computer capability, this study aims at simulations that improve agreement with the available experimental data and give further insights on the influence of the rigid-lid approach on the flow field. Brief insights from the selected simulation with the ‘best’ agreement with the experiments will be discussed.

The present study is partially motivated by the fact that simulations performed with a flat and frictionless rigid-lid implementation may not be valid very close to the recirculation zone and the contracted section because, in the latter, the flow undergoes an acceleration that causes the water surface to drop appreciably (±20%, according to Creëlle et al. (2018) and Rice et al. (2008)). The structure of this paper is as follows. The flow case that will be studied in this paper is presented briefly next. The following section describes the adopted numerical methodology. The results from the simulations with flat and curved rigid-lids, both in terms of time-averaged and instantaneous flow fields, are presented in the next section. Wherever possible, results will be compared with the available experimental data. The assessment of the possibilities and limitations of the rigid-lid approach, and insights about the flow patterns gained from the simulation with better overall agreement with the experiments, are discussed in the final section.

Over the past 40 years, a large number of studies (see a recent literature review in Luo et al. (2018)) have been carried out to investigate the flow structures at a river or a man-made canal confluence. Important parameters that govern these structures are the confluence geometry (confluence angle, width ratio of channels, level of bed discordance between channels, etc.) and the momentum ratio of the merging flows.

Best (1987) developed a conceptual model for such a confluence. The model discerns different flow features that are commonly utilized to schematically represent the complex flow patterns: the flow stagnation zone, the recirculation zone (RZ), the flow contraction zone, the flow recovery zone and the mixing layer (Figure 1).

Table 1 summarizes the flow conditions from the experiment that will be simulated in the present paper.

The numerical modelling in the present contribution is conducted within the computational fluid dynamics (CFD) software OpenFOAM, version 4.0. LES will be adopted to simulate the turbulent, incompressible, viscous and three-dimensional flow in the confluence. As a Subgrid Scale Model (SGS), the standard Smagorinsky model is used, with a constant C_s of 0.158 (Rodi et al. 2013). Although more advanced SGS models exist, the standard Smagorinsky model is selected as it has been successfully applied to confluence flows by Bradbrook et al. (2000) and Schindfessel et al. (2015), among others. Likewise, the influence of the SGS model in hydraulic flows is small, compared to that of the boundary conditions or the wall function (Rodi et al. 2013). This SGS model is used in the simulations in combination with a van Driest damping function, which guarantees that the subgrid-scale eddy viscosity tends towards zero in the vicinity of the wall, as it is physically required (Rodi et al. 2013).

The elevations of the free-surface in the considered flow case, as measured by Weber et al. (2001) at three transects along the main channel (y/W = 0.167, y/W = 0.5 and y/W = 0.833), are shown in Figure 2. As mentioned before, the rigid-lid approach will be adopted in this paper to treat the free-surface. The first three simulations will adopt a flat (F) rigid-lid at different elevations above the bed, which will be further referred to as F_DS, F_US and F_RZ. These levels are defined as follows, making use of the average of the measured free-surface elevations at the three transects:

• (a)

  F_DS is based on the average downstream level at x/W = 8;

• (b)

  F_US is based on the average upstream level at x/W = –1;

• (c)

  F_RZ is based on the average level in the recirculation zone, around x/W = 2.

Note that the latter equation implicitly assumes the hydrostatic pressure law to hold. On both the flat and the curved rigid-lids, zero shear stress and zero normal velocity conditions are imposed. In order to maintain the same discharge ratio, the inlet discharges Q_m and Q_t remain unaltered for all four simulations, but the inlet velocities U_m and U_t are adapted according to the respective inlet cross-sections (Table 2).

In the present LES simulations, the wall boundary layers will not be fully resolved, Instead, the wall function approach validated in Schindfessel et al. (2015) will be adopted. Such an approach requires the first node to be located at y⁺ ≈ 30 − 300 (Rodi et al. 2013). Since an LES resolves a large part of the turbulence on the computational mesh, and for the sake of a realistic boundary condition, the inlet velocity should also be turbulent and fully developed. In the present numerical set-up this is achieved by means of a so-called precursor simulation. In this subdomain, a periodic channel is simulated and its turbulent outlet velocity is used as an inlet condition. In order to prevent the need to store the results of precursor simulations, the precursor channel is enclosed within the main domain of the confluence (Figure 3). In the present simulations, this means that the flow simulated at (x/W = –3; y/W = –3), i.e. 2 W downstream of the main channel and tributary inlets of the computational domain, is emulated to the respective inlets at (x/W = –5; y/W = –5), preserving the specified mass flux.

After the precursor domain, there is still a length of 3 W of inlet channel left. This is considered sufficient to capture the upstream flow disturbances induced by the confluence, as described by Creëlle et al. (2018), without these interfering with the precursor simulation. The cyclic precursor channel provokes a certain periodicity to the simulated flow, the frequency of which equals approximately the bulk velocity in the channel divided by the length of the precursor simulation (2 W). For the pressure variable, a zero gradient boundary condition is applied at the inlets (x/W = –5; y/W = –5), the walls and the rigid-lid. At the outlet (x/W = 8), the pressure is set to a constant value, such that in all simulations the virtual water depth (Equation (4)) at the downstream section is always identical to the h_d value specified in Table 1. The boundary conditions for the subgrid-scale eddy viscosity are zero gradient everywhere, except at the walls where the aforementioned wall function is implemented. Due to the use of the wall function, the near-wall turbulence is not fully resolved. Although the near-wall turbulence can be important when predicting the location of flow separation, the sharp-edged confluence investigated in the present contribution makes sure that the point of separation is unequivocal. Consequently, modelling the near-wall flow is up to the standards (Rodi et al. 2013).

For the three flat rigid-lid simulations (F_DS, F_US, F_RZ), the computational meshes have an identical number of cells in the horizontal direction and lateral directions (Table 3). Grading of the cell size is adopted, yielding a higher resolution in the confluence zone. Also, along the vertical direction the cell sizes of the computational meshes are quite comparable, but the number of cells differs, given the different elevations of the flat rigid-lids (Table 2). Note that the aforementioned meshes have been defined after a mesh sensitivity analysis of the results.

The computational mesh for the curved rigid-lid simulation (C_DS) is obtained by deforming the mesh for the F_DS simulation, as illustrated in Figure 4. Along each vertical grid line (i.e. the z direction), the highest grid point will be shifted to coincide with the curved surface defined by Equation (4), whereas the near-wall point will be kept in place in order to maintain the dimensionless wall-normal distance, z⁺, constant and apply the wall function always under the same circumstances. The grid points in between will be gradually redistributed along the vertical gridline. This approach is similar to the one described in Rameshwaran & Naden (2004).

Given the simplicity of the geometry of the computational domain, fulfilling the usual mesh quality parameters is not a complex challenge for the first three simulations (F_DS, F_US, F_RZ). Due to the limited gradients in the curved surface and the gradual redistribution of the points in the mesh deformation, the resulting C_DS mesh also has cell aspect ratios, skewness and mesh non-orthogonality parameters that reach values within the limits suggested by guidelines for a proper LES simulation (Rodi et al. 2013).

In addition, all meshes meet the criterion for the wall function (see above under ‘Other boundary conditions’) to apply, i.e. the wall-normal distance of the first cells, y⁺, is verified for the sidewalls and bed of the numerical domains (30 < y⁺ < 300, see references in Rodi et al. (2013)). Those conditions are met for all four simulations except for the zones of low flow velocity. In principle, the latter is not critical and, moreover, is expected, especially in the stagnation zone (Schindfessel et al. 2015).

Finally, the LES simulations based on the defined meshes were verified to meet the criterion of Pope (2004) and Rodi et al. (2013), stating that the turbulent kinetic energy in the subgrid-scales should not exceed 20%. Figure 5 confirms that for the majority of the grid cells in each of the meshes, more than 80% of the turbulent kinetic energy is indeed resolved on the mesh.

In the OpenFOAM toolbox, the equations governing the LES, i.e. the spatially-averaged Navier–Stokes and continuity equations, are discretized using the Finite Volume Method. Therefore, the LES equations are written in function of the cell-averaged quantities and fluxes through the cell boundaries (Rodi et al. 2013). The discretized equations are coupled and solved using the PISO algorithm. The selected discretization schemes are second order accurate in time and space.

Each simulation comprised an initial period of 825 seconds before the data collection started. This initialization time corresponds to 50 T, where T (⁠13 W/U_d) is a rough estimate of the flow-through time for the 13 W long main channel. Data collection and time-averaging span an additional 100 T of simulation (2,472 s). Each run was computed on a 4 × 8-core Intel E5-2670 (Sandy Bridge @ 2.6 GHz). The total computational cost of a simulation is approximately 4,600 CPU hours. Since the numerical domain is decomposed in 32 sub-domains, the real computational time, due to the parallel processing capabilities, is about 144 hours for each simulation.

The position of the free-surface is not explicitly calculated in any of the simulations, but a virtual time-averaged free-surface position can be deduced from the predicted time-averaged pressure on the rigid-lid, according to Equation (4). For three lateral positions (described in Figure 3 and corresponding to y/W = 0.167, y/W = 0.5 and y/W = 0.833), Figure 6 gives the longitudinal free-surface profiles along the main channel as predicted by the three flat rigid-lid simulations (F_DS, F_US, F_RZ), in comparison to the experimental data by Weber et al. (2001). Note that all predicted free-surface profiles have an identical downstream water depth, because of the imposed pressure boundary condition at the outlet (see above under ‘Other boundary conditions’).

The general shape of all predicted profiles is very similar and agrees qualitatively well with the experiments. Due to the typical head losses induced by the confluence flow (Rice et al. 2008), the water depths upstream are remarkably higher than the water depths downstream of the confluence (approximately 12% for the experimental case of Weber et al. (2001)). In the flow contraction zone, the free-surface profiles experience a local depression due to the flow acceleration. This depression also extends to the recirculation zone, which is adjacent to the flow contraction zone. More downstream, the flow gradually decelerates in the flow recovery zone, which goes along with a rising free-surface.

Despite these similarities, some qualitative and quantitative differences between the predicted and measured free-surface profiles can be discerned. Upstream of the confluence, the F_DS simulation agrees best with the experimental data. The agreement with the measurements in the zone with the free-surface depression is best for the F_RZ simulation in the flow contraction zone (y/W = 0.833), whereas the F_US simulation is superior in the recirculation zone (y/W = 0.167). None of the flat rigid-lid simulations, however, is capable of predicting the superelevation at mid-channel (x/W = 1.8; y/W = 0.5). In the flow recovery zone, all simulations show a satisfactory agreement with the measurements.

In order to parameterize the key free-surface effects in the confluence, several depth ratios have been expressed throughout the literature. The most documented ones are the up- to downstream water depth ratios for the main channel (h_m/h’_d) and for the tributary (h_t/h’_d). These parameters are a quantitative measure for the head losses associated with the confluence hydrodynamics. Table 4 summarizes these ratios for the present case, comparing the experimental with the numerical results (note that the downstream water depth, h’_d, is defined at x/W = 4.66 in Table 4, whereas h_d in Table 1 is defined at x/W = 8). The performance of each simulation, regarding these ratios, is relatively similar.

Based on the foregoing observations, none of the flat rigid-lid simulations is in every aspect superior to the others, when it comes to predictions of the free-surface. Therefore, the free-surface predicted by the F_DS simulation has been chosen, somewhat arbitrarily, to define the curved rigid-lid C_DS (as was already mentioned above under ‘Other boundary conditions’). For the sake of conciseness, only the results of the F_DS and C_DS simulations will be further discussed in this paper.

Figure 7 shows vertical profiles of the time-averaged streamwise velocity component at five locations along the centerline of the main channel. The predictions of the F_DS and C_DS simulations are compared to the acoustic Doppler velocimetry (ADV) measurements of Weber et al. (2001). The scales adopted in the dimensionless presentation of Figure 7 are the cross-sectionally averaged downstream velocity, U_d, and the vertical distance between the bed and the local top boundary (i.e. the rigid-lid for the simulations and the free-surface for the measurements) at a given location, h_ref.

Note that the C_DS simulation produces slightly higher velocities than the F_DS, which is logical since the discharge Q_d is flowing through a reduced cross-section, because of the water surface depression that occurs in the presented locations.

The simulations capture a reduction of velocity magnitude when approaching the free-surface. One may be satisfied by the predictions of the simulation with the flat rigid-lid, but, regarding the region near the surface, slight improvements are made when a curved rigid-lid is adopted (Figure 7).

From analysis of Figure 7, it is not clear if the deformation of the rigid-lid improves the velocity profiles predictions. Nevertheless, and although the flat rigid-lid simulation also predicts a small curvature near the water surface, the simulation with the curved rigid-lid shows that feature in more velocity profiles. This velocity dip phenomenon is a well-known feature of open-channel flow, caused by the presence of secondary flow structures (Nezu & Nakagawa 1993). The latter structures are addressed below under ‘Turbulent structures and secondary flow’.

For some locations in Figure 7 (i.e. b and c), the shape of the predicted profiles in the near-bed zone does not fully agree with the experiments, which may be due to the choice of the wall function.

Figure 8 depicts the contours of the time-averaged velocity magnitude in a plane at z/W = 0.166 (corresponding to z/h_d = 0.5), showing the typical features of an open-channel confluence, already described in Figure 1. Only slight differences between the two simulations (flat vs. curved rigid-lid) are noticeable.

The maximum velocity magnitudes in Figure 8 occur in the flow contraction zone and depend on the extent of the adjacent recirculation zone. The relevant dimensions of the recirculation zone are, traditionally, the width (W_RZ) and the length (L_RZ), as defined in Figure 1. The determination of these dimensions is neither trivial nor unambiguous, because of the complex flow patterns. Consequently, several methods for delineating the recirculation zone are suggested in the literature (e.g. Đorđević (2013); Schindfessel et al. (2015)), yielding significant differences in dimensions. For the present case, only the isovel method is applied to the time-averaged velocity field to predict the dimensions of the recirculation zone in the plane z/W = 0.166. The isovel method delineates this zone by means of the isoline of zero streamwise velocity (Best 1987). The dimensions of the RZ (Table 5) do not seem to depend strongly on the curvature of the rigid-lid.

As expected and previously described in the literature (Đorđević 2013), numerical models with flat rigid-lids underpredict the TKE budget (Figure 9). That underprediction is particularly severe at the points represented in panels d, e, f, g and h of Figure 9. The use of a curved rigid-lid as the top boundary (CDS) visibly increases the TKE values in the flow contraction zone and improves the agreement with the experiments. However, a lack of agreement with the experiments (acquired with a rather low sampling rate of 10 Hz) is still visible, especially, at the points illustrated in Figures 9(e)–9(g). The increment, especially visible at panels f, g and h (Figure 9), is induced by the shear layer at the interface of the two flows originating from the main channel and the tributary. It is highly visible in both simulations, but the increment is higher in the simulation with the curved rigid-lid.

A very important characteristic of confluences is the generation of a mixing interface between the converging flows and the development of large-scale coherent turbulence structures throughout this interface.

However, it seems that the individual tubes that are being generated in the stagnation zone are longer in the CDS simulation compared to FDS.

Besides the turbulent eddies seen in Figure 10, another type of coherent structure affects mass exchange processes: coherent streamwise-oriented vortical cells of helical motion may form adjacent to the mixing interface (as is explained, for example, by Rice et al. (2008)). Those kind of coherent structures are visualized for the F_DS and C_DS simulations, using the Q-criterion for the time-averaged flow in Figure 11. Those structures, seen in the instantaneous plot (Figure 11, right panel), also seem to persist in time, mainly for the simulation C_DS. At least four structures are visible in C_DS (S1, S2, S3 and S4, Figure 11), meanwhile only two of them seem to persist in the time-averaged field of F_DS (S1′ and S2′, Figure 11).

To gain insight into the secondary flow differences between the use of a flat and a curved rigid-lid, Figures 12 and 13 display the cross-section view x/W = 2 of F_DS and C_DS, respectively. These figures show the secondary-circulation (v w) vectors in the cross-section, together with the experimental and numerically predicted water-surface profiles. Through the v, w vector field, circulatory cells (corresponding to the abovementioned streamwise-oriented coherent structures) are visible when the rigid-lid is curved. In the flat rigid simulation, it appears that some recirculation exists, but not so intense or notable as in the simulation C_DS. The secondary flow changes with those different approaches and every simulation underpredicts the intensity of the secondary flow: the vectors v, w are significantly larger in the experiments. It is worth noting that Đorđević (2013) reported results with the same issue: v, w values were underpredicted by at least 50% in the numerical simulations, speculating that the cause of such a large discrepancy was both the rigid-lid and hydrostatic pressure distribution assumptions. The vector field suggests the existence of coherent structures, as was expected by Figure 11. Looking at the experimental data, in the outer half-cross-section (y/W > 0.5) seems to exist only one secondary cell, whereas the numerical data show three, with their respective centre of rotations (y/W, z/W) at (0.55; 0.09), (0.675, 0.285) and (0.95, 0.05). Note that the measurement grid in the experiments did not allow the capture of the small secondary cell near the outer bank (whose centre of rotation is approximately at y/W = 0.96). Likewise, the upper secondary cell may also not be represented, because the experimental data do not reach the water surface (red circles show the free-surface in the experiments). However, it should be noted that the size of that cell may be slightly exaggerated in the numerical simulations because the curved rigid-lid, calculated based on the pressures, is slightly higher than the experimentally observed one. That difference is very small. Thus the authors believe it does not compromise the finding of an upper secondary cell around that location. For the inner half-section (y/W < 0.5) the agreement is good, although a higher spatial resolution of the ADV measurements would be necessary in order to enable a more rigorous comparison.

Firstly, the influence of the position of the flat rigid-lid is briefly investigated. In a general way, and compared to previous studies, the numerically predicted water depth profiles, derived from the pressure field on the lid, all show good agreement with the experimental data. However, the initial level of the flat rigid-lid does have an influence on the prediction of local water levels (Figure 6). In the past (e.g. Đorđević (2013); Schindfessel et al. (2015)), there has been a tendency to choose the downstream level as the rigid-lid position (i.e. F_DS case of the present study). Nevertheless, Figure 6 shows that if, for example, the goal is to obtain the water levels in the vicinity of the downstream corner, the upstream level of the flow is preferable (F_US). The classic position of the flat rigid-lid (F_DS) fails at the prediction of the minimum water depth, for y/W = 0.167, y/W = 0.5 and y/W = 0.833, but studies also show results obtained by VOF presenting the same problem (e.g. Yang et al. (2013); Luo et al. (2018)). Based on these observations, it follows that it is worth making a more careful choice of the elevation of a flat rigid-lid, depending on the region of interest.

The simulation C_DS in the location of the shown profiles (contracted region) overpredicts the streamwise velocity while F_DS slightly underpredicts due to the fact that the flat rigid-lid is higher than the real water level on that location (note that the use of a rigid-lid does induce a mass error). The C_DS computation can better simulate the typical and expected dip in the streamwise velocity profiles, which later proved to be important for a more accurate reproduction of the secondary flow patterns. The changes in the recirculation zone dimensions are not very significant (<10%) between the flat rigid-lid (F_DS) and the curved rigid-lid (C_DS), although it is important to mention that a certain ambiguity is inherent to the methods of calculating the dimensions of such a feature.

Importantly, the results show that the more realistic top boundary has a major influence on the intensity of the TKE in the flow contraction zone. The simulation with the curved lid shows better agreement with the experimental data, as also already shown, for example, by Kara et al. (2015) for the flow around an abutment. Without the deformation of the rigid-lid, the flow depths (hence the cross-sectional area) in the flat rigid-lid simulations are too high and the velocity gradients are not high enough to reach the same high values of TKE. Also, Weber et al. (2001) show that in the flow contraction zone, the TKE values are higher close to the free-surface than to the bed.

The additional turbulence production within the simulation with the curved lid (C_DS) has a visible impact on vortex characteristics (Figure 10). Under the specific flow conditions of the reproduced case, four large-scale rotating cell structures aligned in the longitudinal direction persist in space and time, acting as a secondary motion to the mean flow and merging themselves further downstream (x/W > 5). Not all these features are unequivocally present in simulations by others focused on the same case and coincidently neither in the flat rigid-lid simulations and/or coarser meshes simulations of the present study. By comparison, the simulation F_DS does not exhibit coherent structures of similar organization in the post-confluence zone. This may suggest that this occurs because of the oversimplification of the free-surface. These observations follow the findings about TKE because secondary cells in the downstream channel are typically associated with higher TKE values (Schindfessel et al. 2015). A numerical method capable of fully resolving the water surface will arguably predict the TKE better, because in such a flow case, there are turbulent fluctuations of the water surface, which are closely related to the TKE budget of the flow. Although those vortical structures were identified, it is not clear how they effectively contribute to the turbulent energy in the flow. The visualizations using the Q-criterion (Figures 10 and 11), from which a number of small vortex structures could be identified, should be further quantified in terms of their contributions to the overall TKE. As the flat rigid-lid is adopted in this range of Froude numbers (Fr_d ≈ 0.37), the results have shown that some important flow features are missing or their presence is weaker (e.g. lower TKE, Figure 9).

As opposed to the case C_DS, the F_DS case with a coarser resolution of the mesh (not shown) did not allow the capture of some of the secondary flow features reproduced in the simulation with the curved rigid-lid and a finer mesh (C_DS). Many recent studies have shown the same problem (as reported in Pouchoulin et al. (2018)). This implies that mesh resolution is still a challenge when it comes to reproducing secondary flows and turbulence. During the sensitivity analysis on the mesh refinement, it is noticed that different mesh resolutions (not shown) may produce the same results in terms of streamwise flow but still produce very different secondary flow patterns. Therefore, the validation of numerical simulations of open-channel flows shall be extended further than the classical comparison of streamwise flow velocities and water depths.

Most previous simulations based on the present case suggest contradictory information regarding the orientation of the rotation and/or the number of secondary cells in the flow contraction and recovery zone (as the authors of the present paper already explained in Pouchoulin et al. (2018)). The experimental data collected in Weber et al. (2001) are not sufficiently dense and/or do not reach enough locations to clearly confirm every secondary cell predicted in fine mesh numerical simulations.

To the best knowledge of the authors, among the many numerical studies reproducing the present case of Weber et al. (2001), the secondary flow patterns found (at x/W = 2) in the simulation C_DS of the present study shows one of the better agreements with the experiments, albeit that there are still important discrepancies. Regrettably, these discrepancies are rarely discussed in the literature and deserve further attention, in view of the practical importance of secondary flow for mixing and transport.

This paper is a call for the use of a more realistic free-surface, independently from the numerical method. One could argue that a level-set or a VOF method implementation together with LES would bring a more conclusive study or even get a better agreement with experimental data. Yet, because computational power plays a major role in our predictive capabilities, the adoption of the present method allows an abovementioned objective to be fulfilled: to have a finer mesh, enabling the computation of flow features of interest, which could be more difficult in demanding, complex methods such as VOF. The implemented curved rigid-lid does not mean that re-gridding is required. In contrast to previous research, which was generally restricted either to coarse measuring points grid and/or the effect of a flat rigid-lid, detailed numerical information on three-dimensional mean and turbulent flow structure is provided in the present study. This allows the quantitative characterization of streamwise oriented vortices and how they have been affected by the implementation of flat rigid-lids and/or coarser numerical meshes. These structures are important for practical applications because these vortices will eventually play an important contribution to mixing and to the shear stresses at the bed and the walls.

It is traditionally accepted that the flat rigid-lid approximation is only strictly applicable to low downstream Froude numbers (Fr_d < 0.5; Constantinescu et al. (2011); Schindfessel et al. (2015), etc.) but the findings of the present study suggest that it may be needed to narrow the range Fr_d < 0.5, at least depending on the main interest of the CFD user: while the flat and curved rigid-lid simulations show performances very close to each other on the water surface and mean velocity profiles, the results for the TKE, instantaneous turbulence and secondary flow patterns are significantly different. While Kara et al. (2015) have already stated that the computation of the free-surface of a flow around an abutment may be problematic while dealing with flow with a downstream Froude number equal to the present study (Fr_d ≈ 0.37), higher local Froude numbers (Fr ≈ 0.78) were reached as a result of constriction imposed on the flow by the abutment. In the present study, the highest local Froude number is even lower (Fr ≈ 0.54 at x/W = 2) compared to the Froude number in Kara et al. (2015), which suggests narrowing the range of applicability of the rigid-lid simulations even more. Also Đorđević (2013) mentioned concerns regarding the rigid-lid approach in the simulation of the case of Weber et al. (2001), suggesting this oversimplification as the cause of underpredicted vertical and lateral velocities and TKE by approximately 50%. Further research is needed to fully understand whether and in which situations in cases of Fr_d < 0.37 the rigid-lid approach is valid.

The results from numerical simulations were used to discuss the advantages and limitations of two different ways of using the rigid-lid approach (flat vs. curved) as a simplification of the free-surface within LES computations of an open-channel confluence.

The numerically predicted water depth ratios do not differ significantly, depending on the elevation of each flat rigid-lid. Thus, globally, there is not a single optimal option regarding such an aspect. Instead, the adequacy of the elevation of a flat rigid-lid depends somewhat on the region of interest (i.e. recirculation zone, stagnation zone, contracted flow).

Investigation of the time-averaged flow properties obtained by both flat and curved rigid-lid simulations showed that the agreement with experimental data regarding the velocities of the longitudinal flow is very similar in both cases. This may suggest that the implementation of a curved rigid-lid or a more sophisticated approach for the free-surface simulation is not needed when the main goal is just to obtain the longitudinal flow velocities.

3D instantaneous structures, and especially secondary flow, are significantly different and they lead to significant differences in parameters such as turbulent kinetic energy and in the secondary flow patterns, namely in the size of the large-scale tube-like vortical structures. Although a perfect agreement with the experimental data is not obtained in the present study (or in others, to the best knowledge of the authors), the more realistic rigid-lid (e.g. curved) induces secondary flow patterns very close to those observed in the experiments. This paper suggests that the inaccuracy in the free-surface modelling and/or low resolution result in poor modelling of secondary current cells. Hence, it also represents an appeal to a more realistic approach in the calculation of the free-surface of open channel flows, to the detriment of the simple flat rigid-lid.

As these results and findings provide motivation for further study of the influence of several parameters (junction angle, bed discordance ratio, discharge ratio, etc.) onto the open-channel confluence hydrodynamics, the present method (deformation of the rigid-lid based on the time-averaged pressure field) can be used by others. The findings and final discussion of the present study may also be interpreted as a call for extra efforts for additional experiments in laboratory flumes aiming at a more detailed characterization of the hydrodynamics of a confluence zone.

The authors are not aware of any biases that might be perceived as affecting the objectivity of this research.

This work was performed using the computational facilities of the HPC infrastructure of Ghent University.

Additional details, contours, plots and videos are available at http://bit.ly/2IXiY5o.