3.1. Laminar boundary layer

The> first validation example considers the well-known case of the two-dimensional, laminar and incompressible boundary layer, where Prandtl (Reference Prandtl1904) showed that some terms of the Navier–Stokes equations are negligible, and reduced the latter to the ‘boundary layer equations’ (BLE). By utilising the streamfunction $\varPsi$ , the BLE can be expressed in a single equation. Blasius (Reference Blasius1908) noted that, for semi-infinite flat plate boundary conditions, the BLE is symmetrical under a dilational transformation, which reduces the problem to an ordinary differential equation for the non-dimensional streamfunction $f(\tilde {y})=\varPsi /\sqrt {\nu U_\infty x}$ , with $\tilde {y} = y \sqrt {U_\infty /(\nu x)}$ . In the above, $U_\infty$ is the free-stream velocity away from the flat plate, $x$ and $y$ are the streamwise and normal-to-the-wall distances, respectively, and $\nu$ is the kinematic viscosity of the fluid. The BLE are thus reduced to the boundary value problem

(3.1a) \begin{align}&\qquad\quad 2 f_{\tilde {y} \tilde {y} \tilde {y}} + \textit{ff}_{\tilde {y} \tilde {y}} =0, \end{align}

(3.1b) \begin{align}& f(0) = 0, f_{\tilde {y}}(0) = 0, f_{\tilde {y}}(\infty ) = 1. \\[12pt] \nonumber \end{align}

Figure 1. Data-driven identification of self-similarity in the Blasius boundary layer. (a) Streamwise velocity field (Blasius solution). The dashed vertical lines denote the nine velocity profile sampling stations. (b) Sampled velocity profiles. (c) Algorithmically collapsed velocity profiles. (d) Algorithmically identified scaling $\alpha$ of the wall-normal coordinate $y$ . (e) Algorithmically identified scaling $\beta$ of the streamwise velocity $u$ .

Knowledge of $\tilde {u} = f_{\tilde {y}}$ (e.g. via numerical integration of the Blasius equation) allows for the calculation of the streamwise velocity distribution $u = U_\infty \tilde {u}$ at any location in the boundary layer. Here, we derive the Blasius similarity directly from data, without prior knowledge of the Navier–Stokes equations, Prandtl’s scaling analysis or Blasius’ similarity arguments. Our dataset consists of $u(x,y)$ data from a flat plate laminar boundary layer flow (figure 1 a), acquired by solving the Blasius ODE with a shooting method (Bakarji et al. Reference Bakarji, Callaham, Brunton and Kutz2022). In particular, we extract velocity profiles $u(y)$ from nine stations at different streamwise coordinates $x$ (figure 1 b). The algorithm commences by seeking similarity under the dilational transformations $\tilde {y} = \alpha (x) y$ , $\tilde {u} = \beta (x) u$ . Figure 1(c) shows that the transformed profiles have collapsed onto a single curve. The identified transformations for the wall-normal coordinate $y$ and the streamwise velocity $u$ are shown in figures 1(d) and 1(e), respectively.

The second step consists of expressing the identified transformations as functions of the governing parameters and independent variables, i.e. a symbolic regression task. The library of candidate variables is composed of the free-stream velocity $U_\infty$ , the viscosity $\nu$ and the streamwise coordinate of the extracted profiles $x$ , i.e. $\alpha = \alpha (U_\infty ,\nu ,x)$ and $\beta =\beta (U_\infty ,\nu ,x)$ . Owing to the nature of the assumed dilational transformations, we look for expressions in the form of monomials, i.e. $\alpha =c_1 U_\infty ^{c_2} \nu ^{c_3} x^{c_4}$ and $\beta =c_5 U_\infty ^{c_6} \nu ^{c_7} x^{c_8}$ , where $c_1$ and $c_5$ are multiplicative constants that reflect the arbitrariness of the dilation factors between the original and similarity variables. Dimensional homogeneity is enforced by assuming $\tilde {y}$ and $\tilde {u}$ to be dimensionless and a non-zero weighting constant $w_D$ in the objective function (2.4). The interpreted expressions are

