2 Bounds on the growth of enstrophy

We consider the incompressible Navier–Stokes system defined on the 3-D unit cube $\unicode[STIX]{x1D6FA}=[0,1]^{3}$ with periodic boundary conditions

(2.1a)$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}p-\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\boldsymbol{u}=0\quad \text{in }\unicode[STIX]{x1D6FA}\times (0,T], & \displaystyle\end{eqnarray}$$

(2.1b)$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0\quad \text{in }\unicode[STIX]{x1D6FA}\times [0,T], & \displaystyle\end{eqnarray}$$

(2.1c)$$\begin{eqnarray}\displaystyle & \boldsymbol{u}(0)=\boldsymbol{u}_{0}, & \displaystyle\end{eqnarray}$$

$\boldsymbol{u}=[u_{1},u_{2},u_{3}]^{\text{T}}$

$p$

$\unicode[STIX]{x1D708}>0$

$\unicode[STIX]{x1D708}=0.01$

$\unicode[STIX]{x1D735}\boldsymbol{u}$

$[\unicode[STIX]{x1D735}\boldsymbol{u}]_{ij}=\unicode[STIX]{x2202}_{j}u_{i}$

$i,j=1,2,3$

$\unicode[STIX]{x1D70C}$

$\unicode[STIX]{x1D70C}=1$

${\mathcal{K}}(\boldsymbol{u}(t))$

(2.2)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{K}}(\boldsymbol{u}(t)):=\frac{1}{2}\int _{\unicode[STIX]{x1D6FA}}|\boldsymbol{u}(t,\boldsymbol{x})|^{2}\,\text{d}\boldsymbol{x}, & \displaystyle\end{eqnarray}$$

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}{\mathcal{K}}(\boldsymbol{u}(t))}{\text{d}t}=-\unicode[STIX]{x1D708}\int _{\unicode[STIX]{x1D6FA}}|\unicode[STIX]{x1D735}\boldsymbol{u}(t,\boldsymbol{x})|^{2}\,\text{d}\boldsymbol{x}, & \displaystyle\end{eqnarray}$$

where ‘$:=$’ means ‘equal to by definition’. The enstrophy ${\mathcal{E}}(\boldsymbol{u}(t))$ and its rate of growth are given by

(2.4)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{E}}(\boldsymbol{u}(t)):=\frac{1}{2}\int _{\unicode[STIX]{x1D6FA}}|\unicode[STIX]{x1D735}\times \boldsymbol{u}(t,\boldsymbol{x})|^{2}\,\text{d}\boldsymbol{x}, & \displaystyle\end{eqnarray}$$

(2.5)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}{\mathcal{E}}(\boldsymbol{u}(t))}{\text{d}t}=-\unicode[STIX]{x1D708}\int _{\unicode[STIX]{x1D6FA}}|\unicode[STIX]{x0394}\boldsymbol{u}|^{2}\,\text{d}\boldsymbol{x}+\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x0394}\boldsymbol{u}\,\text{d}\boldsymbol{x}=:{\mathcal{R}}(\boldsymbol{u}(t)). & \displaystyle\end{eqnarray}$$

(We note that unlike energy, cf. (2.2), enstrophy is often defined without the factor of $1/2$. However, for consistency with earlier studies belonging to this research programme (Ayala & Protas Reference Ayala and Protas2011, Reference Ayala and Protas2014a,Reference Ayala and Protasb, Reference Ayala and Protas2017; Yun & Protas Reference Yun and Protas2018), we choose to retain this factor here.)

For incompressible flows with periodic boundary conditions, we also have the following identity (Doering & Gibbon Reference Doering and Gibbon1995):

(2.6)$$\begin{eqnarray}\int _{\unicode[STIX]{x1D6FA}}|\unicode[STIX]{x1D735}\times \boldsymbol{u}|^{2}\,\text{d}\boldsymbol{x}=\int _{\unicode[STIX]{x1D6FA}}|\unicode[STIX]{x1D735}\boldsymbol{u}|^{2}\,\text{d}\boldsymbol{x}.\end{eqnarray}$$

Hence, combining (2.2)–(2.6), the following system of ordinary differential equations is obtained for energy and enstrophy:

(2.7a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}{\mathcal{K}}(\boldsymbol{u}(t))}{\text{d}t}=-2\unicode[STIX]{x1D708}{\mathcal{E}}(\boldsymbol{u}(t)), & \displaystyle\end{eqnarray}$$

(2.7b)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}{\mathcal{E}}(\boldsymbol{u}(t))}{\text{d}t}={\mathcal{R}}(\boldsymbol{u}(t)). & \displaystyle\end{eqnarray}$$

A standard approach at this point is to try to bound $\text{d}{\mathcal{E}}/\text{d}t$ and using classical techniques of functional analysis, it is possible to obtain the following well-known estimate in terms of ${\mathcal{K}}$ and ${\mathcal{E}}$ (Doering Reference Doering2009):

(2.8)$$\begin{eqnarray}\frac{\text{d}{\mathcal{E}}}{\text{d}t}\leqslant -\unicode[STIX]{x1D708}\frac{{\mathcal{E}}^{2}}{{\mathcal{K}}}+\frac{c}{\unicode[STIX]{x1D708}^{3}}{\mathcal{E}}^{3},\end{eqnarray}$$

where $c$ is a known constant. A related estimate expressed entirely in terms of the enstrophy ${\mathcal{E}}$ is given by

(2.9)$$\begin{eqnarray}\frac{\text{d}{\mathcal{E}}}{\text{d}t}\leqslant \frac{27}{8\,\unicode[STIX]{x03C0}^{4}\,\unicode[STIX]{x1D708}^{3}}{\mathcal{E}}^{3}.\end{eqnarray}$$

By simply integrating the differential inequality in (2.9) with respect to time we obtain the finite-time bound

(2.10)$$\begin{eqnarray}{\mathcal{E}}(\boldsymbol{u}(t))\leqslant \frac{{\mathcal{E}}_{0}}{\sqrt{1-{\displaystyle \frac{27}{4\,\unicode[STIX]{x03C0}^{4}\,\unicode[STIX]{x1D708}^{3}}}\,{\mathcal{E}}_{0}^{2}\,t}},\end{eqnarray}$$

which clearly becomes infinite at time $t_{0}=4\,\unicode[STIX]{x03C0}^{4}\,\unicode[STIX]{x1D708}^{3}/(27\,{\mathcal{E}}_{0}^{2})$. Thus, based on estimate (2.10), it is not possible to establish the boundedness of the enstrophy ${\mathcal{E}}(\boldsymbol{u}(t))$ and, hence, also the regularity of solutions globally in time. Therefore, the question about the finite-time singularity formation can be recast in terms of whether or not there exists initial data $\boldsymbol{u}_{0}$ with enstrophy ${\mathcal{E}}_{0}<\infty$ such that in the resulting flow evolution the enstrophy ${\mathcal{E}}(\boldsymbol{u}(t))$ becomes unbounded in finite time, as allowed by estimate (2.10). A systematic search for such worst-case initial data using variational optimization methods is the main theme of this study. We add that while the analysis presented here was carried out based on the vorticity and enstrophy, an inequality analogous to (2.9) can also be obtained in terms of strain, i.e. the symmetric part of the velocity gradient $\unicode[STIX]{x1D735}\boldsymbol{u}$, resulting in a smaller value of the constant prefactor (Miller Reference Miller2019).

The question of sharpness of the instantaneous bound (2.9) was addressed in the seminal study by Lu & Doering (Reference Lu and Doering2008), see also Lu (Reference Lu2006), and further elaborated by Ayala & Protas (Reference Ayala and Protas2017), who constructed a family of divergence-free velocity fields $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ parameterized by the enstrophy ${\mathcal{E}}_{0}$ which saturate this bound. These fields were obtained by numerically solving the following variational optimization problem.

Problem 2.1. Given ${\mathcal{E}}_{0}\in \mathbb{R}_{+}$ and the objective functional ${\mathcal{R}}(\boldsymbol{u})$ from (2.5), find

$$\begin{eqnarray}\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}=\mathop{\text{argmax}}_{\boldsymbol{u}\in {\mathcal{S}}_{{\mathcal{E}}_{0}}}{\mathcal{R}}(\boldsymbol{u}),\end{eqnarray}$$

where

$$\begin{eqnarray}{\mathcal{S}}_{{\mathcal{E}}_{0}}=\{\boldsymbol{u}\in H^{2}(\unicode[STIX]{x1D6FA})\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,{\mathcal{E}}(\boldsymbol{u})={\mathcal{E}}_{0}\}\end{eqnarray}$$

for the enstrophy ${\mathcal{E}}_{0}$ spanning a broad range of values. Since problem 2.1 is in general non-convex, ‘$\operatorname{argmax}$’ represents a local maximizer, which might also be global. The symbol $H^{2}(\unicode[STIX]{x1D6FA})$ denotes the Sobolev space of functions with square-integrable second derivatives endowed with the inner product (Adams & Fournier Reference Adams and Fournier2005)

(2.11)$$\begin{eqnarray}\forall \;\boldsymbol{z}_{1},\boldsymbol{z}_{2}\in H^{2}(\unicode[STIX]{x1D6FA})\quad \langle \boldsymbol{z}_{1},\boldsymbol{z}_{2}\rangle _{H^{2}(\unicode[STIX]{x1D6FA})}=\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{z}_{1}\boldsymbol{\cdot }\boldsymbol{z}_{2}+\unicode[STIX]{x1D735}\boldsymbol{z}_{1}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{z}_{2}+\unicode[STIX]{x0394}\boldsymbol{z}_{1}\boldsymbol{\cdot }\unicode[STIX]{x0394}\boldsymbol{z}_{2}\,\text{d}\boldsymbol{x},\end{eqnarray}$$

which, for simplicity, was chosen in a somewhat nonstandard form as it does not involve mixed derivatives in the last term. For sufficiently large values of ${\mathcal{E}}_{0}$, the thus obtained instantaneously optimal fields $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ have the form of a pair of colliding vortex rings, cf. figure 1(a). The corresponding rate of growth of enstrophy $\text{d}{\mathcal{E}}/\text{d}t$ was found to be proportional to ${\mathcal{E}}^{3}$, cf. figure 1(b), demonstrating that estimate (2.9) is sharp up to a numerical prefactor. However, the sharpness of the instantaneous estimate alone does not allow us to make a conclusion about the possibility of singularity formation, because for this situation to occur, a sufficiently large enstrophy growth rate would need to be sustained over a finite time window $[0,t_{0})$. At the same time, it was shown by Ayala & Protas (Reference Ayala and Protas2017) that when the extreme vortex states $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ are used as the initial data $\boldsymbol{u}_{0}$ in the Navier–Stokes system (2.1), the initially maximal rate of growth of enstrophy is immediately depleted producing only a very modest increase $\max _{t\geqslant 0}{\mathcal{E}}(\boldsymbol{u}(t))-{\mathcal{E}}_{0}$ during subsequent evolution. We note that in the limit ${\mathcal{E}}_{0}\rightarrow \infty$ the extreme vortex states $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ are given by pairs of axisymmetric vortex rings without swirl and such flows are known to be globally well posed in the classical sense (Feng & Šverák Reference Feng and Šverák2015). The fact that a singularity cannot form in such axisymmetric configurations with no vorticity on the axis can also be deduced from the celebrated Caffarelli–Kohn–Nirenberg theorem (Caffarelli, Kohn & Nirenberg Reference Caffarelli, Kohn and Nirenberg1982).

Figure 1.

