Searching for the log law in open channel flow

Abstract We carry out direct numerical simulations of flow in a plane open channel at friction Reynolds number up to ${{Re}}_{\tau } \approx 6000$. We find solid evidence for the presence of universal large-scale organization in the outer layer, with eddies that are larger and stronger than in the closed channel flow. As a result, velocity fluctuations are found to be stronger than in closed channels, throughout the depth. The inner-layer peak of the streamwise velocity variance is observed to grow logarithmically, as in Townsend's attached-eddy model (Townsend, The Structure of Turbulent Shear Flow, 2nd edn, Cambridge University Press, 1976), but saturation of the growth cannot be discarded based on the present data. Although we do not observe a clear outer peak of the streamwise velocity variance, we present substantial evidence that such a peak should emerge at a Reynolds number barely higher than achieved herein. The most striking feature of the flow is the presence of an extended logarithmic layer, with associated Kármán constant asymptoting to $k \approx 0.375$, in line with observations made in shear-free Couette–Poiseuille flow (Coleman et al., Flow Turbul. Combust., vol. 99, issue 3, 2017, pp. 553–564). The virtual absence of a wake region and of corrective terms to the log law in the present flow leads us to conclude that deviations from the log law observed in internal flows are likely due to the effects of the opposing walls, rather than the presence of a driving pressure gradient.


Introduction
Turbulent flows in open channels, namely in the presence of a (nearly) shear-free, approximately flat interface, are relevant in a large number number of engineering and environmental situations, including rivers, lakes and oceanic flows (Nezu 2005;Chaudhry 2007).During the past decades, significant attention has been given to the study of open channel flow, with special emphasis given to processes of sediment transportation (Soldati & Marchioli 2012), and mass and momentum exchanges at the free surface (Lovecchio, Marchioli & Soldati 2013;Mashayekhpour et al. 2019;Pinelli et al. 2022).Open channel flow is also used as a convenient model for the atmospheric surface layer (Nieuwstadt 2005;Flores & Riley 2011), and for the oceanic bottom layer (Taylor 2008).Studying how turbulence is affected by the presence of a free surface is obviously important to understand the complex interactions that take place near gas-liquid interfaces.From a modelling standpoint, access to first-and second-order statistics of turbulence near a free interface is crucial for the improvement of Reynolds-averaged Navier-Stokes turbulence models.In this respect, classical experimental studies (Nezu & Rodi 1986;Kumar, Gupta & Banerjee 1998) have shown that velocity fluctuations in the direction normal to the interface are damped, whereas fluctuations parallel to the free surface are enhanced.However, experiments in open channel flow suffer from serious difficulties in measuring very close to a free surface owing to finite deformation effects, which results in lack of accuracy.
Direct numerical simulations (DNS) have been shown to be very useful to analyse open channel flow, and important contributions were given for instance by Lam & Banerjee (1992), Komori et al. (1993), Campagne et al. (2009) and Ahmed et al. (2021).Large-eddy simulations have also been used for the purpose of studying open channel flow by Salvetti et al. (1997), Taylor, Sarkar & Armenio (2005), Hinterberger, Fröhlich & Rodi (2008) and Ahmed et al. (2021), and especially to study mass transfer across the shear-free interface (Calmet & Magnaudet 2003;Magnaudet & Calmet 2006).However, those studies achieved relatively modest Reynolds numbers, measured in terms of the friction Reynolds number at the no-slip wall, namely Re τ = u τ h/ν, where u τ = (τ w /ρ) 1/2 is the friction velocity, h is the channel depth, and ν is the fluid kinematic viscosity.An important step forward towards realism has been made recently by Yao, Chen & Hussain (2022), who carried out well-resolved DNS of open channel flows at friction Reynolds numbers up to Re τ = 2000.Those authors found that the mean velocity profiles differ from the case of the closed channel in the outer region, and found the presence of a logarithmic layer with Kármán constant k = 0.363.Very-large-scale motions with streamwise wavelength λ x > 3h, or spanwise wavelength λ z > 0.5h, were found to be stronger than in closed channels, resulting in a slightly higher streamwise velocity variance.
An outstanding issue in the wall turbulence community is to what extent the commonly used paradigms are robust, such as the occurrence of a logarithmic layer for the mean velocity profile, and the validity of wall scaling in general.Regarding the first item, whereas most authors agree that a logarithmic layer should arise in all wall-bounded flows at sufficiently high Reynolds number, universality of the log-law constants is still debated (Nagib & Chauhan 2008).Deviations from the assumed log law were discussed by Luchini (2017), who claimed that discrepancies among flows in different (circular or plane) geometries can be ascribed to the effect of the pressure gradient, and that when this effect is accounted for in the form of a higher-order perturbation, universal agreement emerges.Monkewitz (2021) determined the two-term inner and outer asymptotic expansions for closed channel flow, uncovering a change of the logarithmic slope of the mean velocity at a wall distance of several hundred wall units.The slope change was connected with the flow symmetry, and motivated the hypothesis that the breakpoint between the possibly universal short inner logarithmic region and the actual overlap log-law corresponds to the penetration depth of large-scale turbulent structures originating from the opposite wall.In this respect, we believe that the analysis of open channel flow can yield important insight, in that absence of shear at the free surface implies that turbulence kinetic energy production is also zero as at the centre of channels and pipes, however without influence of other walls.Hence we expect the effect of the free surface on the no-slip wall to be much less than in the case of other canonical internal flows.
The distributions of the velocity variances in wall-bounded flows at high Re also suggest violation of strict wall scaling.In fact, it is now well established (Marusic & Monty 2019) that the streamwise (u) and spanwise (w) velocity variances increase slowly with the Reynolds number, with commonly accepted logarithmic growth as in Townsend's attached-eddy model (Townsend 1976).On the other hand, the wall-normal velocity fluctuations seem to level off to a maximum value approximately 1.30 (Smits et al. 2021).It is remarkable that the general growth of the longitudinal and spanwise fluctuations is more evident in the outer layer, and in fact it has long been argued about the possible occurrence of a secondary peak of u 2 , besides the primary buffer-layer peak.Experiments carried out in the Superpipe (Hultmark et al. 2012) and CICLoPE (Willert et al. 2017) facilities indeed support the occurrence of such peak at Re τ 10 4 .It would be extremely interesting to verify whether these similarities also persist at higher Reynolds number, and whether any hint for the onset of an outer peak of the streamwise velocity variance can also be found in open channel flow.
Given this background, our intent here is pushing DNS of flow in open channels to yet higher Reynolds number than achieved so far, in order to test current notions and infer possible infinite-Re extrapolations.Using the DNS solver to be described in the forthcoming sections, we reach Re τ ≈ 6000 (based on friction at the stationary wall), which is a factor three beyond the current state of the art.Important theoretical and practical inferences include assessing the proper form of the log law of the wall, and establishing similarities and differences with other canonical wall-bounded flows.

