  Coherent structures/motions in turbulence inherently give rise to
intermittent signals with sharp peaks, heavy-skirt, and skewed distributions of
velocity increments, highlighting the non-Gaussian nature of turbulence. That
suggests that the spatial nonlocal interactions cannot be ruled out of the
turbulence physics. Furthermore, filtering the Navier-Stokes equations in the
large eddy simulation of turbulent flows would further enhance the existing
nonlocality, emerging in the corresponding subgrid scale fluid motions. That
urges the development of new nonlocal closure models, which respect the
corresponding non-Gaussian statistics of the subgrid stochastic motions. To
this end and starting from the filtered Boltzmann equation, we model the
corresponding equilibrium distribution function with a \textit{L\'evy}-stable
distribution, leading to the proposed fractional-order modeling of
subgrid-scale stresses. We approximate the filtered equilibrium distribution
function with a power-law term, and derive the corresponding filtered
Navier-Stokes equations. Subsequently in our functional modeling, the
divergence of subgrid-scale stresses emerges as a single-parameter fractional
Laplacian, $(-\Delta)^{\alpha}(\cdot)$, $\alpha \in (0,1]$, of the filtered
velocity field. The only model parameter, i.e., the fractional exponent,
appears to be strictly depending on the filter-width and the flow Reynolds
number. We furthermore explore the main physical and mathematical properties of
the proposed model under a set of mild conditions. Finally, the introduced
model undergoes \textit{a priori} evaluations based on the direct numerical
simulation database of forced and decaying homogeneous isotropic turbulent
flows at relatively high and moderate Reynolds numbers, respectively. Such
analysis provides a comparative study of predictability and performance of the
proposed fractional model.
