Fourier Neural Operators (FNOs) have demonstrated exceptional accuracy in mapping functional spaces by leveraging Fourier transforms to establish a connection with underlying physical principles. However, their opaque inner workings often constitute an obstacle to physical interpretability. This work introduces Hierarchical Fourier Neural Operators (HFNOs), a novel FNO-based architecture tailored for reduced-order modeling of turbulent fluid flows, designed to enhance interpretability by explicitly separating fluid behavior across scales. The proposed architecture processes wavenumber bins in parallel, enabling the approximation of dispersion relations and non-linear interactions. Inputs are lifted to a higher-dimensional space, Fourier-transformed, and partitioned into wavenumber bins. Each bin is processed by a Fully Connected Neural Network (FCNN), with outputs subsequently padded, summed, and inverse-transformed back into physical space. A final transformation refines the output in physical space as a correction model, by means of one of the following architectures: Convolutional Neural Network (CNN) and Echo State Network (ESN). We evaluate the proposed model on a series of increasingly complex dynamical systems: first on the one-dimensional Kuramoto-Sivashinsky equation, then on the two-dimensional Kolmogorov flow, and finally on the prediction of wall shear stress in turbulent channel flow, given the near-wall velocity field. In all test cases, the model demonstrates its ability to decompose turbulent flows across various scales, opening up the possibility of increased interpretability and multiscale modeling of such flows.