Methodology
Numerical simulations of fully developed forced turbulent flow in a plane open channel are carried out, assuming periodic boundary conditions in the streamwise (x) and spanwise (z) directions, no-slip boundary conditions at the bottom wall, and zero shear stress boundary conditions at the top boundary (see figure 1).For that purpose, the wall-normal velocity (v) is set to zero, whereas the coefficients for the discrete evaluation of the derivatives in the viscous terms for the streamwise (u) and spanwise (w) velocity components are redefined to have zero wall-normal derivative.The bulk velocity (u b ) is kept strictly constant during the simulations through the use of a time-varying, spatially uniform body force.The single controlling parameter is the bulk Reynolds number Re b = 2hu b /ν, with h the channel depth, and ν the fluid kinematic viscosity.The computer code used for the DNS is based on the classical second-order marker-and-cell method (Harlow & Welch 1965;Orlandi 2000), whereby pressure is located at the cell centres, whereas the velocity components are located at the cell faces, thus removing odd-even decoupling phenomena and guaranteeing discrete conservation of the total kinetic energy in the inviscid limit.The Poisson equation resulting from enforcement of the divergence-free condition is solved efficiently by double trigonometric expansion in the periodic streamwise and spanwise directions, and inversion of tridiagonal matrices in the wall-normal direction (Kim & Moin 1985).An extensive series of previous studies about wall-bounded flows from this group proved that second-order finite-difference discretization in practical cases of wall-bounded turbulence yields results that are by no means inferior in quality to those  of pseudo-spectral methods (e.g.Pirozzoli, Bernardini & Orlandi 2016).The governing equations are advanced in time by means of a hybrid third-order low-storage Runge-Kutta algorithm, whereby the diffusive terms are handled implicitly, and convective terms are handled explicitly.The code was adapted to run on clusters of graphic accelerators (GPUs), using a combination of CUDA Fortran and OpenACC directives, and relying on the CUFFT libraries for efficient execution of fast Fourier transforms (Ruetsch & Fatica 2014;Pirozzoli et al. 2021).
A list of the main simulations that we have carried out is given in table 1.All DNS are carried out in a box with L x = 6πh, L z = 2πh.Extensive preliminary studies of sensitivity to the computational box size and grid spacing have been carried out in a preliminary stage, for flow case C, which have shown that the maximum effect on the one-point statistics presented hereafter is much less than 1 %.This is well in line with what was reported in previous studies (Yao et al. 2022).The mesh resolution for these cases is designed based on the criteria discussed by Pirozzoli & Orlandi (2021).In particular, the collocation points are distributed in the wall-normal direction so that approximately 30 points are placed within y + ≤ 40, with the first grid point at y + < 0.1, and the mesh is stretched progressively in the outer wall layer in such a way that the mesh spacing is proportional to the local Kolmogorov length scale, which there varies as η + ≈ 0.8y + 1/4 (Jiménez 2018).Mild stretching towards the free surface is also used in order to resolve the thin layer in which vertical velocity damping is effective.Based on experience accumulated in a number of previous studies, the grid resolution in the wall-parallel directions is set to Δx + ≈ 8.5, Δz + ≈ 4.0, for all the flow cases.According to the established practice (Hoyas & Jiménez 2006;Lee & Moser 2015), the time intervals used to collect the flow statistics are reported in terms of eddy turnover times, namely h/u τ .It should be noted that such intervals are longer than those used typically in closed channels, as we have verified slower time convergence in open channel flow.
The sampling errors for some key properties discussed in this paper have been estimated using the method of Russo & Luchini (2017), based on extension of the classical batch means approach.Consistent with what we found in pipe flow (Pirozzoli et al. 2021), we find that the sampling and discretization errors are at most 0.3 % for the friction coefficient, 0.5 % for the mean velocity at the free surface, and 0.6 % for the peak velocity variances.
Hereafter, capital letters will be used to denote flow properties averaged in the homogeneous spatial directions and in time, brackets denote the averaging operator, and lower-case letters denote fluctuations from the mean.Instantaneous values will be denoted with tildes, hence ũ = U + u.

