Simulation of Dynamic Two-Phase Flow During Multistep Air Sparging

Air sparging is an in situ soil/groundwater remediation technology, which involves the injection of pressurized air through air sparging well below the zone of contamination. To investigate the rate-dependent flow properties during multistep air sparging, a rule-based dynamic two-phase flow model was developed and applied to a 3D pore network which is employed to characterize the void structure of porous media. The simulated dynamic two-phase flow at the pore scale or microscale was translated into functional relationships at the continuum-scale of capillary pressure–saturation (Pc–S) and relative permeability—saturation (Kr–S) relationships. A significant contribution from the air injection pressure step and duration time of each air injection pressure on both of the above relationships was observed during the multistep air sparging tests. It is observed from the simulation that at a given matric potential, larger amount of water is retained during transient flow than that during steady flow. Shorter the duration of each air injection pressure step, there is higher fraction of retained water. The relative air/water permeability values are also greatly affected by the pressure step. With large air injection pressure step, the air/water relative permeability is much higher than that with a smaller air injection pressure step at the same water saturation level. However, the impact of pressure step on relative permeability is not consistent for flows with different capillary numbers (Nca). When compared with relative air permeability, relative water permeability has a higher scatter. It was further observed that the dynamic effects on the relative permeability curve are more apparent for networks with larger pore sizes than that with smaller pore sizes. In addition, the effect of pore size on relative water permeability is higher than that on relative air permeability.

sizes than that with smaller pore sizes. In addition, the effect of pore size on relative water permeability is higher than that on relative air permeability. Keywords Pore network model · Dynamic two-phase flow · Air sparging · Capillary pressure-saturation · Relative permeability-saturation

