2. Three-dimensional Kolmogorov flow at low Reynolds number

We consider a Kolmogorov flow on a periodic domain of dimensions $6\pi \times 2\pi \times 2\pi$ . The governing equation is the incompressible Navier–Stokes equation

(2.1a) \begin{align} \partial _t \boldsymbol u + \boldsymbol{u} {\cdot }\boldsymbol{\nabla }\boldsymbol{u} &= - \boldsymbol{\nabla }p + \nu {\nabla} ^2 \boldsymbol u - \bar {\mu } \bar {\boldsymbol u}+ f \sin \left ( k_f y \right )\boldsymbol e_x ,\end{align}

(2.1b) \begin{align} \boldsymbol{\nabla \boldsymbol{\cdot }u} &= 0. \end{align}

Here, $\boldsymbol u$ is the three-dimensional velocity field, $p$ is the (kinematic) pressure and $\nu$ is the kinematic viscosity. The flow is damped by a large-scale drag with damping parameter $\bar {\mu }$ where (Lalescu & Wilczek Reference Lalescu and Wilczek2021)

(2.2) \begin{equation} \bar {\boldsymbol u} = \sum _{k \lt k_f} \hat {\boldsymbol u }(\boldsymbol k) \mathrm{e}^{\mathrm{i}\boldsymbol k \boldsymbol{\cdot }\boldsymbol x} \end{equation}

only contains modes with wavenumbers smaller than the forcing wavenumber $k_f$ . Lalescu & Wilczek (Reference Lalescu and Wilczek2021) showed that the introduction of a drag term that only acts on the large scales could modify the number of vortex pairs. The Kolmogorov flow is driven by forcing with amplitude $f$ , which is constant in time and varies sinusoidally in a spanwise direction with wavenumber $k_f$ . We non-dimensionalise the equations using $L = 2\pi / k_f$ and $T = (L/f)^{1/2}$ as spatial and temporal scales, respectively. We define the forcing-based Reynolds number ${\textit{Re}}= UL / \nu = (2 \pi / k_f)^{3/2} f^{1/2}/ \nu$ , the non-dimensional large-scale drag coefficient $\mu = T \bar {\mu }$ and the aspect ratio $L_x/L_y = \alpha = 1/3$ . We keep the forcing amplitude fixed at $f = 1/2$ and the forcing wavenumber at $k_f = 1$ . By rescaling space as $\boldsymbol x = L \boldsymbol x^\prime$ , time as $t = T t^\prime$ and dropping the primes we obtain

(2.3) \begin{equation} \partial _t \boldsymbol u + \boldsymbol{u \boldsymbol{\cdot }\boldsymbol{\nabla }u} = - \boldsymbol{\nabla } p + \frac {1}{\textit{Re}} {\nabla} ^2 \boldsymbol u - \mu \bar {\boldsymbol u}+ \sin \left ( 2 \pi y\right ) \, \boldsymbol e_x . \end{equation}

For every value of $ Re$ and $\mu$ there exists a laminar solution given by $\boldsymbol u_{\mathrm{L}} = ({\textit{Re}}/{4\pi ^2}) \sin (2\pi y) \boldsymbol e_x,$ and homogeneous pressure $p = p_0$ . The stability of this laminar solution has been characterised several times, e.g. by Meshalkin & Sinai (Reference Meshalkin and Sinai1961), Thess (Reference Thess1992) and Chatterjee & Verma (Reference Chatterjee and Verma2020), to mention a few. We extend the aforementioned linear stability analysis to include the large-scale drag. The Kolmogorov flow with a drag acting on all scales has already been considered by Thess (Reference Thess1992), who showed that the drag shifts the critical Reynolds number at which the laminar solution becomes unstable. To study perturbations to the laminar solution it is convenient to note that, since the laminar solution of (2.3) only depends on the $y$ -coordinate, the most unstable perturbations to such a flows are independent of the $z$ -coordinate, a result known as Squire’s theorem (Squire Reference Squire1933). We determine the stability of the laminar solution by studying the linear evolution of a streamfunction $\psi (x,y,t)$ , from which the velocity is determined via $\boldsymbol u(x,y,t) = \boldsymbol u_{\mathrm{L}}(y) + \boldsymbol{\nabla }\times (\psi (x,y,t) \boldsymbol e_z)$ . By taking the curl of (2.3), we obtain the following equation for the streamfunction:

(2.4) \begin{align} \partial _t {\nabla} ^2 \psi &= \underbrace {\left (\frac {1}{\textit{Re}} {\nabla} ^4 -{\textit{Re}}\sin {( 2\pi y )} \left (1 + \frac {1}{4\pi ^2}{\nabla} ^2 \right )\partial _x \right ) \psi - \mu {\nabla} ^2 \bar {\psi }}_{\mathcal{L}\left (\mu , \boldsymbol{Re}; \psi \right )}\notag\\&\quad -\underbrace { \left ( \partial _y \psi \partial _x - \partial _x \psi \partial _y \right ) {\nabla} ^2 \psi }_{\mathcal{N}(\psi )} ,\end{align}

where we have used $\omega _z = (\boldsymbol{\nabla }\times \boldsymbol{u}) \boldsymbol{\cdot }\boldsymbol{e}_z = -(Re/2 \pi) \cos({2 \pi y})- {\nabla} ^2 \psi$ . We study perturbations of the form

(2.5) \begin{equation} \psi _n(x,y,t) = \hat {\psi }_n(y) \mathrm{e}^{2\pi \mathrm{i} \alpha x + \lambda _n(\boldsymbol{Re} , \mu ) t} \text{,} \end{equation}

with $n \in \mathbb{Z}$ . Each $ \psi _n(x,y,t)$ corresponds to an eigenvector of the operator $\mathcal{L}(\mu , {\textit{Re}}; \psi )$ . A mode $\psi _n$ becomes unstable when the real part of $\lambda _n( Re, \mu )$ crosses $0$ from below. For the aspect ratio ( $\alpha = 1/3$ ) we are considering, only the modes $\psi _1$ and $\psi _2$ can become unstable, namely, for every value of $\mu$ there exist $ {\textit{Re}}_1(\mu )$ and $ {\textit{Re}}_2(\mu )$ such that $\mathrm{\textit{Re}} \left \{\lambda _1( {\textit{Re}}_1(\mu ), \mu )\right \} =\mathrm{\textit{Re}} \left \{\lambda _1( {\textit{Re}}_2(\mu ), \mu )\right \}=0$ . When $\psi _1$ becomes unstable for sufficiently high $ Re$ , a new state with a single vortex pair emerges as a secondary flow. On the other hand, when $\psi _2$ becomes unstable, a state with two vortex pairs emerges instead. The linear stability analysis does not tell us whether these new states are stable solutions of the Navier–Stokes equations (2.3), we will determine the stability of these states using direct numerical simulations (DNS) and a weakly nonlinear analysis. In the frictionless case, perturbations with the smallest streamwise wavenumber are the most unstable ones (Meshalkin & Sinai Reference Meshalkin and Sinai1961), namely $ {\textit{Re}}_1(0)\lt {\textit{Re}}_2(0)$ , therefore the one-vortex-pair configuration is preferred at low $ Re$ . However, for $\mu \gt \mu _c \approx 0.90$ , we have $ {\textit{Re}}_2(\mu )\lt {\textit{Re}}_1(\mu )$ , and the two-vortex-pair configuration is preferred. At $\mu =\mu _c$ we have

(2.6) \begin{align} {\textit{Re}}_2(\mu _c) = {\textit{Re}}_1(\mu _c) = {\textit{Re}}_c ,\end{align}