(a) Extreme vortex state $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ obtained by Ayala & Protas (Reference Ayala and Protas2017), see also Lu & Doering (Reference Lu and Doering2008), as a solution of the instantaneous optimization problem 2.1 for ${\mathcal{E}}_{0}=20$; shades of red correspond to the magnitude of the vorticity $|(\unicode[STIX]{x1D735}\times \widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}})(\boldsymbol{x})|$ (see the colour bar). (b) The maximum rate of growth of enstrophy ${\mathcal{R}}(\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}})$ as a function of the enstrophy ${\mathcal{E}}_{0}$ (the black solid symbol corresponds to the value of ${\mathcal{E}}_{0}$ characterizing the extreme vortex state shown in (a)).

In addition to sharpness under worst-case conditions, another question pertaining to inequality (2.9) is whether the upper bound on its right-hand side can also be realized under generic conditions in turbulent flows. This problem was studied by Schumacher, Eckhardt & Doering (Reference Schumacher, Eckhardt and Doering2010) who demonstrated that in turbulent flows the rate of change of ensemble-averaged squared vorticity grows at most as $(\text{d}/\text{d}t)\langle \unicode[STIX]{x1D74E}^{2}\rangle \sim \langle \unicode[STIX]{x1D74E}^{2}\rangle ^{3/2}$, where $\langle \cdot \rangle$ denotes ensemble averaging. This observation is consistent with the statistical theory of turbulence, more specifically, the Kármán–Howarth equation (Davidson Reference Davidson2004).

The key conclusion from the results recalled above is that if a significant, let alone unbounded, growth of enstrophy is to be achieved in finite time, it must be associated with initial data $\boldsymbol{u}_{0}$ other than the extreme vortex states $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ saturating the upper bound in estimate (2.9) on the instantaneous rate of growth of enstrophy, cf. figure 1. More specifically, assuming the instantaneous rate of growth of enstrophy in the form $\text{d}{\mathcal{E}}/\text{d}t=C{\mathcal{E}}^{\unicode[STIX]{x1D6FC}}$ for some prefactor $C>0$, any exponent $\unicode[STIX]{x1D6FC}>2$ will lead to blow-up of ${\mathcal{E}}(\boldsymbol{u}(t))$ at some finite time $t_{0}=t_{0}(\unicode[STIX]{x1D6FC})$ if this rate of growth is sustained over the interval $[0,t_{0})$. The fact that there is no blow-up when $1<\unicode[STIX]{x1D6FC}\leqslant 2$ follows from the observation that one factor of ${\mathcal{E}}$ in (2.9) can be bounded in terms of the initial energy ${\mathcal{K}}_{0}:={\mathcal{K}}(\boldsymbol{u}_{0})$ using (2.7a) as follows:

(2.12)$$\begin{eqnarray}\int _{0}^{t}{\mathcal{E}}(\boldsymbol{u}(s))\,\text{d}s=\frac{1}{2\unicode[STIX]{x1D708}}[{\mathcal{K}}_{0}-{\mathcal{K}}(\boldsymbol{u}(t))]\leqslant \frac{1}{2\unicode[STIX]{x1D708}}{\mathcal{K}}_{0}\end{eqnarray}$$

which, upon employing Grönwall’s lemma, yields the bound

(2.13)$$\begin{eqnarray}\max _{0\leqslant t\leqslant T}{\mathcal{E}}(\boldsymbol{u}(t))\leqslant {\mathcal{E}}_{0}\exp \left[\int _{0}^{T}{\mathcal{E}}(\boldsymbol{u}(s))\,\text{d}s\right]\leqslant {\mathcal{E}}_{0}\exp \left[\frac{1}{2\unicode[STIX]{x1D708}}{\mathcal{K}}_{0}\right]\end{eqnarray}$$

valid for $0\leqslant \unicode[STIX]{x1D6FC}\leqslant 2$. Evidently, as the rate of growth of enstrophy slows down and $\unicode[STIX]{x1D6FC}\rightarrow 2^{+}$, for blow-up to occur a certain minimum growth rate must be sustained over windows of time with increasing length, i.e. $t_{0}\rightarrow \infty$. To assess the feasibility of such a scenario, in the next section we introduce an optimization approach which will allow us to systematically search for worst-case flow evolutions producing the maximum possible growth of enstrophy in finite time.

3 Maximization problem

In order to probe the upper bound in estimate (2.10) for realizability, our objective in this section is to construct initial data $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ for the Navier–Stokes system (2.1) with the prescribed enstrophy ${\mathcal{E}}_{0}>0$, such that at the given time $T>0$ the corresponding flow evolution will produce the maximum possible value of enstrophy ${\mathcal{E}}(\boldsymbol{u}(T))$ under the assumption that the Navier–Stokes system admits smooth (classical) solutions on the time interval $(0,T]$. Defining the objective function ${\mathcal{E}}_{T}:H^{1}(\unicode[STIX]{x1D6FA})\rightarrow \mathbb{R}_{+}$, where

(3.1)$$\begin{eqnarray}{\mathcal{E}}_{T}(\boldsymbol{u}_{0}):={\mathcal{E}}(\boldsymbol{u}(T))=\frac{1}{2}\int _{\unicode[STIX]{x1D6FA}}|\unicode[STIX]{x1D735}\times \boldsymbol{u}(T,\boldsymbol{x})|^{2}\,\text{d}\boldsymbol{x},\end{eqnarray}$$

we thus arrive at the following.

Problem 3.1. Given ${\mathcal{E}}_{0},T\in \mathbb{R}_{+}$ and the objective functional ${\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ from (3.1), find

$$\begin{eqnarray}\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}=\mathop{\text{argmax}}_{\boldsymbol{u}_{0}\in {\mathcal{Q}}_{{\mathcal{E}}_{0}}}\,{\mathcal{E}}_{T}(\boldsymbol{u}_{0}),\end{eqnarray}$$

where

$$\begin{eqnarray}{\mathcal{Q}}_{{\mathcal{E}}_{0}}=\{\boldsymbol{u}\in H^{1}(\unicode[STIX]{x1D6FA})\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,{\mathcal{E}}(\boldsymbol{u})={\mathcal{E}}_{0}\},\end{eqnarray}$$

where $H^{1}(\unicode[STIX]{x1D6FA})$ denotes the Sobolev space of functions with square-integrable derivatives endowed with the inner product (Adams & Fournier Reference Adams and Fournier2005)

(3.2)$$\begin{eqnarray}\forall \;\boldsymbol{z}_{1},\boldsymbol{z}_{2}\in H^{1}(\unicode[STIX]{x1D6FA})\quad \langle \boldsymbol{z}_{1},\boldsymbol{z}_{2}\rangle _{H^{1}(\unicode[STIX]{x1D6FA})}=\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{z}_{1}\boldsymbol{\cdot }\boldsymbol{z}_{2}+\ell _{1}^{2}\unicode[STIX]{x1D735}\boldsymbol{z}_{1}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{z}_{2}\,\text{d}\boldsymbol{x},\end{eqnarray}$$

where $\ell _{1}\in \mathbb{R}_{+}$ is a parameter with the meaning of a characteristic length scale (the reasons for introducing this parameter in the definition of the inner product will become clear below). The inner product in the space $L^{2}(\unicode[STIX]{x1D6FA})$ is obtained from (3.2) by setting $\ell _{1}=0$. The maximizers in problem 3.1 are constrained to belong to the manifold ${\mathcal{Q}}_{{\mathcal{E}}_{0}}$ which represents the intersection of the subspace of divergence-free vector fields $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0)$ and a nonlinear manifold defined by the enstrophy constraint $({\mathcal{E}}(\boldsymbol{u})={\mathcal{E}}_{0})$ in the Sobolev space $H^{1}(\unicode[STIX]{x1D6FA})$, where the smoothness requirement is necessary to ensure that the initial enstrophy ${\mathcal{E}}(\boldsymbol{u}_{0})$ is well defined. Without loss of generality, we assume that $\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{u}_{0}\,\text{d}\boldsymbol{x}=\mathbf{0}$, a property which is also invariant during the flow evolution. The fact that problem 3.1 admits solutions is a consequence of the assumption that with the given parameters ${\mathcal{E}}_{0}$ and $T$ solutions of the Navier–Stokes system (2.1) are smooth on the time interval $[0,T]$.

The key insight we seek to deduce is how the maximum growth of enstrophy $\max _{T>0}\,{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ obtained for time intervals $(0,T]$ with different lengths $0<T<\infty$ scales with the initial enstrophy ${\mathcal{E}}_{0}$. In order to evaluate this quantity for a given value of the initial enstrophy ${\mathcal{E}}_{0}$, we thus need to solve problem 3.1 for different values of $T\in (0,T_{max}]$, where $T_{max}<\infty$ is the maximum considered length of the time interval, and with fixed ${\mathcal{E}}_{0}$, so that the maximum with respect to $T$ can be evaluated. We note that this approach is justified by the existence of bounds (expressed in terms of norms of the initial data $\boldsymbol{u}_{0}$ and viscosity $\unicode[STIX]{x1D708}$) on the largest time $T_{max}$ when a singularity might occur (Ohkitani Reference Ohkitani2016). Assuming that the optimal initial data $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained as a solution of problem 3.1 has (at least piecewise) continuous dependence on the parameter $T$, we will refer to the mapping $T\longmapsto \widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ with ${\mathcal{E}}_{0}$ fixed as a ‘maximizing branch’.

In order to solve problem 3.1 for given values of ${\mathcal{E}}_{0}$ and $T$ we adopt an ‘optimize-then-discretize’ approach (Gunzburger Reference Gunzburger2003) in which a gradient method is first formulated in the infinite-dimensional (continuous) setting and only then the resulting equations and expressions are discretized for the purpose of numerical solution. An essentially identical approach was used by Ayala & Protas (Reference Ayala and Protas2011) to solve an optimization problem analogous to problem 3.1, but formulated for the 1-D Burgers equation. The maximizer $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ can be found as $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}=\lim _{n\rightarrow \infty }\boldsymbol{u}_{0;{\mathcal{E}}_{0},T}^{(n)}$ using the following iterative procedure representing a discretization of a gradient flow projected on ${\mathcal{Q}}_{{\mathcal{E}}_{0}}$:

(3.3)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{u}_{0;{\mathcal{E}}_{0},T}^{(n+1)}=\mathbb{P}_{{\mathcal{Q}}_{{\mathcal{E}}_{0}}}(\boldsymbol{u}_{0;{\mathcal{E}}_{0},T}^{(n)}+\unicode[STIX]{x1D70F}_{n}\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0;{\mathcal{E}}_{0},T}^{(n)})),\\ \boldsymbol{u}_{0;{\mathcal{E}}_{0},T}^{(1)}=\boldsymbol{u}^{0}.\end{array}\right\}\end{eqnarray}$$

Here $\boldsymbol{u}_{0;{\mathcal{E}}_{0},T}^{(n)}$ is an approximation of the maximizer obtained at the $n$th iteration, $\mathbb{P}_{{\mathcal{Q}}_{{\mathcal{E}}_{0}}}:H^{1}(\unicode[STIX]{x1D6FA})\rightarrow {\mathcal{Q}}_{{\mathcal{E}}_{0}}$ is the projection operator, $\boldsymbol{u}^{0}$ is the initial guess and $\unicode[STIX]{x1D70F}_{n}$ is the length of the step. A key element of the iterative procedure (3.3) is the evaluation of the gradient $\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ of the objective functional ${\mathcal{E}}_{T}(\boldsymbol{u}_{0})$, cf. (3.1), representing its (infinite-dimensional) sensitivity to perturbations of the initial data $\boldsymbol{u}_{0}$ in the governing system (2.1). We emphasize that it is essential for the gradient to possess the required regularity, namely, $\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})\in H^{1}(\unicode[STIX]{x1D6FA})$.

