We study decaying turbulence in the one-dimensional (1-D) Burgers equation (Burgulence) and 3-D Navier–Stokes (NS) turbulence. We first investigate the decay in time
$t$ of the energy
$E(t)$ in Burgulence, for a fractional Brownian initial potential, with Hurst exponent
$H$, and demonstrate rigorously a self-similar time decay of
$E(t)$, previously determined heuristically. This is a consequence of the non-trivial boundedness of the energy for any positive time. We define a spatially forgetful oblivious fractional Brownian motion (OFBM), with Hurst exponent
$H$, and prove that Burgulence, with an OFBM as initial potential
$\varphi _0(x)$, is not only intermittent, but it also displays a, hitherto unanticipated, large-scale bifractality or multifractality; the latter occurs if we combine OFBMs, with a distribution of
$H\hbox{-}$values. This is the first rigorous proof of genuine multifractality for turbulence in a nonlinear hydrodynamical partial differential equation. We then present direct numerical simulations (DNSs) of freely decaying turbulence, capturing some aspects of this multifractality. For Burgulence, we investigate such decay for two cases: (a)
$\varphi _0(x)$ a multifractal random walk that crosses over to a fractional Brownian motion beyond a cross-over scale
$\mathcal{L}$, tuned to go from small- to large-scale multifractality; (b) initial energy spectra
$E_0(k)$, with wavenumber
$k$, having one or more power-law regions, which lead, respectively, to self-similar and non-self-similar energy decay. Our analogous DNSs of the 3-D NS equations also uncover self-similar and non-self-similar energy decay. Challenges confronting the detection of genuine large-scale multifractality, in numerical and experimental studies of NS and Magnetohydrodynamics turbulence, are highlighted.