**6. Materials and Methods**

For solving the incompressible Navier–Stokes equations in channel geometry, we used our in-house code as described in [24,25], which adopts a high-order finite-difference method with a centered nine-point stencil in the wall-normal direction and Fourier-spectral method in the periodic streamwise and spanwise directions. Readers are referred to OPENPIPEFLOW [33] for details about the finite-difference scheme and the parallelization of the code. The Navier–Stokes equations were integrated using the method of Hugues and Randriamampianina [34], which adopts a second-order-accurate backward-differentiation scheme, combined with the Adamas–Bashforth scheme for the nonlinear term, for the temporal discretization and a projection method to impose the incompressibility condition. The time-step size was fixed at Δ*t* = 0.01 for the simulations presented in this paper, which was shown to be sufficiently small for the Reynolds number regime considered [4,6].

We adopted the method proposed by Song and Xiao [25] to generate turbulent bands in large domains. The method firstly derives a body force that is needed to maintain an inflectional velocity profile that bears a sufficiently strong instability. Given a target velocity profile *U*(*y*), the body force is derived as

$$f = -\frac{1}{\mathcal{R}\varepsilon} \nabla^2 \mathcal{U}(y). \tag{2}$$

Then, the body force is multiplied by a localization factor such that the force is localized in the *x*-*z* plane. The size of the localization region should be comparable with the size of the head of a turbulent band and the forcing region is moved at the speed of *cz* = 0.1 (absolute value) and *cx* = 0.85 (see Figure 2). If the profile *U*(*y*) is sufficiently inflectional, the instability can generate sufficiently strong tilted waves (streaks and vortices) and trigger turbulent bands. Once triggered, the length of the band increases, and the force can be switched off after the band has sufficiently developed. The tilt direction of the band can be controlled by the signs of the spanwise component of *U*(*y*) and the moving speed *cz*.

We measured the advection speed of the low-speed streaks inside the bulk using the Structural Similarity Index Measure (SSIM) method proposed by Wang et al. [28], which is commonly used in image processing to measure the similarity between two images. The SSIM index is defined as:

$$SSIM(\mathbf{x}, y) = l(\mathbf{x}, y)^a \cdot \mathbf{c}(\mathbf{x}, y)^\beta \cdot \mathbf{s}(\mathbf{x}, y)^\gamma,\tag{3}$$

where *x* and *y* are one-dimensional vectors containing all the pixel values of the two images to be compared, respectively, and

$$d(x, y) = \frac{2\mu\_x \mu\_y + c1}{\mu\_x^2 + \mu\_y^2 + c1},\tag{4}$$