Results
We begin by inspecting the instantaneous velocity fields in a cross-stream plane in figure 2. The figure highlights well that, as in all wall-bounded flows, the near-wall region with low-speed fluid becomes thinner as the Reynolds number increases, and finer-scale organization of the turbulent eddies is clearly visible.At the same time, large-sized towering eddies are well visible, which originate from the no-slip wall and extend up to the free surface.Despite differences in the details, the organization of the large-scale eddies seems to be similar, as will be confirmed quantitatively from the velocity spectra.
The flow organization in wall-parallel planes is considered in figure 3, where we show streamwise velocity contours near the no-slip wall and at the free surface, at two Reynolds numbers.The flow near the bottom wall is dominated by streaks of alternating high-and low-speed fluid, which are the hallmark of wall-bounded turbulence (Kline et al. 1967).Organization of the streaks on two different scales is quite apparent, with large-scale streaks whose size is similar at the two Reynolds numbers, and smaller-scale streaks whose size scales in wall units, and which hence become smaller at higher Re.The large-scale streaks are the roots of the towering eddies noted previously.The imprinting of their tops is clear in the free-surface planes, which show the presence of streamwise-elongated eddies, whose spanwise size is similar at the two Reynolds numbers under scrutiny.Detailed information about the large-scale organization of the velocity field and the dynamical contributions of the large eddies can be found in the studies of Duan et al. (2020) and Yao et al. (2022).
Quantitative information about the flow organization can be gained from the spectral maps of u, which are depicted in figure 4, where k z = 2π/λ z is the relevant wavenumber for the spanwise direction.At low Reynolds number, a single peak is present, associated with the wall regeneration cycle (Jiménez & Pinelli 1999).At higher Reynolds number, a secondary energetic site emerges, which is associated with outer-layer, large-scale motions (Hutchins & Marusic 2007).In between the two sites, a spectral ridge is present, corresponding to eddies with typical spanwise length scale λ z ∼ y, here encompassing over one decade in flow case F, which can be interpreted as the footprint of a hierarchy of wall-attached eddies as in Townsend's model (Townsend 1976).Although the spanwise spectra may seem noisy in the outer part, their time convergence has been checked carefully.Hence the observed undulations are rather the result of difficulties for the largest 3 z/h y/h . (a-f ) Instantaneous streamwise velocity field (normalized by the mean free-surface velocity U FS ) in a cross-stream plane, for the flow cases A-F, respectively, listed in table 1.Here, y = 0 corresponds to the no-slip wall, and y/h = 1 corresponds to the free surface.Note that only a limited part of the domain span is shown.
structures to be accommodated in the relatively narrow computational box.These effects are commented on in Appendix A.
The spectral maps are compared with those found in closed channel flow (Pirozzoli et al. 2016), at matching Reynolds number (Re τ = 2000) in figure 5. Whereas the inner spectral peak is pretty much the same in the two flows (namely λ + z ≈ 117, at y + ≈ 13), the outer peak is quite different.In fact, it occurs at a wavelength λ z = πh/2, at a distance y = 0.28h in the open channel, and it occurs at λ z = πh/3, at a distance y = 0.18h in the closed channel.This shows that in the open channel, the eddies that emerge from the no-slip wall can penetrate a greater distance, which also implies that the eddies are larger.
Reynolds number effects on the streamwise velocity spectra are scrutinized more closely in figure 6, where we show spectral densities at the inner energy site location (y + = 15), and at the outer site location (y/h = 0.28).Figure 6(a) clearly brings out universality of the spectra in inner units at the small scales.Additional energy appears at large scales as  More refined information on the behaviour of the mean velocity profile can be gained from inspection of the log-law diagnostic function, namely y + dU + /dy + , which is shown in figure 8, and whose constancy would imply the presence of a genuine logarithmic layer.The figure supports universality of the inner-scaled streamwise velocity for Re τ 10 3 , up to y + ≈ 60, where the diagnostic function attains a minimum.The figure also shows clear flattening of this indicator as the Reynolds number increases.At Re τ = 2000, an extended flat region is first observed, for which the extrapolated Kármán constant is in perfect agreement with the results of Yao et al. (2022), and in very good agreement with the data of Coleman et al. (2017) for Couette-Poiseuille flow with a zero-shear wall, from which k ≈ 0.37 can be inferred.At higher Re τ the Kármán constant seems to decrease a little, and at Re τ = 6000 there is a well developed logarithmic layer up to y/h = 0.5, corresponding to k = 0.375.These values of the Kármán constant seem to differ somewhat from the commonly quoted value k ≈ 0.39 for closed channel and pipe flow (Lee & Moser 2015;Pirozzoli et al. 2021;Hoyas et al. 2022).However, in those cases, significant scatter and ambiguity persist on how the fitting is carried out, in the absence of extended regions with flat diagnostic functions.In fact, systematic deviations from a logarithmic distribution are observed in internal flows, which are typically modelled through additional linear terms (Afzal 1982;Jiménez & Moser 2007;Luchini 2017).In particular, the latter study attributed the occurrence of additional linear terms in the overlap layer expansion of the mean velocity profile in internal flows to the effect of the imposed pressure gradient.The case of the open channel flow, in which a pressure gradient is present, but deviations from the log law are not evident, seems to rule out such explanation.On the other hand, the present findings seem to corroborate the interpretation given by Monkewitz (2021), that the late start of the logarithmic layer in internals flows is due to intrusions from the opposite wall, which in fact are not present in the open channel flow.Differences with other canonical wall-bounded flows are made more explicit in figure 8(b), which makes it quite clear that the log-law diagnostic function is much flatter in the present case than in other flows, and even closed channel flow at substantially higher Reynolds number (Hoyas et al. 2022) exhibits a much wider range of variation.
The structure of the wake region of the open channel is examined in detail in figure 9, where the mean velocity profiles are shown in defect form.Full outer-layer similarity is achieved even at the lowest Reynolds numbers under scrutiny, starting at y/h ≈ 0.2.The defect logarithmic profile fits the data very well with k = 0.375, and with additive constant B = 0.128.Deviations from the outer log law seem to be very minor, and concentrated in a narrow layer adjacent to the free surface.
The previous observations regarding the structure of the mean velocity field can be leveraged to deduce that the dependence of the bulk and mean free-surface velocity on the friction Reynolds number should be logarithmic.This is checked in figure 10, where we show the mean velocity at the free surface and the bulk velocity normalized by the friction velocity, as functions of the friction Reynolds number.Consistently with the case of internal flows (e.g.Monkewitz 2021), the data in fact support logarithmic increase with Re τ , according to with data fitting suggesting k FS ≈ k = 0.375, and B b = 1.57,B FS = 4.02.As a result of the observed identity (or very close vicinity) of the Kármán constant for the centreline and for the bulk velocity, figure 10(b) highlights that their ratio approaches unity at large Re, supporting the inference that open channel flow asymptotes to plug flow in the infinite-Reynolds-number limit; see Pullin, Inoue & Saito (2013).Regarding that study, it is worthwhile noticing that one of the assumptions made in the analysis is that the wall-normal location of the onset of the logarithmic region is either finite or increases no faster than Re τ .Interpreting the near-wall minimum of the diagnostic function in figure 8 as the root of the (near) logarithmic layer, our data well support that assumption.However, as figure 10(b) suggests, this trend is extremely slow.Equations (3.2a,b) can be exploited in order to determine friction relationships for flow in smooth open channels, namely whence the friction coefficient (either ) can be evaluated as an implicit function of the bulk Reynolds number, or of the Reynolds number at the free surface (Re FS = U FS h/ν).Re τ The distributions of the velocity variances along the coordinate directions are shown in figure 11, in inner scaling.As is well established (Marusic & Monty 2019), the streamwise (u) and spanwise (w) velocity fluctuations show slow increase with the Reynolds number, with commonly accepted logarithmic growth as in Townsend's attached-eddy model (Townsend 1976).On the other hand, the wall-normal velocity fluctuations seem to level off to maximum value approximately 1.30 (Smits et al. 2021).It is remarkable that the general growth of the longitudinal and spanwise fluctuations is more evident in the outer layer, and in fact it has long been argued about the possible occurrence of a secondary peak of u 2 , besides the primary buffer-layer peak.Experiments carried out in the Superpipe (Hultmark et al. 2012) and CICLoPE (Willert et al. 2017) facilities indeed support the occurrence of such a peak at Re τ 10 4 .Whereas DNS data do not reach sufficiently high Re τ to show this secondary peak, it appears that in flow case F, the streamwise velocity variance has attained a nearly horizontal inflectional point at y + ≈ 140.Comparison with the data from closed channels shows generally higher values of all fluctuating velocity components, which is consistent with additional energy resulting from the more intense outer-layer eddies.This also yields higher probability that the outer-layer peaks of u 2 , if found, would occur earlier in the open channel than in the closed channel.Additional differences are related to the behaviour of the velocity variances near the free surface, where of course v 2 must go to zero as a consequence of the rigid-lid assumption, and where the horizontal velocity variances show some pile-up and become very similar in magnitude (Nezu 2005).Comparison with the experimental data of Duan et al. (2020) in figures 11(a,b) shows generally good agreement, especially in the outer part of the flow.Differences observed towards the wall, and especially overestimation of the vertical velocity variance in the experiment, are the likely result of well-known difficulties in measuring in the vicinity of walls.
Figure 11(d) displays the trend of the streamwise velocity variance with Re τ .This quantity is the subject of a lot of theoretical interest, and in particular its behaviour in the infinite-Re limit.In fact, the attached-eddy model (Townsend 1976;Marusic & Monty 2019) predicts that it should be increasing logarithmically with Re τ , on account of influence from distant, outer-layer eddies, which would imply a violation of the strict wall scaling.On the other hand, arguments have been made recently that this peak should stay finite, as a result of saturation of growth of the turbulence kinetic energy dissipation rate at the wall (Chen & Sreenivasan 2021).Although difference between slow logarithmic growth and attainment of an asymptotic value is quite tenuous in practice, the theoretical interest is high as in the latter case universality of wall scaling would be restored eventually.The data shown in figure 11(d) generally seem to support logarithmic growth, with slope approximately 0.63, very close to the case of closed channel (Lee & Moser 2015), and as suggested from a collection of DNS and experiments on boundary layers (Marusic, Baars & Hutchins 2017).At the same time, the theoretical prediction of Chen & Sreenivasan (2021) (see the dotted purple line of figure 11d) seems to conform well with DNS data at the highest Reynolds number achieved.Discriminating between these two trends would probably require much higher Reynolds numbers than are currently feasible in DNS of wall-bounded flows.
Whereas our DNS data cannot be used to evaluate directly the theoretical predictions owing to limited achievable Reynolds number, they can be used to better scrutinize the foundations of the theoretical arguments.The main argument made by Chen & Sreenivasan (2021), although not strictly justified mathematically, was that since turbulence kinetic energy production is bounded, the wall dissipation must also stay bounded.Letting P = − uv dU/dy be the turbulence kinetic energy production rate, and 11 = ν |∇u| 2 be the dissipation rate of the streamwise velocity variance, those authors argued that the wall limiting value of 11 should scale as with β a suitable constant.In figure 12(a), we explore deviations of 11w and of the peak turbulence kinetic energy production, say P PK , from their asymptotic value, namely 1/4.According to analytical constraints (Pope 2000), we find that production tends to its asymptotic value quite rapidly, as approximately 1/Re τ .Consistent with (3.4), the wall dissipation rate also tends to 1/4, more or less at the predicted rate, thus validating the first assumption empirically.The next argument advocated by Chen & Sreenivasan (2021) is that wall balance between viscous diffusion and dissipation and Taylor series expansion of the streamwise velocity variance near the wall yields whence, from assumed invariance of the peak location of u 2 (say, y + PK ), saturation of growth of the peak velocity variance would follow.Figure 12(b) reports the position of the streamwise velocity variance peak position as a function of Re τ , as inferred from cubic spline interpolation of the DNS data.Also accounting for uncertainty associated with data interpolation, the figure suggests that this second assumption may not be valid, as the (inner-scaled) peak position shifts away from the wall at increasing Re τ , with non-negligible effect on the peak variance as it appears in squared form in (3.5).As a consequence, logarithmic growth of the peak velocity variance would still be possible even in the presence of saturation of the wall dissipation rate, at least based on DNS data in the currently accessible range of Reynolds numbers.
It is well established that the intensity of wall pressure fluctuations in boundary layers and channels increases logarithmically with the Reynolds number, when scaled by the mean wall shear stress (Willmarth 1975;Farabee & Casarella 1991;Hu, Morfey & Sandham 2006;Tsuji et al. 2007;Jiménez & Hoyas 2008;Sillero, Jiménez & Moser 2013).Tsuji et al. (2007) found a similar trend for the peak pressure variance, which occurs at approximately y + = 30.Jiménez & Hoyas (2008) suggested that this logarithmic increment should be attributed to a growing hierarchy of self-similar wall-attached eddies (Townsend 1976;Hwang 2015;Marusic & Monty 2019).Pressure, like the wall-parallel velocity components, can indeed be regarded as a wall-attached quantity (Bradshaw 1967), which explains why profiles of the pressure variance also follow a logarithmic trend with respect to the wall distance.The same behaviour can, however, also be predicted based on inner/outer layer overlap arguments (Panton, Lee & Moser 2017).Mehrez et al. (2019) explored the DNS database of closed channel flow at friction Reynolds number up to 4000, and found that the logarithmic trend shows up only for Re τ 500.
The distributions of the pressure variances obtained from the present DNS are displayed in figure 13(a).These are found to remain nearly constant within the viscous sublayer (say, y + ≤ 5), then increase to attain a peak value at y + ≈ 30, and decrease thereafter.A region with distinctly logarithmic decrease is found from y + ≈ 100 basically all the way to the free surface, with slight deviations.As the Reynolds number increases, the pressure variances at the wall and their peak values increase systematically, whereas the value at the free surface is not much affected.Compared with the variances in closed channels (dashed lines), we note that pressure fluctuations in open channels tend to be a bit less in the buffer layer, whereas they are virtually identical towards the wall and in the outer layer.When plotted against outer wall coordinates, the logarithmic portions of the pressure variance distribution collapse for Re τ 500, consistent with the case of closed channel flow (Mehrez et al. 2019).Expressing the outer distribution as thus supporting theoretical inferences quite well.The logarithmic growth rate seems to be identical to that for closed channel flow.