(3.2a) \begin{align}&\, \tilde {y} \propto y U_\infty ^{0.5018} \nu ^{-0.5018} x^{-0.4981} ,\end{align}

(3.2b) \begin{align}& \tilde {u} \propto u U_\infty ^{-0.9995} \nu ^{-0.0005} x^{0.0005} ,\end{align}

which are in close agreement with the theoretically derived similarity variables of Blasius ( $\tilde {y} = y \sqrt {U_\infty /(\nu x)}$ , $\tilde {u} = u / U_\infty$ ). As previously noted, the identified and theoretical scalings can vary by an arbitrary multiplicative constant (i.e. the offset in figure 1 d), depending as to whether and how the user normalises the collapsed profiles. Appendix B shows results from a sensitivity and robustness analysis with respect to the input data, parameters and noise.

3.2. Burgers’ equation

The second validation example considers Burgers’ equation, which is a partial differential equation that finds wide application in fluid dynamics, nonlinear acoustics and traffic flow, among others. It acts as a prototype for a variety of phenomena including shock wave formation, rarefaction waves and turbulence (Zel’Dovich & Raizer Reference Zel’Dovich and Raizer1967; Whitham Reference Whitham2011). In the absence of diffusion, it is the simplest model for gas dynamics. Its mathematical expression in one dimension, along with its initial condition is

(3.3) \begin{equation} u_t + u u_x = 0, \quad u(x,0) = f(x) ,\end{equation}

where $u(x,t)$ is the fluid velocity and $x$ and $t$ the spatial and temporal coordinates, respectively. The subscript denotes partial differentiation. Using the method of characteristics, the solution can be expressed as $u(x,t) = f(\xi ) = f(x-ut) = f(x+\alpha )$ , i.e. the problem described by (3.3) is invariant under a non-uniform translational transformation. In this work, we consider a sinusoidal initial condition, $u(x,0) = \sin {(2 \pi x)}$ , with $x \in [0,1 ]$ . Figures 2(a) and 2(b) show the spatio-temporal evolution of velocity and profiles extracted at seven different time instants, respectively, illustrating the steepening of the solution as time evolves. The times shown are non-dimensionalised with the time of shock formation $t_c$ , $t^*=t/t_c=2 \pi t$ . The algorithm is provided with the data shown with markers in figure 2(b).

Figure 2. Data-driven identification of self-similarity in Burgers’ equation. (a) Spatio-temporal evolution of the velocity. (b) Extracted profiles at different time instants. Markers show the data given to the algorithm. (c) Algorithmically collapsed profiles.

Figure 2(c) shows the algorithmically collapsed velocity profiles on the transformed coordinates $ (\tilde {x},t )$ , i.e. the output of step 1 of the algorithm, following a transformation of the form $\tilde {x} = x+ \alpha (x,t)$ (simpler transformations do not succeed at collapsing the data). Since this is an initial value problem, we only consider the distance of the profiles with respect to the initial condition, i.e. the objective function in (2.2) is simplified to $\sum _{i=1,\ldots ,n_t} \| u(\tilde {x},t_i) - u(x,0) \|_2^2$ .

We use PySR to interpret the algorithmically computed transformation $\alpha$ . The library of candidate variables is composed of the problem’s independent and dependent variables, i.e. $\alpha = \alpha (u,x,t)$ , and the library of operators includes the four basic mathematical operations ( $+,-,\times ,\div$ ). $\texttt {PySR}$ interprets the algorithmically computed transformation as $\alpha = - u t$ , which is identical to the analytical similarity transformation. Close inspection of figure 2(c) shows that the profiles are not perfectly collapsed (for example, at $t^*=0.936$ ). The range-normalised mean absolute error between the algorithmically computed discrete values of the transformation and the analytical ones is $2.2\,\%$ . However, in step 2 of the workflow (symbolic regression), the algorithm attempts to fit expressions that balance complexity and accuracy, allowing it to recover the exact analytical transformation $\alpha = -ut$ .

3.3. Free turbulent flow

