Modeling of nonaqueous phase liquid (NAPL) migration in heterogeneous saturated media: Effects of hysteresis and fluid immobility in constitutive relations

Abstract

The confidence in model predictions for nonaqueous phase liquid (NAPL) transport in stochastically heterogeneous systems is limited. The fundamental approaches as well as the constitutive models have not been sufficiently validated, mainly because of the lack of appropriate experimental data. Recently, Fagerlund et al. (2007a, 2007b) presented a set of well‐controlled laboratory data that are used here (1) to analyze the overall performance of the continuum‐based approach for predicting two‐phase NAPL‐water flow in stochastically heterogeneous media and (2) to compare the predictions from different constitutive models. The five models tested were the nonhysteretic Brooks‐Corey‐Burdine (BCB) and van Genuchten‐Mualem (VGM) models, the hysteretic versions of these models (HBCB and HVGM), and the HVGBCB model, a model combining the hysteretic van Genuchten (HVG) Pc‐S relation and the hysteretic Brooks‐Corey‐Burdine kr‐S relation. Two cases of NAPL migration were considered: a layered system of two homogeneous sands separated by a dipping interface and a system where one of the layers was stochastically heterogeneous. The results showed that the best models could indeed capture the main characteristics of the spreading and immobilization well, demonstrating the validity of the continuum‐based approach for this level of stochastic heterogeneity. Implementation of hysteresis was necessary for correct prediction of the observed speed of NAPL migration as well as the amount of immobilized NAPL. The three hysteretic models were similar in their overall prediction error‐based performance. The HVGM model produced less overestimation of NAPL saturations but instead underestimated the entrapment at capillary barriers in comparison to the HBCB and HVGBCB models. The HVGM model also overestimated the speed of NAPL migration, which is attributed to its closed‐form kr‐S function, for which the VG parameter m has to be fitted under the constraint m = 1 − 1/n. The HVG (and VG) Pc‐S function, in contrast, used a different set of VG parameters produced with no constraint on m, which better represented the pore size distributions of the sands. A relation for partial nonwetting phase immobility during drainage is also presented.