This study develops a new numerical model adopting a generic relation between the nonaqueous phase liquid (NAPL) mass and aqueous-phase NAPL concentration to simulate the relationship between NAPL contaminant mass discharge and contaminant mass reduction in the source zone, which plays a critical role to assist with site management decisions on contaminated zone remediation. The model can accommodate any contaminant mass and concentration relations and applicable to the situations when groundwater flowrate in the NAPL source zone varies in any form temporally. Therefore, the combined effects of mass-concentration relation and groundwater flowrate variations can be examined. It is hypothesized that the NAPL mass-concentration relations reflect the spatial variability of porous media in the subsurface. The developed model is compared with results from field monitoring sites and found to exhibit high flexibility and capability in capturing the observed complex NAPL source zone dynamics. Using six contaminant mass and concentration functions of varying shapes, we show that contaminant mass and concentration relation has pronounced effects on contaminant mass discharge dynamics in addition to the groundwater flowrate temporal variations. In general, the coupled mass-concentration relation and groundwater flowrate variations demonstrate stronger capability in capturing the NAPL source zone dynamics under a wide range of field porous medium conditions reported in the literature than the models that only consider groundwater flux variations. In particular, the concave mass-concentration models can be used in less heterogeneous porous media, while the convex mass-concentration models are more appropriate in more heterogeneous site conditions.