As a third example, we consider the statistically stationary flow past a slender bluff body at high Reynolds numbers, which is a characteristic example of a free (i.e. unconfined) turbulent shear flow. The term slender signifies an effectively infinite body length in one cross-flow direction, along which the flow statistics can be considered homogeneous. Such flows can be modelled using a boundary layer approximation, similar to the one employed by Prandtl for laminar boundary layers, applied to the time-averaged flow (Townsend Reference Townsend1976), i.e.

(3.4a) \begin{align}&\qquad\quad \bar {u}_x + \bar {v}_y = 0 , \end{align}

(3.4b) \begin{align}& \bar {u} \bar {u}_x + \bar {v} \bar {u}_y = \nu \bar {u}_{\textit{yy}} - (\overline {u'v'})_y , \end{align}

where the overbar denotes time averaging and $\nu$ is the kinematic viscosity of the fluid. Here, $u$ and $v$ are the flow velocities along the streamwise $x$ and cross-wise $y$ directions, which are decomposed into time-averaged $(\bar {u},\bar {v})$ and fluctuating $(u',v')$ components. The formulation leading to (3.4) introduces the Reynolds stress term $\overline {u'v'}$ into the flow governing equations, i.e. an additional unknown which cannot be calculated implicitly (this is the well-known closure problem of turbulence). The identification of similarity in this case is, therefore, not possible from the equations themselves, unless an assumption is made regarding the relation of the unknown Reynolds stress terms to the mean velocity distribution (i.e. turbulence modelling). However, it is customary to hypothesise that the turbulent BLE accept self-similar solutions far from initial conditions, a fact also supported by experimental observations (Townsend Reference Townsend1976; Cantwell Reference Cantwell2002). This approximate self-similarity is the origin behind various scaling laws for the evolution of free shear flows, widely used in many applications of the energy and transportation sectors, such as wind farm planning (Bempedelis & Steiros Reference Bempedelis and Steiros2022) and jet engine noise prediction (Tam Reference Tam2019). Returning to the particular example of the slender bluff body wake, we may consider the time-averaged velocity deficit $\zeta (x,y) = U_\infty - \bar {u}(x,y)$ that the body produces far (i.e. tens of characteristic body lengths) downstream. In the above, $U_\infty$ is the constant free-stream velocity. The flow statistics are homogeneous along the spanwise direction $z$ . Wake self-similarity is then known to assume the form (Townsend Reference Townsend1976; Pope Reference Pope2000)

(3.5) \begin{equation} \zeta (x,y) = \tilde {\zeta }(\xi ) (U_\infty -\bar {u}_{\textit{cntr}}(x)), \quad \text{with } \xi = \frac {y}{y_{1/2}(x)}\,, \end{equation}

Figure 3. Data-driven identification of self-similarity in the wake of a porous plate. (a) Schematic of the experimental apparatus showing the flume, porous plate and PIV configuration. (b) Mean streamwise velocity profiles at different locations downstream of the plate. (c) Algorithmically collapsed velocity profiles.

where $\bar {u}_{\textit{cntr}}=\bar {u}(x,0)$ denotes the centreline velocity and $y_{1/2}(x)$ is the wake half-width defined such that $\bar {u}(x,\pm y_{1/2}) = 0.5 (U_\infty +u_{\textit{cntr}} )$ . We attempt to extract the self-similar relation (3.5) directly from experimental data. To this end, we measured the turbulent wake of a plate of 53 % porosity immersed in a water flume normal to the flow (see figure 3 a) at a Reynolds number based on the free-stream velocity and plate width $Re \approx 6000$ . The velocity fields at various positions downstream of the plate were measured using Particle Image Velocimetry (PIV) (Phantom 4MP camera at 50 Hz acquisition frequency) as shown in figure 3(a). Each mean streamwise velocity profile was then calculated by averaging 3000 vector fields. More information regarding the experimental procedure can be found in Bekoglu et al. (Reference Bekoglu, Bempedelis and Steiros2024, Reference Bekoglu, Bempedelis and Steiros2025). We consider the mean velocities at five stations in the wake of the plate (figure 3 b). Figure 3(c) shows the collapse that is obtained via the proposed method, assuming transformations of the form $\tilde {y} = \alpha (x) y$ and $\tilde {u} = \beta (x) \bar {u} + \gamma (x)$ (simpler transformations cannot collapse the data). The identified transformations are regressed as functions of the characteristic scales and variables of the problem using $\texttt {PySR}$ , with $\alpha = \alpha (x, y_{w}, y_{1/2})$ (with $y_w(x)$ defined as $\bar {u}(x,\pm y_{w}) = 0.99 U_\infty$ ), $\beta =\beta (U_\infty , \bar {u}_{\textit{cntr}}, x, \nu )$ and $\gamma =\gamma (U_\infty , \bar {u}_{\textit{cntr}}, x, \nu )$ . The library of operators consists of the four basic mathematical operations ( $+,-,\times ,\div$ ). The interpreted similarity transformations are