The first step to determine the gradient $\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ is to consider the Gâteaux (directional) differential ${\mathcal{E}}_{T}^{\prime }(\boldsymbol{u}_{0};\cdot ):H^{1}(\unicode[STIX]{x1D6FA})\rightarrow \mathbb{R}$ of the objective functional ${\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ defined as ${\mathcal{E}}_{T}^{\prime }(\boldsymbol{u}_{0};\boldsymbol{u}_{0}^{\prime })=\lim _{\unicode[STIX]{x1D716}\rightarrow 0}\unicode[STIX]{x1D716}^{-1}[{\mathcal{E}}_{T}(\boldsymbol{u}_{0}+\unicode[STIX]{x1D716}\boldsymbol{u}_{0}^{\prime })-{\mathcal{E}}_{T}(\boldsymbol{u}_{0})]$ for some arbitrary perturbation $\boldsymbol{u}_{0}^{\prime }\in H^{1}(\unicode[STIX]{x1D6FA})$. The gradient $\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ can then be extracted from the Gâteaux differential ${\mathcal{E}}_{T}^{\prime }(\boldsymbol{u}_{0};\boldsymbol{u}_{0}^{\prime })$ recognizing that, when viewed as a function of its second argument, this differential is a bounded linear functional on the space $H^{1}(\unicode[STIX]{x1D6FA})$ and we can therefore invoke the Riesz representation theorem (Luenberger Reference Luenberger1969)

(3.4)$$\begin{eqnarray}{\mathcal{E}}_{T}^{\prime }(\boldsymbol{u}_{0};\boldsymbol{u}_{0}^{\prime })=\langle \unicode[STIX]{x1D6FB}^{L_{2}}{\mathcal{E}}_{T}(\boldsymbol{u}_{0}),\boldsymbol{u}_{0}^{\prime }\rangle _{L^{2}(\unicode[STIX]{x1D6FA})}=\langle \unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0}),\boldsymbol{u}_{0}^{\prime }\rangle _{H^{1}(\unicode[STIX]{x1D6FA})},\end{eqnarray}$$

where the gradient $\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ is the Riesz representer in the function space $H^{1}(\unicode[STIX]{x1D6FA})$. In (3.4) we also formally defined the gradient $\unicode[STIX]{x1D6FB}^{L_{2}}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ computed with respect to the $L^{2}$ topology as it will be useful in subsequent computations. Given the definition of the objective functional in (3.1), its Gâteaux differential can be expressed as

(3.5)$$\begin{eqnarray}{\mathcal{E}}_{T}^{\prime }(\boldsymbol{u}_{0};\boldsymbol{u}_{0}^{\prime })=\int _{\unicode[STIX]{x1D6FA}}(\unicode[STIX]{x1D735}\times \boldsymbol{u}(T,\boldsymbol{x}))\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\times \boldsymbol{u}^{\prime }(T,\boldsymbol{x}))\,\text{d}\boldsymbol{x}=\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D71F}\boldsymbol{u}(T,\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{u}^{\prime }(T,\boldsymbol{x})\,\text{d}\boldsymbol{x},\end{eqnarray}$$

where the last equality follows from integration by parts and the vector identity $\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D735}\times \boldsymbol{z})=\unicode[STIX]{x1D735}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{z})-\unicode[STIX]{x1D71F}\boldsymbol{z}$, whereas the perturbation field $\boldsymbol{u}^{\prime }=\boldsymbol{u}^{\prime }(t,\boldsymbol{x})$ is a solution of the Navier–Stokes system linearized around the trajectory corresponding to the initial data $\boldsymbol{u}_{0}$ (Gunzburger Reference Gunzburger2003), i.e.

(3.6a)$$\begin{eqnarray}\displaystyle & {\mathcal{L}}\left[\begin{array}{@{}c@{}}\boldsymbol{u}^{\prime }\\ p^{\prime }\end{array}\right]:=\left[\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{\prime }+\boldsymbol{u}^{\prime }\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime }+\unicode[STIX]{x1D735}p^{\prime }-\unicode[STIX]{x1D708}\unicode[STIX]{x1D71F}\boldsymbol{u}^{\prime }\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\prime }\end{array}\right]=\left[\begin{array}{@{}c@{}}\mathbf{0}\\ 0\end{array}\right], & \displaystyle\end{eqnarray}$$

(3.6b)$$\begin{eqnarray}\displaystyle & \boldsymbol{u}^{\prime }(0)=\boldsymbol{u}_{0}^{\prime }, & \displaystyle\end{eqnarray}$$

$p^{\prime }$

We note that expression (3.5) for the Gâteaux differential is not consistent with the Riesz form (3.4), because the perturbation $\boldsymbol{u}_{0}^{\prime }$ of the initial data does not appear in it explicitly as a factor, but is instead hidden as the initial condition in the linearized problem, cf. (3.6b). In order to transform (3.5) to the Riesz form, we introduce the adjoint states $\boldsymbol{u}^{\ast }:[0,T]\times \unicode[STIX]{x1D6FA}\rightarrow \mathbb{R}^{3}$ and $p^{\ast }:[0,T]\times \unicode[STIX]{x1D6FA}\rightarrow \mathbb{R}$, and the following duality-pairing relation:

(3.7)$$\begin{eqnarray}\displaystyle \left({\mathcal{L}}\left[\begin{array}{@{}c@{}}\boldsymbol{u}^{\prime }\\ p^{\prime }\end{array}\right],\left[\begin{array}{@{}c@{}}\boldsymbol{u}^{\ast }\\ p^{\ast }\end{array}\right]\right) & := & \displaystyle \int _{0}^{T}\int _{\unicode[STIX]{x1D6FA}}{\mathcal{L}}\left[\begin{array}{@{}c@{}}\boldsymbol{u}^{\prime }\\ p^{\prime }\end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}c@{}}\boldsymbol{u}^{\ast }\\ p^{\ast }\end{array}\right]\,\text{d}\boldsymbol{x}\,\text{d}t=\left(\left[\begin{array}{@{}c@{}}\boldsymbol{u}^{\prime }\\ p^{\prime }\end{array}\right],{\mathcal{L}}^{\ast }\left[\begin{array}{@{}c@{}}\boldsymbol{u}^{\ast }\\ p^{\ast }\end{array}\right]\right)\nonumber\\ \displaystyle & & \displaystyle +\,\underbrace{\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{u}^{\prime }(T,\boldsymbol{ x})\boldsymbol{\cdot }\boldsymbol{u}^{\ast }(T,\boldsymbol{ x})\,\text{d}\boldsymbol{ x}}_{{\mathcal{E}}_{T}^{\prime }(\boldsymbol{u}_{0};\boldsymbol{ u}_{0}^{\prime })}-\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{ u}^{\prime }(0,\boldsymbol{ x})\boldsymbol{\cdot }\boldsymbol{u}^{\ast }(0,\boldsymbol{ x})\,\text{d}\boldsymbol{ x}=0.\qquad\end{eqnarray}$$

Here ‘$\boldsymbol{\cdot }$’ in the first integrand expression denotes the Euclidean dot product evaluated at $(t,\boldsymbol{x})$. Performing integration by parts with respect to both space and time then allows us to define the adjoint system as

(3.8a)$$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{L}}^{\ast }\left[\begin{array}{@{}c@{}}\boldsymbol{u}^{\ast }\\ p^{\ast }\end{array}\right]:=\left[\begin{array}{@{}c@{}}-\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{\ast }-[\unicode[STIX]{x1D735}\boldsymbol{u}^{\ast }+(\unicode[STIX]{x1D735}\boldsymbol{u}^{\ast })^{\text{T}}]\boldsymbol{u}-\unicode[STIX]{x1D735}p^{\ast }-\unicode[STIX]{x1D708}\unicode[STIX]{x1D71F}\boldsymbol{u}^{\ast }\\ -\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}^{\ast }\end{array}\right]=\left[\begin{array}{@{}c@{}}\mathbf{0}\\ 0\end{array}\right], & \displaystyle\end{eqnarray}$$

(3.8b)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{\ast }(T)=\unicode[STIX]{x1D71F}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

$\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{u}^{\prime }(T,\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{u}^{\ast }(T,\boldsymbol{x})\,\text{d}\boldsymbol{x}$

(3.9)$$\begin{eqnarray}{\mathcal{E}}_{T}^{\prime }(\boldsymbol{u}_{0};\boldsymbol{u}_{0}^{\prime })=\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{u}_{0}^{\prime }(\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{u}^{\ast }(0,\boldsymbol{x})\,\text{d}\boldsymbol{x}.\end{eqnarray}$$

Applying the first equality in Riesz relations (3.4) to (3.9), we obtain the $L^{2}$ gradient as

(3.10)$$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{L^{2}}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})=\boldsymbol{u}^{\ast }(0).\end{eqnarray}$$

In order to obtain the required Sobolev $H^{1}$ gradient $\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$, we identify the Gâteaux differential in (3.9) with the $H^{1}$ inner product, cf. (3.2). Then, recognizing that the perturbations $\boldsymbol{u}_{0}^{\prime }$ are arbitrary, we obtain the elliptic boundary-value problem (Protas, Bewley & Hagen Reference Protas, Bewley and Hagen2004)

(3.11)$$\begin{eqnarray}[\operatorname{Id}-\ell _{1}^{2}\unicode[STIX]{x1D71F}]\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})=\unicode[STIX]{x1D6FB}^{L_{2}}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})\quad \text{in }\unicode[STIX]{x1D6FA}\end{eqnarray}$$

subject to the periodic boundary conditions, which must be solved to determine $\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$. The gradient fields $\unicode[STIX]{x1D6FB}^{L_{2}}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ and $\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ can be interpreted as infinite-dimensional sensitivities of the objective functional ${\mathcal{E}}_{T}(\boldsymbol{u}_{0})$, cf. (3.1), with respect to perturbations of the initial data $\boldsymbol{u}_{0}$. While these two gradients may point towards the same local maximizer, they represent distinct ‘directions’, since they are defined with respect to different norms ($L^{2}$ versus $H^{1}$). As shown by Protas et al. (Reference Protas, Bewley and Hagen2004), extraction of gradients in spaces of smoother functions such as $H^{1}(\unicode[STIX]{x1D6FA})$ can be interpreted as low-pass filtering of the $L^{2}$ gradients with the parameter $\ell _{1}$ acting as the cut-off length scale. Although Sobolev gradients obtained with different $0<\ell _{1}<\infty$ are equivalent, in the precise sense of norm equivalence (Berger Reference Berger1977), the value of $\ell _{1}$ tends to have a significant effect on the rate of convergence of gradient iterations (3.3) (Protas et al. Reference Protas, Bewley and Hagen2004) and the choice of its numerical value will be discussed in § 4. We emphasize that, while the $H^{1}$ gradient is used exclusively in the actual computations, cf. (3.3), the $L^{2}$ gradient is computed first merely as an intermediate step.

Evaluation of the $L^{2}$ gradient at a given iteration via (3.10) requires solution of the Navier–Stokes system (2.1) followed by solution of the adjoint system (3.8). We note that this system is a linear problem with coefficients and the terminal condition determined by the solution of the Navier–Stokes system obtained earlier during the iteration. The adjoint system (3.8) is a terminal value problem, implying that it must be integrated backwards in time from $t=T$ to $t=0$ (since the term with the time derivative has a negative sign, this problem is well posed). Once the $L^{2}$ gradient is determined using (3.10), the corresponding Sobolev $H^{1}$ gradient can be obtained by solving problem (3.11). We add that the thus computed gradient satisfies the divergence-free condition by construction, i.e. $\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0}))=0$.