which means that both states can emerge when ${\textit{Re}}\gt {\textit{Re}}_c$ . Analytical expressions for $\psi _1(x,y;\mu )$ and $\psi _2(x,y;\mu )$ are approximately given by

(2.7a) \begin{align} \psi _1(x,y;\mu ) &\approx \left ( 2 \mathrm{i}\sin {(2\pi y)} + \frac {16 \pi ^3\left (1 + \alpha ^2 \right )^2}{\alpha {\textit{Re}}_{1}^2(\mu )\left ( 1 - \alpha ^2 \right )} \right ) \mathrm{e}^{2 \pi \mathrm{i}\alpha (x -x_1) } , \end{align}

(2.7b) \begin{align} \psi _2(x,y;\mu ) &\approx \left ( 2 \mathrm{i}\sin {(2\pi y)} + \frac {8 \pi ^3\left (1 + 4\alpha ^2 \right )^2}{\alpha {\textit{Re}}_{2}^2(\mu )\left ( 1 - 4\alpha ^2 \right )} \right ) \mathrm{e}^{4 \pi \mathrm{i}\alpha (x -x_2)} , \end{align}

Figure 1.

(a) Regime map of the three-dimensional Kolmogorov flow at low Reynolds numbers, along with visualisations of the possible states. The blue curve with squares ( $ {\textit{Re}}_1^*(\mu )$ ) corresponds to the lower stability threshold of the solution with one vortex pair, while the orange curve with squares ( $ {\textit{Re}}_2^*(\mu )$ ) corresponds to the lower stability threshold of the solution with two vortex pairs. Note that for $\mu \lt \mu _c \approx 0.90$ , $ {\textit{Re}}_1^*(\mu )$ coincides with $ {\textit{Re}}_1(\mu )$ plotted in a dotted blue curve corresponding to the $ Re$ value where the state with one-vortex-pair emerges in a supercritical bifurcation, while for $\mu \gt \mu _c$ , $ {\textit{Re}}_2^*(\mu )$ coincides with $ {\textit{Re}}_2(\mu )$ , plotted in a dotted orange curve. Vertical black lines correspond to the section where we show the corresponding bifurcation diagrams in figure 2. (b) Visual representation of the possible flow configuration. The vertical component of the vorticity $\omega _z$ is coded in the colour map superimposed by a visualisation of the velocity field.

Figure 1. Long description

The image contains two main sections. The left section features a line graph showing the regime map of the three-dimensional Kolmogorov flow at low Reynolds numbers. The x-axis represents the parameter mu, ranging from 0 to 2, and the y-axis represents the Reynolds number, ranging from 20 to 45. The graph includes two curves: a blue curve with squares indicating the lower stability threshold of the solution with one vortex pair, and an orange curve with squares indicating the lower stability threshold of the solution with two vortex pairs. The graph also highlights a bistable zone and regions of stable laminar flow and one vortex pair. The right section features visual representations of the possible flow configurations. The top image shows no vortex pairs, the middle image shows one vortex pair, and the bottom image shows two vortex pairs. Each image depicts the vertical component of the vorticity with a color map and a visualization of the velocity field. The visualizations illustrate the flow patterns and the emergence of coherent structures in the Kolmogorov flow.

where $x_1$ and $x_2$ characterise the position along the $x$ -axis of the emerging unstable solutions. The expressions are approximated in the sense that the number of Fourier modes in the $y$ -direction has been truncated, see Appendix B for a detailed derivation of these expressions along with the computation of $\mu _c$ and $ {\textit{Re}}_c$ . Combining these results with DNS allows us to obtain a regime map of this system at a low Reynolds number, which we display in figure 1 along with snapshots of the possible flow configurations in this region of the parameter space $( \mu ,{\textit{Re}})$ . We construct this regime map in the following way: we run two sets of simulations with different initial conditions, one starting in the one-vortex-pair configuration, and the other in the two-vortex-pair configuration. We determine the streamfunction $\psi (x,y)$ , project it onto $\psi _1(x,y)$ and $\psi _2(x,y)$ and define the following amplitudes:

(2.8a) \begin{align} a_1 &= \frac {\int _{0}^{1} \int _{0}^{3} \psi _1(x,y)^* \psi (x,y) \mathrm{d}x \mathrm{d}y}{\int _{0}^{1} \int _{0}^{3}|\psi _1(x,y)|^2 \mathrm{d}x \mathrm{d}y} , \end{align}

(2.8b) \begin{align} a_2 &= \frac {\int _{0}^{1} \int _{0}^{3} \psi _2(x,y)^* \psi (x,y) \mathrm{d}x \mathrm{d}y}{\int _{0}^{1} \int _{0}^{3}|\psi _2(x,y)|^2 \mathrm{d}x \mathrm{d}y} . \end{align}

Figure 2.

(a) Bifurcation diagram of the flow at $\mu = 0.25$ . Solid lines represent stable states while dashed lines represent unstable states. At $ {\textit{Re}}_1$ the flow undergoes a pitchfork bifurcation which gives rise to a configuration with one vortex pair. At $ {\textit{Re}}_2$ the flow undergoes another pitchfork bifurcation which leads to a configuration with two vortex pairs, however, this state is unstable. At $ {\textit{Re}}_2^*$ , there is a saddle-node bifurcation that gives rise to a mixed state where $a_1 \neq 0$ and $a_2 \neq 0$ . After this bifurcation, the two-vortex-pair configuration becomes stable while the mixed state remains unstable. For $ Re\gt {\textit{Re}}_2^*$ the flow exhibits either one or two vortex pairs depending on the initial condition. (b) Bifurcation diagram of the flow at $\mu = 1.5$ , the situation is analogous to the plot on the left-hand side, but now for $ {\textit{Re}}_1 \lt{\textit{Re}}\lt {\textit{Re}}_1^*$ the one-vortex-pair configuration is unstable and a saddle-node bifurcation gives rise to an unstable mixed state at $ {\textit{Re}}_1^*$ .

Figure 2. Long description

The image contains two bifurcation diagrams side by side. The left diagram shows the flow at a specific condition, with solid lines representing stable states and dashed lines representing unstable states. At a Reynolds number of approximately 25, the flow undergoes a pitchfork bifurcation, leading to a configuration with one vortex pair. At a Reynolds number slightly below 30, another pitchfork bifurcation occurs, resulting in a configuration with two vortex pairs, which is unstable. At a Reynolds number of approximately 32, a saddle-node bifurcation gives rise to a mixed state where both configurations coexist, with the two-vortex-pair configuration becoming stable. The right diagram shows the flow at a different condition, where the one-vortex-pair configuration is unstable, and a saddle-node bifurcation gives rise to an unstable mixed state at a Reynolds number of approximately 40. The diagrams illustrate how different flow states emerge and stabilize at and values of the large-scale drag.