(3.6a) \begin{align}&\qquad\qquad\qquad\quad \tilde {y}= \alpha y = \frac {1.0605}{y_{1/2}(x)} y , \end{align}

(3.6b) \begin{align}& \tilde {u}(x,\tilde {y}) = \beta \bar {u} + \gamma = \frac {1}{U_\infty - \bar {u}_{\textit{cntr}}} \bar {u} - \frac {\bar {u}_{\textit{cntr}}}{U_\infty - \bar {u}_{\textit{cntr}}} = 1 - \tilde {\zeta } , \end{align}

which match (ignoring the arbitrary multiplicative constant) the self-similar expressions for turbulent wakes found in the literature (Tennekes & Lumley Reference Tennekes and Lumley1972; Townsend Reference Townsend1976; Pope Reference Pope2000). Besides our own experiments, we also test our method by using experimental data available in the literature (Cimbala, Nagib & Roshko Reference Cimbala, Nagib and Roshko1988), and again retrieve the self-similarity relations for turbulent wakes (see Appendix C).

3.4. Cavity collapse

The bursting of bubbles at the sea’s surface plays a crucial role in the exchanges between oceans and the atmosphere, thereby being important for climate and weather (Deike Reference Deike2022). When a bubble bursts at the sea’s surface, aerosol is produced via two mechanisms: the rupture of the bubble’s cap film, and the formation of a vertical jet that breaks into droplets following cavity collapse (Deike Reference Deike2022). The jet and aerosol properties are a consequence of the nonlinear fluid dynamics near the points where the topology changes (Eggers Reference Eggers1997; Eggers & Fontelos Reference Eggers and Fontelos2015). In the case of a collapsing cavity, several studies have found that, near the collapse of the travelling capillary waves at the axis of symmetry, the liquid–gas interface evolves in a self-similar manner, enabling the derivation of scaling laws for the bubble and jet dynamics (Keller & Miksis Reference Keller and Miksis1983; Zeff et al. Reference Zeff, Kleber, Fineberg and Lathrop2000; Duchemin et al. Reference Duchemin, Popinet, Josserand and Zaleski2002; Ghabache et al. Reference Ghabache, Antkowiak, Josserand and Séon2014; Gañán-Calvo Reference Gañán-Calvo2017; Deike et al. Reference Deike, Ghabache, Liger-Belair, Das, Zaleski, Popinet and Séon2018).

