High-resolution shock-capturing numerical simulations of three-phase immiscible fluids from the unsaturated to the saturated zone

Numerical modeling of immiscible contaminant fluid flow in unsaturated and saturated porous aquifers is of great importance in many scientific fields to properly manage groundwater resources. We present a high-resolution numerical model that simulates three-phase immiscible fluid flow in both unsaturated and saturated zone in a porous aquifer. We use coupled conserved mass equations for each phase and study the dynamics of a multiphase fluid flow as a function of saturation, capillary pressure, permeability, and porosity of the different phases, initial and boundary conditions. To deal with the sharp front originated from the partial differential equations’ nonlinearity and accurately propagate the sharp front of the fluid component, we use a high-resolution shock-capturing method to treat discontinuities due to capillary pressure and permeabilities that depend on the saturation of the three different phases. The main approach to the problem’s numerical solution is based on (full) explicit evolution of the discretized (in-space) variables. Since explicit methods require the time step to be sufficiently small, this condition is very restrictive, particularly for long-time integrations. With the increased computational speed and capacity of today’s multicore computer, it is possible to simulate in detail contaminants’ fate flow using high-performance computing.


Flow governing equations and mathematical setup
The equation that describes a multiphase fluid flow in a porous medium is calculated using the conservation equation for the mass and momentum for each fluid phase. For three-phase fluid flow in terms of nonaqueous (n), water (w), and air (a) is given by where α = (n, w, a) and, Darcy's velocity for each phase is given by, , φ is the porosity, S α is the dimensionless volumetric saturation which is the fraction of the total available volume (inside the porous media) occupied by that component. These equations include 16 dependent variables in the general case, ρ w , ρ n , ρ a , S w , S n , S a , φ, k w , k n , k a , p w , p n , p a , µ w , µ n , µ a , and the three terms, q w , q n , q a . As a consequence of this, we need 13 additional relationships to obtain a solution for the system (3). One relationship corresponds to the sum of the volumetric saturations equal to one, www.nature.com/scientificreports/ For a three-phase fluid flow, we need two capillary pressures. In our formulation, we used the capillary pressure for the air-water phase and the capillary pressure for the air-nonaqueous phase, respectively indeed, p a , p w and p n are not independent and we have that the capillary pressure of the nonaqueous-water phase is given by, p cnw = (p n − p w ) = (p caw − p can ).
In this paper, we developed a general formulation for three-phases immiscible fluids flow . We consider nonzero pressure gradients (see Ref. 25,26 where the air gradient pressure is assumed negligible). In general, the capillary pressure is a function of saturation. From Eq. (5) we have two relationships. The porosity φ is a function of the pressure (and is one relationship). Densities and viscosities are funcions of phase pressures (and correspond to six relationships). The effective permeability tensor can be written as k ij w = k rw k ij , k ij n = k rn k ij , k ij a = k ra k ij , where k rw , k rn , and k ra are the dimensionless relative permeabilities for the phases w, n, a respectively. k ij is the absolute permeability, that depends on the properties of the porous medium. They are function of saturations (and are three relationships).
Finally, we can rewrite Eq. (3) using the capillary pressures (5) and the relative permeabilities, in terms of the variables p a , S w , S n , and S a .
The system of Eqs. (4,6) need to specify the five functions, k rα = k rα (S a , S n , S w )) , p can = p can (S a , S n , S w ) , and p caw = p caw (S a , S n , S w ) . We will be back to this point at the end of this section. But it is worth noticing that choices of solutions correspond to the possible different porous medium. But the numerical solution and the method we apply here is not affected by any particular option.
To complete the mathematical description, we need the expressions for q α , the boundary conditions, and the initial conditions. Before doing that, let us write the rock compressibility c R as a function of the porosity and pressure, and consider the Taylor expansion up to order one in c R where we get a linear approximation for the porosity Let us define the product of the porosity φ and the saturation for each phase as, σ w ≡ φS w , σ n ≡ φS n , σ a ≡ φS a , then Eqs. (4) and (8) can be written as, where φ 0 is the porosity at the reference pressure p 0 , which we consider to be the atmospheric pressure. Now we can write the left-hand side of Eq. (8) as a function of the pressure, where p will be associated to p a and the system of PDEs (6) becomes and ( α = w, n, a): where the do not depend on the spatial derivative of the saturation and the (4) S w + S n + S a = 1.
The PDEs system to be numerically resolved is composed by (11,12) and variables p, σ w , σ n , σ a , together with the functional form of the relative permeabilities (15) and capillary pressures (19) (see next subsection). We separated the right-hand side of the system (12) in advection (hyperbolic) PDEs, the one that also depends on the spatial derivative of the pressure p and the gravity g [Eq. (13)], and parabolic PDEs, the one proportional to the capillary pressures p c [Eq. (14)] in the variable saturation. These two pieces will be treated differently when we used numerical methods. The hyperbolic PDEs is responsible for the shock formation when the flow passes through a discontinuity and has to be treated using a mass-conservative numerical method.
Permeabilities and capillary pressures model. Now that we have the complete system of three-phase flow equations, we need to define the functional form for the relative permeabilities and capillary pressures. The relative permeabilities for three phases has been extended from the two-phase expressions 27 and are given by, where the total effective liquid saturation, S et , is defined as and S wir is the irreducible wetting phase saturation. For the capillary pressure, we use the van Genuchten model 19 where the effective saturation, S e , has the following form, and α and n are model parameters. Defining m = 1 − 1 n , from which we get, n = 1 1−m , and resolving for p c we obtain where p c0 = α −1 is the capillary pressure at S e = 0 . The water content S w is entirely determined by the capillary pressure between the NAPL and water, p cnw = p cnw (S w ) . The total liquid content S t = S w + S n (or alternatively, the air saturation S a = 1 − S t ) is completely determined by the capillary pressure between air and NAPL, p can = p can (S t ) . From the definition of capillary pressure, it follows that the capillary pressure between air and water is not an independent quantity but is given by, p caw = p can + p cnw . As a consequence of these wettability assumptions, the local saturated distribution of a three-phase water-NAPL-air system is determined by twophase capillary pressure-saturation relationships for air-water and air-NAPL. The capillary pressures are given by,