Since the modes $\psi _1(x,y)$ and $\psi _2(x,y)$ are made from different Fourier modes, they are orthogonal to each other under this inner product, namely $\int _{0}^{1} \int _{0}^{3} \psi _1(x,y)^* \psi _2(x,y) \mathrm{d}x \mathrm{d}y = 0$ . If $|a_1|\neq 0$ and $|a_2| = 0$ , we say that the flow is in the one-vortex-pair configuration, while if $|a_2|\neq 0$ and $|a_1| = 0$ , the flow is in the two-vortex-pair configuration. In principle, mixed states where $|a_2|\neq 0$ and $|a_1| \neq 0$ would also be possible; however, numerical simulations reveal that there are no such states that are stable. We then start each set of simulation at fixed value of $\mu$ in a state with $|a_1|\neq 0$ and $|a_2| = 0$ ( $|a_2|\neq 0$ and $|a_1| = 0$ ), decrease $ Re$ and determine $ {\textit{Re}}_1^*(\mu )$ ( $ {\textit{Re}}_2^*(\mu )$ ) when $|a_1| = 0$ ( $|a_2|\neq 0$ ). We present these results in figure 1. We note that when $\mu \lt \mu _c$ , the two-vortex-pair configuration becomes unstable at $ {\textit{Re}}_2^*(\mu ) \gt {\textit{Re}}_2(\mu )$ , namely, at a higher Reynolds number than the one predicted by the linear theory, while for $\mu \gt \mu _c$ , the two curves coincide, namely $ {\textit{Re}}_2^*(\mu ) = {\textit{Re}}_2(\mu )$ . Analogously, for $\mu \gt \mu _c$ , we have $ {\textit{Re}}_1^*(\mu ) \gt {\textit{Re}}_1(\mu )$ instead, while $ {\textit{Re}}_1^*(\mu ) = {\textit{Re}}_1(\mu )$ for $\mu \lt \mu _c$ . For $ Re\gt {\textit{Re}}_1^*(\mu ), {\textit{Re}}_2^*(\mu )$ both states are equilibrium states of the flow. Whether the flow displays one or two vortex pairs depends on the initial conditions, and thus we label this region the bistable zone in figure 1. To clarify the type of bifurcations the flow undergoes when changing $ Re$ , in figure 2 we show the bifurcation diagram at two values of $\mu$ . On the left plot, we consider $\mu = 0.25 \lt \mu _c$ . Solid lines represent stable equilibrium values of $|a_1|$ and $|a_2|$ , while the dashed ones represent unstable ones. In this region of the parameter space, we have

(2.9) \begin{align} {\textit{Re}}_1 \lt {\textit{Re}}_2 \lt {\textit{Re}}_2^* . \end{align}

When ${\textit{Re}}\gt {\textit{Re}}_1$ , the configuration with one vortex pair is always a stable equilibrium; for $ {\textit{Re}} \in ( {\textit{Re}}_2, {\textit{Re}}_2^*)$ , the two-vortex-pair solution exists but is only an unstable equilibrium; At ${\textit{Re}}= {\textit{Re}}_2^*$ , the two-vortex-pair configuration becomes stable and an unstable mixed state, namely a state where $|a_1|, |a_2| \neq 0$ emerges. In the right plot we consider $\mu = 1.5 \gt \mu _c$ , where we have

(2.10) \begin{align} {\textit{Re}}_2 \lt {\textit{Re}}_1 \lt {\textit{Re}}_1^* . \end{align}

The situation is analogous to one we just described: For ${\textit{Re}}\gt {\textit{Re}}_2$ , the two-vortex-pair configuration is always a stable equilibrium; for ${\textit{Re}}\in ( {\textit{Re}}_1, {\textit{Re}}_1^*)$ , the one-vortex-pair configuration is an unstable equilibrium; at ${\textit{Re}}= {\textit{Re}}_1^*$ , the one-vortex-pair configuration becomes stable, and again, an unstable mixed state emerges.

The bifurcation diagrams were obtained in the following way: we start by finding the vorticity of the two equilibrium configurations for a given $ Re\gt {\textit{Re}}_1^*, Re^*_2$ . Let $\omega _1( Re)$ and $\omega _2( Re)$ be these vorticity configurations. To find the unstable saddle point, namely the unstable mixed state, we conduct multiple simulations using

(2.11) \begin{align} \omega _z(x,y) = \alpha \omega _1( Re) + (1 -\alpha )\omega _2( Re) \quad \text{with}\quad \alpha \in [0,1] \end{align}

as an initial condition for the vorticity. We then solve the Navier–Stokes equation in the vorticity formulation and keep track of whether the flow converges back to $\omega _1( Re)$ or to $\omega _2( Re)$ for each value of $\alpha$ . This way, we identify a point in the separatrix between the two configurations in the $(|a_1|, |a_2|)$ plane. We finally set an initial condition lying in the vicinity of the separatrix and track the trajectory $(a_1(t), a_2(t))$ . We identify the saddle point as the slowest point of this trajectory, i.e. the point where $\sqrt {(\mathrm{d}a_1/\mathrm{d}t )^2 + (\mathrm{d}a_2/\mathrm{d}t )^2}$ is minimised. We call this mixed state a saddle point because the trajectory $(a_1(t),a_2(t))$ approaches it, slows down and then is repelled so that $(a_1(t),a_2(t))$ ends up in one of the stable vortex configurations.

At $\mu = \mu _c$ the curves $ Re(\mu )$ and $ {\textit{Re}}_2(\mu )$ intersect. At this point, both stable solutions emerge in supercritical bifurcations. We call the point $(\mu _c, {\textit{Re}}_c)$ the onset of bistability, because the laminar solution becomes unstable simultaneously to $\psi _1$ and $\psi _2$ in a supercritical bifurcation. At this point we can study how the modes $\psi _1$ and $\psi _2$ interact with each other by choosing $\mu = \mu _c + \varepsilon \mu _0$ and ${\textit{Re}}= {\textit{Re}}_{c}( 1 +\varepsilon r )$ , where $\varepsilon \ll 1$ is a scaling parameter and $\mu _0$ and $r$ are arbitrary values of order $1$ . Under these conditions, we introduce the following ansatz into (2.4):

(2.12) \begin{equation} \psi = \varepsilon ^{1/2} \left ( a_1(t)\psi _1 + a_2(t)\psi _2 + \mathrm{c.c.} \right ) + w(a_1, a_2) ,\end{equation}

where $a_1(t),a_2(t)$ are complex amplitudes of the two modes, c.c. stands for complex conjugate, and $w(a_1, a_2)\sim \mathcal{O}(\varepsilon )$ is a nonlinear function of the amplitudes. In Appendix C, we employ normal form theory (Carr Reference Carr1981; Elphick et al. Reference Elphick, Tirapegui, Brachet, Coullet and Iooss1987; Haragus & Iooss Reference Haragus and Iooss2011) to derive a set of gradient equations describing the evolution of amplitudes

(2.13a) \begin{align} \frac {\mathrm{d} a_1}{\mathrm{d}t^*} &= \delta _1 a_1 - \gamma _1 |a_1|^ 2a_1 - \eta _1 |a_2|^ 2a_1 = - \eta _1 \frac {\partial \mathcal{F}}{\partial a_1^*} ,\end{align}

(2.13b) \begin{align} \frac {\mathrm{d} a_2}{\mathrm{d}t^*} &= \delta _2 a_2 - \gamma _2 |a_2|^ 2a_2 - \eta _2 |a_1|^ 2a_2 = - \eta _2 \frac {\partial \mathcal{F}}{\partial a_2^*}\mathrm{.} \end{align}

Here, $t^* = \varepsilon t$ is a slow time. All the coefficients are real numbers implying that these two equations minimise a functional $\mathcal{F}$ given by

(2.14) \begin{equation} \mathcal{F} = - \delta _1 |a_1|^2/\eta _1 - \delta _2 |a_2|^2/\eta _2 + \gamma _1 |a_1|^4/2\eta _1 + \gamma _2 |a_2|^4/2\eta _2 + |a_1|^2 |a_2|^2/2 ,\end{equation}