We simulate the collapse of a cavity ( ${Bo}=\rho g R_0^2/\gamma = 10^{-3}$ , ${{La}}=\rho \gamma R_0/ \mu ^2 = 2500$ ) by solving the two-phase incompressible axisymmetric Navier–Stokes equations using Basilisk (www.basilisk.fr), an open-source flow solver that has been previously used in computational studies of bursting bubbles (Lai, Eggers & Deike Reference Lai, Eggers and Deike2018; Berny et al. Reference Berny, Deike, Séon and Popinet2020; Sanjay, Lohse & Jalaal Reference Sanjay, Lohse and Jalaal2021). In the above, $\rho$ and $\mu$ are the liquid density and viscosity, respectively, $\gamma$ is the interfacial tension, $R_0$ is the initial bubble radius and $g$ is the gravitational acceleration. The evolution of the liquid–gas interface is shown in figures 4(a)–4(d). We extract profiles of the interface $h(r,t)$ as it approaches the axis of symmetry, at $t_* = (t_0 - t )/t_c = [10, 9, 8, 7, 6 ]$ , where $t_c$ is the characteristic time of the horizontal capillary wave (Lai et al. Reference Lai, Eggers and Deike2018), and $t_0$ is the moment where the capillary waves meet at the axis of symmetry, $r=0$ . The extracted profiles are shown in figure 4(e). Duchemin et al. (Reference Duchemin, Popinet, Josserand and Zaleski2002) and Lai et al. (Reference Lai, Eggers and Deike2018) observed that, in this time window, the interface approximately follows the $ (t_0-t )^{2/3}$ scaling of inviscid theory (Keller & Miksis Reference Keller and Miksis1983). The algorithm successfully collapses the extracted profiles into a single curve (figure 4 f) for transformations of the form $\tilde {r} = \alpha (t) r$ for the radial coordinate and $\tilde {h} = \beta (t) h + \gamma (t)$ for the liquid–gas interface. By regressing the identified transformations $\alpha$ , $\beta$ and $\gamma$ as power-law functions of the non-dimensional time $t_*$ (figure 4 g–4i), we retrieve the expressions

(3.7a) \begin{align}&\qquad\,\, \tilde {r} = c_1r t_*^{-0.61} ,\end{align}

(3.7b) \begin{align}& \tilde {h} = c_2 h t_*^{-0.69} - c_3 h_0 t_*^{-0.66} , \end{align}

which are close (ignoring the arbitrary multiplicative constants $c_1$ , $c_2$ and $c_3$ ) to the theoretically derived ones (Keller & Miksis Reference Keller and Miksis1983; Zeff et al. Reference Zeff, Kleber, Fineberg and Lathrop2000; Duchemin et al. Reference Duchemin, Popinet, Josserand and Zaleski2002; Ghabache et al. Reference Ghabache, Antkowiak, Josserand and Séon2014; Lai et al. Reference Lai, Eggers and Deike2018).

Figure 4. Data-driven identification of self-similarity in a collapsing cavity. Three-dimensional visualisation of the liquid–gas interface at (a) $t/\tau =0$ , (b) $t/\tau =0.25$ and (c) $t/\tau =0.5$ , where $\tau =\sqrt {\rho R_0^3 / \gamma }$ is the inertio-capillary time scale. (d) Liquid–gas interface time evolution, $t/\tau = (0, 0.05, \ldots , 0.5 )$ . (e) Interface profiles near cavity collapse. ( f) Algorithmically collapsed interface profiles near cavity collapse. (g–i) Identified transformations $\alpha$ , $\beta$ and $\gamma$ . Comparison with theoretical scaling laws.

3.5. Decaying turbulence

We consider the case of homogeneous decaying turbulence, experimentally realised by passing a stream of fluid through a uniformly spaced grid inside a wind tunnel (figure 5 a). We analyse measurements of velocity time series, obtained via hot-wire anemometry at seven positions downstream of the grid (the details of the experimental campaign can be found in Steiros Reference Steiros2022a ). The measured turbulence is fully developed, approximately homogeneous and yields the −5/3 law for the energy spectrum at intermediately sized eddies, as predicted by Kolmogorov’s K41 framework (Kolmogorov Reference Kolmogorov1941a , Reference Kolmogorovb , Reference Kolmogorovc ).

Figure 5. Data-driven identification of self-similarity in decaying turbulence. (a) Schematic of the experimental set-up. (b) Experimentally measured power spectral densities. The dashed vertical line delineates the range of the spectrum that is used. (c) Measured spectrum normalised by the inertial scales. (d) Measured spectrum normalised by the Kolmogorov scales. (e) Measured spectrum normalised by the algorithmically identified expression (3.9).

The theoretical derivation of the −5/3 law assumes that, at sufficiently high Reynolds numbers, an intermediate self-similar region forms in the cascade, which is independent of large- and small-scale effects (Batchelor Reference Batchelor1953; Frisch Reference Frisch1995; Pope Reference Pope2000). However, the above description is part of a more complex flow picture, as the widely studied paradigm of grid turbulence sufficiently far from initial conditions shows. The turbulence cascade in that case is generally accepted to be characterised by two self-similarities, both present at the same time and at different eddy sizes: One at large scales where viscosity is negligible (Steiros Reference Steiros2022a ; Lundgren Reference Lundgren2003), and one at small scales where non-equilibrium effects are negligible (Pope Reference Pope2000; Lundgren Reference Lundgren2003). The energy spectrum thus accepts the following general expression (Pope Reference Pope2000):