Concluding comments
The study of flow in an open channel has the merit of isolating the effect of wall-imposed shear from other complicating effects, for instance due to the presence of opposing walls (as in closed channel flow) or geometrical confinement in general (as in pipe flow).Hence, despite limited practical impact, as by far most instances of open channel flow occur over rough beds, we believe that the theoretical interest of studying open channel flow is quite large.We find that flow in an open channel is a very good candidate to observe the log law of the wall in its 'pure' form, and in fact a clear logarithmic layer is present even at 'modest' Reynolds number (Yao et al. 2022), and in the present DNS there is a clear evidence that a genuine logarithmic layer exists, at least up to half of the channel depth.The estimated Kármán constant is found to increase a little with the Reynolds number, starting from k = 0.363 at Re τ = 2000 (in perfect agreement with Yao et al. 2022), to arrive at k = 0.375 at Re τ = 6000.A certain variability of the Kármán constant was in fact noticed in previous systematic experimental studies (Nagib & Chauhan 2008), which is possibly related to finite-Re effects.In any case, the scenario is much clearer than in internal flows such as closed channels and pipes.In the case of a closed channel, convincing evidence for a (small) region with flat behaviour of the log-law diagnostic function has come only recently from DNS at Re τ = 10 000 (Hoyas et al. 2022).Otherwise, data at lower Re seem to feature obvious deviations from the expected logarithmic behaviour that are traditionally modelled as additional linear terms (Jiménez & Moser 2007;Luchini 2017), and the Kármán constant is extrapolated based on fitting those compound distributions to the DNS data (e.g.Pirozzoli et al. 2021).The resulting value of the Kármán constant is then k = 0.38-0.39,hence a bit higher than found in the open channel.It is also noteworthy that Coleman et al. (2017) carried out numerical simulation of the Couette-Poiseuille flow (namely, channel flow with combined relative motion of the walls and imposed pressure gradient) under conditions such that one of the walls has zero shear, and found very similar behaviour, namely an extended logarithmic layer with Kármán constant k = 0.37.In our opinion, universality of wall turbulence in the near-wall region requires universality of the Kármán constant, hence we believe that variations observed in all studies so far are related to limited Reynolds number, or to numerical and experimental uncertainties, or effects of geometrical confinement, or all of these.Study of flow in open channels has the merit of removing at least the latter contaminating effects.Regarding deviations from the log law in internal flows, Luchini (2017) argued that additive linear terms should be present on account of the driving pressure gradient.However, a driving force is present in the open channel flow, but clear deviations from logarithmic behaviour are not observed.Hence it seems more reasonable to accept that deviations or late occurrence of a logarithmic behaviour in internal flows is rather due to geometric confinement effects, and specifically to influence of the opposing wall in closed channels, as argued by Monkewitz (2021).Study of flow in open channels also bears some interesting implications regarding interactions between the inner and outer wall layers, and their universality.Our study corroborates previous findings (Duan et al. 2020;Yao et al. 2022) that large outer-layer eddies form in open channels, which are larger and more energetic than in closed channels.Here, we find evidence that these eddies also become universal at sufficiently high Reynolds number, when lengths are expressed in outer units.This would support universality of the outer-layer dynamics, in addition to the well-recognized universality of the inner-layer dynamics, as proposed by Dennis & Sogaro (2014).The signature of wall-attached eddies is also evident in the spectral footprints, as well as in the statistics of the spanwise velocities and pressure fluctuations, which exhibit extended logarithmic layers.The imprinting of the energetic outer-layer eddies manifests itself with higher values of the velocity variances.This in turn translates to earlier occurrence of the (possible) outer-layer peak, which we believe can be observed in open channel flow at Re τ ≈ 10 4 .Regarding the growth of the inner-layer peak of the streamwise velocity variance, available evidence seems to support logarithmic growth with Re τ , although firm asymptotic trends are difficult to discern.Specifically, we find empirical evidence that the turbulence kinetic energy dissipation rate at the wall reaches an asymptote, as argued by Chen & Sreenivasan (2021), but we also find that the peak location at which the streamwise velocity variance is maximum slowly varies with Re τ , which would allow for logarithmic growth of the peak value even if the growth of the wall dissipation was to saturate.
In summary, we find that, besides its intrinsic interest, open channel flow can be exploited usefully as a testbed to explore the behaviour of wall turbulence at extreme Reynolds number, and especially to disentangle various complicating effects.Follow-up studies will also include analysis of the process of exchange of mass and momentum at the free surface at computationally high Reynolds number.(A1), however, with j b = 15, c η = 1.03.The two distributions are smoothly connected at y/h = 0.73, as shown in figure 14.The wall-normal resolution is checked a posteriori in figure 15, where we show the grid spacing expressed in local Kolmogorov units, which is found to be no larger than 2.6, throughout the channel.Adequate resolution of the free-surface layer is also evident in figure 15, which shows that it is resolved with ten collocation points at least.
Sensitivity of the computed results to the computational box size and grid resolution has also been verified by carrying out additional simulations at Re b = 20 000, as listed in table 2. As one could argue from the instantaneous flow visualizations of figure 2, spanwise confinement effects are expected, especially towards the free surface, where structures are wider.This is why flow cases C-M, C-W, C-WW were carried out, corresponding to box width up to L x = 5πh.Some confinement effect is in fact noticeable in the spectral maps (see figure 16), which show some variability of the outer spectral peak position and wavelength.This is made quantitative in table 2, where the outer peak wavelength is reported.Apparently, having a wider domain allows the outer eddies to become a little larger.The effect of rearrangement of the outer part of the spectral maps is, however, quite small when it comes to the one-point statistics, which we show in figure 17.In fact, scatter in the mean and first-order statistics, shown in figure 17(a), is much less than 1 %, hence within the statistical sampling error.This error is magnified when the diagnostic function is reported, as in figure 17(b).In this case, the scatter is larger, but not to the extent of modifying the results qualitatively.

