It　is　now　well　known　that　the　magnetic　resonance　imaging　(MRI)　signal　decay　deviates　from　the　classical　mono-exponential　relaxation.　This　deviation　is　referred　to　in　the　literature　as　anomalous　relaxation.　The　modelling　of　this　anomalous　relaxation　can　provide　a　better　understanding　of　MRI　magnetization.　The　purpose　of　this　work　is　to　investigate　the　utility　of　the　distributed-order　time　fractional　Bloch　equations　to　describe　anomalous　relaxation　processes　in　human　brain　MRI　data.　Two　choices　of　continuous　distribution　weight　functions,　which　are　parameterised　by　their　mean　μ　and　standard　deviation　σ,　are　studied　to　investigate　their　impact　on　the　model　solution　behaviour.　An　implicit　numerical　method　implemented　on　a　graded　mesh　is　proposed　to　solve　the　model　and　the　stability　and　convergence　analysis　are　presented.　We　also　derive　semi-analytical　solutions　of　the　fully　coupled　Bloch　equations　using　the　Laplace　transform　technique　to　assess　the　accuracy　of　the　numerical　scheme.　Furthermore,　three　different　voxel　models　of　continuous　distribution　weight　functions,　namely　a　single　continuous　probability　distribution　(model　1),　two　distinct　continuous　probability　distributions　(model　2)　and　a　mixture　of　two　continuous　probability　distributions　(model　3),　are　applied　to　the　in　vivo　human　brain　MRI　data,　and　a　feasible　and　reliable　parameter　estimation　method　based　on　a　modified　hybrid　Nelder-Mead　simplex　search　and　particle　swarm　optimization　is　presented　to　perform　the　voxel-level　temporal　fitting　of　the　MRI　data.　The　application　of　these　distributed-order　time　fractional　Bloch　models　highlights　the　validity　of　the　proposed　models,　and　based　on　the　mean　square　error　we　conclude　that　models　2　and　3　might　be　more　suitable　than　model　1　to　characterize　anomalous　relaxation　processes　in　human　brain　MRI　data.　©　2022　Elsevier　Inc.