<p>This paper examines the design and study of a computational technique to address a class of 2D nonlinear time-fractional variable-order advection-reaction-diffusion equations. In the temporal direction, the Caputo fractional variable-order derivative is discretized as a linear B-spline basis function. The spatial variables are then discretized and analyzed using a modified Bi-cubic B-spline basis methodology on a piecewise uniform mesh. It is shown that the resultant discrete scheme exhibits unconditional stability and convergence with an order of convergence <InlineEquation ID="IEq1"> <EquationSource Format="TEX">\((\Delta \varsigma ^{2- \varrho (\bar{\varvec{x}},\varsigma )}+h^{2})\)</EquationSource> <EquationSource Format="MATHML"><math> <mrow> <mo stretchy="false">(</mo> <mi mathvariant="normal">Δ</mi> <msup> <mi>ς</mi> <mrow> <mn>2</mn> <mo>-</mo> <mi>ϱ</mi> <mo stretchy="false">(</mo> <mover accent="true"> <mrow> <mrow> <mi mathvariant="bold-italic">x</mi> </mrow> </mrow> <mrow> <mo stretchy="false">¯</mo> </mrow> </mover> <mo>,</mo> <mi>ς</mi> <mo stretchy="false">)</mo> </mrow> </msup> <mo>+</mo> <msup> <mi>h</mi> <mn>2</mn> </msup> <mo stretchy="false">)</mo> </mrow> </math></EquationSource> </InlineEquation>, where <i>h</i> is the maximum value of <InlineEquation ID="IEq2"> <EquationSource Format="TEX">\((h_{\texttt {x}},h_{\texttt {y}})\)</EquationSource> <EquationSource Format="MATHML"><math> <mrow> <mo stretchy="false">(</mo> <msub> <mi>h</mi> <mi mathvariant="monospace">x</mi> </msub> <mo>,</mo> <msub> <mi>h</mi> <mi mathvariant="monospace">y</mi> </msub> <mo stretchy="false">)</mo> </mrow> </math></EquationSource> </InlineEquation>. In order to confirm the theoretical conclusions and illustrate the efficiency of the approach, certain examples that have been solved are presented. In summary, the numerical findings illustrate that the suggested approach is straightforward, effective, adaptable, and reliable, and also validate the precision of the error estimates reported in the study.</p>