Flow Statistics in the Transitional Regime of Plane Channel Flow

The transitional regime of plane channel flow is investigated above the transitional point below which turbulence is not sustained, using direct numerical simulation in large domains. Statistics of laminar-turbulent spatio-temporal intermittency are reported. The geometry of the pattern is first characterized, including statistics for the angles of the laminar-turbulent stripes observed in this regime, with a comparison to experiments. High-order statistics of the local and instantaneous bulk velocity, wall shear stress and turbulent kinetic energy are then provided. The distributions of the two former quantities have non-trivial shapes, characterized by a large kurtosis and/or skewness. Interestingly, we observe a strong linear correlation between their kurtosis and their skewness squared, which is usually reported at much higher Reynolds number in the fully turbulent regime.


Introduction
Laminar and turbulent flows are two different regimes encountered sometimes at the same parameters for a given geometry. In many flows they are in competition from the point of view of the state space. Shear flows next to solid walls however show this surprisingly robust property that both laminar and turbulent regions coexist spatially on very long time scales, when the laminar state is locally stable. This phenomenon, called 'laminar-turbulent intermittency' is well known in circular pipe flow since the days of O. Reynolds [1] and has lead recently to a burst of interest, a review of which is provided in Reference [2]. Such laminar-turbulent flows have been identified and partly characterized in Taylor-Couette flow [3,4] and in plane Couette flow [4][5][6]. They also have been identified in other set-ups involving curvature [7][8][9] or stabilizing effects [10]. The transitional regimes of plane Poiseuille flow, the flow between two fixed parallel plates driven by a fixed pressure gradient, have not received as much attention although this flow is the archetype of wall-bounded turbulent flows. Although this flow is frequently cited as an example of flow developing a linear instability (under the form of Tollmien-Schlichtling waves) [11], coherent structures typical of laminar-turbulent coexistence have been frequently reported in channel flow well below the linear instability threshold and a series of experimental and cutting-edge numerical studies in the 1980s and 1990s have focused on the development of spots [12][13][14][15][16]. Sustained intermittent regimes have not been identified as such before the mid-2000s, when Tsukahara [17] reported large-scale coherent structures from numerics in larger numerical domains. Like their counterpart in Couette flows, these structures display obliqueness with respect to the mean flow direction and a complicated long-time dynamics. The dynamics at onset in particular have remained mysterious [18] and, although this is currently debated, could follow a scenario different from the directed percolation one proposed for Couette flow. [9,19,20]. In recent years, the so-called transitional regime of plane channel flow has attracted renewed attention after new experimental studies. Although the works in Refs [21][22][23] focused on the minimal transition amplitude for spot development, other studies [24][25][26][27][28][29] focused on the sustained intermittent regimes and their statistical quantification.
Experimentally the finite length of the channel sets a limitation to most statistical approaches. Numerical simulation in large domains combined with periodic boundary conditions is a well-established way to overcome such limitations. Surprisingly, despite a large number of numerical studies of transitional channel flow, investigation of spatio-temporal intermittency in large enough domains has not been possible before the availability of massive computational resources. Owing to recent numerical studies [30][31][32], there is currently a good consensus about a few facts concerning the transitional regime: laminar-turbulent bands with competing orientations emerge progressively as the Reynolds number is reduced below Re τ « 100, and their mean wavelength increases as the Reynolds number is decreased. At even lower flow rate the bands turn into isolated spots with ballistic dynamics rather than forming a seemingly robust stripe pattern [33][34][35]. The global centerline Reynolds number for the disappearance of the stripes is close to 660 [18,27]. However, many questions remain open. The most sensible theoretical issues revolve around the (still open) question of the universality class of the transition process (see Reference [18]), the role of the large-scale flows [23,25,36,37] in the sustainment of the stripes, or the mutual way different stripes interact together.
There is also a lack of quantitative data about the patterning regime itself. The present special issue is an opportunity to document the geometric characteristics of the stripe patterns in unconstrained settings. Moreover, there is an ongoing philosophical question about whether traces of spatio-temporal intermittency can be found in the fully turbulent regimes commonly reported at higher Reynolds numbers. In the present paper, using numerical simulation in large domains, we focus on three specific points hitherto undocumented: the angular distribution of turbulent stripes, the statistics of the laminar gaps between them, and high-order statistics of the local and instantaneous bulk velocity, wall shear stress and turbulent kinetic energy. The outline of the paper is as follows: Section 2 introduces the numerical methodology with the relevant definitions. The geometrical statistics of the stripe angles are presented in Section 3.1. The statistics of a few global quantities are presented in Sections 3.2-3.4. A discussion of the results is made in Section 4 with the conclusions and outlooks in Section 5.