(3.8) \begin{equation} E_{11}(\kappa ,x) = \epsilon (x)^{2/3} \kappa ^{-5/3} f(\kappa L) g(\kappa \eta )\,, \end{equation}

where $x$ , $\kappa$ , $E_{11}(\kappa ,x)$ and $\epsilon (x)$ represent the distance from the grid, wavenumber, energy spectrum and dissipation rate, respectively. Here, $L$ and $\eta$ are the integral and Kolmogorov scales, characteristic of the large- and small-scale self-similarities, respectively. To retrieve the −5/3 law, $E_{11}(k,x) = \epsilon (x)^{2/3} \kappa ^{-5/3}$ , one needs to be asymptotically far from large scales (i.e. $\kappa L \to \infty$ ) and far from small scales ( $\kappa \eta \to 0$ ) at the same time, as in that case K41 assures that both $f$ and $g$ tend to unity. It is of interest to see how our algorithm fares in this two-scale problem.

Figure 5(b) plots the measured power spectral densities versus the wavenumbers at different measurement stations. Figure 5(c,d) shows the collapse of the power spectral densities using large- and small-scale similarity variables, respectively. Using a library that consists of $(L, \eta )$ for the $x$ axis and $(\epsilon , L, \eta , k)$ for the $y$ axis, symbolic regression of the algorithmically identified dilational transformations ( $\tilde {\kappa } = \alpha (x) \kappa$ , $\tilde {E_{11}} = \beta (x) E_{11}$ ) yields the following expressions, plotted in figure 5(e):

(3.9) \begin{equation} \tilde {\kappa } \propto \kappa L^{0.368} \eta ^{0.632}, \qquad \tilde {E_{11}} \propto E_{11} \epsilon ^{0.695} L^{0.709} \eta ^{-1.015} k^{-2.042}. \end{equation}

As expected, the algorithm cannot derive the two-scale similarity expression (3.8), as the assumed dilational transformations effectively reduce it to a single-scale algorithm. Further insight may be obtained by using the dissipation scaling $\epsilon \propto k^{3/2}/L$ , where $k$ is the turbulence kinetic energy, which is known to characterise homogeneous decaying turbulence far from initial conditions (Steiros Reference Steiros2022a , Reference Steirosb ; Vassilicos Reference Vassilicos2015). Expression (3.9) then becomes

(3.10) \begin{equation} \tilde {\kappa } \propto \kappa L^{0.368} \eta ^{0.632}, \qquad \tilde {E_{11}} \propto E_{11} \epsilon ^{-0.667} \big(L^{0.391} \eta ^{0.609} \big)^{-5/3}. \end{equation}

Inspection of the above expression reveals that our algorithm has, seemingly, collapsed the totality of the cascade dynamics based on a single empirical length scale, $l \approx L^{0.38} \eta ^{0.62}$ , which is a combination of $L$ and $\eta$ . There have been several attempts to develop single length scale theories of turbulence, derived by assuming self-similarity of the turbulence statistics, i.e. see von Karman & Howarth (Reference von Karman and Howarth1938) and Sedov (Reference Sedov2018), and the later theories of Barenblatt & Gavrilov (Reference Barenblatt and Gavrilov1974) and George (Reference George1992). The last two studies predict the Taylor microscale as the appropriate similarity variable, in close agreement with our data-driven methodology. In homogeneous decaying turbulence the Taylor microscale becomes $\lambda = L^{1/3} \eta ^{2/3}$ (Pope Reference Pope2000), which is close to our empirically derived length scale.