As regards the fixed-enstrophy constraint ${\mathcal{E}}(\boldsymbol{u}_{0})={\mathcal{E}}_{0}$, this property is enforced in iterations (3.3) using the projection operator $\mathbb{P}_{{\mathcal{Q}}_{{\mathcal{E}}_{0}}}$ defined as the normalization

(3.12)$$\begin{eqnarray}\mathbb{P}_{{\mathcal{Q}}_{{\mathcal{E}}_{0}}}(\boldsymbol{u}_{0}):=\sqrt{\frac{{\mathcal{E}}_{0}}{{\mathcal{E}}(\boldsymbol{u}_{0})}}\,\boldsymbol{u}_{0},\end{eqnarray}$$

which clearly preserves the divergence-free property of the argument.

The step size $\unicode[STIX]{x1D70F}_{n}$ in algorithm (3.3) is computed as

(3.13)$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{n}=\mathop{\text{argmax}}_{\unicode[STIX]{x1D70F}>0}\{{\mathcal{E}}_{T}[\mathbb{P}_{{\mathcal{Q}}_{{\mathcal{E}}_{0}}}(\boldsymbol{u}_{0;{\mathcal{E}}_{0},T}^{(n)}+\unicode[STIX]{x1D70F}\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0;{\mathcal{E}}_{0},T}^{(n)}))]\},\end{eqnarray}$$

which is done using a suitable derivative-free line-search approach, such as a variant of Brent’s algorithm (Press et al. Reference Press, Flannery, Teukolsky and Vetterling1986; Nocedal & Wright Reference Nocedal and Wright1999). Equation (3.13) can be interpreted as a modification of the standard line-search problem where optimization is performed following an arc (a geodesic in the limit of infinitesimal step sizes) lying on the constraint manifold ${\mathcal{S}}_{{\mathcal{E}}_{0}}$, rather than along a straight line. This approach was already successfully employed to solve similar problems in Ayala & Protas (Reference Ayala and Protas2011, Reference Ayala and Protas2014a, Reference Ayala and Protas2017).

Maximizing branches are computed using a continuation approach by fixing one parameter, e.g. ${\mathcal{E}}_{0}$, and then solving problem 3.1 with procedure (3.3) repeatedly for increasing values of $T$. In this process the maximizer $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained for some ${\mathcal{E}}_{0}$ and $T$ is employed as the initial guess $\boldsymbol{u}^{0}$ in (3.3) to compute the maximizer $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T+\unicode[STIX]{x0394}T}$ on a larger time interval $[0,T+\unicode[STIX]{x0394}T]$, or $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0}+\unicode[STIX]{x0394}{\mathcal{E}}_{0},T}$ for a larger initial enstrophy ${\mathcal{E}}_{0}+\unicode[STIX]{x0394}{\mathcal{E}}_{0}$, for some sufficiently small $\unicode[STIX]{x0394}T$ or $\unicode[STIX]{x0394}{\mathcal{E}}_{0}$. Since in the limit $T\rightarrow 0$ solutions of the finite-time optimization problem 3.1 coincide with the solutions of the instantaneous optimization problem 2.1, for small initial enstrophy values ${\mathcal{E}}_{0}$ the instantaneous maximizers $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$, cf. figure 1, are used as ‘seeds’ to initiate the computation of the maximizing branch for the given value of ${\mathcal{E}}_{0}$, i.e. as the initial guess for $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\unicode[STIX]{x0394}T}$. The procedure outlined above is summarized as algorithm 1. For larger initial enstrophy values, it is also convenient to perform continuation with respect to ${\mathcal{E}}_{0}$ with $T$ fixed and in such cases essentially the same procedure applies, except that the order of the two outermost ‘repeat’ loops in algorithm 1 is reversed. While there exist alternatives to the continuation approach, provided $\unicode[STIX]{x0394}{\mathcal{E}}_{0}$ and $\unicode[STIX]{x0394}T$ are sufficiently small, this technique in fact results in the fastest convergence of iterations (3.3) and also ensures that computed optimal initial data lie on a maximizing branch.

Finally, we remark that the ‘optimize-then-discretize’ approach adopted here has the key advantage that in such a continuous formulation the expressions representing the sensitivity of the objective functional ${\mathcal{E}}_{T}(\boldsymbol{u}_{0})$, i.e. the gradients $\unicode[STIX]{x1D6FB}^{L_{2}}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$ and $\unicode[STIX]{x1D735}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$, are independent of the specific discretization approach chosen to evaluate them. This should be contrasted with the discrete (‘discretize-then-optimize’) formulation, where a change of the discretization method would require rederivation of the gradient expressions. In addition, the continuous formulation allows us to strictly enforce the regularity of maximizers required in problem 3.1. Last, but not least, the continuous formulation is also arguably better adapted to study problems motivated by questions in mathematical analysis such as the problems considered here. Our strategy for the numerical implementation of the different elements of algorithm 1 is presented in the next section.

4 Numerical approach

In this section we briefly describe the key elements of the numerical approach used to implement algorithm 1 and also comment on the validation strategies we employed. Evaluation of the objective functional (3.1) requires solution of the Navier–Stokes system (2.1) on the time interval $[0,T]$ with the given initial data $\boldsymbol{u}_{0}$. This system is solved numerically with an approach combining a pseudo-spectral approximation of spatial derivatives with a fourth-order semi-implicit Runge–Kutta method (Bewley Reference Bewley2009) used to discretize the problem in time. In the evaluation of the nonlinear term dealiasing is performed using the Gaussian filtering approach proposed by Hou & Li (Reference Hou and Li2007). Massively parallel implementation based on MPI (message passing interface) and using the software library FFTW (Frigo & Johnson Reference Frigo and Johnson2003) to perform Fourier transforms allowed us to employ resolutions varying from $128^{3}$ to $512^{3}$ in the low-enstrophy and high-enstrophy cases, respectively. Solution of optimization problem 3.1 for a large initial enstrophy ${\mathcal{E}}_{0}\lessapprox 1000$ and an intermediate length $T$ of the time interval typically required a computational time of $\mathit{O}(10^{2})$ hours on $\mathit{O}(10^{2})$ CPU cores. A number of different diagnostics were checked to ensure that all flow solutions discussed in § 5 are well resolved. We refer the reader to the dissertation by Ayala (Reference Ayala2014) for additional details and a validation of this approach.

In addition to the Navier–Stokes system (2.1), evaluation of the $L^{2}$ gradient (3.10) also requires solution of the adjoint system (3.8). This problem is solved using essentially the same numerical approach as used to solve the Navier–Stokes system (2.1). The velocity field $\boldsymbol{u}=\boldsymbol{u}(t,\boldsymbol{x})$ needed to evaluate the coefficients and the terminal condition in the adjoint system is saved at discrete time levels during solution of the Navier–Stokes system (since the considered time intervals $[0,T]$ are not very long, data for entire flow evolutions could be stored with our available temporary storage resources). The accuracy of the evaluation of the gradient $\unicode[STIX]{x1D6FB}^{L^{2}}{\mathcal{E}}_{T}(\boldsymbol{u}_{0})$, determined as described in § 3, is verified by examining the quantity $\unicode[STIX]{x1D705}(\unicode[STIX]{x1D716}):=\unicode[STIX]{x1D716}^{-1}[{\mathcal{E}}_{T}(\boldsymbol{u}_{0}+\unicode[STIX]{x1D716}\boldsymbol{u}_{0}^{\prime })-{\mathcal{E}}_{T}(\boldsymbol{u}_{0})]/\langle \unicode[STIX]{x1D6FB}^{L_{2}}{\mathcal{E}}_{T}(\boldsymbol{u}_{0}),\boldsymbol{u}_{0}^{\prime }\rangle _{L^{2}(\unicode[STIX]{x1D6FA})}$ which represents the ratio of a forward finite-difference approximation of the Gâteaux differential ${\mathcal{E}}_{T}^{\prime }(\boldsymbol{u}_{0};\boldsymbol{u}_{0}^{\prime })$ and its expression given in terms of the Riesz formula (3.4) with the gradient given in (3.10). The evaluation of gradients is validated when $\unicode[STIX]{x1D705}(\unicode[STIX]{x1D716})\approx 1$ for intermediate values of $\unicode[STIX]{x1D716}$ and different choices of $\boldsymbol{u}_{0}$ and $\boldsymbol{u}_{0}^{\prime }$. For such intermediate values of $\unicode[STIX]{x1D716}$, we observe that, as expected, $|\unicode[STIX]{x1D705}(\unicode[STIX]{x1D716})-1|$ is reduced as the numerical discretization parameters are refined (when $\unicode[STIX]{x1D716}\rightarrow 0$, $|\unicode[STIX]{x1D705}(\unicode[STIX]{x1D716})|$ becomes unbounded as a result of round-off errors due to subtractive cancellation, whereas when $\unicode[STIX]{x1D716}$ is large, $\unicode[STIX]{x1D705}(\unicode[STIX]{x1D716})$ deviates from unity because of truncation errors in the finite-difference approximation of the Gâteaux differential, both of which are well-known effects (Ayala Reference Ayala2014)).

As regards the computation of the Sobolev $H^{1}$ gradients, cf. (3.11), the parameter $\ell _{1}$ is adjusted during the optimization iterations and is chosen so that $\ell _{1}\in [\ell _{min},\ell _{max}]$, where $\ell _{min}$ is the length scale associated with the spatial resolution $N=128,256,512$ used for computations and $\ell _{max}$ is the characteristic length scale of the domain $\unicode[STIX]{x1D6FA}$, that is, $\ell _{min}\sim \mathit{O}(1/N)$ and $\ell _{max}\sim \mathit{O}(1)$. We remark that, given the equivalence of the Sobolev inner products (3.2) corresponding to different values of $\ell _{1}$ (as long as $\ell _{1}\neq 0$), these choices do not affect the maximizers found, but only how rapidly they are approached by iterations (3.3). The computational results presented in the next section have been thoroughly validated to ensure they are converged with respect to refinement of the different numerical parameters discussed above.

5 Computational results

In this section we proceed to present our computational results obtained by solving the finite-time optimization problem 3.1 for a broad range of initial enstrophies ${\mathcal{E}}_{0}$ and lengths $T$ of the time interval. The ultimate goal is to understand what is the largest growth of enstrophy possible under the 3-D Navier–Stokes dynamics, in particular, whether this growth could saturate estimate (2.10) and become unbounded in finite time, then how this maximum growth depends on the initial enstrophy ${\mathcal{E}}_{0}$ and, finally, what is the flow mechanism realizing this extreme behaviour. To address these questions we compute branches of maximizing solutions corresponding to different fixed values of the initial enstrophy ${\mathcal{E}}_{0}$ and increasing lengths $T$ of the time window. For smaller values of ${\mathcal{E}}_{0}$, this can be done using algorithm 1 in which solutions of the instantaneous maximization problem 2.1 are used to initialize the computation of the maximizing branches for short optimization times $T$. Since for large initial enstrophy values (${\mathcal{E}}_{0}>100$), the instantaneous problem 2.1 is harder to solve because of increased resolution requirements (more on this below), in such cases it is more convenient to initiate computation of the new branches corresponding to increased values of ${\mathcal{E}}_{0}$ by performing continuation with respect to ${\mathcal{E}}_{0}$ for some intermediate values of $T$. Only then can we again perform continuation with respect to $T$ at a new fixed value of ${\mathcal{E}}_{0}$.

Figure 2.