where $\delta _1 = \delta _1(\mu _0, r)$ and $\delta _2 = \delta _2(\mu _0, r)$ . Expressions for both functions, as well as expressions for $\gamma _1, \gamma _2,\eta _1$ and $\eta _2$ , as functions of $\mu _c$ and $ {\textit{Re}}_c$ can be found in Appendix C. Such amplitude equations can also be obtained purely by symmetry arguments. The Navier–Stokes equation (2.1a ) remains invariant under translations along the $x$ -axis. This means that vortex pairs can arise at any position along this axis. By looking at the unstable modes in (2.7), it becomes clear that displacing the vortex pairs along the $x$ -axis is equivalent to a phase shift of the complex exponential. These additional phases can then be absorbed into $a_1$ and $a_2$ . This implies that their respective evolution equations must be phase invariant so that all positions along the $x$ -axis are equally accessible. The phase invariance demands that only terms like $|a_1|^n a_1$ , $|a_2|^n a_2$ , $|a_1|^m |a_2|^n a_1$ and $|a_1|^n |a_2|^m a_1$ can appear in the amplitude equations, which is consistent with our result. The symmetry approach provides a basic structure for the amplitude equation, but gives no information about the values of the coefficients. We find that all the coefficients in (2.13) are real and positive, which restricts the range of behaviours that such an equation can exhibit. To summarise the implications of (2.13):

1. (i) When $\delta _1 \lt 0$ and $\delta _2 \lt 0$ , (2.13) have only a one stable equilibrium point: the origin $(|a_1|, |a_2|) = (0,0)$ representing the laminar state with no vortex pairs.

2. (ii) When $\delta _1$ becomes positive, the equilibrium lying on the $|a_1|$ -axis emerges in a supercritical bifurcation, namely at $(|a_1|,|a_2|) =(\sqrt {({\delta _1}/{\gamma _1})}, 0)$ representing the one-vortex-pair configuration. This emerging equilibrium is stable if $( {\gamma _1}/{\eta _2}) \lt ( {\delta _1}/{\delta _2})$ or if $\delta _2 \lt 0$ , and unstable otherwise.

3. (iii) When $\delta _2$ becomes positive, the equilibrium lying on the $|a_2|$ -axis emerges in a supercritical bifurcation, namely at $(|a_1|,|a_2|) = (0 ,\sqrt {( {\delta _2}/{\gamma _2})})$ representing the two-vortex-pair configuration. This emerging equilibrium is stable if $({\delta _1}/{\delta _2}) \lt ( {\eta _1}/{\gamma _2})$ or if $\delta _1 \lt 0$ , and unstable otherwise.

4. (iv) When $( {\gamma _1}/{\eta _2}) \lt ({\delta _1}/{\delta _2}) \lt ({\eta _1}/{\gamma _2})$ , there is an additional equilibrium point at $( {\delta _2 \eta _1 - \delta _1 \eta _2})/({\eta _1 \eta _2 - \gamma _1 \gamma _2}), ({\delta _1 \eta _2 - \delta _2 \eta _1})/({\eta _1 \eta _2 - \gamma _1 \gamma _2})$ . This equilibrium corresponds to the mixed state and is a saddle point, i.e. it is an unstable equilibrium.

These analytical results are in agreement with the results we present in the regime map in figure 1 and the bifurcation diagrams in figure 2. The system evolves until reaching a stationary state corresponding to a minimum of the functional. Even though amplitude equations are, in principle, valid only around the point $(\mu _c, {\textit{Re}}_c)$ , our numerical study suggests that the structure of the equation remains valid further from this point and even when turbulence develops.

3. Connecting turbulent superstructures with the dynamics at low Reynolds number

For large Reynolds numbers, the flow is fully three-dimensional, and the results of the linear stability analysis are no longer valid. However, we observe persistent large-scale vortex pairs (Sarris et al. Reference Sarris, Jeanmart, Carati and Winckelmans2007; Lalescu & Wilczek Reference Lalescu and Wilczek2021). The large-scale vortices can be isolated using a spatial averaging procedure. Following Lalescu & Wilczek (Reference Lalescu and Wilczek2021), we define the vertically averaged velocity field as $\tilde {\boldsymbol u} = \int _0^{1} \boldsymbol u(x,y,z)\,\mathrm{d} z$ . We generate numerical data using the code TurTLE (Lalescu et al. Reference Lalescu, Bramas, Rampp and Wilczek2022) for two values of the Reynolds number ( ${\textit{Re}}= 2851, 7184$ ) and for a range of values of the non-dimensional friction coefficient ( $\mu \in [0.0, 1.24]$ ), see Appendix A for details. To illustrate the different large-scale structures and how they coexist with a fully three-dimensional flow, we show snapshots of two different instants of a simulation of the three-dimensional Kolmogorov flow at $\mu = 0.89$ and ${\textit{Re}}= 2851$ in figure 3. In figure 3 (a), we show a volume rendering of the squared vorticity $ \omega ^2 = (\boldsymbol {\nabla }\times \boldsymbol u)^2$ to illustrate the complexity of the three-dimensional flow. In figure 3 (b), we display the streamlines of the averaged horizontal velocity $\tilde {\boldsymbol u }_\bot = \tilde {u}_x \boldsymbol e_x + \tilde {u}_y \boldsymbol e_y$ superimposed with a colour map of the averaged $z$ -component of the vorticity $\tilde {\omega }_z = \partial _x \tilde {u}_y - \partial _y \tilde {u}_x$ which makes large-scale vortex pairs readily visible. We observe two types of large-scale structures, one with a single pair of vortices and another displaying two vortex pairs. These states show a strong resemblance with the attractors that emerge after the laminar solution loses stability at a low Reynolds number, as discussed in § 2 (see figure 1). Such structures have also been reported in higher-resolution simulations with a Reynolds number up to ${\textit{Re}}= 18\,100$ (Lalescu & Wilczek Reference Lalescu and Wilczek2021). To quantify the resemblance between the dynamics at high and low Reynolds number we define the following velocity fields: $\boldsymbol u_0 = \sin (2 \pi y) \boldsymbol e_x$ , $\boldsymbol u_{1} = \boldsymbol{\nabla }\times (\psi _1 \boldsymbol e_z)$ and $\boldsymbol u_{2} = \boldsymbol{\nabla }\times (\psi _2 \boldsymbol e_z)$ which correspond to velocity fields parallel to the forcing, the complex velocity generated by the streamfunction $\psi _1$ and the complex velocity field generated by $\psi _2$ , respectively. In addition, we introduce an inner product for two-component vector fields $\langle \boldsymbol F, \boldsymbol G\rangle = 1/3\int _0^{1} \int _0^{3} (F_x G_x^* + F_y G_y^* )\mathrm{d}x \, \mathrm{d}y$ . The velocity fields $\boldsymbol u_0, \boldsymbol u_1$ , and $\boldsymbol u_2$ are orthogonal to each other with respect to this inner product. We use these large-scale velocity fields to decompose the average velocity

(3.1) \begin{align} \tilde {\boldsymbol u} = \underbrace {A_0(t)\boldsymbol u_0 + \left [A_1(t)\boldsymbol u_1+ A_2(t)\boldsymbol u_2 + \mathrm{c.c.}\right ]}_{\boldsymbol u_{\textit{TS}}(t,x,y)}\,+\,{\boldsymbol u}_{{R}}(t,x,y) + \tilde {u}_z(t,x,y) \boldsymbol e_z ,\end{align}

