Abstract A better understanding of the distribution of nonaqueous phase liquid (NAPL) plumes is of great importance to groundwater pollution remediation and control. However, the efficiency of surrogate models in simulating the transport is still not well addressed. Selecting a leakage problem as an example, 50 sets of random permeability distributions are generated using the Monte Carlo method, and a numerical model is used to obtain benchmark data of NAPL transport. Four machine learning methods are employed to simulate dense NAPL transport under point leakage sources across spatiotemporal scales. The validation of the models demonstrate that the random forest model can also effectively capture the spatial-temporal distribution of the plume in heterogeneous aquifers, with a maximum mean absolute error and root mean square error smaller than 5.55 × 10−4 and 5.88 × 10−5, respectively. Meanwhile, the multiple phase outcome from the random forest model fits well with the numerical results under the scenarios of linear leakage sources and light NAPL transport. The total time consumed in the computation is reduced by over 150 times after using the surrogate models. The results suggest that surrogate models can provide a promising way to understand NAPL transport in heterogeneous aquifers.