The identification of dense non-aqueous phase liquid (DNAPL) distribution in the subsurface is important for optimizing remediation strategies. However, the DNAPL distribution is highly sensitive to subsurface heterogeneity and is controlled by multiphase physics. Typical multiphase models are computationally too expensive to be applied for inverse modeling. Therefore, we first propose a sequential inversion strategy to estimate the hydraulic conductivity (K) and DNAPL saturation (S_N) fields separately, which can implicitly consider the interdependences between K and S_N by using indirect conditioning data. Nevertheless, when the observation data is limited and less informative, ignoring the underlying multiphase physics may result in poor estimation. To improve DNAPL imaging, we further propose a machine learning-based method to capture the underlying multiphase physics between K-S_N. In conclusion, we suggest that: (1) With enough and informative measurements, one can use the sequential strategy to produce accurate DNAPL estimation. (2) With limited data, one can use the machine learning-based method to capture the multiphase physics in a computationally efficient way.