High performance computing and CactusHydro
Most of the multiphase/three-phase groundwater flow equation solvers employ either FDM or finite element method (FEM) in the three coordinate space. Hybrid combinations are also possible 9 . Other used methods are: transversal methods of lines for the numerical modeling of vertical infiltration into the vadose zone 28,29 , or a WENO based method of lines 30 . Although there were very few attempts to use a supercomputer to simulate subsurface solute transport, the situation has changed in the last ten years.. New codes use HPC, for example, PFLOTRAN 31,32 , Parflow Hydrologic model 33,34 , RichardsFOAM 35 , CATHY (CATchment HYdrology) 36 . But it should be noticed that none of them systematically use explicit forward in time evolutions and conservative highresolution shock capturing methods to efficiently study the dynamics of multicomponent flow in porous soil.
Most of these codes are written in Fortran language (as, for example, PFLOTRAN), while few are written in C++/C. Both programming languages are efficient for large computation problems and are well supported compared with other languages. Some codes, as for example, RichardsFOAM, only treat Richard's equation and study the water fluxes at the watershed scale. Most of the codes used FDM, FEM, and less often, FVM. Some codes are commercial software such as, COMSOL Multiphysics 37 , FEFLOW (Finite Element subsurface FLOW system) 38 . It is expected that, if properly applied, different spatial discretization methods should not change results too much when other numerical methods are used.
Our interest is focused on the application of software and codes to large-scale hydrological modeling, similar to CATHY, ParFlow, PFLOTRAN, etc., and use sufficiently finer mesh size, especially in the vertical direction or refinement in the area where the contaminant or the principal problem is evolved. This paper presents a new code (CactusHydro) that performs high-resolution numerical simulations in multiphase flow from the vadose zone to the saturated one in a porous medium. CactusHydro is a code written in C/Python language and is  20,23,39 , an open-source software framework for developing parallel highperformance simulations codes. Data are evolved on a cartesian mesh with six refinement level using Carpet 40,41 . We discretize the domain problem using the FDM and implement the explicit (forward) Euler method for the space-time evolution of the system.