Introduction
Air sparging is a well-developed method to remediate saturated soils and groundwater contaminated with volatile organic compounds (VOCs). Over the past two decades, the use of in situ air sparging (IAS) in conjunction with soil vapor extraction systems (SVE) has been proved to be an efficient, fast, and relatively economical system to remediate saturated soil deposits contaminated with VOCs (Leonard and Brown 1992). During air sparging, pressurized air/oxygen is injected below the region of contamination, and while migrating through VOCs-contaminated zone, it partitions the contaminants into the vapor phase and transports them to the unsaturated zone. The contaminant-laden air can be captured by SVE system, and then the extracted contaminated air is cleaned using conventional methods, such as carbon filters and/or combustion. During air sparging mass transport is the primary mechanism controlling contaminant removal (Bausmith et al. 1996); hence, the efficiency of in situ sparging system is mainly controlled by the extent of contact between injected air and contaminated soils and/or pore fluid. Therefore, the determination of size and shape of zone of influence (ZOI) and the understanding of air flow dynamics and spatial air distribution within ZOI are helpful in design of in situ air sparging systems (Hu et al. 2010. The two important functional relationships for modeling the two-phase flow during air sparging are the correlation between relative permeability-saturation and the correlation between capillary force-saturation. Usually these correlations are developed based on experimental data gathered under static equilibrium or steady state conditions. Those relationships are applied to analyses of both steady and transient flows, which implies that those relationships between water content, matric potential, during monotonous draining or imbibition are not affected by the dynamics of water flow. However, the relative permeability depends upon the microscopic distribution of the fluids within the pore space, and under different air injection rates with different capillary numbers, the distribution of liquids at the pore scale and flow patterns are quite different depending on the level of saturation, the saturation history, and rate of change of saturation (Jerauld and Salter 1990). As a result, the relationship between relative permeability and saturation is not unique during air sparging under different air injection pressures. One of the objectives in this study is to investigate the impact of air injection pressure on relative water/air permeabilities. Similarly, the other functional relationships of the correlation between capillary force and saturation, expressed in the form P n − P w = P c = f (S), is usually obtained experimentally under equilibrium conditions. P n and P w are the average pressures of non-wetting and wetting phases, respectively. P c is capillary pressure and S is the saturation for wetting phase. To obtain a drainage (or imbibition) curve, one starts with a wet (or dry) soil sample, then the capillary pressure is increased (or decreased) incrementally. At each step the water content is measured after equilibrium is reached. The typical time needed to construct a complete capillary pressure-saturation curve is couple of weeks or months. However, during multistep air sparging operation, each step of injection pressure may not be maintained for sufficiently long time to reach the equilibrium or steady state. Hence, the capillary pressure-saturation relationship curves obtained under equilibrium conditions cannot adequately describe the relationship between capillary force and saturation within hours. In fact, Topp et al. (1967), Davidson et al. (1966), Smiles et al. (1971), and Vachaud et al. (1972) showed that both the retention characteristic and the hydraulic conductivity characteristic are different under the static equilibrium/steady state and unsteady state cases, and they depend on both the history and the change in rate of saturation. Hence, if dynamic non-equilibrium occurs between the water content and the water potential during transient water flow, the water retention characteristic measured in the laboratory under equilibrium cannot be reliably applied to that under transient flow. Hence, the other objective of this study is to study the dependence of capillary curves on the change in rate of saturation.
To investigate the air flow dynamics on both functional relationships during air sparging, a 3D pore network can be employed to characterize the void space in the porous media. In the pore network, the pore spaces of porous media can be represented by a 3D network of interconnected pores and throats. Then, a dynamic two-phase flow model can be developed to capture the essential features of flow through natural porous media based on solution of flow equations and interface settings in the 3D pore network. A brief review of the current research on both characterizing void space in porous media and simulating the two-phase flow through porous media is presented in the following sections.
The most critical part of characterizing the void space in a porous medium is the identification of the geometrical and topological properties of the voids within the porous media. The predictive value of pore network models depends on the accuracy with which the network captures the complex geometrical and topological properties of pores. Blunt et al. (2002) suggested that if complex pore geometries can be adequately represented, then pore scale models can be developed into predictive models. Arns et al. (2004) who performed a comprehensive study of effect of network topology on relative permeability showed that matching the average coordination number of a network is not sufficient to predict relative permeability. However, reasonable relative permeability values can be obtain by matching coordination number distribution. Gao et al. (2011) proposed and validated one method of generating a 3D equivalent pore network of porous media to capture the pore fabric. This pore network can match both main geometrical and topological characteristics of the porous media including pore size distribution, coordination number distribution, and spatial correlation between pore size and coordination number. In addition, Gao et al. (2011) introduced the Biconical Abscissa Asymmetric CONcentric (BACON) bond ( Fig. 1) as a replacement for traditional cylindrical bond with constant cross-sectional dimension to describe the connection between two adjacent voids, and the size and position of pore throat was correlated to its two neighboring pore bodies.
After characterizing the void space of porous media by a 3D pore network, another main task of this research is to simulate the dynamic air-water two-phase flow in the characterized pore space. There are two main types of models to study the displacement of one phase by the other in porous media. One is percolation model while the other is diffusion-limited aggregation (DLA) model. In the limiting case of capillary-dominated flow, the invasion percolation model (Chandler et al. 1982;Mohanty and Salter 1982;Wilkinson and Willemsen 1983;Lenormand and Zarcone 1984) adequately describes two-phase flow through porous media using concepts based on percolation theory. In the other limiting extreme of viscous-dominated flow, the DLA model describes the displacement of high-viscosity fluid by a low-viscosity fluid, thereby giving rise to pore scale viscous fingering (Meakin 1983a,b;Meakin and Stanley 1984;Witten and Sander 1983;Chen and Wilkinson 1985). Both theories described the porous media as a set of interconnected pores. These theories expand the capillary tube models by allowing virtually any interconnectivity and pore size distribution and associated pore space. In summary, percolation theory was successful in simulating flows dominated by capillary pressures, while DLA theory was successful in simulating the fingering of a highly mobile fluid into a very viscous fluid. However, these models do not contain any physical time for the front evolution and they cannot describe the crossover between the two major flow regimes.
In general, two-phase flow situations during air sparging both viscous and capillary forces are present. Hence, a suitable combination of rule-based model is needed to adequately describe two-phase flow. For this, one can use a rule-based model to update the fluid configuration and then to determine the pressure field while immobilizing the interfaces and such models are known as dynamic network models. The essential difference between the static models and dynamic models lies in the computation of the pressure field. In static network models (e.g., percolation theory), the pressure field is assumed to be uniform in each phase (in the absence of gravity). However, the dynamic network models estimate the pressure field at each time step and make rule-based interface movements based on the computed pressure drops in both fluids.
Two main goals of this research are to develop a rule-based dynamic two-phase flow model that can simulate the dynamic of two-phase flow during air sparging and to investigate dynamic effects of two-phase flow on functional relationships of two-phase flow during multistep air sparging under different boundary conditions using the developed dynamic two-phase flow model.

Network Description
The equivalent pore network model of porous medium is employed to characterize the void space of porous media. The procedure for constructing an equivalent network is presented in Gao et al. (2011). The main geometric and topological parameters that can be considered in the network include total number of pores in the network, average pore size, maximum and minimum pore size, standard deviation of pore size, pore spacing, parameters for characterizing the BACON bonds connecting two neighboring pores, average coordination number, standard deviation of coordination number. All these parameters can be extracted directly from packed spheres . Gao et al. (2011) showed that both the pore size distribution and the coordination number distribution approximately obey normal distribution for the random close packing of particles. There is no direct proof that the pore bodies are spatially correlated in size; however, there is weak spatial correlation between the pore body size and coordination number. As a result, in this research, the spatial pore space correlation is not considered and all the pore sizes are randomly assigned to the nodes in the network. While assigning a generated coordination number distribution to each node, the larger coordination number is assigned to a bigger pore body.

Model Assumptions and Initial & Boundary Conditions
In the dynamic two-phase flow model, the injected fluid is non-wetting fluid-air while the defending fluid is the wetting fluid-water. The non-wetting fluid can be assumed to be Newtonian, incompressible and immiscible. It is assumed that the wetting phase is "perfect wettable" with no contact angle hysteresis. The contribution from gravity of water and air during the flow are neglected in this microscale model. The two types of fluids can coexist both in bond and nodes. Pore body and throat conductance are computed by assuming stokes or creeping flow at pore level. It is also reasonable to assume that all pressure drops occur in the bonds while the pressure drop in nodes are neglected.
The medium is initially saturated with the wetting phase and the non-wetting phase is injected at the bottom of the network as in the case of air sparging. During multistep air sparging simulation, a constant air injection pressure is maintained at each pressure step at the bottom of the network, while at the top atmospheric pressure is maintained. Boundary conditions in the other two directions are no flow boundaries. During the simulation of twophase flow all the bottom nodes are air saturated and the air-water interfaces in the bonds connecting to the bottom nodes are introduced at the beginning of the simulation.

Rule-Based Dynamic Two-Phase Flow Model
A rule based on dynamic air-water two-phase flow model was developed to simulate the twophase flow in the pore network during air sparging. Different from the percolation model and the DLA model, the rule-based dynamic two-phase flow model can consider both the viscous and the capillary forces during two-phase flow. Under very low air injection pressures, the capillary fingering can be simulated while under very high air injection pressures, the viscous fingering can be simulated. The flow mechanisms and fluid configurations are briefly described below.

Flow Through a Bond
In a given bond, there may be either one or two phases present, corresponding to single-phase or two-phase flows. In both the cases, linearly of the Poiseuille Law relates the pressure drop and flow rate in each bond as expressed by equation below: where subscripts i and j represent pore bodies i and j, respectively, q i j is the volumetric flow rate in the bond connecting pores i and j, G i j is the conductance of the bond connecting pores i and j, P nw i is non-wetting phase pressure in pore i, P w j is wetting phase pressure in pore j, and P ci j is the capillary force at the interface between wetting and non-wetting phases, which depends on the position of meniscus inside a bond.
According to Gao et al. (2011), the flow rate in a BACON bond i-j as shown in Fig. 1 connecting pores i and j for two-phase flow can be expressed as In Eq. (2), an effective viscosity μ eff is introduced to consider the coexistence of the two fluids inside a bond and depends on the volume fraction of each phase in the bond. μ eff can be calculated as P i , P j and P c are the nodal pressure of nodes i and j and the capillary pressure in the bond connecting nodes i and j. Other parameters y i , y j , n, b s , b e , b t , l are the parameters describing the configuration of the BACON bond as shown in Fig. 1.

Movement of Meniscus
The moving direction of a meniscus depends on the pressure drop in the bond. The menisci can both invade into or retreat from a bond. This interface movement is different from that used in static multiphase flow models (Blunt et al. 2002), where the air-water interface is frozen by capillary force if P i − P c ≺ P j , as a result, there is no flow cross a bond until P i increases.

Flow Configurations
Due to the movement of menisci in bonds and pores, a bond can be occupied by water or air or both fluids. Based on the types of fluids in a bond and also in the nodes on the both sides of the bond, there can be three different flow configurations. They are bonds with one air-water interface, bonds with two air-water interfaces, and bonds with no fluid interfaces. It should be noted here that, the fluid configuration of a bond with more than two air-water interfaces are not considered in this research. A pore body has two main types of configurations: the first one is called two-phase node while the second one is called one-phase node.
Two-phase node: For this type of nodes, both air and water phase can coexist in the node and can flow in/out of the node simultaneously. When there are two types of fluids in a node, each type of fluid should have at least one bond connected to it and filled with the same fluid. The typical configurations of two-phase nodes are shown in Fig. 2.
One-phase node: For this type of nodes, there can be either two types of fluids or single type of fluid in a node. When there are two types of fluids in a node, only one type of fluid is continuous in the pore network, while the other one is isolated in the pore network without any bond connected to it with the same fluid. In this case, the saturation of the isolated phase in the node is called residual saturation of the node. During the simulation, it is assumed that the isolated phase in the node does not move until the node is filled with the isolated phase. The typical configurations of one-phase nodes are shown in Fig. 3.

Determination of the Flow Field
As fluids are assumed incompressible and immiscible, the volume flux is conserved everywhere in the network. Hence, the pressure field of the pore network during two-phase flow simulation can be obtained by solving conservation of flow equations at each node. At each node, the following equation should be maintained: In Eq. (4), q i j is the volumetric flow rate from pore j to i. The summation of j runs over the nearest neighbor nodes to the ith pore where i can be all nodes that do not belong to the top or bottom planes. Equations 2 and 4 constitute a set of linear equations which are to be solved for the nodal pressure with the constraints that the pressures at the nodes belonging to the upper and lower rows are kept fixed. Equation 4 can be solved using Gauss-Seidel iterative method or conjugate gradient method.

Selection of Time Step
For the bond i-j which connects nodes i and j, a time step t i j is chosen such that every meniscus is allowed to travel at most a maximum step length x max during that time step. The calculation of time step depends on the positions of interfaces and different fluid configurations in bonds. For the configuration of one interface in a bond, the maximum step length x max is chosen such that the air-water interface can travel to the throat of a bond if it is behind the throat or can travel to the end of a bond if it is beyond the throat. If the air-water interface is already at the end of bond (at the entrance to the connected pore body), the time step is chosen such that the pore body can be filled during the selected time step. When calculating the time step for this case, attention should be paid to the special situation where more than one bond may be filled by one node simultaneously. Hence, all the bonds connected to the pore should be scanned to check the number of bonds filling pores and the time step can be selected for filling the pore by all such bonds.
Once the time step of each bond is determined, the time step for the whole pore network under current configuration can be determined which is the minimum time step of all the bonds in the network.

Determination of the Duration of Each Air Injection Pressure
Step and Updating of Fluid Configurations

Determination of the Duration of Each Air Injection Pressure
During multistep air sparging simulation, the pressure at the bottom boundary is increased by a specified pressure step each time and is maintained for a certain time until spatial distribution of fluid reaches an equilibrium configuration or until the steady state is reached.
To determine the duration of each air injection pressure step, criterion for determining the equilibrium condition and steady state should be first defined. In this research, steady state is determined by the changing rate of air saturation of the pore network. During the simulation of two-phase (air-water) flow, the saturation of the whole network and the moving average of the change in saturation per time (| s w / t|) is successively calculated for a window of 50 time steps. The equilibrium condition or the steady state is set to reach when | s w / t| ≺ 10 −3 . It should be noted here that while calculating the saturation of the pore network, both the volume of the pore bodies and bonds are considered. As the pore bodies at the bottom are always filled with air, the nodes in this layer are excluded in the saturation calculation.
The saturation calculation of each layer is defined by Eq. 5.
where n z is the total number of layers in the pore network, n p , n b are the total number of pore bodies and bonds in each layer, respectively, V p , V b are the volume of pore bodies and bonds, respectively. Once the saturation of each layer of the network is obtained, the saturation of the whole network can also be calculated based on the total volume of each layer and the corresponding saturation of each layer.

Fluid Configuration Update
Once the time step for each bond in the pore network is determined, the minimum time step t min of all the calculated time steps of each bond can be chosen as the time step for the whole network with the current fluid configurations. Then, the fluid configurations including number of interfaces in the bond, the positions of interfaces, saturation of each bond and nodes can be updated based on the selected time step and the flow rate in the bond. Then, the parameters for simulation of two-phase flow through the pore network including the capillary force, effective viscosity of fluids in each bond can be updated. When updating configurations in pore bodies, there are two special situations: one, where a pore body is initially fully saturated with water but gets completely air saturated at the end of the time step (the pore body is called a drained node) and the second situation, where a pore body is initially filled with air but imbibed by water and becomes water saturated at the end of the time step (the pore body is called an imbibed node). For both drained and imbibed nodes, while there is a phase transition in the node, it leads to an increase/decrease of number of interfaces in the bonds connected to it as shown in Fig. 4.
The flow chart of simulating the two-phase flow dynamic model is presented in Fig. 5.

The Capillary Number and Determination of Pore Network Model Parameters
To make the simulated two-phase flow results in pore networks with different geometrical properties under different boundary conditions comparable, a permeability-based capillary number is defined as where μ is the viscosity of the injected fluid, K is the single-phase permeability of the network, P is the pressure drop across the pore network according to different boundary conditions, σ is the interfacial tension, and L is the axial length of the pore network. This definition of capillary number was first proposed by Singh and Mohanty (2003). One main advantage of using the permeability-based capillary number is that the capillary number can be predetermined before performing numerical experiments by knowing the boundary conditions. As a result, it is very helpful in designing scheme of experiments. For example, to compare the two-phase flow with same capillary number in different networks, once the capillary number is given, the air injection pressures for the different networks can be determined based on the given capillary number and the permeability of each pore network.
Another important parameter need to be defined before the simulation of two-phase flow is the network dimensions. The two-phase flow simulation is time consuming for large 3D pore network; however, if the size of the pore network model is too small, it cannot produce reproducible results. To obtain the appropriate size of the network, two-phase flow in three networks with same pore size distributions and pore spacing but with different network  sizes (total number of pores) were studied and compared. The main parameters of the three networks are shown in Table 1, where ther , r min , r max , and σ r are the average pore size, the minimum pore size, the maximum pore size, and the standard deviation of pore size distribution respectively, while n and σ n are the average coordination number and standard deviation of coordination number distribution. In addition, curvature constant is a parameter for describing the configuration of a bond as shown in Fig. 1. The two-phase flow was simulated in three networks with the permeability-based capillary number of 2.2E−3. The relative air and water permeability values obtained from three networks are compared. Here, relative permeability is defined as the ratio of the effective permeability of the porous medium to the absolute permeability of the porous medium. The effective permeability of a permeable medium is a measure of the ability of the material to conduct one fluid phase of a multiphase fluid. When air-water two phases are present in the network, if flow rates are computed using the same pressure drop for single-phase flow, then the relative permeability is where μ nw is the viscosity of air, Q mw is the total water flow rate across the network, which is obtained by summing water flow over all bonds connected to the outlet, and Q mnw is the total air flow rate across the network for two-phase flow, which is obtained by summing air flow over all bonds connected to the outlet. Please note that in both formulae for calculating relative permeability values, the flow rates of each phase are based on outflows of each phase at the outlet of the network; hence, the flow of the isolated clusters of each phase in the network during dynamic flow cannot be considered. The P inlet and P outlet are pressures at the inlet and outlet of the network, respectively. The K is the absolute permeability of the network which can be calculated from the Darcy's law: where μ w is the viscosity of the water, Q sw is the total flow rate across the network for flow of water (this is the flow summed over all bonds connected to the inlet or outlet), A is the cross-sectional area of the network model, and L is the axial size of the network. The simulated relative air/water permeability values from three networks are compared in Fig. 6. In Fig. 6, the solid symbols represent relative air permeability values while the hollow ones represent relative water permeability values. It can be observed that the predicted Fig. 7 The applied three multistep air injection pressures relative permeability values from the three networks are similar. To some extent, the relative permeability values can reflect the connection between the developed air flow channels in the pore network. As the relative permeability curves in the three pore networks are quite similar, it can be concluded that the air channel development in the three pore networks are almost same. Based on the above observation, it can be concluded that the network with size of 12-12-12 can generate reasonable results as those from the network with size 15-15-15 and the network with size of 17-17-17. As a result, in this research, the pore network size of 12-12-12 was used to simulate remaining numerical experiments of two-phase flow, and this network is denoted as Network I in this manuscript.

Simulation of Multistep Air Injection Pressure
In this section, the dynamic two-phase flow during air sparging with multistep air injection pressures was investigated. The constitutive relationship governing multiphase flow which can be expressed as functional relationships of capillary pressure, saturation, relative permeability of coexisting phases in both steady state and transient state are compared. The parameters of the Network I with the dimension of 12-12-12 used in the simulation are presented in Table 1.

Unsteady State Capillary Pressure Curve
To study the dynamic effects of two-phase flow, a 20-step pressure with the maximum air injection pressure of 20 kPa was applied to the Network I at different time intervals. Here, the time interval is defined as the elapse time or duration of each air injection pressure. Three time intervals were employed in this numerical simulation: 0.125, 0.25, and 0.5 ms. The applied air injection pressure with each time interval is shown in Fig. 7.
During the numerical simulation, both the air injection pressure and the saturation of the whole network were recorded for each set of injection pressures. The dynamic capillary pressure curves for all the three numerical tests are presented in Fig. 8. It can be observed from Fig. 8 that the dynamic capillary pressure curves are always above the steady state capillary pressure curve. This is consistent with experimental observations of Topp et al. (1967), who compared the correlation between capillary force and water saturation during steady and unsteady flows, observed that during unsteady state, the dynamic capillary pressure is significantly higher than that performed under equilibrium or steady state. Topp et al. (1967) Fig. 8 Comparison between dynamic and steady capillary pressure curve stated that the presence of disconnected pendular water rings or the accumulation of traces of contaminants in the air-water interfaces within the sample were possible contributors for this deviation. Smiles et al. (1971) also observed that when a large drainage step was imposed, the measured capillary force for a given water content was higher than that with static drainage. In addition, Smiles et al. (1971) observed that the size of the dynamic effect depended on the rate at which water content changes or the size of the change in the boundary suction head.
From Fig. 8, it can be further observed that at very low water saturation values, there is no obvious difference between dynamic capillary curve and steady state capillary curve. This can be explained by the movement of isolated water clusters at lower water saturation values. At lower water saturation, which is corresponding to high capillary pressure, the spatial distribution of water in different tests may be different; however, as most of air flow channels have been developed at this stage, there are only a few of isolated water clusters left in different simulations. As a result, there are no obvious differences between the two capillary curves at lower water saturation. Based on this explanation, it can be predicted that there will be significant differences between dynamic capillary curve and steady state capillary curves when simulating air sparging with lower air injection pressures in networks with wide range of pore size distribution. Under such situations two phases of fluids will be non-uniformly distributed leading to additional isolated clusters during the two-phase flow.

Dynamic Relative Permeability Curve
The contribution of pressure steps with different magnitudes of dynamic relative permeability values was studied by conducting three multistep air sparging tests using Network I. During the three tests, the maximum air injection pressure was 20 kPa with the corresponding capillary number of N ca = 6.1E−3. Three different sets of multistep pressures with steps 1, 5, and 10 were employed in the three tests with the corresponding pressure steps of 20, 4, and 2 kPa, respectively. As the transport of injected air during the two-phase flow is important for air sparging, the relative air permeability and relative water permeability are discussed separately.
The relative air permeabilities during the different simulations with varied pressure steps are compared and presented in Fig. 9. From Fig. 9, it can be clearly observed that the simulated relative air permeability curve can be divided into three parts. Part I is for the water saturation greater than 0.75. This part corresponds to the initial stage of air injection. In this Fig. 9 Air relative permeability in Network I with N ca = 6.1E−3

Fig. 10
Front dynamics for the 1-step air injection pressure part, relative air permeability is almost zero, which means there is no continuous air flow channels developed in the soil for all the three multistep pressure sets.
Part II is for the water saturation between 0.35 and 0.75. In this part, the air-water twophase flow does not reach an equilibrium state and air flows for all the three pressure steps are in transient state. In this part, it can be clearly observed that the relative air permeability for smaller pressure steps is much lower than that with higher pressure steps. The main reason accounting for this is that for the flow with small pressure step, the capillary number is very small; as a result, the flow pattern is similar to capillary fingering. Hence, for the same water saturation, there are more isolated air clusters in the network than that with larger pressure step. To confirm this, the front dynamics of flow in the networks during the 1-and 10-step flows are compared at water saturation values of 0. 40,0.42,0.44,0.46,0.48,and 0.50 and presented in Figs. 10 and 11,respectively. From the two figures, it can be concluded that the water distribution during the 1-step pressure simulation is uniformly distributed when compared with that for 10-step pressure simulation. For example, at the average water saturation level of 0.4, the highest difference in water saturation in different layers for 1 step pressure is less than 0.15; however, the maximum difference of water saturation for all network layers for the 10-step pressure is around 0.5. This non-uniform spatial distribution of air-water phases accounts for the lower relative permeability value for smaller pressure step. However, by comparing the simulation results of 1-and 5-step pressures in Fig. 9, it can be further concluded that the dynamic effect on the relative permeability is also weakened for the flow with higher capillary number.
The third part of relative air permeability curve shown in Fig. 9 is for the water saturation less than 0.3 which corresponds to the steady state of two-phase flow during air sparging. At this stage, it can be found that relative air permeability values of the three numerical tests are almost identical to each other. To verify this observation, another similar numerical experiment was conducted. In this simulation, the maximum air injection pressure was maintained at 7.5 kPa with the corresponding capillary number of N ca = 2.2E−3. Three experiments were performed with 1-, 7-, and 15-step pressures. The variation between relative air permeability and water saturation for the three experiments are presented in Fig. 12. From Fig. 12, it can be found that relative air permeability curve can be also divided into three parts: initial stage before breakthrough, transient stage, and steady state, which is consistent with observations obtained from Fig. 9. As to the third part of relative air permeability, there are some differences between the results shown in Figs. 9 and 12. Figure 9 shows that the relative permeabilities of the three tests overlapped, whereas Fig. 12 shows that the relative permeability values of the tests with 7-and 15-step pressures are much lower than that with 1-step pressure. This difference between the two sets of simulations indicates that the effect of pressure step on relative permeability is not consistent for flow with all different capillary numbers. The main reason accounting for this observation is that for the flow with higher capillary number, the isolated air clusters formed can still be driven out during the two-phase flow, which is not the case for that with lower capillary number. For the second set of the tests, the maximum air injection pressure was maintained at 7.5 kPa and this is much lower than that of the first set of tests with maximum air injection pressure of 20 kPa. Hence, for the second set of simulations with lower capillary number, the isolated clusters developed during the two-phase flow were not able to move like that in the first set of simulations.
As the relative air/water permeability is a function of spatial distribution of non-wetting and wetting phases, it should also depend on geometrical and topological properties of pore network as well as on the capillary number. To check the contribution of geometrical properties of pore network to relative permeability, a two-phase flow simulation was performed using a network with smaller pore sizes or Network II. The parameters of the Network II are presented in Table 2. Two sets of tests were performed with this pore network. For the first set of tests, the maximum air injection pressure was kept at 78 kPa with the corresponding capillary number of N ca = 6.1E−3, which is same as that with simulation using Network I with the maximum air injection of 20 kPa. For the second set of tests, the maximum air injection pressure was 27 kPa with the corresponding capillary number of N ca = 2.2E−3, which is same as that with simulation using Network I with the maximum air injection of 7.5 kPa. During the two sets of tests two types of simulations with 1-and 10-step pressures were performed, and Figs. 13 and 14 show the comparison of relative air permeability values.
As shown in Fig. 13, the air relative permeability of 1-step pressure test and that of the 10-step pressure test are almost identical during the first set of tests with the maximum capillary number of N ca = 6.1E−3; however, there are obvious differences between the relative permeabilities of the two sets of tests with the maximum capillary number of N ca = 2.2E−3 as presented in Fig. 14. From this it can be further concluded that the effect of pressure step on relative air permeability is not consistent but varies with the capillary number. The effect of pressure step on relative air permeability is noticeable for the flow with lower capillary numbers than that with higher capillary numbers. By comparing Figs. 9 and 13, it can be observed that for the flow with same maximum capillary number, the difference between the relative permeability of 1-step test and that of 10-step test is smaller for Network II than that for Network I. The main reason accounting for this result is that comparing the flow with same capillary number in the both networks, it is much easier to develop continuous air channels in a small-pore network than in a large-pore network. As a result, the effect of pressure step on the air permeability is less evident in the small-pore networks than that in large-pore networks. Figure 15 shows both relative air and water permeability values from three experiments using Network I with the maximum air injection pressure of 20 kPa. From Fig. 15, it can be observed that there is higher scatter of relative water permeability than that of relative air permeability. It can be also found from Fig. 15 that relative water permeability for the flow with smaller number of steps (high pressure step) is smaller than that with larger number of steps (low pressure step). The main reason accounting for this is that for the flow with lower pressure step, it is much easier for water phase to get trapped and isolated which leads to lower water permeability.
To check the impact of network geometrical properties on relative water permeability, the results of the simulation performed using Network II shown in Fig. 16 can be compared with that shown in Fig. 15. From Fig. 16, the same conclusion as that in Fig. 15 about the effects of pressure steps on water relative permeability can be drawn. As shown in Fig. 16, relative water permeability in 1-step simulation is much higher than that in 10-step simulation. By  Fig. 16 Air-water relative permeability in Network II with N ca = 6.1E−3 comparing Figs. 15 and 16, it can be further observed that the difference between the 1-and 10-step tests in Network II is much smaller than that in Network I. In addition, as can be observed from both Figs. 15 and 16 that, the difference between relative water permeability values are much larger than those between air permeability values. As shown in Fig. 16, the two relative air permeability curves almost overlapped; however, the two relative water permeability curves were quite different. The main reason accounting for this is the different flow pattern transition between air and water. During the air-water two-phase flow, the air is transferred from discrete phase to continuous phase while the water is changed from continuous to discrete phase.

Summary and Conclusions
In this study, a dynamic two-phase flow model was developed to investigate the transport properties including relative permeability and capillary pressure of two-phase flow during air sparging. The void space of the porous medium was characterized using a 3D pore network consisting of pore bodies and bonds. Based on the 3D pore network, a dynamic two-phase flow model was developed to study temporal evolution of air pressure, air flow rate, and time dependence of spatial air distribution in porous media during air sparging. Then, the developed dynamic air-water two-phase flow model was applied to study dynamic flow during air sparging tests with multistep air injection pressures.
During multistep air injection pressure test, both the dynamic capillary pressure curve and relative air/water permeability curves were obtained and studied. As to the dynamic capillary pressure curve, it was always above the steady capillary pressure curve. This conclusion is consistent with experimental observations of other researchers. As to relative air permeability, relative air permeability with saturation curve can be divided into three parts. Part I is for water saturation greater than 0.75. This part corresponds to the initial stage of air injection without air breakthroughs. Part II is the part with water saturation between 0.35 and 0.75. In this part, the air-water two-phase flow is in transient state and it can be clearly observed that relative air permeability for flows with smaller pressure steps are much lower than that with larger pressure steps. The third part is the part with water saturation less than 0.35 which corresponds to steady state flow. In this part, the differences between relative permeabilities of the tests with different pressure steps become much smaller than that in the second part of the relative permeability curve. It was concluded from the above simulations that relative air permeability is affected by the air injection pressure step. Relative air permeability is higher for the flow with larger pressure steps. However, the effect of pressure step on relative air permeability is not consistent for flow with all different capillary numbers and it is weakened with high capillary number flows. Compared with the relative air permeability values, the relative water permeability values have higher scatter. It was also observed that relative water permeability for flows with larger pressure steps is higher than that with smaller pressure steps.
It was further observed that the dynamic effects on the relative permeability curve are more apparent in pore networks with larger pore sizes than those with smaller pore sizes. In addition, the effect of pore size on relative water permeability is more noticeable than that on relative air permeability.
Based on the proposed approach of studying the dynamic flow properties during twophase flow, more accurate constitutive equations of simulating air sparging at macroscale including capillary pressure-saturation (P c -S) and relative permeability-saturation (K r -S) relationships can be obtained for different soil structures and for flow with different capillary numbers. This research provides substantial insight to improve the accuracy of simulating multistep air sparging tests in granular soils. However, the void structure of soil is characterized by a 3D network of connected spherical pores, the proposed approach and the observations are not applicable to clay soils with non-spherical particles.