Mathematical model-based accurate evaluation of the remediation process at organic pollution sites serves as an efficient approach to the management and remediation of contaminant source zones. Numerical and upscaled analytical solution models are effective mathematical methods for reproducing the Dense Nonaqueous Phase Liquid (DNAPL) remediation process. However, in the current design of pollutant removal schemes, effective mass transfer models for characterizing the elution behaviors of contaminants remain lacking. In this study, two mathematical methods integrated with improved mass transfer models were employed to simulate the multi-stage contaminant elution behaviors under two distinct scenarios: A mixed-source region subjected to continuous water flushing and a residual DNAPL source treated with shorter-duration pulse flushing of the ethanol solution. Both the improved numerical model and upscaled analytical solution model demonstrated enhanced accuracy, which was attributed to the incorporation of solubilization mechanisms into mass transfer processes and the adoption of a multi-source region division method. The Mean Absolute Errors (MAE) of the numerical simulation for the two scenarios were 20.68 mg/L and 6.93 mg/L, respectively, whereas those of the upscaled model were 33.29 mg/L and 8.60 mg/L, respectively. Comparing the two improved models, the numerical model exhibited higher accuracy, while the upscaled model was characterized by faster computation speed and fewer input parameters.