Numerical methods
The strong nonlinearities in the previous system of PDEs are represented by the relationship between the permeability-saturation, which is responsible for creating the shocks, and capillary pressure-saturation of each phase. This situation has to be described accurately in time (and space) and needs high-resolution and refinement at a small scale (vertical and horizontal) up to a large scale. The crucial part to notice is the fact that the dominant piece of the multiphase flow PDEs that has to do with the water and nonaqueous phase is dominated by the hyperbolic part (the one proportional to gravity and the gradient of the pressure, rather than the elliptic part of the equation (proportional to the capillary pressure). We will use explicit methods to resolve the part of the PDEs equation that has a hyperbolic structure and resolve it in an explicit way to eliminate the oscillations/shocks that are caused due to the term proportional to the gravity constant. When introducing the formal properties of numerical solutions methods of PDEs associated with conservation laws, two fundamental theorems are underlining the importance of using a conservative formulation. The first one is due to Lax, and Wendorff 42 and the second by Hou and LeFloch 43 . Namely, that conservative numerical schemes, if convergent, do converge to the weak solution of the problem. The second theorem states that nonconservative numerical schemes do not converge to the correct solution if a shock wave (or discontinuity) is present in the flow. These two theorems state that if a conservative formulation is used, then we are guaranteed that the numerical solution will converge to the correct one, while if a conservative formulation is not used, we are guaranteed to converge to the incorrect solution in the likely event in which the flow develops a discontinuity.
The main approach to the numerical solution of the problem is based on a class of semi-discrete methods (Ordinary differential equation (ODE) in time and discretized in space) conservative flux methods of the type discussed in Kurganov and Tadmor (KT) 22 using upwind fluxes. This approach is based on (full) explicit in the time evolution of the discretized equation using the so-called methods of lines (MoL) for the time evolution of the discretized (in space) variables. This technique allows to overcome the major limitation of the solution methods used before and to ensure mass conservation of the various components easily and accurately propagate the sharp front of the fluid component. This has the main drawback that the size of the time step that must be used to ensure convergence is much smaller than the one that can be used with an implicit scheme or semi-implicit one, like the one used in IMPES simulations 21 . In particular, to deal with the sharp front originated from the nonlinearity of the equations, we use a high-resolution shock-capturing (HRSC) method to treat discontinuities due to the capillary pressure (air-nonaqueous), (water-nonaqueous) and the permeabilities that depend on the saturation of the three different phases.
The conservation form (second-order scheme) in the 1D-case (straightforwardly generalized to 2D and 3D) is given by 22 , with the numerical flux, and the discretized variable u j (t) are approximated TVD limited piecewise linear (discontinuous at the cell interface) that assume the values u − j+1/2 and u + j+1/2 at the boundary between the cell j and j + 1 and a j+1/2 (t) is the maximum of the spectral radius of F ′ (u) . We can now use the crucial simplification (that it is now accurate at firstorder in space) that holds true when F(u) does not change sign (as it is for our equations) for different values of u(t) at the interface. Namely we assume: The HRSC methods possess the following properties: sharp resolution of discontinuities without considerable smearing; at least second-order of accuracy on smooth parts of the solution; absence of spurious oscillations in the solution; convergence to the "weak" solution as the grid is refined; no use of artificial-viscosity terms. The method belongs to the class of Monotonic Upstream-centered Scheme for Conservation Law (MUSCL) suggested by van Leer in 1973, and the KT scheme is second-order accurate in space. To achieve second-order accuracy, one needs to use analytical information on the form of the flux. However, as in this work, the use of the first-order upwind formula for the fluxes, and the minmod flux limiter, only the point values of the flux are required. This is of great advantage since it can use tabulated values for permeabilities. The only penalty is that it is first-order not only at the discontinuities but also over the whole simulation grid. But one has to keep in mind that any methods will be first-order at the physical discontinuities, and the problems we target to study involve following real-discontinuities.

Results
In this study, some tests are performed to verify the accuracy and reliability of the HRSC method and the Cac-tusHydro code. The results are compared with analytical models and experimental data results.

The inviscid Burgers' equation.
To show the effectiveness of the method, we consider the inviscid Burgers' hyperbolic equation, which represents a model for nonlinear wave propagation and how the use of a (simplified) HRSC method is capable of dealing with sharp fronts (or discontinuities) and to evolve them correctly as well as to follow the shock-formation and the rarefaction. Figure 1 shows (in black) the initial conditions u(x, 0) = cos(x) for |x| < π/4 and its evolution, moving from the left to the right, using the employed HRSC method at different values of t (colored lines). One may observe the formation of the jump discontinuity (shock). It is also observed the formation of a rarefaction when the propagation velocity ( v(x, t) = 2u(x, t) ) grows with respect to x. When it decreases, there is the formation of a shock. When the shock is completely formed (purple color), the method correctly evolves the shock front. This technique allows overcoming the major limitation of the solution methods described in the previous section, easily ensuring mass-conservation of the various components, and accurately propagating the sharp front of the fluid component. We also considered the Courant-Friedrichs-Lewy (CFL) restrictions in our simulations. In particular, at any fixed spatial resolution, one has to select a time step size so that the evolution is convergent in time. We checked that the used time step size at any resolution fulfills the CFL conditions. For example, in this test case of the Burgers' equation, we set, to accurately follow the formation of the shock discontinuity the CFL factor to be 0.001.

Linear waterflood.
To verify that the numerical code is properly solving the governing equations and dealing with sharp front coming from discontinuities, it is important to compare our numerical results with analytical results from simplified examples. An analytical solution for a one-dimensional problem is given by the Buckley-Leverett model 44 which estimates the advance of a fluid displacement front in an immiscible displacement process in the horizontal direction (zero gravity force) and zero capillary pressure.
The one-dimensional test is performed using full 3D-code evolution and imposing constant values of all variables along with two directions with constant boundary conditions achieving an effective 1D evolution. Along this direction is keeping fixed pressure difference and the evolved direction in x is L = 25 m long with a cell resolution of dx = 0.125 m . We consider a pipe in which the initial pressure was set to be a fixed gradient. This is not a profile of a stationary two-phase fluid flow with constant flux. Nevertheless, in a few second it stabilizes to the right initial conditions. The data for this problem is given in Table 1. Relative permeabilities were taken from Eq. (15) as a function of water saturation. Figure 2, shows the behavior of the product of saturation and the porosity ( σ alpha ) as a function of the distance for different values of the time, for a two-phase nonaqueous-water flood (with σ a = 0 ), and using CactusHydro code. Notice that σ a + σ w + σ n = φ which is equal to 0.30 (see Table 1). We also show the comparison with a double resolution cell, dx = 0.250 m which perfectly match with the one corresponsing to dx = 0.125 m . This also shows very good convergence for the results of the numerical simulations.We also compare the analytical result coming from the Buckley-Leverett model using the same parameters in the theoretical expression. Notice how the numerical results perfectly match the theoretical expectation even in the discontinuity zone. This is www.nature.com/scientificreports/ due to the use of a mass-conservative method implemented in the code. Relative permeabilities were chosen as a function of water saturation, as appears in Table 1, where the zero gravity and zero capillary pressure were imposed. From this Fig. 2   Three-dimensional unsaturated three-phase fluid flow with gravity. An additional numerical test code, before showing a 3D example, is provided in the Supplementary Information. Our code results are compared with a two-dimensional numerical unsaturated-saturated water flow model 9 , together with the sand tank experimental data conducted by Koichi et al. 24 . See Table S1 for a list of parameters used in the 2D example, and Figure S1 for a comparison between the water table elevations measured in a sand tank experiment 24 and the water table elevations calculated with CactusHydro as a function of the distance at different times.. Now that we have shown some tests, we can apply CactusHydro to a specific 3D example. In this example, we show the impact of an undetected leak of NAPL phase liquid onto a surface and, as a consequence of the gravity, moves from the unsaturated to the saturated zone. We consider a three-phase fluid case in which the contaminant is a light nonaqueous phase liquid (LNAPL) for which the density is less than the water. See Table 2 for details on the parameters. It provides a summary of the physical properties of the porous medium and fluids, initial conditions, and numerical model discretization data. The relative permeabilities and capillary pressure were taken from Eqs. (15,19) which consider three-phase fluid flow relations. Figure 4 shows an example result of a numerical simulation of three-phase fluid flow (water, LNAPL, air) performed using a spatial grid resolution of dx = 0.5 m and a total grid dimension of 80 m × 30 m and 20 m in the z-x axis. The simulation has been successfully performed using CactusHydro on a parallel cluster dedicating to HPC at Parma University (Parma).
Initially, there is a constant leak of NAPL situated on top of this figure (left-hand side), in the unsaturated zone composed by a porous medium and a contaminant density situated at z = 10 meters. The water table is at z = 0 coordinate, and there is a gradient of gravity of 15 degrees. The LNAPL contaminant flows downward saturated zone due to the gravitational force. It then arrives at the water table, where a capillary pressure between air-contaminant and contaminant-water is present. The contaminant has a density lighter than the water. It remains around the water table while part of it goes in the left direction due to the pressure gradient (bottom, left-hand side). On the other column at the right-hand side, the same situation is viewed in the z-y plane. Notice the effect of the capillary pressure between contaminant and water, which causes a movement of the water flow in their vicinity. The right-hand side of this figure shows the z-y axis where there is a zero-gravity component effect and thus, no privileged direction in the y axes. The transient numerical simulation shows the behavior of this contaminant for a period of time of around five days, although it is possible to go further in the simulation.

Conclusions
In this paper we have presented high-resolution 3D numerical results of three-phase immiscible fluids flow from the unsaturated to the saturated zone in a porous medium. They were obtained using a new code, CactusHydro, based on the Cactus toolkit (an open-source that provides a range of parallel computational capabilities). The governing equations are coupled general Richards' equations for each phase (water, NAPL, and air). We investigated the three-phase immiscible fluid flow dynamics as a function of the saturation, capillary pressure, relative permeability of the different phases, using different initial and boundary conditions.
To deal with the sharp front originated from the partial differential equations' nonlinearity and accurately propagate the fluid component sharp front, we use a high-resolution shock-capturing (HRSC) method to treat discontinuities due to capillary pressures and permeabilities that depend on the saturation of the three different phases. The innovative method considered in CactusHydro is the (first-order accurate) upwind fluxes in the conservative-flux methods discussed in 22 . This approach is based on (full) explicit in the time evolution of the discretized equation using the so-called methods-of-lines (MoL) for the time evolution of the discretized Table 2. Data used in the 3D transient three-phase fluid flow, with an LNAPL released in the unsaturated porous medium zone that goes downward directed to the saturated aquifer. Data results from numerical simulations using CactusHydro.

Parameter Value
Absolute permeability, k 4.14 × 10 −10 m 2 We compare our numerical results with known analytical models such as the Burgers' equation and the Buckley-Leverett model. An additional numerical test code is the comparison with a two-dimensional unsaturated-saturated water flow model, together with the sand tank experimental data conducted by Koichi et al. 24 (see Supplementary Information). Lastly, we consider a simple three-phase fluid flow case in which a contaminant www.nature.com/scientificreports/ (LNAPL) leak onto the surface and, as a consequence of the gravity, moves from the unsaturated to the saturated zone. These techniques allow overcoming the previous section's solution methods' major limitation and ensuring the various components' mass-conservation. Also accurately propagating the sharp front of the fluid component. This has a drawback: the size of the time step that must be used to ensure convergence is much smaller than the one used with an implicit scheme or semi-implicit one like the one used in IMPES. The numerical results are very encouraging and show the robust capabilities of this code (and methods). We plan to apply it to more realistic cases in future works.