Enstrophy ${\mathcal{E}}(\boldsymbol{u}(t))$ as a function of time $t$ obtained from the solution of the Navier–Stokes system (2.1) with the initial condition $\boldsymbol{u}_{0}$ given by the maximizer $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ of the instantaneous optimization problem 2.1 (blue dashed line) and the maximizers $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ of the finite-time optimization problem 3.1 (red solid lines) for the indicated values of ${\mathcal{E}}_{0}$ and $T$ (the curves corresponding to the optimal lengths of the time window, $\widetilde{T}_{{\mathcal{E}}_{0}}=0.3$ for ${\mathcal{E}}_{0}=50$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.23$ for ${\mathcal{E}}_{0}=200$, are marked with thick lines, whereas the insets represent magnifications of the initial stages of evolution).

A typical time evolution of enstrophy ${\mathcal{E}}(\boldsymbol{u}(t))$ is presented in figure 2, where we show the results produced by solving the Navier–Stokes system (2.1) with the initial data $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ and $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained as the solutions of the instantaneous and finite-time optimization problems 2.1 and 3.1 for ${\mathcal{E}}_{0}=50$ and ${\mathcal{E}}_{0}=200$, and different $T$. We note that when the instantaneously optimal initial data $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ is used, the enstrophy grows very rapidly for short times which is followed by an immediate depletion of its growth, as already analysed by Ayala & Protas (Reference Ayala and Protas2017). On the other hand, when the optimal initial data $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained as a solution of the finite-time optimization problem 3.1 with ‘long’ time windows $T$ is used, the enstrophy grows very slowly at first (or even decreases when $T$ is sufficiently large), but eventually a much larger growth is achieved at the end of the time window, i.e. at $t=T$. As is also evident from figure 2, when $T$ is sufficiently large, the enstrophy achieves a well-defined maximum as a function of time and there exists a value of $T$ for which the maximum enstrophy is achieved precisely at the end of the optimization window, i.e. for $t=T$. Clearly, this value, which we shall denote $\widetilde{T}_{{\mathcal{E}}_{0}}:=\text{argmax}_{T>0}\,{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$, defines the optimal length of the optimization interval in problem 3.1 for which the largest maximum enstrophy is achieved for a given ${\mathcal{E}}_{0}$. Then, ${\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}})=\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ is the largest enstrophy achievable using an initial condition with enstrophy ${\mathcal{E}}_{0}$. Approximating the optimal length of the time interval $\widetilde{T}_{{\mathcal{E}}_{0}}$ for a given value of the initial enstrophy ${\mathcal{E}}_{0}$ requires solution of the finite-time optimization problem 3.1 for several different values of $T$.

Figure 3.

Evolution of (black thick solid lines) the total enstrophy ${\mathcal{E}}(\boldsymbol{u}(t))$ and its components (red thin solid lines) ${\mathcal{E}}_{1}(\boldsymbol{u}(t))$, (blue dotted lines) ${\mathcal{E}}_{2}(\boldsymbol{u}(t))$ and (pink dashed lines) ${\mathcal{E}}_{3}(\boldsymbol{u}(t))$ in the solution of the Navier–Stokes system (2.1) with the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 for the indicated values of ${\mathcal{E}}_{0}$ and $T$. Panels (a) and (b–d) correspond to the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ from the symmetric and asymmetric branch, respectively.

Closer inspection of the solutions to the Navier–Stokes system (2.1) corresponding to the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving problem 3.1 for ${\mathcal{E}}_{0}\lessapprox 100$ and ${\mathcal{E}}_{0}\gtrapprox 100$ indicates that these solutions have in fact distinct properties in terms of symmetry. To analyse this issue, in figure 3 we show the time histories of the enstrophy components ${\mathcal{E}}_{i}(\boldsymbol{u}(t))$, $i=1,2,3$, associated with the three coordinate directions defined as

(5.1)$$\begin{eqnarray}{\mathcal{E}}_{i}(\boldsymbol{u}(t)):=\int _{\unicode[STIX]{x1D6FA}}|(\unicode[STIX]{x1D735}\times \boldsymbol{u}(t))\boldsymbol{\cdot }\boldsymbol{e}_{i}|^{2}\,\text{d}\boldsymbol{x},\quad i=1,2,3,\end{eqnarray}$$

where $\boldsymbol{e}_{1}$, $\boldsymbol{e}_{2}$, $\boldsymbol{e}_{3}$ are the unit vectors of the Cartesian coordinate system and we have the obvious identity $\forall t\;{\mathcal{E}}(\boldsymbol{u}(t))=\sum _{i=1}^{3}{\mathcal{E}}_{i}(\boldsymbol{u}(t))$. In solutions corresponding to optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained for ${\mathcal{E}}_{0}\lessapprox 100$ we have $\forall t\in [0,\widetilde{T}_{{\mathcal{E}}_{0}}]\quad {\mathcal{E}}_{1}(\boldsymbol{u}(t))={\mathcal{E}}_{2}(\boldsymbol{u}(t))={\mathcal{E}}_{3}(\boldsymbol{u}(t))=(1/3){\mathcal{E}}(\boldsymbol{u}(t))$, indicating that the enstrophy is equipartitioned among the three coordinate directions, cf. figure 3(a) for the case ${\mathcal{E}}_{0}=100$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.27$. On the other hand, this special property is absent in the solutions corresponding to optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained for ${\mathcal{E}}_{0}\gtrapprox 100$, cf. figure 3(b–d), in which the three enstrophy components all contain different fractions of the total enstrophy. Moreover, in these cases the ordering of ${\mathcal{E}}_{1}(\boldsymbol{u}(t))$, ${\mathcal{E}}_{2}(\boldsymbol{u}(t))$ and ${\mathcal{E}}_{3}(\boldsymbol{u}(t))$ changes during the flow evolution. For example, for ${\mathcal{E}}_{0}=1000$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.12$, cf. figure 3(d), while we have ${\mathcal{E}}_{3}(\boldsymbol{u}(0))<{\mathcal{E}}_{1}(\boldsymbol{u}(0))<{\mathcal{E}}_{2}(\boldsymbol{u}(0))$ and ${\mathcal{E}}_{3}(\boldsymbol{u}(\widetilde{T}_{{\mathcal{E}}_{0}}))<{\mathcal{E}}_{1}(\boldsymbol{u}(\widetilde{T}_{{\mathcal{E}}_{0}}))<{\mathcal{E}}_{2}(\boldsymbol{u}(\widetilde{T}_{{\mathcal{E}}_{0}}))$ at the initial and final instant of time, this ordering changes six times and is completely reversed during the flow evolution. Such transfer of enstrophy between different vorticity components is known to signal the phenomenon of reconnection of vortex lines (Hussain & Duraisamy Reference Hussain and Duraisamy2011; Jaque & Fuentes Reference Jaque and Fuentes2017). We add that in all the flow evolutions discussed above the helicity ${\mathcal{H}}(\boldsymbol{u}(t)):=\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{u}(t)\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\times \boldsymbol{u}(t))\,\text{d}\boldsymbol{x}$ remains identically equal to zero, i.e. $\forall t\in [0,T]$${\mathcal{H}}(\boldsymbol{u}(t))=0$.

Figure 4.

Maximum attained enstrophy ${\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ as a function of the length $T$ of the time interval over which maximization is performed in problem 3.1 for initial enstrophies (a) $0<{\mathcal{E}}_{0}\leqslant 100$ and (b) $100<{\mathcal{E}}_{0}\leqslant 1000$. The blue dashed and red solid curves correspond to flow evolutions starting from, respectively, symmetric and asymmetric optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ with the same value of the initial enstrophy and different $T$ (the trends with the increase of ${\mathcal{E}}_{0}$ are indicated with arrows, whereas solid symbols represent the values of ${\mathcal{E}}_{0}$ and $T$ for which problem 3.1 was solved).

We thus have evidence for the existence of maximizers $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ of two distinct types: those for which the flow evolution is characterized by the properties $\forall t\in [0,\widetilde{T}_{{\mathcal{E}}_{0}}]{\mathcal{E}}_{1}(\boldsymbol{u}(t))={\mathcal{E}}_{2}(\boldsymbol{u}(t))={\mathcal{E}}_{3}(\boldsymbol{u}(t))=(1/3){\mathcal{E}}(\boldsymbol{u}(t))$ and ${\mathcal{H}}(\boldsymbol{u}(t))=0$, and those for which these properties do not hold. We will refer to these families of maximizers as ‘symmetric’ and ‘asymmetric’, respectively. The symmetric maximizers are normally found when the finite-time optimization problem 3.1 is solved for ${\mathcal{E}}_{0}<100$, whereas the asymmetric maximizers are found when this problem is solved for ${\mathcal{E}}_{0}>100$. However, maximizers of both types were found for ${\mathcal{E}}_{0}\approx 100$. These results are illustrated in figure 4(a,b) where we show the dependence of the maximum attained enstrophy ${\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ on the length $T$ of the time interval over which optimization was performed in problem 3.1 for initial enstrophies ${\mathcal{E}}_{0}\leqslant 100$ and ${\mathcal{E}}>100$, respectively. In figure 4(a) we see that indeed two distinct branches of symmetric and asymmetric maximizers are simultaneously present when ${\mathcal{E}}_{0}=100$, thus demonstrating that these maximizers are nonunique. There are indications that symmetric maximizers also exist for ${\mathcal{E}}_{0}>100$ and asymmetric maximizers exist for ${\mathcal{E}}_{0}<100$; however, since they correspond to suboptimal local maxima, it is very difficult to compute them with the gradient approach (3.3). It is possible to capture both the symmetric and asymmetric branch for ${\mathcal{E}}_{0}=100$, because for this value of the initial enstrophy the maximum attained enstrophy ${\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ is comparable for both branches. The asymmetry of the optimal initial data $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ on the asymmetric branch vanishes as $T\rightarrow 0.13$, which indicates that, for ${\mathcal{E}}_{0}=100$, the asymmetric branch bifurcates off the symmetric branch at $T\approx 0.13$. In figure 4(b) a significant, albeit finite, increase of the enstrophy growth is observed as a function of ${\mathcal{E}}_{0}$ when asymmetric maximizers $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ are used as the initial data for system (2.1), and this figure also confirms that, for each value of the initial enstrophy ${\mathcal{E}}_{0}$, the maximum enstrophy ${\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ achieves a well-defined maximum for a certain optimal time $\widetilde{T}_{{\mathcal{E}}_{0}}$. We note that the accuracy with which $\widetilde{T}_{{\mathcal{E}}_{0}}$ can be determined for given ${\mathcal{E}}_{0}$ is limited by the number of time intervals $[0,T]$ for which the finite-time optimization problem 3.1 can be solved, which is a matter of the computational cost. However, we use spline interpolation to improve the accuracy with which the quantities $\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ and $\widetilde{T}_{{\mathcal{E}}_{0}}$ are determined for a given initial enstrophy ${\mathcal{E}}_{0}$.

Figure 5.

Dependence of the maximum enstrophy growth $\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ on the initial enstrophy using (a) linear scaling for small ${\mathcal{E}}_{0}$ and (b) logarithmic scaling for large ${\mathcal{E}}_{0}$. The blue squares and red circles correspond to, respectively, the symmetric and asymmetric branches of solutions of the finite-time optimization problem 3.1 for different ${\mathcal{E}}_{0}$. The black solid line in (b) represents the least-squares fit (5.2).

Figure 6.

Dependence of the time $\widetilde{T}_{{\mathcal{E}}_{0}}$ when the maximum enstrophy $\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ is achieved (with asymmetric maximizers) on the initial enstrophy ${\mathcal{E}}_{0}$. The red circles represent the interpolated data corresponding to solutions of the finite-time optimization problem 3.1 for different ${\mathcal{E}}_{0}$, whereas the black solid line represents the least-squares fit (5.3).

As the central result of this study, the maximum growth of enstrophy $\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ achieved using initial conditions with enstrophy ${\mathcal{E}}_{0}$ is shown in figure 5 as a function of ${\mathcal{E}}_{0}$. In figure 5(a), which corresponds to small values of ${\mathcal{E}}_{0}$, we see that the dependence of $\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ on ${\mathcal{E}}_{0}$ appears slightly different for the symmetric and asymmetric branches, i.e. for ${\mathcal{E}}_{0}<100$ and ${\mathcal{E}}_{0}>100$. As is evident from figure 5(b), for asymmetric maximizers this growth follows a well-defined power–law relation

(5.2)$$\begin{eqnarray}\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})\sim (0.224\pm 0.006){\mathcal{E}}_{0}^{1.490\pm 0.004}\end{eqnarray}$$

obtained by performing a least-squares fit of the relation between $\log _{10}\! [\!\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})\!]$ and $\log _{10}{\mathcal{E}}_{0}$ for $100\leqslant {\mathcal{E}}_{0}\leqslant 1000$. The fit is performed in the logarithmic coordinates in order for the least-squares error not to be dominated by contributions from the data corresponding to large values of ${\mathcal{E}}_{0}$. As regards the times when the maxima are achieved for different values of ${\mathcal{E}}_{0}$, when asymmetric maximizers are considered and the initial enstrophy is sufficiently large $({\mathcal{E}}_{0}\gtrapprox 400)$, figure 6 shows that $\widetilde{T}_{{\mathcal{E}}_{0}}$ also exhibits a power–law dependence on ${\mathcal{E}}_{0}$, described by the relation

