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.