where $A_i(t) = \langle \boldsymbol u_i, \boldsymbol u \rangle / \langle \boldsymbol u_i, \boldsymbol u_i \rangle$ for $i=0,1,2$ . We name $\boldsymbol u_{\textit{TS}}$ the component of the average velocity field that belongs to the subspace generated by $\boldsymbol u_0, \boldsymbol u_1$ and $\boldsymbol u_2$ , where TS stands for turbulent superstructures. That means we effectively identify the vortices at low Reynolds numbers with the turbulent superstructures at high Reynolds numbers. By construction, we have that the remainder velocity field $\boldsymbol u_{{R}}$ is orthogonal to both $\boldsymbol u_{\textit{TS}}$ and $\tilde {u}_z \boldsymbol e_z$ . We use the decomposition (3.1) to study the instantaneous two-dimensional energy decomposition, which is given by

(3.2) \begin{equation} E_{2{D}}(t) = \frac {1}{3}\int _{0}^{1} \int _{0}^{3} \frac {1}{2}\left ( \boldsymbol u_{\textit{TS}}^2 + \boldsymbol u_{{R}}^2 + \tilde {u}_z^2\right ) \mathrm{d}x \, \mathrm{d}y = E_{\textit{TS}}(t) + E_{{R}}(t) + E_z(t) . \end{equation}

In figure 4 we display the time evolution of the instantaneous ratio between $E_{\textit{TS}}, E_{{R}}, E_z$ and $E_{{2D}}$ for the same simulation shown in figure 3. We observe that, on average, around $76\,\%$ of the two-dimensional energy is contained within the large-scale velocity modes. We consider the total energy decomposition as well

(3.3) \begin{equation} E_{3{D}}(t) = \frac {1}{3}\int _{0}^{1} \int _{0}^{1} \int _{0}^{3} \frac {1}{2}\left ( \tilde {\boldsymbol u}^2 + \boldsymbol u^{\prime 2}\right ) \mathrm{d}x\, \mathrm{d}y\,\mathrm{d}z = E_{2{D}}(t) + E^\prime (t), \end{equation}

where $\boldsymbol{u}^\prime = \boldsymbol{u} - \tilde{\boldsymbol{u}}$ . We find that energy in the turbulent superstructures makes up most of the two-dimensional energy for a wide range of $\mu$ and, therefore, a significant fraction of the total energy; we present these results in figure 5.

Figure 3.

(a) Volume rendering of the squared vorticity $\omega ^2$ at two different instants for the simulation at ${\textit{Re}}= 2851$ and $\mu = 0.89$ (see Appendix A). (b) Corresponding streamlines of the averaged horizontal velocity $\tilde {\boldsymbol u}_\bot (x,y)$ superimposed on the averaged $z$ -component of the vorticity to the volume renderings above.

Figure 3. Long description

The image contains two sets of visualizations related to turbulent Kolmogorov flow. The top row features volume renderings of the squared vorticity at two different time instants for a specific simulation. The color scale on the right indicates the intensity of the squared vorticity, ranging from 0 to 1. The bottom row displays the corresponding streamlines of the averaged horizontal velocity superimposed on the averaged z-component of the vorticity. The color scale on the right of the bottom row indicates the intensity of the z-component of the vorticity, ranging from −2.5 to 2.5. The streamlines illustrate the flow patterns and the vorticity component provides insights into the turbulent structures within the flow. The visualizations highlight the complex spatio-temporal dynamics and the emergence of large-scale flow structures in turbulent Kolmogorov flow.

Figure 4.

Time evolution of the energy decomposition of the averaged two-dimensional flow for a simulation at ${\textit{Re}}= 2851$ and $\mu = 0.89$ . Each time series represents the share of each term in (3.2) of the total two-dimensional energy $E_{2{D}}$ . Here, $E_{\textit{TS}}$ is the energy contained in the turbulent superstructures, $E_{{R}}$ is the energy contained in the remaining modes of the averaged horizontal velocity $\tilde {\boldsymbol u}_\bot$ which are orthogonal to $\boldsymbol{u}_{\textit{TS}}$ and $E_z$ is the energy contained in the averaged vertical velocity $u_z$ . The horizontal lines correspond to the time averages of the different fractions.

Figure 4. Long description

The line graph displays the time evolution of energy decomposition in a flow simulation. The x-axis represents time normalized by a time scale T, ranging from 0 to 1500. The y-axis represents the energy fractions normalized by the total two-dimensional energy E2D. Three data lines are present: a blue line representing the energy contained in turbulent superstructures (ETS/E2D), an orange line representing the energy contained in the remaining modes of the averaged horizontal velocity (ER/E2D), and a green line representing the energy contained in the averaged vertical velocity (Ez/E2D). The blue line fluctuates around a time average of 0.76, the orange line fluctuates around a time average of 0.19, and the green line fluctuates around a time average of 0.05.

Figure 5.

Time-averaged energy ratios of turbulent superstructures. Shown are the fractions of energy contained in the large-scale modes relative to the total energy $E_{3D}$ and to the two-dimensional average $E_{2D}$ .

Figure 5. Long description

The image contains two line graphs side by side, each representing the time-averaged energy ratios of turbulent superstructures at different Reynolds numbers. The left graph corresponds to a Reynolds number of 2851, while the right graph corresponds to a Reynolds number of 7184. Both graphs plot the fractions of energy contained in the large-scale modes relative to the total energy and to the two-dimensional average. The x-axis in both graphs represents the variable mu, ranging from 0 to 1.25. The y-axis represents the energy ratios, ranging from 0 to 1. In each graph, two data series are plotted: one with blue stars and another with orange circles. The blue stars represent the ratio of energy in turbulent superstructures to the three-dimensional total energy, while the orange circles represent the ratio of energy in turbulent superstructures to the two-dimensional average energy. In the left graph, the orange circles show an increasing trend as mu increases, starting around 0.6 and reaching approximately 0.8. The blue stars show a relatively stable trend around 0.4. In the right graph, the orange circles also show an increasing trend, starting around 0.6 and reaching approximately 0.85. The blue stars again show a relatively stable trend around 0.4. The graphs illustrate how the energy ratios vary with different values of mu and Reynolds numbers, highlighting the presence of turbulent superstructures in the flow.

Figure 6.

(a) Time evolution of the amplitude $A_0(t)$ , with the best-fit linear regression from (3.4) superimposed for a simulation at ${\textit{Re}}= 2851$ and $\mu = 0.89$ . (b) The $10T$ -averaged signal $\langle A_0(t) \rangle _{10T}$ compared with the averaged regression.

Figure 6. Long description

Two line graphs. The top graph shows the time evolution of the amplitude with the best-fit linear regression superimposed. The x-axis represents time normalized by T, ranging from 0 to 1500, and the y-axis represents the amplitude A0(t), ranging from 0.5 to 2.5. The orange line represents the best-fit linear regression equation 1.74 minus 6.74 times the square of the absolute value of A 1 minus 10.79 times the square of the absolute value of A2. The blue line represents the actual data points. The bottom graph shows the signal averaged over a time window of length 10T with the averaged regression. The x-axis represents time normalized by T, ranging from 0 to 1500, and the y-axis represents the T-averaged amplitude A0(t), ranging from 0.5 to 2.5. The orange line represents the averaged regression equation 1.74 minus the average of 6.74 times the square of the absolute value of A1 averaged over a time window of 10T plus 10.79 times the square of the absolute value of averaged over a time window of 10T. The blue line represents the actual averaged data points.

Figure 7.

