Numerical investigation on three-phase flow in unsaturated soil considering porosity effects on soil hydraulic properties

. Ground contamination by mineral oil and hydrocarbons is a serious problem worldwide. These non–aqueous phase liquid (NAPL) represents long-term contamination sources to the soil, groundwater, and ecological environment. In the existing numerical simulation models utilized in most simulation programs, the relation between porosity and the hydraulic properties of unsaturated soils is highly simplified as nonlinear. However, the nonlinear porosity effect on soil permeability and the capillary pressures cannot be well considered in this way. In this study, a new mathematical model has been developed to incorporate porosity effects into the hydraulic properties of soils, including the retention behaviour and permeability function of NAPL and water. This model has been implemented in MATLAB using the finite difference method and verified by a centrifuge test about NAPL infiltration. Then, a series of parametric studies will be carried out to investigate the multi-phase flow under a continuously leaking NAPL source. The influence of porosity magnitude and distribution will be investigated. Special attention will be paid to the influence of porosity distribution by soil consolidation on NAPL flow, which is often ignored by the existing simulations using constant porosity values.


Introduction
The presence of non-aqueous phase liquid (NAPL) in soil and groundwater is a common problem at industrial sites such as petrol stations, workshops, and oil refineries. To properly assess and manage this contamination, it is important to determine the extent and distribution of NAPL in unsaturated soil [1][2][3]. Numerical predictions of multi-phase flow and distribution can be useful for optimizing NAPL remediation and recovery efforts and have been found to be effective in complex ground conditions [4].
There have been many numerical models developed for simulating multi-phase flow in unsaturated soil, such as NAPL Simulator [5], STOMP [6], TMVOC [7], and PFLOTRAN [8]. These models all require an accurate representation of the hydraulic properties of unsaturated soils, including the retention behaviour and permeability function of water and NAPL, which can be challenging to model due to the complex interactions between NAPL, air, water, and soil. Previous studies often simplified or overlooked the effects of porosity on soil properties. While some researchers assumed a linear relationship between porosity and permeability [9], others [10] incorporated effective porosity, which is the ratio of the volume of conducting pores to the total volume. However, these approaches did not consider the influence of porosity on retention behaviour, which is described by constitutive models such as those of Brooks and Corey [11] and van Genuchten [12]. Experimental studies have shown that porosity can significantly affect NAPL retention behaviour and the * Corresponding author: c.zhou@polyu.edu.hk shape of water retention curves [13][14][15]. In order to more accurately predict multi-phase flow, it is important to take a unified approach in considering the effects of porosity on hydraulic properties, rather than simply adjusting model parameters for different soil materials. This is especially relevant in field situations where porosity may change due to factors such as soil deformation under self-weight or external loads. Incorporating the influence of porosity on soil properties into models can help to more accurately predict the distribution and flow of NAPL and other soil behaviours.
In this study, a numerical model was applied to simulate the flow and distribution of NAPL in various ground conditions, taking into account porosity effects on soil hydraulic properties. The model was implemented in MATLAB using the finite difference method, and was verified using a centrifuge test reported by Pasha, Hu [16]. The results of the model indicated that porosity can significantly affect NAPL flow by directly influencing the fluid retention equation.  Volumetric fractions of these four phases can be described using the following variables:

