<p>The main objective of this article is to investigate approximate solutions for the two-dimensional space/multi-time fractional Bloch-Torrey model involving the Riesz fractional operator. To achieve this, an effective hybrid numerical approach is developed and employed, enabling accurate and efficient approximation of the model’s complex fractional dynamics. A numerical approach based on the L1-2 approximation is employed for the discretization of the Caputo time-fractional operator with multiple orders. The classic Galerkin finite element approach is employed to approximate the Riesz fractional operator in the spatial directions. The convergence order for the proposed semi-discrete numerical scheme is <InlineEquation ID="IEq1"> <EquationSource Format="TEX">\(\big (\big (\Delta t\big )^{3-\sigma }\big )\)</EquationSource> <EquationSource Format="MATHML"><math> <mrow> <mrow> <mo maxsize="1.2em" minsize="1.2em" stretchy="true">(</mo> </mrow> <mrow> <mo maxsize="1.2em" minsize="1.2em" stretchy="true">(</mo> </mrow> <mi mathvariant="normal">Δ</mi> <mi>t</mi> <msup> <mrow> <mo maxsize="1.2em" minsize="1.2em" stretchy="true">)</mo> </mrow> <mrow> <mn>3</mn> <mo>-</mo> <mi>σ</mi> </mrow> </msup> <mrow> <mo maxsize="1.2em" minsize="1.2em" stretchy="true">)</mo> </mrow> </mrow> </math></EquationSource> </InlineEquation>. In this article, we provide a detailed error and stability analysis of the developed numerical method. To validate the effectiveness and accuracy of the fully discrete scheme, several numerical examples are presented. The corresponding results are illustrated through figures and tables, highlighting the performance of the method under various conditions.</p>