Time series of the evolution of the moduli of the complex amplitude $A_1(t)$ and $A_2(t)$ for a simulation at ${\textit{Re}}= 2851$ and $\mu = 0.89$ . The vertical lines correspond to the time stamps at where we extracted the visualisations in figure 3.

Figure 7. Long description

The two line graph displays the evolution of the moduli of the complex amplitudes over time in a simulation. The x-axis represents time normalized by T, ranging from 0 to 750. The y-axis represents the moduli values, ranging from 0 to 0.4. Two data series are plotted: one in blue representing the modulus of the complex amplitude A1(t), and the other in orange representing the modulus of the complex amplitude A2(t). The graph shows fluctuations in both data series over time, with notable interactions and variations. Vertical dashed lines at specific time stamps indicate where visualizations were extracted for further analysis.

To simplify the problem even further, we conjecture that $A_0(t)$ is not an independent degree of freedom but instead is a function of the instantaneous values $A_1(t)$ and $A_2(t)$ . We base this assumption on the fact that, at low $ Re$ , disturbances to the laminar solution are proportional to $|a_1|^2$ and $|a_2|^2$ , as shown in Appendix C. We propose an analogous decomposition

(3.4) \begin{equation} A_0(t) = c_0 + c_1 |A_1(t)|^2 + c_2 |A_2(t)|^2 + \mathrm{h.o.t.} \end{equation}

Here, the set of coefficients $\{ c_i \}$ is determined using a linear regression. In figure 6, we keep track of the time evolution of the amplitude $A_0(t)$ , superimposed with the best linear regression. We find a good agreement between the linear regression and the instantaneous value of $A_0(t)$ up to short-term fluctuations. This agreement becomes clearer after eliminating the fluctuations by averaging the data over a time window of length $10T$

(3.5) \begin{equation} \langle A_0(t)\rangle _{10T} = \frac {1}{10T} \int _{-5T}^{5T} A_0(t - t^\prime )\mathrm{d}t^\prime .\end{equation}

The time-averaged data and linear regression results are displayed in figure 6(b). From these results, we conclude that the amplitude $A_0$ is not the independent degree of freedom but instead enslaved to the evolution of $A_1$ and $A_2$ . We focus on the time evolution of these amplitudes; in figure 7, we show the temporal evolution of the moduli $|A_1|$ and $|A_2|$ , which exhibit two distinct time scales. On the slowest time scale of the order of $100T$ , there are time intervals where $|A_1| \gt |A_2| \sim 0$ and vice versa. The times $t_1$ and $t_2$ are the time stamps corresponding to the snapshots of the flow we display in figure 3. There is a correspondence between the number of large-scale vortices and the dominating amplitude. The system spends a finite time around one state before switching into the other, a phenomenon referred to as dynamical switching, which has already been observed by Lalescu & Wilczek (Reference Lalescu and Wilczek2021). We observe fast fluctuations around the equilibrium states on a shorter time scale than the switching scale. We proceed to characterise the different regions of the $(A_1, A_2)$ space according to the number of large-scale vortices; for this, we take a sample of long time series of $A_1(t)$ and $A_2(t)$ and use the classification algorithm KMeans (MacQueen Reference MacQueen1967) to partition the space into two disjoint regions $\varOmega _1$ and $\varOmega _2$ . The result of the partition is displayed in figure 8: the upper row shows scatter plots of the real and imaginary parts of $A_1$ and $A_2$ , respectively. The blue points correspond to the sample points belonging to $\varOmega _1$ , which we identify as the instants where the flow is in a one-vortex-pair configuration, while the orange stars correspond to the samples belonging to $\varOmega _2$ , which we associate with the two-vortex-pair configuration. In the lower row of figure 8, we present histograms of the moduli $|A_1|$ and $|A_2|$ , the blue bars correspond to the samples in $\varOmega _1$ , while the orange ones to samples in $\varOmega _2$ . These histograms emphasise the switching between the dominant mode. Using conditional averages over these two sets, we can extract the modulus of equilibria around which the amplitudes oscillate, for $\varOmega _1$ we find that $A_1^{\textit{mean}} = 0.31 \pm 0.04$ and for $\varOmega _2$ we find $A_2^{\textit{mean}} = 0.14 \pm 0.05$ .

Figure 8.

(a) Partition of $(A_1,A_2)$ using the KMeans algorithm from the data of the simulation at ${\textit{Re}}= 2851$ and $\mu = 0.89$ , we observe that the region $\varOmega _1$ is characterised by $A_1$ oscillating around a ring of radius $A_1^{\textit{mean}}$ while $A_2$ oscillates around the origin of the complex plane. In the region $\varOmega _2$ , we observe that $A_2$ oscillates around a ring of radius $A_2^{\textit{mean}}$ and it is $A_1$ that remains close to the origin. (b) Histograms of $|A_1|$ and $|A_2|$ conditioned on $\varOmega _1$ and $\varOmega _2$ .

Figure 8. Long description

The image contains two scatter plots and two histograms. The top row features two scatter plots. The left scatter plot shows data in blue points and orange stars. The x position of each data point corresponds the value of the real part of A1, while the y position corresponds to imaginary part of A1. Blue points correspond to instants where the flow is in the one-vortex configuration while the orange stars correspond to the instants when the flow is in the two-vortex configuration. The dashed circle marks A1mean, the average distance of the blue dots to the origin. The right scatter plot also displays data in blue points and orange stars. The x position of each data point corresponds now to the value of the real part of A2, while y position to its imaginary part. The dashed circle marks A2mean, the average distance of the orange stars to the origin. The bottom row contains two histograms. The left histogram shows the counts of the values of |A_1|, while the right one counts of |A_2|. The blue bars correspond to the data from the region \Sigma_1, meaning the flow is the one-vortex configuration, while the orange bars correspond to the data from the region \Sigma_2, namely, when the flow is in two-vortex configuration. The scatter plots and histograms illustrate the partitioning of data using the KMeans algorithm, highlighting regions where variables scatter around specific values.

The switching between states and the fluctuations resemble the dynamics of a stochastic process governed by a Langevin equation whose deterministic part corresponds to a double-well potential, such as that found in (2.14). Such analogies are not a novelty; for example, Benzi (Reference Benzi2005) argues how a background of turbulence can be regarded as a noise term on an effective equation for the large scales and uses this as a model to understand spontaneous flow reversals as a result of a noise-induced transition. An analogy with a low-dimensional stochastic model is also presented in Bouchet & Simonnet (Reference Bouchet and Simonnet2009) to model transitions between large-scale jets and box-size vortices on a stochastically forced two-dimensional Navier–Stokes equation. In the next section, we use this observation to construct a model for the amplitudes $A_1$ and $A_2$ .

4. Reduced-order model

4.1. Effective amplitude equations

We have shown that only two complex degrees of freedom can describe the turbulent superstructures fairly well. We observe three main qualitative features in the fully three-dimensional flow: first, the dynamics of $A_1(t)$ and $A_2(t)$ exhibits two time scales. Second, on average, a significant fraction of the total two-dimensional energy is contained within the large-scale modes. Third, the dynamics of $A_1(t)$ and $A_2(t)$ resembles that of a bistable system. The bistable structure we characterised at a low Reynolds number remains in this regime. Considering the translation symmetry in the $x$ -direction of the original Navier–Stokes equation (2.1a ), we conclude that the evolution equations for $A_1$ and $A_2$ must be phase invariant. Therefore, if they were the only relevant degrees of freedom, then their equations must be of the form

(4.1a) \begin{align} \frac {\mathrm{d} A_1}{\mathrm{d}t} &= \sum _{j,k=0}^{\infty } C_{j,k} |A_1|^{2j} |A_2|^{2k} A_1 ,\end{align}