(5.3)$$\begin{eqnarray}\widetilde{T}_{{\mathcal{E}}_{0}}\sim (4.0\pm 1.1){\mathcal{E}}_{0}^{-0.51\pm 0.04}\end{eqnarray}$$

obtained from a least-squares fit of the relation between ${\mathcal{E}}_{0}$ and $\widetilde{T}_{{\mathcal{E}}_{0}}$. We observe that for the maximum enstrophy $\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})$ the power–law dependence on ${\mathcal{E}}_{0}$ becomes evident starting from smaller values of ${\mathcal{E}}_{0}$ and the fit is characterized by smaller error bars than for $\widetilde{T}_{{\mathcal{E}}_{0}}$, cf. figures 5(b) and 6. We also remark that the dependence of $\widetilde{T}_{{\mathcal{E}}_{0}}$ on ${\mathcal{E}}_{0}$, given in (5.3), appears different from the form of the upper bounds on the maximum blow-up time obtained by Ohkitani (Reference Ohkitani2016). We add that the times $t_{0}$ when the upper bound in estimate (2.10) becomes unbounded are about 8–10 orders of magnitude shorter than the times $\widetilde{T}_{{\mathcal{E}}_{0}}$ shown in figure 6.

Figure 7.

Flow trajectories corresponding to the optimal initial data $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ with different ${\mathcal{E}}_{0}\in [100,1000]$, cf. figure 3(b,c), shown using the coordinates $\{{\mathcal{E}},\text{d}{\mathcal{E}}/\text{d}t\}$ (blue solid lines with the arrow indicating the trend with the increase of ${\mathcal{E}}_{0}$). The thick red line represents the relation $\text{d}{\mathcal{E}}/\text{d}t=1.72\times 10^{-3}{\mathcal{E}}^{3}$, cf. figure 1, whereas the dashed black line represents the relation $\text{d}{\mathcal{E}}/\text{d}t=10^{2}{\mathcal{E}}^{0.85}$.

In order to understand how close the flow evolutions corresponding to the optimal initial data $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ come to saturating a priori bounds on the rate of growth of enstrophy, cf. (2.9), in figure 7 we plot the corresponding trajectories using the coordinates $\{{\mathcal{E}},\text{d}{\mathcal{E}}/\text{d}t\}$, such that each trajectory is parameterized by time $t$ (since the logarithmic scale is used, initial parts of the trajectories when $\text{d}{\mathcal{E}}/\text{d}t\lessapprox 0$ are not shown). The slope of the tangent to each of the curves thus represents the exponent $\unicode[STIX]{x1D6FC}$ characterizing the instantaneous rate of enstrophy production $\text{d}{\mathcal{E}}/\text{d}t\sim {\mathcal{E}}^{\unicode[STIX]{x1D6FC}}$. In figure 7 we also indicate the relation $\text{d}{\mathcal{E}}/\text{d}t=1.72\times 10^{-3}{\mathcal{E}}^{3}$ describing the maximum rate of enstrophy growth realized by the solutions of the instantaneous optimization problem 2.1 (Ayala & Protas Reference Ayala and Protas2017), cf. figure 1 (we note that the prefactor in this relation is approximately seven orders of magnitude smaller than the prefactor in estimate (2.9)). We observe that the rate of growth of enstrophy realized by the trajectories corresponding to the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ is at all times and for all values of ${\mathcal{E}}_{0}$ several orders of magnitude smaller than the maximum rate of growth achieved by the instantaneous maximizers $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$, cf. figure 1. On the other hand, in figure 7 we also note that the exponent $\unicode[STIX]{x1D6FC}$ characterizing the growth of enstrophy can be much larger than three at intermediate stages along these trajectories (we recall from the discussion in § 2 that the rate of enstrophy production $\text{d}{\mathcal{E}}/\text{d}t\sim {\mathcal{E}}^{\unicode[STIX]{x1D6FC}}$ must be sustained with $\unicode[STIX]{x1D6FC}>2$ over a sufficiently long time interval for enstrophy to become unbounded in finite time). We add that at the final stages of the flow evolutions before the enstrophy maximum is reached at $t=\widetilde{T}_{{\mathcal{E}}_{0}}$ the enstrophy is amplified at the approximate rate $\text{d}{\mathcal{E}}/\text{d}t\sim {\mathcal{E}}^{0.85}$.

Figure 8.

Optimal initial conditions (a) $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ obtained by solving the instantaneous optimization problem 2.1, and (b–h) $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 for the initial enstrophy ${\mathcal{E}}_{0}=100$ and indicated lengths $T$ of the time interval. Optimal initial conditions from both the symmetric and asymmetric branch are shown. Shades of red correspond to the magnitude of the vorticity $|(\unicode[STIX]{x1D735}\times \widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}})(\boldsymbol{x})|$ or $|(\unicode[STIX]{x1D735}\times \widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})(\boldsymbol{x})|$ (see the colour bars), whereas white curves represent vortex lines chosen to pass through regions with strong vorticity.

Figure 9.

Energy spectra of the optimal initial conditions $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ and $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the instantaneous optimization problem 2.1 (red circles) and the finite-time optimization problem 3.1 for the initial enstrophy ${\mathcal{E}}_{0}=100$ and $T=0.1$ (pink squares), $T=0.2$ (blue diamonds) and $T=0.3$ (black triangles). In the latter case the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ come from the asymmetric branch.

We now move on to characterize the structure of the optimal initial conditions which give rise to the maximum enstrophy growth reported in figures 2, 4 and 5. First, in figure 8 we analyse the structure of both the symmetric and asymmetric initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 for a fixed ${\mathcal{E}}_{0}=100$ and different $T$, and also show the instantaneously optimal initial condition $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ obtained as a solution of problem 2.1 for the same value of ${\mathcal{E}}_{0}$ (Ayala & Protas Reference Ayala and Protas2017). These initial conditions are presented in figure 8 in terms of their vorticity $\unicode[STIX]{x1D735}\times \widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ and $\unicode[STIX]{x1D735}\times \widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$, whose magnitude is visualized via volume rendering and, in addition, we also show a number of vortex lines (i.e. lines everywhere tangent to the vorticity field) chosen to pass through regions with strong vorticity. Such an approach allows us to simultaneously assess both the intensity and the structure of the vorticity field. In figure 8 we see that as $T$ increases the structure of the optimal initial condition gradually changes from two colliding vortex rings characterizing the instantaneous maximizers $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ (Lu & Doering Reference Lu and Doering2008; Ayala & Protas Reference Ayala and Protas2017) to a more complex vorticity distribution filling the entire flow domain. There are also evident differences between the optimal initial conditions belonging to the symmetric and asymmetric branches, cf. figures 8(b,c,e,g) versus 8(d,f,h) – while the former retain the structure of deformed rings, the latter develop a tubular form. Since the asymmetric optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ give rise to a much larger growth of enstrophy, cf. figure 5, hereafter we will exclusively focus on this case.

The fact that the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 are less localized than the maximizers of the instantaneous optimization problem 2.1 is also evident from figure 9 showing their energy spectra $e(|\boldsymbol{k}|,0)$ defined as

(5.4)$$\begin{eqnarray}e(|\boldsymbol{k}|,t):=\int _{S_{|\boldsymbol{k}|}}|\boldsymbol{k}|^{2}|[\widehat{\boldsymbol{u}}(t)]_{\boldsymbol{ k}}|^{2}\,\text{d}\unicode[STIX]{x1D70E},\end{eqnarray}$$

where $[\widehat{\boldsymbol{u}}(t)]_{\boldsymbol{k}}$ are the Fourier coefficients of the velocity field $\boldsymbol{u}(t)$, $\unicode[STIX]{x1D70E}$ is the solid angle in the wavenumber space and $S_{|\boldsymbol{k}|}$ denotes the sphere of radius $|\boldsymbol{k}|$ in this space, such that ${\mathcal{K}}(\boldsymbol{u}(t))=\int _{0}^{\infty }e(k,t)\,\text{d}k$ (with some abuse of notation justified by simplicity, here we have treated the wavevector $\boldsymbol{k}$ as a continuous rather than discrete variable). In figure 9 we see that while all the fields $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ remain real-analytic (their energy spectra vanish exponentially fast for high wavenumbers $|\boldsymbol{k}|$), interestingly, the ones corresponding to intermediate lengths $T\approx 0.2$ of the optimization interval are in fact the smoothest in the sense that their energy spectra $e(|\boldsymbol{k}|,0)$ decay fastest with $|\boldsymbol{k}|$. In particular, the energy spectrum of the instantaneously optimal initial condition $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ decays much slower. These observations indicate that the finite-time optimization problems 3.1 corresponding to intermediate and longer time intervals $T$ are in fact less demanding in terms of space resolution $N$, which is why for larger values of the initial enstrophy ${\mathcal{E}}_{0}$ it is easier to perform continuation with respect to ${\mathcal{E}}_{0}$ at some fixed intermediate $T$, rather than the other way round starting from the instantaneous maximizers $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ which are hard to compute when ${\mathcal{E}}_{0}>100$.

Figure 10.

Asymmetric optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ obtained by solving the finite-time optimization problem 3.1 for different indicated values of the initial enstrophy ${\mathcal{E}}_{0}$ and the corresponding optimal lengths $\widetilde{T}_{{\mathcal{E}}_{0}}$ of the time interval (cf. figure 6). Shades of red correspond to the magnitude of the vorticity $|(\unicode[STIX]{x1D735}\times \widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}})(\boldsymbol{x})|$ (see the colour bars), whereas white curves represent vortex lines chosen to pass through regions with strong vorticity.

Figure 11.

Vorticity components of the asymmetric optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ obtained by solving the finite-time optimization problem 3.1 for the initial enstrophy ${\mathcal{E}}_{0}=500$ and the corresponding optimal length $\widetilde{T}_{{\mathcal{E}}_{0}}=0.17$ of the time interval, cf. figure 10(e).

