Modeling of wall-bounded turbulent flows is still an open problem in classical physics, with relatively slow progress in the last few decades beyond the log law, which only describes the intermediate region in wall-bounded turbulence, i.e., 30–50 y+ to 0.1–0.2 R+ in a pipe of radius R. Here, we propose a fundamentally new approach based on fractional calculus to model the entire mean velocity profile from the wall to the centerline of the pipe. Specifically, we represent the Reynolds stresses with a non-local fractional derivative of variable-order that decays with the distance from the wall. Surprisingly, we find that this variable fractional order has a universal form for all Reynolds numbers and for three different flow types, i.e., channel flow, Couette flow, and pipe flow. We first use existing databases from direct numerical simulations (DNSs) to lean the variable-order function and subsequently we test it against other DNS data and experimental measurements, including the Princeton superpipe experiments. Taken together, our findings reveal the continuous change in rate of turbulent diffusion from the wall as well as the strong nonlocality of turbulent interactions that intensify away from the wall. Moreover, we propose alternative formulations, including a divergence variable fractional (two-sided) model for turbulent flows. The total shear stress is represented by a two-sided symmetric variable fractional derivative. The numerical results show that this formulation can lead to smooth fractional-order profiles in the whole domain. This new model improves the one-sided model, which is considered in the half domain (wall to centerline) only. We use a finite difference method for solving the inverse problem, but we also introduce the fractional physics-informed neural network (fPINN) for solving the inverse and forward problems much more efficiently. In addition to the aforementioned fully-developed flows, we model turbulent boundary layers and discuss how the streamwise variation affects the universal curve.