Abstract: To address the challenges of complex grouting channels in goafs, nonlinear grout migration, and the high computational cost of multi-parameter analysis,a surrogate model for the grout front advancement velocity is established by combining the material point method (MPM) with non-intrusive polynomial chaos expansion (PCE). Grout density, bulk modulus,plastic viscosity, yield stress, solid volume fraction, and grouting pressure are selected as input variables, and two cases—a single goaf and a complex branched goaf—are considered. MPM is used to calculate the grout velocity field. A total of 168 modeling samples are used to fit the PCE coefficients, and independent samples are used to compare the surrogate predictions with the MPM results. The results show that in the single goaf,high-velocity zones extend from the vertical grout-conducting channel to the horizontal channel, whereas in the complex goaf, the velocity distributions among branches differ markedly. For the complex case, the in-sample CVRMSE values of the first-, second-, and third-order PCE models are 24. 91%, 11. 37%, and 3. 05%, respectively, with the third-order model yielding the smallest fitting error. The coefficient grouping results show that grouting pressure has the highest total-effect index in both cases,while in the complex case,the total-effect ranking of grout density rises from third to second. In a field fluidized backfill trial,a fly ash ∶ gangue ∶ water ratio of 2 ∶ 5 ∶ 3 was adopted;continuous pumping was achieved,and no obvious pipeline blockage or grout leakage was observed. The combination of MPM and PCE can be used to analyze the grout velocity response within the given parameter ranges;the field trial indicates that the mix proportion and backfilling process used meet the pumping requirements of the trial.