A new model for the coupled water and NAPL flow with effects of soil porosity
where n is the porosity; SW and SN are the water and NAPL saturation, respectively, and their sum St is the total liquid saturation, so the air saturation is thus equal to (1-St); V with the subscripts are the volumes of the above four phases. In addition, the following two variables are defined and used: where W S and t S are the effective water saturation and effective total liquid saturation, respectively; Sm is the residual degree of saturation, which is also referred to the apparent minimum or "irreducible" wetting fluid saturation in the literature [17].
When multiple types of fluid exist in the soil, they interact with each other and cause curved interfaces. The pressure difference between the fluids is known as capillary pressure, with water having the lowest pressure due to its highest wettability, and air having the highest pressure due to its lowest wettability [18]. The difference between air pressure (PA) and water pressure (PW) is defined as the air-water capillary pressure (PcAW): Similarly, the air-NAPL capillary pressure (PcAN) and NAPL-water capillary pressure (PcNW) are defined as follows: where PN is the NAPL pressure. The relationship between fluid saturations and the capillary pressures is the key governing equation in multi-phase modelling in unsaturated soil. The water retention curve (WRC) model of van Genuchten [12] was modified and used by previous researchers [17] for this purpose: where PcNW is the NAPL-water capillary pressure; α is a model parameter; βNW is a fluid type-dependent scaling factor, defined as the ratio of air-water interfacial tension to NAPL-water interfacial tension. It has been verified by some experimental results [19], but it has an obvious limitation (i.e., without consideration of porosity effects). This limitation has caused several issues to arise from both theoretical and practical perspectives. Some theoretical modellers [20][21][22] summarised the data in the literature and found that with a decrease in porosity, the displacement pressure becomes larger associated with a smaller α, and the value of m can be approximated by a constant value. Based on these findings and the theoretical model of Tarantino [21], the current study adopts the following equation to describe porosity effects on α: where m is the parameter in the van Genuchten model (see Eq. (10)) and a is a new parameter associated with soil type. Eq. (11) suggests that as soil porosity decreases, α decreases and hence displacement pressure and water retention ability increase.
The newly proposed equation for the WRC of NAPL-contaminated unsaturated soils incorporates porosity effects by introducing a nonlinear α -n relation (Eq. (11)). This equation differs from existing equations in the literature.
Leverett [24] proposed an NAPL retention model, in which the total liquid saturation is a function of the air-NAPL capillary pressure. Based on this model, the following equations can be derived [17]: where PcAN is the air-NAPL capillary pressure; t S is the effective total liquid saturation; m and α are the same soil parameters as those in Eq. (10); βAN is the scaling factor defined as the ratio of air-NAPL and air-water interfacial tension. These equations were verified by Busby, Lenhard [19]. Substituting the α -n relation (i.e., Eq. (11)) into Eq. (13), it is obtained that . 3. NAPL retention curves with porosity effects. Fig. 3 demonstrates how the NAPL saturation is affected by two capillary pressures (i.e., PcAN and PcNW) and porosity in two different soils. It is clear that porosity has a significant effect on the NAPL retention behaviour.
The fluid flow is described using Darcy's law: where k is the effective permeability of the fluid (air, water and NAPL, m 2 ); μ is the dynamic viscosity of the fluid (Ns/m 2 ); P  is the gradient of the fluid pressure; g is the gravitational acceleration which is equal to 9.8 m/s 2 ; z is the elevation (m).
The effective permeability in Eq. (15) is a measure of soil's liquid conductivity. It is a function of the intrinsic permeability and relative permeability [25,26]: where K is the intrinsic permeability (m 2 ); kr is the relative permeability, which is the ratio of the permeability at a given saturation to the permeability at the saturated condition. These two variables are cribbedbed below: where C1 (m 2 ) and C2 are model parameters. Parameter C1 mainly describes the influence of soil type, while parameter C2 governs the sensitivity of permeability to porosity. When C2 is equal to 1, the equation is similar to the formulation of Leverett (1941) and predicts a linear relationship between permeability and porosity. At other C2 values, this equation predicts a nonlinear permeability-porosity relation. The relative permeability in Eq. (16) is a scaling factor for considering saturation effects on permeability. In this model, equations derived by Mualem [25] and Parker [17] to describe the relative permeability of water and NAPL, respectively:

Numerical implementation and model verification
This study implements mathematical equations in MATLAB based on the finite difference method (FDM).
To accomplish this, a simple orthogonal grid system is used for spatial discretization, and an explicit FDM method for temporal discretization. After the spatial and temporal discretisation, nonlinear algebraic equations were obtained and solved. A very small time increment was used for each step to minimize numerical errors. Then, a range of spatial and temporal discretization is tested to ensure stable and convergent results. The primary goal of this study is to investigate the influence of soil porosity with the help of a newly developed mathematical model. In addition, there are ten independent model parameters. Parameters ρW, ρN, μW, μN, TNW and TAN are fluid parameters, and they can be determined readily based on the types of fluid. Parameters m, a, C1 and C2 are soil parameters and all of them can be calibrated from experimental results. The first three parameters are obtained by fitting the water retention curves at different densities, while the remaining two parameters are determined from the relationship between water permeability and porosity.
The performance of the model was evaluated by comparing it to the results of a centrifuge test reported by Pasha, Hu [16]. This test was designed to study the 2-D flow of NAPL in a fine-grained soil and was conducted at 50 g-level. The prototype scale test lasted for 27.2 days, ending at 33.2 days. The numerical simulation in this study used the same soil and fluid parameters, initial and boundary conditions as the centrifugal test. The wetting fronts of NAPL at different times are summarized in Fig. 4a (centrifugal test) and Fig. 4b (numerical simulation). The average error between the measured and computed wetting front depths is 5.9%. Fig. 4c shows the profile of NAPL saturation at the center line of the plume after 27.2 days of leakage. The differences between measured and computed results are generally less than 10%. After this verification, the new model is applied to carry out parametric studies, as shown in the following section. Table 1 summarizes the results of three series of parametric studies conducted to investigate how soil porosity condition influence NAPL flow. These parameters were chosen because they are known to have a major impact on the hydraulic properties and NAPL flow, and because soil porosity and soil type can be easily determined in the field. The values of the parameters were determined based on experimental data from the literature.  3 ) from 0.408 to 0.519 porosity are uniformly distributed along the depth of each model. Series P2, using the same clay as P1, applies a descending porosity with depth, which is common in the field. By comparing the results of P1 and P2, the influence of porosity distribution on NAPL flow can be revealed. To eliminate the impact of descent porosity, In both two series, the water table is at a depth of 12 m and water pressure follows the hydrostatic distribution, with the corresponding water content distribution dependent on the water retention behaviour of soils.  Table 2 shows the properties of clay used in the parametric studies. It should be noted that the clay used in this study has low plasticity. Its water retention curves at various porosities were measured by Zhang, Wang [23]. These measurements indicate that the displacement pressures range from 5 to 25 kPa, depending on the porosity. This translates to a capillary rise of approximately 0.5 -2.5 m, which is much smaller than the thickness of the soil layer above the water table.