Materials and Methods
The present section is devoted to the methodology used for the numerical simulation of pressure-driven plane channel flow. The flow is governed by the incompressible Navier Stokes equations. Channel flow is described here using the Cartesian coordinates x,y,z, respectively the streamwise, wall-normal and spanwise coordinates. The velocity field upx, y, z, tq is decomposed into the steady laminar base flow solution Upyq " pU x , 0, 0q and a perturbation field u 1 px, y, z, tq. Similarly, the pressure field is decomposed as ppx, y, z, tq " xG`p 1 px, y, z, tq. The equation governing the steady base flow for an incompressible fluid with constant density ρ and kinematic viscosity ν is given by with G a constant. Together with the no-slip condition at the walls Equation (1) yields the analytic Poiseuille solution U x 91´py{hq 2 . The equation governing the perturbation field involves the base flow and reads The channel geometry is formally infinitely extended, yet in the numerical representation it is given by its extent L xˆ2 hˆL z as in Figure 1, with stationary walls at y "˘h and periodic boundary conditions in x and z. The flow is driven by the imposed pressure gradient G assumed negative. The spanwise pressure gradient is explicitly constrained to be null. The centerline velocity u cl of the laminar base profile with the same pressure gradient is chosen as the velocity scale (U) and the half gap h of the channel is chosen as the lengthscale used for non-dimensionalization. Time is hence expressed in units of h{U. In these units the laminar velocity profile is given by Ux py˚q " 1´y 2 . From Chapter 3 onwards only dimensionless quantities will be used and the˚notation will be dropped from there on. Primed quantities denote perturbations to the base flow while non-primed quantities involve the full velocity field, including the laminar base flow.
In the following we shall consider, both locally and temporally fluctuating quantities, as well as their time and space averages. We denote by ‚ the space (x, z) average and s ‚-the time average. Space-time averages are indicated by Ď ¨ . More explicitly the space-average operator is defined as the discrete average over the grid points, and the time average is the discrete average sum over the total number of snapshots in the steady regime.
Different velocity scales characterize the flow. One such scale is the centerline velocity u cl of the corresponding laminar flow with the same value of G. Another one is the total streamwise flow through the channel, is the so-called local bulk flow. Finally, the friction velocity is defined as U τ " p Ď τ {ρq 1 2 , where τ " pτ t`τb q {2 ą 0, with τ t and τ b the net shear stress on the top and the bottom wall, respectively given by: where µ " ρν is the dynamic viscosity of the fluid. The three Reynolds numbers arising from these velocity scales are Re cl " u cl h{ν, Re b " U b h{ν and Re τ " U τ h{ν. For the laminar base flow, they are inter-related as Re 2 τ " 3Re b " 2Re cl . Imposing a pressure gradient G < 0 translates into a fixed average shear stress Ď τ on the walls which sets an imposed value of Re τ " Re G τ to stress that this is the control parameter.
Direct numerical simulation (DNS) of Equation (2) is carried out using the open source, parallel solver called Channelflow [38,39] written in C++. It is based on a Fourier-Chebychev discretization in space and a 3rd order semi-implicit backward difference scheme for timestepping. It makes use of the 2{3 dealiasing rule for the nonlinear terms. An influence matrix method is used to ensure the no-slip boundary condition at the walls. The numerical resolution is specified in terms of the spatial grid points pN x , N y , N z q which translates into a maximum of pN x {2`1, N z {2`1q Fourier wavenumbers and N y Chebychev modes. Please note that the definitions of N x and N z take into account the aliasing modes. The domain sizes used in this study, expressed in units of h, are L x " 2L z " 250 for 55 ă Re G τ ď 100 and L x " 2L z " 500 for 39 ď Re G τ ď 55. The local numerical resolution used is N x {L x " N z {p2L z q " 4.096 and N y " 65, comparable to that used in Reference [34]. The simulation follows an "adiabatic descent": a first simulation is carried out at sufficiently high value of Re G τ , known to display space-filling turbulence. After the stationary turbulent regime is reached, Re G τ is lowered and the simulation advanced further in time. This step-by-step reduction has been performed down to Re G τ " 39. The initial condition for the first simulation is a random distribution of localized seeds of the kind described in Reference [40]. The time required T to reach a stationary regime gradually increases as Re G τ is decreased. As an order of magnitude, for Re G τ " 100, T « 1500, while for Re G τ " 50, T « 3000. Statistics are computed, after excluding such transients, from time series of lengths up to 2ˆ10 4 time units.

Results
The entire adiabatic descent is shown using a space-time diagram of the crossflow energy shown in Figure 2a evaluated at an arbitrary value of z (here z " L z {2). The space variable is expressed in a frame moving in the streamwise direction with the mean bulk velocity U b pGq for that particular value of Re G τ . Since Re G τ is lowered over the course of time, this allows one to capture the different flow regimes preceding full relaminarization. The intensity of turbulence, measured here by the value of E c f , is seen to gradually increase as Re G τ is lowered. At high Re G τ , the so-called featureless turbulence occupies the full domain, as shown in Figure 2b at Re G τ " 100 using isocontours of τ 1 px, zq " τpx, zq´τ lam . As Re G τ is lowered, turbulence self-organizes into the recognizable pattern regime [17] shown in Figure 2c for Re G τ " 80. As Re G τ is further reduced the turbulent zones become sparser (see Figure 2d for Re G τ " 60). The spatially localized turbulent regions emerge as narrow stripes throughout the process of decreasing Re G τ while the gaps between them constantly increase in size. The emerging patterns never feature an array of strictly parallel stripes like in former computational approaches [19,31,41], instead they feature competing orientations as in pCf [4], see Figure 2b-d. In this regime the pattern travels with a streamwise convection velocity slightly slower than U b pGq. Within the quasi-laminar gaps, E c f reaches very low values, at least an order of magnitude less than in the core of the turbulent stripes. The lower Re G τ , the lower these values. Below Re G τ " 50 the stripe pattern eventually breaks up to form independent turbulent bands of finite length, all parallel to each other [34], as shown in Figure 2e for Re G τ " 40. The new resulting pattern as a whole shows negligible spanwise advection, while it propagates in x with a velocity close to Ě u b [42]. The independent turbulent bands show enhanced motility in both directions x and z. This motion relative to the frame of reference causes the tilt of the stripes seen in Figure 2a for Re G τ ą 50 as well as the apparent increase of thickness. In pipe flow it was noted recently [43] that the emergence of spatial localization does not imply the proximity to the transitional point (below which turbulence is not sustained) as long as the statistics about the size of the laminar gaps fail at displaying power-laws tails. The laminar gaps are estimated as the streamwise distance l x between local maxima of τ (values lower than τ `σpτq, with σ the standard deviation, have been discarded). The cumulative distribution (CDF) of the laminar gap size is shown in Figure 3 in lin-log coordinates. For all values of Re τ shown, it shows an exponential tails and no algebraic part. Exponential distributions are a hallmark of spatio-temporal intermittency, unlike critical phenomena which are characterized by algebraic/power law related to the scale invariance property. The entire regime of channel flow for 39 ď Re G τ ď 100 can be described as being spatiotemporally intermittent, and is hence far above any critical point. Please note that the critical point of pPf is estimated to approximately Re cl " 660 [18] i.e., Re G τ « 36 and falls outside the range of parameters investigated here.

Angular Statistics of Turbulent Bands
The self-organization of turbulence into long band-like structures, oriented with an angle with respect to the streamwise direction, is depicted in Figure 2. The (signed) angle is computed using two different methodologies. As in Duguet et al. [36] in the case of pCf, the local y-integrated velocity field is found to be parallel to the bands. The same holds for pPf, as is visible in Figure 4a,c for Re G τ " 60 and 40, respectively. Please note that unlike Couette flow, pPf features advection with a non-zero mean bulk velocity. Hence the local velocity field is here computed by removing this mean advection velocity. A first estimation of the local and instantaneous band angle is therefore computed following Equation (6): The second estimation is obtained from Fourier analysis and computed from Equation (7), following Reference [44] : where λ " 2π{k, with k being the leading non-zero wavenumber identified from the power spectra (excluding the k x " k z " 0 mode). The Fourier spectrum is computed for the quantity τpx, z, tq, but similar results have been observed for other observables such as E c f px, z, tq and E v " p1{2q ş u 2 y dy. The angles can be read directly from the Fourier spectra in polar coordinates, see Figure 4b,d for the same values of Re G τ " 60 and 40, respectively. The mean angles Ě θ L and s θ F are then computed by respectively space-time-averaging and time averaging the data obtained from Equation (6) and (7).
The variation of the mean (signed) angles with Re G τ , computed using the two methods, is shown in Figure 5a, where the indices 1, 2 stand for the two band orientations. Both methods provide identical results. The variation of the (unsigned) angle of the band denoted by θ, computed as θ " Ě |θ F | is shown in Figure 5b. It is found that the mean angle θ of the bands remains approximately constant with θ " 25˝˘2.5˝in the range of values 60 ď Re G τ ď 90 and increases for lower value of Re G τ ă 60. In the patterning regime, i.e., for Re G τ ě 50, the angle of the bands is found to be distributed symmetrically with respect to zero, as a consequence of the natural symmetry z Ð´z of the flow. For lower Re G τ these quasi-regular patterns break down into individual localized structures analogous to individual puffs in cylindrical pipe flow. As the pattern dissolves, one single band orientation ends up dominating the dynamics as shown by Shimizu and Manneville [34] for a similar domain size. The angle θ further increases as the regular pattern deteriorates, with θ max « 40 at Re G τ " 39. Previous studies [18,27] have documented that the angle of the bands approach 45˝close to the onset of transition. The present investigation agrees well with these studies (Figure 5b) while covering a wider range in Reynolds number, highlighting the difference between the puff regime for which θ « 40´45˝, and the patterning regime for which θ is almost half this value (see also Figure 2).  Figure 5. (a) Variation of the mean (signed) angle of the turbulent bands with Re G τ , computed from the Fourier spectra ( Ď θ F 1 , Ď θ F 2 ) and the mean (signed) angle of the local velocity ( Ę θ L 1 , Ę θ L 2 ) (b) Variation of the mean unsigned band angle θ along with the data from Reference [27,30]. Figure 4c shows that across a band, the local large-scale velocity changes orientation [41]. This property is used to sort out the local maxima of τ (higher than τ `σpτq) as belonging to one band with a particular inclination. This allows one to define the respective streamwise and spanwise inter-stripe distances l x and l z between bands of the same orientation. Figure 6a,b displays Ě l x and Ě l z for orientations 1 and 2, respectively, as a function of Re G τ . Both increase when decreasing Re G τ . They vary in parallel in the patterning regime, hence the quasi-constant angle θ of the bands. When only one band orientation survives, one observes that the increase in θ amounts to the saturation of Ě l x 1 , while Ě l z 1 keeps increasing.

Global Variables: Moody Diagram
The mean velocity profile Ę ux is defined as the average of u x over x,z and t, expressed in units of u τ . It is shown in Figure 7 as a function of y`" yu τ {ν and compared with the classical DNS data by Kim et al. [45] obtained at higher Re G τ " 180. The whole figure is similar to figures 3 and 10 in Reference [17,46], respectively. As expected for the present low values of Re G τ , the velocity field matches the linearized profile ux " y`next to the wall but does not develop a logarithmic dependence with respect to y`. At a global level of description, the laminar and turbulent flow are traditionally represented in the classical Moody diagram in which the Fanning friction factor C f defined as the ratio between the pressure drop along the channel length and the kinetic energy per unit volume based on the mean bulk velocity U b " Ě u b , is traditionally plotted versus Re b as shown with plain symbols in Figure 8. Another way to express C f is to use inner units, in which case C f " 2{pub q 2 , with ub " u b {u τ . C f is then linked only to the integral of the mean profile displayed in Figure 7.
For the laminar flow, the dependence of C f vs. Re b is analytically given by C f lam " 6{Re b (blue continuous line). In the featureless turbulent regime, it is known empirically as the Blasius' friction law scaling Ě Re b´1 {4 (red continuous line). For intermediate values of Re b , C f clearly deviates from the turbulent branch, and remains far from the laminar value [47]. Here we notice, in agreement with [30] and [34] that C f « 0.01 remains essentially constant in this transitional regime. What is remarkable is that this regime of constant C f coincides with the patterning regime observed for 50 ď Re G τ ď 90, corresponding to 690 ď Re b ď 1225, as if the respective amount of turbulent and laminar domains was precisely ensuring C f " cst. As the pattern fractures, C f increases and approaches the laminar curve. We note that the observation of this property requires large computational domains to be observed, which explains why it had not been noticed until recently, even in experiments. Given the complex spatio-temporal dynamics in the transitional regime, the bulk velocity u b is expected to strongly fluctuate both in space and time. We also report in Figure 8, how these fluctuations would translate on Re b and C f , if the latter were computed using the locally fluctuating field u b instead of its mean value U b . These fluctuations are significant (up to 10-15%) and suggest to further explore them, which is the topic of the next section and the main focus of the present work.

Joint Probability Distribution of Re τ and Re b
Reynolds numbers such as Re τ and Re b are traditionally seen as global parameters characterizing the flow. They are defined based on velocity scales obtained from space-time average. It is straightforward to extend these definitions to the local fields Re b px, z, tq " u b px, z, tqh{ν and Re τ " u τ px, z, tqh{ν, with u τ px, z, tq " pτpx, z, tq{ρq 1{2 . Please note that with this definition, Ğ Re τ is not strictly equal to the imposed Re G τ , because of the nonlinear relation between Re τ and τ. Investigation of the entire transitional regime is provided through a two-dimensional state portrait (Re b´R e τ ) constructed from this local definition of the Reynolds number. The joint probability density distribution is constructed in this state space with the space-time data for different Re G τ . The state space for Re G τ " 100, 80, 60, 40 is shown in Figure 9. The continuous blue and red lines again correspond to the scalings known analytically for the laminar flow, and empirically for featureless turbulent flows for high enough Reynolds numbers. As expected the most probable values of Re b and Re τ , follow the same trend as their global counterpart: they match the continuous curve in the featureless turbulent regime, and progressively depart from it to move towards the laminar branch at the lowest Re G τ explored here. More interesting are the distributions. First, we observe that the relative fluctuations are significantly larger for Re τ than for Re b , the difference being larger for the larger Re G τ . Secondly the distributions are not simple Gaussians. Even in the featureless turbulent regime, the marginal distribution of Re τ is already relatively skewed (Figure 9(a 3 )). . pa 1 q pb 1 q pc 1 q pd 1 q Joint probability distribution of the quantities Re b and Re τ for Re G τ " 100, 80, 60, 40 together with their marginal distribution shown in lin-log scale for Re b in pa 2 q pb 2 q pc 2 q pd 2 q and for Re τ in pa 3 q pb 3 q pc 3 q pd 3 q with the mean value indicated by a vertical/horizontal black line.
As Re G τ is reduced, the overall width of the distribution decreases, but the shape of the marginal distributions of Re τ differs more and more from a Gaussian. More specifically, although the distribution remains unimodal, we note that the marginal distribution of Re τ is more and more skewed. We also note that the right wing of the distribution is not convex anymore. To further quantify these observations, a systematic analysis of the moments of this distribution is conducted in the next section.

Higher-Order Statistics
The higher-order statistics of Re τ , Re b and E c f are presented in this section. For any field A " Apx, z, tq, we compute the spatio-temporal average m " Ě A , the variance σ 2 " Ğ pA´mq 2 and the k th standardized higher-order moment Ğ pA´mq k {σ k (for k ě 3). Their mean values of Re b and Re τ (Figure 10a) simply follow the trends described above for the most probable value of the distribution, connecting the turbulent and the laminar branch, when Re G τ decreases. Away from the turbulent and laminar branches Re τ is linearly related to Re b , in agreement with the observation of a constant C f . The standard deviation σ (Figure 10b) for Re τ and Re b decrease together with Re G τ . This decreasing trend agrees well with the experimental wall shear stress data reported in Reference [29]. The standard deviation for E c f is found to increase with decreasing Re G τ , matching the trend reported in Reference [34].
The variation of the 3rd and 4th moments m 3 and m 4 , i.e., the Skewness (S) and Kurtosis (K), versus Re G τ for the observable Re τ and E c f is shown in Figure 10c. These moments exhibit a strongly increasing trend with reducing Re G τ for both quantities. This similarity in behavior leads to K9S 2 as shown in Figure 10e. This correlation between the third and fourth statistical moments was first noted in Reference [48] for the fluctuating velocity in turbulent boundary layers at high Reynolds number. In the transitional regime, the same relationship has been found to hold in the experiments of Agrawal et al. [29] from wall shear stress data. We therefore confirm this yet-to-be-understood extension of a high Reynolds number scaling down to the spatio-temporal intermittent regime. Furthermore, we observe that the same scaling also holds for the turbulent kinetic energy E c f (Figure 10e). In contrast it does not apply to Re b (inset of Figure 10e). The reason is that while the Kurtosis follows the same trend as for the two other observables, (Figure 10d), the skewness shows a markedly different behavior: it is non-monotonous, changes sign twice and exhibit a maximum in the core of the spatio-temporal intermittent regime.

Discussion
The present simulations of the transitional regime of pPf confirm and extend previously documented knowledge, such as the constancy of C f in the patterning regime and the variation of the band orientations close to the transition point.
The statistical analysis of the distribution of laminar gaps reveals that the distributions are exponentially tailed over the entire parameter range 39 ď Re G τ ď 100, demonstrating that even the value Re G τ " 39 remains away from any sort of critical regime, which would be marked by algebraic distributions. This is consistent with the existing estimation of the location of the transitional critical point Re cl « 660 [18,27], which translates to Re G τ « 36. The entire patterning regime should thus be seen as bona fide spatio-temporal intermittency, with the critical behavior and transition point being relegated to values of Re G τ ă 39. Exploring the statistics of the flow closer to the critical point would require even larger domains and longer observation times. Such an investigation is outside the scope of the current study.
The orientation of the bands in the patterning regime for 60 ď Re G τ ď 90 (1800 ď Re cl ď 4050) is essentially constant, with an angle θ " 25˝˘2.5˝. This validates the choice of θ " 24˝as a suitable value in the numerical approach of Tuckerman et al. [5,31,32], where slender computational domains are tilted at a chosen value of the angle. However, this angle of 24˝no longer fits the mean orientation of the independent turbulent bands in the lower range Re G τ ď 60 (Re cl ď 1800), where the orientation of the bands increases by a factor close to two, with θ « 40˝for Re G τ " 39. We confirm the observation of a constant C f in the patterning regime, which also implies Ğ Re τ " Ğ Re b , as reflected in Figure 10a. This constant value of C f in the transitional regimes further enforces the long lasting analogy with first order phase transitions [49], for which the thermodynamic parameter conjugated to the order parameter remains constant while the system evolves from one homogeneous phase to the other, when a suitable control parameter is varied. At the mean-field level, a trademark of phase coexistence, is then the presence of a bimodal distribution of the order parameter in the coexistence regime. Capturing this bi-modality is however known as being a challenge, even in simulations of standard equilibrium systems: first, not all protocols allow for observing the phase coexistence; second, the order parameter must be coarse-grained on appropriate length-scales as compared to the correlation lengths such that non-mean-field effect do not dominate [50]. More than often, the bi-modality of the order parameter distribution is replaced by a mere concavity and a large kurtosis. If the two phases have very different fluctuations, as is the case here, one also expects a strong skewness of the distribution. Our observations extend the analogy, already reported at the level of the mean observable, to their fluctuations. However, a lot remain to be done to further exploit this analogy, in particular by making more precise what the relevant order and control parameters are. Let us stress that whether the analogy with a first order transition is valid or not, it does not preclude the dynamics at the spinodals from obeying a critical scenario, such as directed percolation close to the laminar phase spinodal [51] and a modulated instability of the turbulent flow close to the turbulent one [4].
Finally, the statistical moments showcased here demonstrate a correlation between the skewness and the kurtosis of both Re τ and E c f . Such a correlation, observed in both the transitional regime and higher Reynolds number turbulence but originally developed for the latter only [48], suggests a universal turbulent character, beyond the mere distinction transitional/featureless.

Conclusions
The transitional regime of pPf has been investigated numerically in large periodic domains. The transitional regime is composed of two sub-regimes each demarcated by a distinct behavior. The patterning regime is characterized, for 50 ď Re G τ ď 90, by a constant value of C f « 0.01 and by a propagation downstream at approximately the mean bulk velocity ă u b ą. For lower Re G τ all the way down to the critical point close to Re G τ ď 36, independent turbulent bands define a regime analogous to the puff regime of cylindrical pipe flow. The patterns are shown to exhibit a near constant angle of inclination θ " 25˝˘2.5˝for 60 ď Re G τ ď 90, which increases with reducing Re G τ . Both sub-regimes can be classified as spatiotemporally intermittent, as demonstrated by the exponential tails of the distribution of laminar gaps. The statistics of the local fields τ and u b reinforce the feeling that a fruitful analogy with first order phase transitions could be developed, but the later remains to be made more precise and exploited.