The optimal asymmetric initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ obtained for the optimal time intervals $[0,\widetilde{T}_{{\mathcal{E}}_{0}}]$ are shown in figure 10 for increasing values of the initial enstrophy ${\mathcal{E}}_{0}$. We note that the structure of these fields does not much change with ${\mathcal{E}}_{0}$, except that, as expected, the vorticity magnitude $|(\unicode[STIX]{x1D735}\times \widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}})(\boldsymbol{x})|$ in the vortex regions increases. In order to shed additional light on the structure of these optimal initial conditions, the field $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ corresponding to the initial enstrophy ${\mathcal{E}}_{0}=500$, cf. figure 10(e), is further analysed in figure 11 in terms of its different vorticity components $[\widetilde{\unicode[STIX]{x1D714}}_{1}(t),\widetilde{\unicode[STIX]{x1D714}}_{2}(t),\widetilde{\unicode[STIX]{x1D714}}_{3}(t)]^{\text{T}}:=\unicode[STIX]{x1D735}\times \boldsymbol{u}(t)$. We observe that the optimal initial condition has in fact the form of three perpendicular pairs of anti-parallel vortex tubes perturbed near the regions where they intersect (these perturbations allow the vorticity field to satisfy the divergence-free condition $\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\times \widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}})=0$). The relative magnitudes of the three vorticity components are given by the approximate relations $0.76:1:0.63$ representing the ratios $\widetilde{\unicode[STIX]{x1D714}}_{1}:\widetilde{\unicode[STIX]{x1D714}}_{2}:\widetilde{\unicode[STIX]{x1D714}}_{3}$.

Figure 12.

Time evolution of the energy spectrum in the solution of the Navier–Stokes system (2.1) with the optimal asymmetric initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 for the initial enstrophy ${\mathcal{E}}_{0}=500$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.17$. The contour lines represent the level sets of the function $\log _{10}e(|\boldsymbol{k}|,t)$ with the distance between two nearby lines corresponding to one order of magnitude.

Next, we analyse details of the flow evolution corresponding to the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$. In figure 12 we illustrate the evolution of the energy spectrum $e(|\boldsymbol{k}|,t)$, cf. (5.4), in time in the solution of the Navier–Stokes system (2.1) with the optimal asymmetric initial condition $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 with ${\mathcal{E}}_{0}=500$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.17$, cf. figure 11 (qualitatively similar behaviour is also observed for other values of ${\mathcal{E}}_{0}>100$). We note that at early stages of the flow evolution the energy spectrum recedes which corresponds to the initially very slow growth of enstrophy already noted in figure 2. Then, the energy begins to flow towards larger wavenumbers such that the energy spectrum becomes the most developed close to the final time $t=\widetilde{T}_{{\mathcal{E}}_{0}}$ when the enstrophy maximum is achieved. As regards this second stage, there are clearly two phases ending at instances of time corresponding approximately to the times when the ordering of ${\mathcal{E}}_{1}(\boldsymbol{u}(t))$, ${\mathcal{E}}_{2}(\boldsymbol{u}(t))$ and ${\mathcal{E}}_{3}(\boldsymbol{u}(t))$ changes, cf. figure 3(c). This indicates that there are in fact two reconnection events occurring during the flow evolution: one approximately in the middle of the time interval $[0,\widetilde{T}_{{\mathcal{E}}_{0}}]$ and another one close to its end.

Figure 13.

Vorticity components $\widetilde{\unicode[STIX]{x1D714}}_{1}$ (blue), $\widetilde{\unicode[STIX]{x1D714}}_{2}$ (red) and $\widetilde{\unicode[STIX]{x1D714}}_{3}$ (green) at the indicated time instances in the solution of the Navier–Stokes system (2.1) with the optimal asymmetric initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 for the initial enstrophy ${\mathcal{E}}_{0}=500$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.17$. Supporting movie 1 is available online at https://doi.org/10.1017/jfm.2020.204.

Figure 14.

Reconnection event occurring at the time $t=0.09$ in the solution of the Navier–Stokes system (2.1) with the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 for ${\mathcal{E}}_{0}=500$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.17$, cf. figure 13(c). White lines represent the vortex lines in the neighbourhood of the reconnection point, whereas the blue regions correspond to locations where $|(\unicode[STIX]{x1D735}\times \boldsymbol{u})(\boldsymbol{x})|\approx 0$.

Finally, we discuss the flow evolution in the physical space and the three vorticity components $\widetilde{\unicode[STIX]{x1D714}}_{1}$, $\widetilde{\unicode[STIX]{x1D714}}_{2}$ and $\widetilde{\unicode[STIX]{x1D714}}_{3}$ are shown in figure 13 at various times $t\in [0,\widetilde{T}_{{\mathcal{E}}_{0}}]$. Since the same colour scale is used in all panels in this figure, some features of the optimal initial condition $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ evident in figure 11 are not visible in figure 13(a). In agreement with figure 12, we observe that at the early stages of the flow evolution the vorticity field becomes less localized, cf. figure 13(b), which is followed by a rapid development of small scales, cf. figure 13(d,e). The final state $\unicode[STIX]{x1D735}\times \boldsymbol{u}(\widetilde{T}_{{\mathcal{E}}_{0}})$ is characterized by a fairly complicated structure with turbulent-like spatial complexity, cf. figure 13(f). In particular, it features a combination of tube-shaped and pancake-shaped regions of concentrated vorticity (animations corresponding to the vorticity fields shown in figure 13 and to optimal flow evolutions obtained for other values of ${\mathcal{E}}_{0}$ and $\widetilde{T}_{{\mathcal{E}}_{0}}$ are available online as movies 1, 2 and 3). To provide more insight into what drives this flow evolution, in figure 14 we illustrate a typical reconnection event occurring approximately at the time when the ordering of ${\mathcal{E}}_{1}(\boldsymbol{u}(t))$, ${\mathcal{E}}_{2}(\boldsymbol{u}(t))$ and ${\mathcal{E}}_{3}(\boldsymbol{u}(t))$ changes, cf. figures 3(c) and 13(c). We see that the bundle of vortex lines aligned with the direction $x_{2}$ and visible in the middle of the figure is split into two parts going in opposite directions, one of which attaches to a perpendicular bundle of vortex lines. In figure 14 we also visualize regions where $|(\unicode[STIX]{x1D735}\times \boldsymbol{u})(\boldsymbol{x})|\approx 0$, which is a necessary condition for vortex reconnection to occur (Hussain & Duraisamy Reference Hussain and Duraisamy2011; Jaque & Fuentes Reference Jaque and Fuentes2017). These regions are therefore potential loci of other reconnection events.

Figure 15.

Evolution of (a) the norm $\Vert \boldsymbol{u}(t)\Vert _{L^{3}}$ and (b) the seminorm $\Vert \boldsymbol{u}(t)\Vert _{{\dot{H}}^{1/2}}$ of the solutions of the Navier–Stokes system (2.1) with the optimal asymmetric initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 for ${\mathcal{E}}_{0}=100$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.27$, ${\mathcal{E}}_{0}=500$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.17$ and ${\mathcal{E}}_{0}=1000$ and $\widetilde{T}_{{\mathcal{E}}_{0}}=0.13$. Arrows indicates the trends with the increase of ${\mathcal{E}}_{0}$.

To close this section, in figure 15 we analyse the time evolution of the $L^{3}$ norm and the $H^{1/2}$ seminorm of the solutions of the Navier–Stokes system (2.1) with the optimal asymmetric initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ obtained by solving the finite-time optimization problem 3.1 for ${\mathcal{E}}_{0}=100,500,1000$, and with the corresponding optimal optimization time intervals $\widetilde{T}_{{\mathcal{E}}_{0}}$. These quantities, defined as

(5.5a)$$\begin{eqnarray}\displaystyle & \displaystyle \Vert \boldsymbol{u}(t)\Vert _{L^{3}}:=\left(\int _{\unicode[STIX]{x1D6FA}}|\boldsymbol{u}(t,\boldsymbol{x})|^{3}\,\text{d}\boldsymbol{x}\right)^{1/3}, & \displaystyle\end{eqnarray}$$

(5.5b)$$\begin{eqnarray}\displaystyle & \displaystyle \Vert \boldsymbol{u}(t)\Vert _{{\dot{H}}^{1/2}}:=\left[\mathop{\sum }_{\boldsymbol{k}\in \mathbb{Z}^{3}}|\boldsymbol{k}||[\widehat{\boldsymbol{u}}(t)]_{\boldsymbol{k}}|^{2}\right]^{1/2}, & \displaystyle\end{eqnarray}$$

$L^{3}$

$H^{1/2}$

$0<t<\widetilde{T}_{{\mathcal{E}}_{0}}$

$\widetilde{T}_{{\mathcal{E}}_{0}}$

6 Discussion and conclusions

In this study we have considered the question of the largest possible growth of enstrophy in finite time in 3-D Navier–Stokes flows starting from initial data with enstrophy ${\mathcal{E}}_{0}$. This problem is motivated by the open question concerning the global-in-time existence of classical (smooth) solutions of the 3-D Navier–Stokes system (Doering Reference Doering2009). One of the landmark results in this context is the conditional regularity result due to Foias & Temam (Reference Foias and Temam1989) asserting that solutions to the Navier–Stokes system (2.1) remain smooth and satisfy this system in the classical sense as long as the enstrophy remains finite. To probe this condition we have considered a family of optimization problems for the Navier–Stokes system (2.1), cf. problem 3.1, in which initial data $\boldsymbol{u}_{0}$ with prescribed enstrophy ${\mathcal{E}}_{0}$ is sought such that the enstrophy at the given time $T$ is maximized. Thus, in each case the problem is solved on a fixed time domain and there are two parameters, ${\mathcal{E}}_{0}$ and $T$. For each value of ${\mathcal{E}}_{0}$, the optimal time window $\widetilde{T}_{{\mathcal{E}}_{0}}$ is then determined by solving problem 3.1 for several different values of $T$, cf. figure 4. In principle, this last step could be eliminated by finding the optimal time window $\widetilde{T}_{{\mathcal{E}}_{0}}$ as a part of the solution of a suitably modified optimization problem 3.1 in which the final enstrophy would be simultaneously maximized with respect to the initial condition $\boldsymbol{u}_{0}$ and the length $T$ of the time window. However, such a modified optimization problem would then be of the free-boundary type, because its time domain $[0,T]$ would not be a priori known, leading to significant complications in its numerical solution. Therefore, it is preferable to solve problem 3.1 for different values of the two parameters ${\mathcal{E}}_{0}$ and $T$, which is facilitated by the continuation approach. For given values of ${\mathcal{E}}_{0}$ and $T$, problem 3.1 is solved computationally with a state-of-the-art adjoint-based gradient ascent technique. Our approach is formulated in the continuous (‘optimize-then-discretize’) setting with gradients defined in a suitable Sobolev space, cf. (3.4) and (3.11), which allows us to ensure that the optimal initial conditions possess the required minimum regularity. The governing and adjoint systems, (2.1) and (3.8), are solved efficiently with a massively parallel pseudo-spectral approach. Being based on the first-order optimality conditions, this gradient optimization approach makes it possible to find local maximizers only and it is not possible to guarantee that these maximizers are global. In order to partially address this issue, in addition to the continuation approach described in § 3, cf. algorithm 1, we have also undertaken an extensive search for other maximizing branches using suitably randomized initial guesses $\boldsymbol{u}^{0}$ in (3.3) which, however, did not yield any maximizers distinct from the two branches, a symmetric and an asymmetric one, already reported in § 5. Nevertheless, the existence of other branches of maximizers cannot be ruled out.