$$\mathbf{c}(\mathbf{x}, \mathbf{y}) = \frac{2\sigma\_{\mathbf{x}}\sigma\_{\mathbf{y}} + c\mathbf{2}}{\sigma\_{\mathbf{x}}^2 + \sigma\_{\mathbf{y}}^2 + c\mathbf{2}'},\tag{5}$$

and

$$s(x, y) = \frac{\sigma\_{xy} + c\mathfrak{A}}{\sigma\_{\mathfrak{X}}\sigma\_{\mathfrak{Y}} + c\mathfrak{A}},\tag{6}$$

measure the luminance, contrast and structural similarity, respectively. The exponents *α* > 0, *β* > 0 and *γ* > 0 are used to tune the relative weight of respective factor, and here we set all of them to 1 according to the suggestion of Wang et al. [28]. In Equation (4), *μ<sup>x</sup>* and *μ<sup>y</sup>* denote the mean of *x* and *y*, respectively. In Equation (5), *σ<sup>x</sup>* and *σ<sup>y</sup>* denote the standard deviation of *x* and *y*, respectively. In Equation (6), *σxy* is the covariance of *x* and *y*. Parameters *c*1 = (*k*1*L*)<sup>2</sup> and *c*2 = (*k*2*L*)2, where *k*1 and *k*2 are set to 0.01 and 0.03, respectively, and *L* is the maximum of the pixel value, which is set to 255 for unit8 data and 1 for floating point data. In our calculation, the flow velocities, which are floating point data, were taken as the pixel value *x* and *y*. The parameter *c*<sup>3</sup> is set such that *c*<sup>3</sup> = *c*2/2 in practice according to the suggestion of Wang et al. [28]. Thus, we have

$$SSIM(\mathbf{x}, y) = \frac{(2u\_{\mathbf{x}}u\_{\mathbf{y}} + c\_1)(2\sigma\_{\mathbf{x}\mathbf{y}} + c\_2)}{(u\_{\mathbf{x}}^2 + u\_{\mathbf{y}}^2 + c\_1)(\sigma\_{\mathbf{x}}^2 + \sigma\_{\mathbf{y}}^2 + c\_2)}.\tag{7}$$

The result is a value between −1 and 1, and the larger is the result, the higher is the similarity.

Firstly, we take the streamwise velocities in the cut plane *y* = −0.5 from two different snapshots *s*<sup>1</sup> and *s*<sup>2</sup> that are separated in time by *δt*, after the tilt angle of the band has stopped changing considerably due to the initial transients. In the frame of reference co-moving with the head, i.e., moving with a streamwise speed of 0.85 and a spanwise speed of −0.1 (we considered a right-going band), the bulk of the band is located in a nearly fixed area (see Figure 10). Therefore, we set a rectangular area in which the data inside were considered for calculating the SSIM index. We set the data outside this area to zero so that we eliminated the influence of the data outside this area. Further, to highlight the low-speed streaks, only the streamwise velocities in this area that satisfies *ux* < 0 and *u*<sup>2</sup> *<sup>x</sup>* > 0.002 were retained. Secondly, we shifted the data from *s*<sup>2</sup> inside the rectangular over the time separation *δt* with a streamwise speed *cx* and a spanwise speed *cz*. The original data from *s*<sup>1</sup> and the shifted data from *s*<sup>2</sup> were used to calculate the SSIM index. Thus, for a given speed pair (*cx*, *cz*), there is a corresponding SSIM index. By varying the speed pair, the SSIM index will maximize with certain speeds, which we considered as the mean advection speeds of the streaks. The contours of the SSIM index in the *cx* and *cz* for *Re* = 750 are shown in Figure 11.

Note that, in practice, we set −*cx* and *cz* to be between 0.1 and 0.4 (the band we considered is a right-going one; therefore, *cx* < 0 and *cz* > 0) because the actual speeds were estimated by eye to exist in this range, and note that the shift speeds are relative to the propagation of the head. Obtaining the contours of the SSIM index, we could estimate the advection speed of the streaks to be *cx* = −0.185 and *cz* = 0.18, i.e., the location of the local peak at the left-bottom corner in Figure 11. It can be seen that there

is another local peak at the right-top corner, which shows a lower SSIM index. That peak was reached when the *s*<sup>2</sup> data were shifted by more than one wave-length associated with the pattern of the low-speed streaks. The lower SSIM index of the top-right peak, i.e., lower similarity, indicates that the streaky pattern slowly change as it is advected in the bulk.

Note that the time separation *δt* between *s*<sup>1</sup> and *s*<sup>2</sup> cannot be too small, otherwise the streaks would have moved too little over the time separation and the speed measurement would be inaccurate. Likewise, it cannot be too large in which case the streaks would have moved by multiple wavelengths, which would also affect the speed calculation. In practice, estimated by eyes, a value between *δt* = 10 and 15 is a good choice, and *δt* = 10 in Figures 10 and 11. In the end, by varying the time instant of *s*1, we can obtain the average advection speed as a function of time and calculate the temporal average, which is plotted in Figure 7 (the speed of the head is added back in that figure).

**Figure 10.** The selection of low-speed streaks for the advection speed calculation. The tilt angle of this rectangle is 45◦, which is close to the tilt angle of the band. The area of the rectangle should be large enough to contain sufficient streaks and meanwhile reduce the influence of the tilt angle of the rectangle. (**a**) The contours of streamwise speedat *y* = −0.5 at *t* = 600. The advection speed of the streaks in the region enclosed by the rectangle is calculated; (**b**,**c**) The filtered low-speed streaks at *t* = 600 and *t* = 610.

**Figure 11.** The contours of the SSIM index in the *cx*-*cz* plane for the case shown in Figure 10.

**Author Contributions:** B.S. designed the research and X.X. performed the simulations. X.X. and B.S. analyzed the data and wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China under grant numbers 91852105 and 91752113 and by Tianjin University under grant number 2018XRX-0027.

**Acknowledgments:** The simulations were performed on Tianhe-2 supercomputer at the National Supercomputer Center in Guangzhou. Discussions with Dwight Barkley and Jianjun Tao are highly appreciated.

**Conflicts of Interest:** The authors declare no conflict of interest.