(4.1b) \begin{align} \frac {\mathrm{d} A_2}{\mathrm{d}t} &= \sum _{j, k =0}^{\infty } D_{j,k} |A_1|^{2j} |A_2|^{2k} A_2. \end{align}

Even allowing $C_{j,k}$ and $D_{j,k}$ to be complex, a system of equations like this cannot exhibit a complex dynamics like the one we observe when measuring these amplitudes on the numerical simulations. By separating these equations into modulus and phase, we find that the moduli evolve independently from the phases, rendering the system effectively a two-component system. Such a system can, at most, allow periodic orbits in the $(|A_1|, |A_2|)$ space. One way of inducing a complex permanent dynamics is to couple the evolution of the amplitudes with the small scales of the three-dimensional flow. Instead of introducing these degrees of freedom as additional amplitudes, we model the effect of the small scales in the form of explicit time-dependent functions in the reduced model. The following Langevin equations give the simplest of these models that can exhibit bistability:

(4.2a) \begin{align} \frac {\mathrm{d}A_1}{\mathrm{d}t} &= -\eta ^\prime _1 \frac {\partial F}{\partial A_{1}^*} + \sigma _1^\prime (t) ,\end{align}

(4.2b) \begin{align} \frac {\mathrm{d}A_2}{\mathrm{d}t} &= -\eta ^\prime _2\frac {\partial F}{\partial A_{2}^*} + \sigma ^\prime _2(t) ,\end{align}

where $F$ is a functional analogous to $\mathcal{F}$ in § 2

(4.3) \begin{equation} F = - \delta _1^\prime |A_1|^2/\eta _1^\prime - \delta _2^\prime |A_2|^2/\eta _2^\prime + \gamma ^\prime _1 |A_1|^4/2\eta _1^\prime + \gamma ^\prime _2 |A_2|^4/2\eta _2^\prime + |A_1|^2 |A_2|^2/2 \, . \end{equation}

We assume $\sigma _1^\prime (t)$ and $\sigma _2^\prime (t)$ in (4.2) to be stochastic processes. When the parameters of the deterministic part of this model are all positive, the system has two two-dimensional attractors in the $(A_1, A_2)$ space: the set of complex pairs $(|A_1^{\textit{eq}}|\mathrm{e}^{\mathrm{i}\varphi _1}, 0)$ and the set of pairs $(0, |A_2^{\textit{eq}}|\mathrm{e}^{\mathrm{i}\varphi _2})$ , with $|A_1^{\textit{eq}}| = \sqrt {\delta _1^\prime /\gamma _1^\prime }$ , $|A_2^{\textit{eq}}| =\sqrt {\delta _2^\prime /\gamma _2^\prime }$ and $\varphi _1, \varphi _2 \in [0, 2\pi )$ . Regarding the stochastic processes in the model, for simplicity we restrict our model to just additive Gaussian noise

(4.4a) \begin{align} &\sigma ^\prime _j(t) = \sigma _j \left ( X_j(t) +\mathrm{i}Y_j(t)\right ), \quad j = 1,2 ,\end{align}

(4.4b) \begin{align} &\langle X_j(t) X_k(t^\prime ) \rangle = \langle Y_j(t) Y_k(t^\prime ) \rangle = \delta _{jk} \delta (t - t^\prime ) ,\end{align}

(4.4c) \begin{align} &\langle X_j(t) Y_k(t^\prime ) \rangle = 0 , \end{align}

where $\sigma _1$ and $\sigma _2$ are positive real quantities. This system of equations minimises the functional $F$ . At the same time, the rapidly varying fluctuations may induce noise-induced transitions and, therefore, account for the switching between the two states with different numbers of vortex pairs. The set of Langevin equations (4.2) along with the noise terms (4.4) allow us to derive a Fokker–Planck equation for the probability density function (PDF) $P_{\textit{RM}}$ of finding the system at $(A_1, A_2)$

(4.5) \begin{equation} \partial _t P_{\textit{RM}} =\frac {\partial }{\partial A_1}\left (\eta _1^\prime P_{\textit{RM}} \frac {\partial F}{\partial A_1^*} + \sigma _1^2 \frac {\partial P_{\textit{RM}}}{\partial A_1^*} \right ) + \frac {\partial }{\partial A_2}\left ( \eta _2^\prime P_{\textit{RM}} \frac {\partial F}{\partial A_2^*} + \sigma _2^2 \frac {\partial P_{\textit{RM}}}{\partial A_2^*} \right ) + \mathrm{c.c.} \end{equation}

This PDF allows us to compare the single-time statistics of the amplitudes measured from the three-dimensional Kolmogorov flow with the ones of the model (4.2).

We can rationalise how we can obtain such a model as (4.2) by studying the evolution equation for the averaged velocity. After taking the spatial average of the Navier–Stokes equation (2.1a )–(2.1b ) we obtain

(4.6a) \begin{align} \partial _t\tilde {\boldsymbol u} + \boldsymbol{\tilde {u} \boldsymbol{\cdot }\boldsymbol{\nabla }\tilde {u}} &= -\boldsymbol{\nabla }\tilde {p} + \frac {1}{\textit{Re}} {\nabla} ^2 \boldsymbol{\tilde {u}} - \mu \boldsymbol{\bar {u}} + \sin {2\pi y} \boldsymbol e_x - \boldsymbol{\nabla \boldsymbol{\cdot }} \underbrace {\left (\widetilde {\boldsymbol u \boldsymbol u} - \tilde {\boldsymbol u} \tilde {\boldsymbol u}\right )}_{ \unicode{x1D64F}} , \end{align}

(4.6b) \begin{align} \partial _x \tilde {u}_x + \partial _y \tilde {u}_y &= 0 . \end{align}

Here, $ \unicode{x1D64F}$ defines the symmetric stress tensor due to the unresolved fluid motion. This tensor depends on the statistical moments of the fluctuating velocity field $\boldsymbol u^\prime = \boldsymbol u - \tilde {\boldsymbol u}$ and, in general, cannot be expressed as a closed function of $\tilde {\boldsymbol u}$ . In the turbulent regime, energy is transported from the large scales, i.e. the scales of the turbulent superstructures, into the smaller scales where the energy is dissipated. The stress tensor $ \unicode{x1D64F}$ is responsible for this process in (4.6a ) and couples the modes of the superstructures to all other small-scale modes. The small scales are chaotic and fast evolving. For this reason, we hypothesise that (4.6a ) can be effectively regarded as a stochastic two-dimensional Kolmogorov flow.

4.2. Optimising the parameters of the reduced-order model

To complete our model, we must determine the values of the parameters $\boldsymbol \varPi = (\delta '_1, \delta '_2, \gamma '_1, \gamma '_2, \eta '_1, \eta '_2, \sigma _1, \sigma _2)$ . We present the detailed estimation procedure in as a separate subsection in Appendix D; a summary of our approach goes as follows:

1. (i) We employ a kernel density estimation from the DNS data (Rosenblatt Reference Rosenblatt1956; Parzen Reference Parzen1962) to obtain an estimation of $P_{{3D}}(A_1, A_2)$ , the PDF of finding the system at a given point of the $(A_1, A_2)$ space.

2. (ii) We approximate the stationary solution of the Fokker–Planck equation (4.5), $P_{\textit{RM}}(A_1, A_2)$ , near the equilibrium points of (4.2). These approximations are expressed as functions of the parameters $\boldsymbol \varPi$ .