Figure 1 .
Figure 1.Definition of the computational set-up for DNS of open channel flow.Here, x, y, z are the streamwise, wall-normal and spanwise directions, respectively, L x and L z are the box sizes along the streamwise and spanwise directions, h is the channel depth, and u b is the bulk velocity.

Figure 5 .Figure 6 .Figure 7 .
Figure 4. (a-f ) Variation of pre-multiplied spanwise spectral densities of fluctuating streamwise velocity with wall distance, for cases A-F, respectively.Wall distances and wavelengths are reported both in inner units (bottom, left axes) and in outer units (top, right axes).The dashed line denotes the channel free surface.The crosses denote the positions of the inner and outer energetic sites.Contour levels from 0.2 to 2.0 are shown, in intervals of 0.2.

Figure 10
Figure 10.(a) Mean free surface velocity (U FS ) and bulk velocity, and (b) their ratio, as functions of the friction Reynolds number.DNS data are shown as circle symbols, and the corresponding logarithmic fits are shown as dashed lines.Black is used for U FS , and purple for u b .

Figure 11 .
Figure 11.Distribution of velocity variances (a-c) and trend of peak streamwise velocity variance with Reynolds number (d).The blue circles in (a,b) denotes experimental data of Duan et al. (2020), at Re τ = 1903.The circles in (d) denote the present open channel data, and the squares the closed channel data of Lee & Moser (2015), at Re τ = 550, 1000, 2000, 5200.The dashed and dot-dashed data denote the corresponding logarithmic fits, and the dotted line the trend predicted by Chen & Sreenivasan (2021).