It might be tempting to interpret our empirical results as being in support of a single-scale picture for homogeneous turbulence, contrary to the commonly accepted view of turbulence as an intrinsically multi-length scale phenomenon (Batchelor Reference Batchelor1953; Tennekes & Lumley Reference Tennekes and Lumley1972; Townsend Reference Townsend1976; Pope Reference Pope2000). However, the results shown in figure 5 (and in most studies of homogeneous decaying turbulence) consider the far region of a turbulence grid where the integral length Reynolds number evolution is relatively slow, leading to a similarly slow evolution of $L/\eta$ (Vassilicos Reference Vassilicos2015). As a result, our data show that all length scales in question (integral, Kolmogorov, Taylor) produce a (visually) ‘good-enough’ collapse, as shown in figure 5(c–e), and inspection of any of these graphs might give the impression of a single length scale self-similar evolution.

To better appreciate the single or multi-length scale nature of the turbulence cascade, we need data for which the integral length scale Reynolds number varies drastically. To achieve this, we analyse the periodic box Direct Numerical Simulation (DNS) data of decaying turbulence taken from Goto & Vassilicos (Reference Goto and Vassilicos2016) (the details of the simulations can be found in the reference). Similar to the grid case, our goal is to collapse energy spectra taken at different decay times $\hat {t}$ (see figure 6 a). Here, however, we consider datasets from two simulations which are identical in every aspect, apart from their size ( $N=1024^3$ and $N=2048^3$ ) and the kinematic viscosity of the fluid, i.e. the Reynolds number of the flow varies drastically between the two cases, making the collapse of the curves less arbitrary.

Figure 6(a) presents the dimensional energy spectra from the two simulations, whereas figures 6(b) and 6(c) show the same data with the axes rescaled according to inertial and Kolmogorov similarity variables, respectively. The collapse of the large (small) scales is adequate only when the axes are normalised using the integral (Kolmogorov) scales, in agreement with a multi-scale view of the cascade. However, data from the same simulation size retain a ‘good-enough’ collapse independent of the choice of normalisation length, similar to our grid measurements, a fact which, if viewed in isolation, might give a false impression of a single length scale process. Figure 6(d) shows the algorithmically identified collapse of the data, that is, using the empirical transformations

(3.11a) \begin{align}&\qquad\qquad\qquad\qquad\qquad\quad \tilde {\kappa } \propto \kappa L^{0.342} \eta ^{0.658} , \end{align}

(3.11b) \begin{align}& \tilde {E_{11}} \propto E_{11} \epsilon ^{-0.655} L^{-0.561} \eta ^{-1.094} k^{-0.017} \approx E_{11} \epsilon ^{-0.667} \big(L^{0.344} \eta ^{0.656} \big)^{-5/3}. \end{align}

The algorithm again identifies, approximately, the Taylor microscale ( $\lambda = L^{1/3}\eta ^{2/3})$ as the appropriate length for the collapse of the turbulence dynamics, but this time it cannot be claimed that the Taylor scale collapses both the large- and small-scale dynamics – only that it is an (algorithmically) optimal compromise when a single length scale collapse is enforced. We conjecture that the reason behind this could be explained from Lundgren’s two-scale turbulence theory which is based on the technique of matched asymptotic expansions (Lundgren Reference Lundgren2002, Reference Lundgren2003). The repercussions of this theory are discussed in Obligado & Vassilicos (Reference Obligado and Vassilicos2019), Meldi & Vassilicos (Reference Meldi and Vassilicos2021): in particular, it is suggested that the Taylor microscale is the length scale at which the cascade dynamics is optimally distanced from both non-equilibrium effects (large scales) and dissipative effects (small scales) at the same time (in a matched asymptotic manner). The algorithmically identified collapse (figure 6 d) produces best results in the intermediate, inertial range of the cascade where both above effects tend to become negligible. Based on this argument, the algorithm identifies the Taylor microscale scale because it is the optimal compromise of the two governing scales (inertial, Kolmogorov) of the cascade.

Figure 6. Data-driven identification of self-similarity in decaying turbulence, using the DNS data of Goto & Vassilicos (Reference Goto and Vassilicos2016). (a) Power spectral densities. (b) Spectrum normalised by the inertial scales. (c) Spectrum normalised by the Kolmogorov scales. (d) Spectrum normalised by the algorithmically identified expression (3.11).