3. (iii) We characterise the time scale at which the amplitudes evolve by computing the correlation times when the flow is in each state from DNS. We match these times with the ones predicted by the reduced model giving us additional relations between the parameters $\boldsymbol \varPi$ .

4. (iv) Finally, we use $P_{{3D}}(A_1, A_2)$ and the correlation times to fit the parameters $\boldsymbol \varPi$ .

Figure 9.

(a) The estimated PDF of the amplitude moduli of the three-dimensional Kolmogorov flow obtained using a kernel density estimation at $\mu = 0.89$ and ${\textit{Re}}= 2851$ and histogram of the reduced model (4.2) obtained with a simulation using the optimised parameters. Contour lines in both plots are located at the same values. (b) Time series of $|A_1|$ (blue line) and $|A_2|$ (orange line) from the three-dimensional data and a simulation of the reduced model.

Figure 9. Long description

The image contains two contour plots and two line graphs. The top row features two contour plots side by side. The left contour plot, labeled ‘PDF from DNS,’ shows the estimated PDF of the amplitude moduli of the three-dimensional Kolmogorov flow. The right contour plot, labeled ‘PDF from reduced model,’ displays a histogram of the reduced model obtained with a simulation using optimized parameters. Both plots use the same contour lines, indicating identical values. The bottom row presents two line graphs side by side. The left line graph shows the time series of the three-dimensional data, with a blue line representing one variable and an orange line representing another. The right line graph depicts the time series from a simulation of the reduced model, again with a blue line for the modulus of A1 and an orange line for the modulus of A2. The x-axis for both line graphs is labeled 't/T,' and the y-axis is labeled with the respective variables. The contour plots and line graphs together illustrate the comparison between the PDFs and time series data of the three-dimensional Kolmogorov flow and the reduced model simulation.

In figure 9 (top), we display a comparison of the PDF from the DNS data with the one of the reduced model obtained with the optimised parameters. Our constructed model can reproduce, up to a fair agreement, the single-time statistics of the three-dimensional system. A comparison of the time series of $|A_1|$ and $|A_2|$ obtained in the three-dimensional simulation and the reduced model is shown in figure 9 (b). We also observe a qualitative agreement of the time scales on which the DNS amplitudes and the reduced model evolve. More quantitatively, in figure 10 we display the following time correlation functions:

(4.7a) \begin{align} C_{1}^{3D}(t) = \frac {\langle A_1(t + \tau ), A_1^*(\tau ) \rangle _{\varOmega _2}}{\langle |A_1|^2(\tau ) \rangle _{\varOmega _2}} ,\end{align}

(4.7b) \begin{align} C_{2}^{3D}(t) = \frac {\langle A_2(t + \tau ), A_2^*(\tau ) \rangle _{\varOmega _1}}{\langle |A_2|^2(\tau ) \rangle _{\varOmega _1}} \text{,} \end{align}

where $\langle \boldsymbol{\cdot }\rangle _{\varOmega _2}$ and $\langle \boldsymbol{\cdot }\rangle _{\varOmega _1}$ correspond to averages when the flow is the two-vortex-pair or one-vortex-pair configuration, respectively, computed from the DNS data, along with their respective counterparts $C_{1}^{RM}(t)$ and $C_{2}^{RM}(t)$ computed from the reduced model. We observe a good agreement between the model and DNS. Even though we have ignored nonlinear terms and assumed that the noise is simply additive, these results show that this type of modelling is a good candidate for describing the dynamics of turbulent superstructures. In table 1, we display the values of the optimised parameters.

Table 1.

Adjusted valued for the parameters of the reduced model (4.2).

Figure 10.

Real parts of the correlation functions $C_1^{3D}(t)$ and $C_2^{3D}(t)$ obtained from DNS data from the Kolmogorov flow at ${\textit{Re}}= 2851$ and $\mu = 0.89$ , superimposed with the real parts of the correlation functions $C_1^{RM}(t)$ and $C_2^{RM}(t)$ obtained from a realisation of the reduced model (4.2).

Figure 10. Long description

Two line graphs compare correlation functions from DNS data and a reduced model in Kolmogorov flow. The left graph shows the real part of the correlation functions C1RM(t) and C13D(t) with an exponential fit exp(−t/tau1). The right graph shows the real part of correlation functions C2RM(t) and C23D(t) with an exponential fit exp(−t/tau2). Both graphs plot t/T on the x-axis and correlation values on the y-axis. The blue lines represent the reduced model, the orange lines represent DNS data, and the black dashed lines represent the exponential fits. The graphs illustrate the decay of correlation functions over time.

4.3. Comparison with additional simulations at higher Reynolds numbers

Lalescu & Wilczek (Reference Lalescu and Wilczek2021) conducted an extensive numerical study of flows with Reynolds numbers up to ${\textit{Re}}= 18\,100$ . They showed that these structures and the dynamical switching between different numbers of large-scale states remain persistent. To test whether our modelling remains plausible at higher Reynolds numbers, we conduct an ensemble of simulations at $\mu = 0.89$ and $\mathrm{\textit{Re}}=7184$ . Similar to the results we already presented, we observe dynamical switching between the two large-scale states. We repeat the same analysis on the data and fitted the coefficients of the reduced-order model. We find that the coefficients at this higher $ Re$ remain fairly close to the ones we computed at a lower $ Re$ , reflecting that the superstructures evolve on a similar time scale and spend a similar time in each state at both values of $\mathrm{\textit{Re}}$ . In figure 11(a), we show the comparisons between the resulting estimated stationary PDF from the three-dimensional Kolmogorov flow data and the one obtained after fitting the model. In the lower row, we show the correlation functions obtained from the DNS and those computed from the reduced-order model. Despite the simplicity of the model, the agreement is fair. Gathering enough data to fit the coefficients at even higher $ Re$ becomes computationally expensive as the flow remains in each state for several dozen integral times. The duration of each simulation from the ensemble at ${\textit{Re}}= 7184$ was roughly two hundred integral times, which amounted to roughly $2 \times 10^5$ CPU-hours, totalling roughly $2 \times 10^6$ CPU-hours for the entire ensemble.

Figure 11.

(a) The estimated PDF of the amplitude moduli of the three-dimensional Kolmogorov flow obtained using a kernel density estimation at $\mu = 0.89$ and ${\textit{Re}}= 7184$ and histogram of the reduced model (4.2) obtained with a simulation using the optimised parameters. Contour lines in both plots are located at the same values. (b) Real parts of the correlation functions $C_1^{3D}(t)$ and $C_2^{3D}(t)$ obtained from DNS data from the Kolmogorov flow at ${\textit{Re}}= 7184$ and $\mu = 0.89$ , superimposed by the real parts of the correlation functions $C_1^{RM}(t)$ and $C_2^{RM}(t)$ obtained from a realisation of the reduced model (4.2).

Figure 11. Long description

The image contains two contour plots and two line graphs. The top contour plots compare the estimated PDF of the amplitude moduli of the three-dimensional Kolmogorov flow obtained using a kernel density estimation and a histogram of the reduced model. The contour lines in both plots are located at the same values. The bottom line graphs show correlation functions obtained from DNS data from the Kolmogorov flow, superimposed by the correlation functions obtained from a realization of the reduced model. The x-axis represents normalized time, and the y-axis represents correlation values. The orange lines represent the reduced model, the blue lines represent the three-dimensional data, and the dashed black lines represent exponential decay. The left panel shows the real parts of C1Rm(t) and C13D(t) and the right one shows the real parts of C2Rm(t) and C23D(t).