For every considered value of the initial enstrophy ${\mathcal{E}}_{0}$ on both the symmetric and asymmetric branch there exists a well-defined time $\widetilde{T}_{{\mathcal{E}}_{0}}$ such that the optimized growth of the enstrophy at this time is maximal. The structure of the optimal initial condition producing this growth is quite different on the symmetric and the asymmetric branch, and changes as the initial enstrophy increases, cf. figure 10. In the limit of large ${\mathcal{E}}_{0}$ the asymmetric branch dominates in terms of the maximum enstrophy growth, and the corresponding initial conditions have the form of three perpendicular pairs of anti-parallel vortex tubes, cf. figure 11. Thus, our subsequent analysis has focused on this family of initial data. Interestingly, as is evident from figure 9, these optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ are smoother than the solutions $\widetilde{\boldsymbol{u}}_{{\mathcal{E}}_{0}}$ of the instantaneous optimization problem 2.1 obtained for the same value of ${\mathcal{E}}_{0}$ (to be precise, the velocity fields are real-analytic in both cases, but in the latter case the energy spectra decay at a slower, though still exponential, rate). Consequently, solution of the finite-time optimization problem 3.1 requires a lower numerical resolution than the corresponding instantaneous problem 2.1 and, as a result, we have been able to solve the finite-time problem for values of the initial enstrophy ${\mathcal{E}}_{0}$ an order of magnitude higher than it was possible for the instantaneous maximization problem in Ayala & Protas (Reference Ayala and Protas2017). Further on this note, we also observe that the extremal flow trajectories originating from the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T}$ involve velocity fields which become even smoother at early stages of their evolution, cf. figures 12 and 13, before developing small-scale components at later stages when the maximum enstrophy is achieved. In other words, in these optimal flow evolutions the energy flux is initially from small to large scales, and then in the opposite direction near the end of the optimization interval $[0,\widetilde{T}_{{\mathcal{E}}_{0}}]$. This behaviour can be understood by considering system (2.7) from which it is evident that large enstrophy implies rapid energy dissipation. At the same time it is known that a 3-D Navier–Stokes flow may not generate significant enstrophy unless its energy is sufficiently large (Ayala & Protas Reference Ayala and Protas2017). Thus, in order to maximize its enstrophy at some large time $T$, an optimal flow seeks to reduce its production at early stages in order to conserve energy while rearranging itself strategically such that a sequence of bursts of energy to small scales will produce maximum enstrophy at the final time $T$.

The main finding of the present study is that the worst-case enstrophy-maximizing Navier–Stokes flow evolutions constructed by solving a suitable optimization problem always lead to a finite only growth of enstrophy. This suggests that unbounded growth of enstrophy allowed by estimate (2.10) cannot, in fact, be realized and, in light of the conditional regularity result (1.1), even under the extremal flow evolution there is no evidence for singularity formation in finite time. This also suggests the possibility of improving estimate (2.10) with the power–law relation (5.2) serving as a potential target. As a matter of course, the validity of these findings is restricted by the possibility that we may not have found true global maximizers of problem 3.1. To be more specific, the maximum enstrophy produced in flows with initial data with enstrophy ${\mathcal{E}}_{0}$ was found to scale in proportion to ${\mathcal{E}}_{0}^{3/2}$, cf. (5.2), which is the same power–law relation as observed by Ayala & Protas (Reference Ayala and Protas2011) in solutions of an analogous finite-time optimization problem for the 1-D Burgers equation. This scaling can be in fact predicted based on dimensional analysis. Noting that the largest flow structures have a characteristic dimension comparable to the size of the domain $L$, cf. figures 10 and 13, there are three physical parameters defining the problem, namely, ${\mathcal{E}}_{0}$, $L$ and the kinematic viscosity $\unicode[STIX]{x1D708}$. The only nontrivial way to combine them gives ${\mathcal{E}}_{T}\sim (L/\unicode[STIX]{x1D708}^{2})^{1/2}{\mathcal{E}}_{0}^{3/2}$. In addition, we can deduce that $T\sim L^{3/2}{\mathcal{E}}_{0}^{-1/2},$ which represents the scaling for $\widetilde{T}_{{\mathcal{E}}_{0}}$ observed in figure 6, cf. (5.3). In this context it would be interesting to know what form the optimal initial conditions would take if problem 3.1 were solved on an unbounded domain $\unicode[STIX]{x1D6FA}=\mathbb{R}^{3}$. Since for a fixed initial enstrophy ${\mathcal{E}}_{0}$ the optimal initial conditions discussed in § 5 would vanish in the limit $L\rightarrow \infty$, the solutions of this problem would likely have a completely different form. We add that in the instantaneous optimization problem in which the characteristic dimension of the flow structures does not scale with $L$, dimensional analysis also correctly predicts that $\text{d}{\mathcal{E}}/\text{d}t\sim {\mathcal{E}}^{3}/\unicode[STIX]{x1D708}^{3}$, cf. (2.9) (Lu & Doering Reference Lu and Doering2008).

Two distinct branches of maximizing initial conditions were found with the one consisting of asymmetric states dominating in terms of the maximum growth of enstrophy in the limit of large ${\mathcal{E}}_{0}$. The corresponding initial conditions have the form of three perpendicular pairs of anti-parallel vortex tubes and the flow evolution resulting from this initial data is quite complex, cf. figure 13, involving an initial flux of energy from small to large scales followed by a sequence of reconnection events, cf. figure 14. At its final stages the flow evolution has the spatio-temporal complexity of turbulent flows, although there is no evidence of a Kolmogorov-type $-5/3$ energy spectrum in figure 9. This can be explained by the fact that we consider an initial-value problem for system (2.1) and no statistically steady state is attained during the evolution. More work is still needed to understand in detail the physical mechanisms responsible for the extreme flow evolutions illustrated in figure 13, for example, in terms of changes of the topology of vortex lines. In particular, the fact that helicity remains zero during all flow evolutions starting from the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ implies the presence of some symmetries in the velocity field. We add that no evidence for reconnection was found in the flow evolutions starting from the optimal initial conditions on the symmetric branch.

In involving pairs of anti-parallel vortex tubes the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ we found resemble the configurations investigated in several earlier studies of extreme flow behaviour (Melander & Hussain Reference Melander and Hussain1989; Kerr Reference Kerr1993, Reference Kerr2013a,Reference Kerrb, Reference Kerr2018). In particular, Orlandi et al. (Reference Orlandi, Pirozzoli and Carnevale2012) used two perpendicular pairs of anti-parallel vortex tubes (referred to as ‘dipoles’ in that study). Although in the case of Navier–Stokes flows none of these studies showed evidence of singularity formation in finite time, they all reported various forms of viscous vortex reconnection. More precisely, while the enstrophy would increase during the reconnection events, it would in all cases remain bounded. Based on three values of the Reynolds number defined as $Re:=({\mathcal{K}}(\boldsymbol{u}_{0}){\mathcal{E}}_{0})^{1/4}/\unicode[STIX]{x1D708}$, the results reported by Orlandi et al. (Reference Orlandi, Pirozzoli and Carnevale2012) exhibit a linear dependence of the maximum attained enstrophy on the Reynolds number, i.e. $\max _{t}{\mathcal{E}}(t)\sim Re$. (We note that this definition of the Reynolds number is different from the one used by Orlandi et al. (Reference Orlandi, Pirozzoli and Carnevale2012).) For comparison, noting that in our case ${\mathcal{K}}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}})\sim {\mathcal{E}}_{0}$ and therefore $Re\sim {\mathcal{E}}_{0}^{1/2}$, the power–law relation (5.2) takes the form $\max _{T>0}{\mathcal{E}}_{T}(\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},T})\sim Re^{3}$ when expressed in terms of the Reynolds number defined above. Thus, even though the range of the Reynolds numbers explored in our study was approximately two orders of magnitude lower than in Orlandi et al. (Reference Orlandi, Pirozzoli and Carnevale2012), in the flows studied here the maximum attained enstrophy grows much more rapidly when the Reynolds number is increased. We add that the ‘trefoil’ initial condition which recently received attention (Kerr Reference Kerr2018) is in fact fundamentally different from the optimal initial data studied here, because it is meant to be defined on an unbounded, rather than periodic, domain. Moreover, unlike the optimal initial conditions $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ and the corresponding flows, it is characterized by non-zero helicity.

Our results indicate that even in the worst-case scenario the nonlinear mechanisms of vortex stretching are significantly depleted, This observation is consistent with the findings of Donzis et al. (Reference Donzis, Gibbon, Gupta, Kerr, Pandit and Vincenzi2013), Gibbon et al. (Reference Gibbon, Donzis, Gupta, Kerr, Pandit and Vincenzi2014) who studied the behaviour of suitably scaled higher-order Lebesgue norms of vorticity in a number of different numerical simulations of turbulent Navier–Stokes flows, both forced and decaying. They provided evidence for a significant depletion of the rate of growth of enstrophy as compared to (2.8). More precisely, they demonstrated that in their simulations enstrophy amplification tends to occur at the rate $\text{d}{\mathcal{E}}/\text{d}t\sim {\mathcal{E}}^{\unicode[STIX]{x1D6FC}}$ with $\unicode[STIX]{x1D6FC}\in [1.575,1.75]$ (we recall from our discussion in § 2 that $\unicode[STIX]{x1D6FC}<2$ implies a regular behaviour of the flow). These findings are also consistent with our observations, cf. figure 7, indicting that while in the flow evolutions corresponding to the optimal initial data $\widetilde{\boldsymbol{u}}_{0;{\mathcal{E}}_{0},\widetilde{T}_{{\mathcal{E}}_{0}}}$ the enstrophy may occasionally be amplified at the much higher rate, the sustained rate is proportional to ${\mathcal{E}}^{\unicode[STIX]{x1D6FC}}$ with $0<\unicode[STIX]{x1D6FC}<2$.

The results summarized above represent a final stage of the research programme outlined in table 1 which aimed to characterize the largest possible growth of enstrophy and enstrophy-like quantities in 1-D Burgers and 2-D and 3-D Navier–Stokes flows. In particular, we sought to understand whether the sharpest estimates on the growth of these quantities obtained using energy-type methods can be realized in actual flows. Somewhat paradoxically, the situation concerning 1-D Burgers flows, where the best finite-time estimate obtained using energy methods was found not to be sharp (Ayala & Protas Reference Ayala and Protas2011), is less satisfactory than in the case of 2-D Navier–Stokes flows where the estimates for both the instantaneous and finite-time growth of palinstrophy were shown to be sharp (and, interestingly, they are realized by the same field (Ayala & Protas Reference Ayala and Protas2014a; Ayala et al. Reference Ayala, Doering and Simon2018)). Moving forward, a natural next step is to consider other conditional regularity results complementary to condition (1.1). More specifically, we will probe the family of the Ladyzhenskaya–Prodi–Serrin conditions asserting that Navier–Stokes flows $\boldsymbol{u}(t)$ are smooth and satisfy system (2.1) in the classical sense provided that (Kiselev & Ladyzhenskaya Reference Kiselev and Ladyzhenskaya1957; Prodi Reference Prodi1959; Serrin Reference Serrin1962)

(6.1)$$\begin{eqnarray}\boldsymbol{u}\in L^{p}([0,T];L^{q}(\unicode[STIX]{x1D6FA})),\quad 2/p+3/q\leqslant 1,q>3.\end{eqnarray}$$

These conditions were recently generalized by Gibbon (Reference Gibbon2018) to include norms of the derivatives of the velocity field. The main technical difficulty in testing conditions (6.1) is that some of the corresponding variational optimization problems will be formulated on function spaces without the Hilbert structure and will therefore require more specialized approaches (Protas Reference Protas2008). Moreover, some of these optimization problems may also be non-smooth. In addition, there are questions concerning potential finite-time singularity formation in inviscid Euler flows which can also be framed in terms of the extreme growth of enstrophy-like quantities. We intend to investigate both these questions in the near future.