Figure 12
Figure 12.(a) Distributions of peak turbulence production (P PK , squares), and wall dissipation of streamwise velocity variance ( 11w , circles).(b) Peak position of streamwise velocity variance.In (a), the dot-dashed and dotted lines denote fits of P PK and 11w in their tendency to the respective assumed asymptotic values.

Figure 13
Figure 13.(a) Wall-normal distribution of pressure variance, and (b) wall pressure variance as a function of Reynolds number.In (a), the dashed lines correspond to closed channel flow DNS data (Lee & Moser 2015) at Re τ = 180, 550, 1000, 2000, 5200, and the grey line denotes the logarithmic fit given in (3.6).In (b), circles are used for the open channel data, and squares for the closed channel data; the dashed lines denote the logarithmic data fit given in (3.7).Refer to table 1 for colour codes.
.6) curve fitting within the range y = 0.02h-0.5hfor flow case F yields the values of parameters A p = −2.40 and B p = 1.23, with error 0.4 %.We further report the wall pressure variances as a function of Re τ in figure 13(b).With exclusion of the lowest Reynolds number case at Re τ ≈ 180, the wall pressure variances can be well characterized in terms of the log law p 2 w + =2.22 log Re τ − 8.99, (3.7)

Figure 17 .
Figure 17.Box and grid sensitivity study for one-point statistics: (a) inner-scaled mean velocity profiles and streamwise velocity variances, and (b) log-law diagnostics function.The dashed line in (a) refers to the logarithmic fit U + = log y + /0.375 + 4.21.The dashed horizontal line in (b) denotes the inverse Kármán constant 1/0.375.The line style is as in table 2.

Table 1 .
Flow parameters for DNS of open channel flow.Here, N x , N y and N z are the numbers of grid points in the streamwise, wall-normal and spanwise directions, respectively, Re b = 2hu b /ν is the bulk Reynolds number, Re τ = u τ h/ν is the friction Reynolds number, and ETT is the time interval used to collect the flow statistics, in units of the eddy turnover time h/u τ .

Table 2 .
Computational parameters for box and grid sensitivity study.