Numerical simulation program and procedures of the parametric studies
The nodes of the simulation had varying sizes from 0.05 m to 0.2 m, and the initial water distribution was theoretically calculated to reach a steady state condition. The boundary conditions for NAPL leakage were applied to the left side of the ground surface, where the NAPL pressure head was maintained at 0.3 m and the air and water fluxes were set to zero. On the right side of the ground surface, a boundary to atmosphere was applied, which imposed zero flux for water and NAPL and zero pressure for air. Along the left and right boundaries, zero flux was applied for all air, water and NAPL. The groundwater table was maintained at the bottom boundary, with zero water pressure and zero flux for air and NAPL.

Influence of porosity on NAPL flow in the uniform clay ground
The total volume of NAPL leaked into the ground, as computed from Series P1 (Fig. 5), is affected by the porosity of the soil. Five porosities ranging from 0.408 to 0.519 were considered. In the first stage, the differences induced by a change in soil porosity were found to be negligible due to the low NAPL saturation and its relative permeability being close to zero. This meant that the computed NAPL volumes at different porosities were similar. At the second stage, the NAPL saturation and its relative permeability increased, making the intrinsic permeability of the soil the dominant factor in the predicted results and making

Series
Objective The corresponding ground condition V To verify the numerical model Same soil properties and ground condition as the centrifugal test [16] P1 To investigate the influence of porosity (uniform) n = 0.408 n = 0.444 n = 0.482 n = 0.500 n = 0.519 [23] P2 To investigate the influence of porosity distribution n decreases with depth from 0.482 n decreases with depth from 0.519 porosity effects more obvious. For example, the NAPL volume accumulated within 36 months was found to be 43 times greater for a porosity of 0.519 than for a porosity of 0.408. Furthermore, the existence of two stages was more obvious in the higher porosity cases due to a lower initial water saturation.
The maximum vertical front depth of the NAPL plume in Series P1 is shown in Fig. 6. It can be seen that the vertical distance of the plume from the leakage source increases from 0.70 m to 5.35 m after 36 months as the porosity increases from 0.408 to 0.519. This is due to two effects of porosity on the hydraulic properties of soils. Firstly, the intrinsic permeability is much higher at a higher porosity. Secondly, a larger porosity results in a lower water retention capacity and hence a smaller initial water saturation (i.e., the condition prior to the leakage of NAPL). The initial water saturation at the ground surface decreased from 53.9% to 33.3% as the porosity increased from 0.408 to 0.519. This reduction in water saturation led to a higher NAPL saturation and then a larger relative permeability for NAPL. Based on the above analysis, it can be concluded that soil porosity has a significant impact on NAPL flow, particularly in the long term. This finding can be used to help engineers predict or minimize NAPL contamination in practical applications, such as when designing petrol and gas stations and industrial sites. To improve the accuracy of numerical simulations, a porosity close to the in-situ condition should be used, as it can vary significantly depending on the geological history and loading path. Fig. 7 shows the influence of porosity on the total NAPL volume leaked into the clay ground. Two cases are compared in Series P2, including uniform porosity and descent porosity along with the depth. For each case, two different porosities at the ground surface (n = 0.482, 0.519) are used. When the porosity distribution is changed from uniform to descent condition, the accumulation rate of NAPL volume keeps approximately constant at the early stage, but then significantly reduces. Furthermore, this trend is more significant when the surface porosity is higher (i.e., n = 0.519). Compared to the case of descent porosity, the total NAPL volume at 36 months in the case of uniform porosity is about 3.5 times and 2.5 times larger when the surface porosity is 0.519 and 0.482, respectively.    8 shows the maximum vertical front depths of the NAPL plume in Series P2. Compared to the case of descent porosity, the vertical front depth at 36 months in the case of uniform porosity is about 3.1 times and 2.4 times deeper when the surface porosity is 0.519 and 0.482, respectively. In the range of higher porosity, the total NAPL volume is more sensitive to the change in soil porosity. In addition, the descent porosity resulted in a higher water retention capacity and a higher initial water saturation than the uniform condition, which caused a lower accumulation rate of NAPL infiltration.

Summary and conclusions
A new theoretical model was developed to simulate NAPL flow in the vadose zone, considering porosity effects on soil hydraulic properties. The finite difference code was verified using a centrifuge test and parametric studies were carried out to investigate the flow of NAPL upon its leakage at the ground surface. Results showed that the wetting front and leakage volume of NAPL are greatly affected by soil porosity, and that heavy compaction can be applied to minimize the NAPL contamination. The leakage volume increases nonlinearly with an increase in porosity and is much larger if the porosity is uniform rather than descending along with the depth. The relationship between the leakage volume of NAPL and time, suggesting that it is more efficient to control NAPL flow at the first stage.