A numerical study of two-phase flow models with dynamic capillary pressure and hysteresis

Saturation overshoot and pressure overshoot are studied by incorporating dynamic capillary pressure, capillary pressure hysteresis and hysteretic dynamic coefficient with a traditional fractional flow equation in one dimension. Using the method of lines, the discretizations are constructed by applying Castillo-Grone's mimetic operators in the space direction and explicit trapezoidal integrator in the time direction. Convergence tests and conservation property of the schemes are presented. Computed profiles capture both the saturation overshoot and pressure overshoot phenomena. Comparisons between numerical results and experiments illustrate the effectiveness and different features of the models.


Introduction
Water infiltrating into initially dry sandy porous media has been shown to produce saturation overshoot and pressure overshoot in Selker et al. (1992); Shiozawa and Fujimaki (2004); DiCarlo (2004DiCarlo ( , 2007. Eliassi and Glass (2001) and  have demonstrated that the traditional Richards equation is unable to describe saturation overshoot. To describe the non-monotonic behaviour, various extensions to the Richards equation have been investigated. Eliassi and Glass (2001) studied three additional forms referred to as hypodiffusive form, hyperbolic form and mixed form, saturation overshoot are obtained by using the hypodiffusive form in Eliassi and Glass (2003). DiCarlo et al. (2008) studied a non-monotonic capillary pressure-saturation relationship and a second order hyperbolic term, but they mentioned these extensions need a regularization term to produce a unique solution. Cueto-Felgueroso and Juanes (2009) explained the formation of gravity fingers during water infiltration in soil by introducing a fourth-order term to Richards equation. The tip and tail saturations after their phase model have good agreement with the experiments in DiCarlo (2004).  obtained non-monotonic saturation profiles by supplementing the Richards equation with a non-equilibrium capillary pressure-saturation relationship, as well as including hysteretic effects. Later, Chapwanya and Stockie (2010) studied gravity-driven fingering instabilities based on the work of , their results demonstrate that the non-equilibrium Richards equation is capable of reproducing realistic fingering flows for a wide range of physically relevant parameters.
Besides extensions to the Richards equation, other approaches to characterize the saturation overshoot have also been proposed. Such as the generalized theory by introducing percolating and non-percolating fluid phases into the traditional mathematical model Hilfer and Besserer (2000); Hilfer et al. (2012); Doster et al. (2010), fractional flow approach DiCarlo et al. (2012) and moment analysis Xiong et al. (2012).
Among all the proposed theories and models, the dynamic (or non-equilibrium) capillary pressure relationship proposed by Stauffer (1978), Hassanizadeh and Gray (1993), Kalaydjian et al. (1992) has received much attention. Results on travelling wave solutions, global existence, phase plane analysis and uniqueness of weak solutions are given in Cuesta et al. (2000); Van Duijn et al. (2007; Mikelić (2010); Spayd and Shearer (2011); Cao and Pop (2015). In order to provide accurate simulations, several numerical methods have been proposed in literature, including a finite difference method with minmod slope limiter in Van Duijn et al. (2007), a cell-centered finite difference method and a locally conservative Eulerian-Lagrangian method in Peszynska and Yi (2008), Godunov-type staggered central schemes in Wang and Kao (2013), two semi-implicit schemes based on equivalent reformulations in , an adaptive moving mesh method in Zegeling (2015) and the fast explicit operator splitting method in Kao et al. (2015).
In the previous studies of saturation overshoot, the dynamic capillary pressure model with a constant dynamic coefficient usually brought oscillations behind the drainage front DiCarlo (2005); Sander et al. (2008); van . The hysteretic non-equilibrium model proposed in Beliaev and Hassanizadeh (2001) postulates that the dynamic capillary effects are significant only outside the main hysteresis loop. Following this idea,  adopted a saturation and pressure dependent dynamic coefficient τ = τ 0 s P ′ w (s)(p 0 − p) γ + . In this treatment, the Mualem hysteresis model was restricted only to the two-stage wetting-drainage process: trajectories for the wetting stage was located within H w (domain above the main wetting curve), while trajectories for the drainage stage was limited to H 0 (main hysteresis loop region). In the wetting stage τ 0 s = τ w and in the drainage stage τ 0 s = 0 (in the numerical simulation the value was a small constant = 10 −3 ). Recently, the hysteretic dynamic capillary pressure effect has been reported in Sakaki et al. (2010). Mirzaei and Das (2013) shows that the dynamic effect in the relationship between capillary pressure and saturation is hysteretic in nature. In this contribution, we consider the hysteresis effects in capillary pressure and study the saturation overshoot and pressure overshoot by adding dynamic capillary pressure, capillary pressure hysteresis and hysteretic dynamic coefficient to a traditional fractional flow equation.
Two-phase flow models in porous media usually consist of coupled, nonlinear partial differential equations. Many reliable discretizations have been proposed for the Richards equation, for instance, the Galerkin finite elements in Arbogast et al. (1993), the multipoint flux approximation in Klausen et al. (2008) and so on. The space and time discretizations usually lead to a large system of nonlinear equations, which makes the numerical simulation of two-phase flow a challenging task. Some popular methods used for the Richards equation are the Newton method Radu et al. (2006), Picard method Zarba et al. (1990) or L-scheme Radu et al. (2015), for an extensive review we refer to List and Radu (2016). The appearance of the dynamic capillary pressure term in the two-phase flow model adds additional difficulty to the numerical treatment. In Sander et al. (2008), a reformulation of the non-equilibrium two-phase flow equation, which consists of an elliptic equation and an ODE, is shown to be effective for numerical simulations,  also presented a reliable and efficient semi-implicit scheme for a similar form. In this paper, we present our schemes based on this reformulation. Because of the conservation property and easy implementation of the Castillo Grone's mimetic (CGM) operators, we apply the mimetic finite difference method to the elliptic equation in the space direction, then integrate the system in the time direction with the implicit trapezoidal rule.
The rest of the paper is organized as follows. In section 2, we first derive the traditional equation for onedimensional two-phase flow, and then we present the extended models by incorporating the dynamic capillary pressure term and hysteresis effects. Section 3 is devoted to presenting the CGM operators. In section 4, we apply the CGM operators in the space direction and an implicit trapezoidal integrator in the time direction to discretize the system. In section 5, numerical experiments are carried out to show the effectiveness and reliability of the extended models. Section 6 summarizes the conclusions.

Mathematical Models
Two-phase flow in porous media can be characterized by the saturation and pressure in each phase. The saturation in each phase is defined as the fraction of the pore volume occupied by the phase and is denoted as S α where α = w, n is an index for wetting and non-wetting phases. For the derivation of the fractional flow equation we refer to Hilfer and Steinle (2014). Let the gravity act in the positive x-direction, for each phase, the mass conservation law is represented by the equation where φ is the porosity of the porous medium, ρ α , v α and F α are the density, volumetric velocity and source of each phase. In Darcy scale, the balance of momentum of each phase is given by the Darcy's law where K is the intrinsic permeability of the porous medium, g is the gravitational acceleration constant, k rα , µ α , λ α = krαK µα and p α are the relative permeability function, viscosity, mobility and pressure of phase α, respectively.
For the two-phase system the following constitutive relation holds: Then k rn and λ n can be written as functions of S w . The total velocity is given by Using Eqs.
(2), (3) and (4) the velocity of the wetting phase can be expressed in terms of the phase mobilities, total velocity and phases pressure difference as Assuming φ and temperature are constant, the phases are incompressible, by substituting Eq. (5) into Eq. (1) for the wetting phase, we obtain a nonlinear equation for the wetting phase In DiCarlo (2004) the experiments were conducted in thin tubes, and water was injected into the tubes with different flux rates. Since the medium is homogeneous and the initial saturation is constant, the experiments are viewed as one-dimensional. Reminding that our aim is to simulate these experiments, we set the source term F w = 0 and consider a flux boundary condition at x L and a Dirichlet boundary condition at x R . Let f (S w ) = λw λw +λn = λw λT , then we have with boundary conditions where q is the flux value used in DiCarlo (2004DiCarlo ( , 2007, S R w is the initial water saturation at x R . Integrating Eq.
in Section 4 we will test the conservation property of the numerical schemes using this formula.

Dynamic capillary pressure model
Under equilibrium conditions, traditional models suggest the difference in phases pressure is equal to the capillary pressure. In the microscale, the capillary pressure is defined as the interfacial tension between two phases and in Darcy scale, it is usually given as a function of the wetting phase saturation: For non-equilibrium conditions, Stauffer (1978), Hassanizadeh and Gray (1990), Kalaydjian et al. (1992) proposed that the phases pressure difference can be written as a function of the capillary pressure under equilibrium condition minus the product of the saturation rate with a dynamic coefficient τ [Pa s]: The parameter τ is also known as damping coefficient and may still be a function of saturation Joekar-Niasar et al. (2010). Adding the dynamic capillary pressure term to Eq. (7) we obtain In the following we mark this model as Model 1.

Play-type capillary pressure hysteresis model
Many studies Morrow et al. (1965); Jerauld and Salter (1990) in recent decades have shown non-uniqueness in the relationship between capillary pressure and saturation, which can depend both on the history of flow displacement and the rate of change of saturation. The dependency of p c -S w on the history of flow is known as capillary pressure hysteresis. Displacement of flow differentiates between drainage and imbibition. The process of drainage describes when the nonwetting phase displaces the wetting phase. Vice versa, imbibition describes the process when the wetting phase displaces the nonwetting phase. In general, for a given saturation S w , p c can lie anywhere within the primary drainage curve P dr c and the primary imbibition curve P im c , depending on the saturation history. Some typical plots of hysteretic capillary pressure curves are presented in Fig. 1. In Parlange (1976); Beliaev and Hassanizadeh (2001); Brokate et al. (2012) different kinds of hysteresis models have been discussed for two-phase flows, in our work we adopt the play-type hysteresis model presented in Brokate et al. (2012).
Assume p c depends only on S w and this relationship is described by the hysteresis operator Note that P hyst c operates on S w as a function of time. In the drainage process, when S w decreases, p c follows the drainage pressure-saturation curve P dr c (S w ). In the imbibition process, when S w increases, p c follows the  Table 2; blue, dashed blue and dash-dotted blue lines illustrate an imbibition hysteresis curve, a drainage hysteresis curve and a play-type hysteresis curve.
imbibition pressure-saturation curve P im c (S w ). In this hysteresis model, between the drainage and imbibition curves, p c and S w evolve as where β is the opposite slope of the hysteresis curve with dimension [Pa]. β is usually chosen to be very large, which means the hysteresis curve is very steep and the hysteresis pressure can move quickly from one equilibrium curve to another when the saturation direction changes. This hysteretic capillary pressure curve is illustrated in Fig. 1 (dashed-dotted blue line).
Since P hyst c acts on the history of S w , it is not possible to compute p c at a given time from S w at that time alone. Consider the system after time discretization, denote p c and S w from the previous time step as p n−1 c and S n−1 w , respectively. The discrete form of Eq. (14) is The algorithm for computing p c is as follows In the numerical simulations, we denote the above algorithm as p n c = P hyst c (S n w ). Combine capillary pressure hysteresis with the dynamic capillary pressure (11) we obtain Substituting Eq. (16) into Eq. (7) we get a model with dynamic capillary pressure effect and play-type capillary pressure hysteresis, in the following this model is marked as Model 2.
2.3 Hysteretic dynamic capillary pressure model with play-type capillary pressure hysteresis Sakaki et al. (2010) conducted a series of experiments to measure values of the dynamic capillary coefficient τ for a porous medium. Their result suggests that τ is hysteretic. Mirzaei and Das (2013) investigated the hysteretic behaviour between τ and S w . The experiments demonstrate that the value of τ for imbibition is generally larger as compared to the τ value for drainage at the same saturation. Thus it is reasonable to introduce hysteresis in the τ -S w relationship. We assume in the imbibition process τ = τ im , in the hysteresis process τ decreases from τ im to τ dr and at the tail the dynamic coefficient is τ dr . To our best knowledge, the ratio τ dr /τ im has not been investigated in the literature. In this work, we follow Van Duijn et al. (2007) to show the influence of τ by using the travelling ansatz. Eq. (12) can be rewritten as where the flux F (S w ) and the capillary induced diffusion Cuesta et al. (2006) H(S w ) are given by In order to find a traveling wave solution for Eq. (17), we introduce the new variable η = x−st. Substituting S w (η) into (17) results in a third order ordinary differential equation (ODE) where prime denotes differentiation with respect to η. This equation is to be solved subject to the boundary conditions at infinities, Integrating Eq. (17) over (η, ∞) and assuming yields the second-order ODE: with s determined by the Rankine-Hugoniot condition When gravity is included into the flux function F (S w ), with different values of v T , F (S w ) may be non-monotone. For simplicity, we only consider (S L w , S R w ) pairs that satisfy s > 0. Next we write Eq. (22) as a first order system of ODEs: Let S α w be the unique root of the equation where S 0 w is the initial water saturation ahead of the wetting front. The Jacobian of (24) reads and has eigenvalues For the case where saturation plateau appears, consider a traveling wave connecting S L w = S B w (equilibrium boundary saturation) and S R w =S P w (plateau saturation) with wave speed Then we can prove (S B w , 0) is an equilibrium of system (24). Using (27) When only saturation overshoot appears, consider a travelling wave connecting From the analysis above, we can conclude that when τ < τ s , saturation oscillation will not appear at the drainage front, which means ∂S w /∂t → 0 − . Therefore, the phase pressure difference p n − p w will tend to P dr c at equilibrium.
Denoting the hysteretic dynamic coefficient as τ hyst , since τ hyst may possibly due to the hysteresis in the retention curve Sakaki et al. (2010), for simplicity, we introduce τ hyst as the linear interpolation between τ dr and τ im by utilizing the capillary pressure hysteresis, Substituting Eqs. (16) and (29) into Eq. (7) will result in a model with hysteretic dynamic capillary pressure and capillary pressure hysteresis. This model is marked as Model 3.

Numerical scheme
In this section we present the numerical scheme based on a reformulation of the non-equilibrium equation, the method of lines is then applied to this reformulation. Denoting p = p n − p w , Eq. (12) can be rewritten as Since Model 2 and Model 3 are incorporated with the capillary pressure hysteresis, we replace P c (S w ) in Eq. (30) by P hyst c (S w ) when solving these two models and replace τ by τ hyst when solving Model 3.

Castillo-Grone mimetic operators
To discretize Eq. (30) in the space direction we adopt the mimetic finite difference method which satisfies the discrete version of continuum conservation law. Castillo and Grone (2003) developed a set of mimetic operators knows as Castillo-Grone's mimetic (CGM) operators. CGM operators have been used in many fields, such as seismic studies Rojas et al. (2008), electrodynamics Runyan (2011) and image processing Bazan et al. (2011). Numerical results in these fields validate the high efficiency and reliability of the CGM operators.
The main features of CGM operators are that they preserve symmetry properties of the continuum and have overall high order accuracy. The CGM operators can be implemented as efficient as the standard finite difference schemes. Here, we briefly describe the CGM operators in one dimension as applied in this work. The CGM 2-D operators can be obtained by the Kronecker products of block matrices. For more details, see Castillo and Miranda (2013) and references therein.
In the one-dimensional situation the Green-Gauss-Stokes theorem reads The cell centers can be indexed as The cell nodes x i and cell centers x i+1/2 build up the uniform staggered grid. We discretize v at the cell nodes: and f at the cell centers and the boundary nodes: LetD (N +2)×(N +1) denote the CGM divergence operator and G (N +1)×(N +2) denote the CGM gradient operator, then the discrete version of the conservation law (31) reads where < x, y > A = y T Ax is the inner product, I is the (N + 2) × (N + 2) identity matrix, Q and P are the weight matrices forD and G respectively. The matrixB = QD + G T P embodies the global conservation requirement for the discrete conservation law and is called the boundary operator. Castillo and Yasuda (2005) presented the second-order divergence mimetic operator aŝ where For this second-order mimetic divergence matrixD, the weights matrix Q is the (N + 2) × (N + 2) identity matrix. The second order CGM gradient operator reads as with weight matrix 3/8 0 · · · · · · · · · · · · 0 . . . . . . 9/8 0 0 · · · · · · · · · · · · 0 3/8 Applying the discrete conservation law (34) with matrices D and G, the boundary operatorB can be written asB 0 0 · · · · · · 0 −1/8 1/8 0 0 0 · · · · · · 0 1/8 −1/8 0 0 0 · · · · · · 0 0 1 3.2 Mimetic discretizations for the two-phase flow equations Using the uniform staggered grid presented in Section 3.1, we discretize S w and p at the centers asS Then we use linear interpolation to get S w at the nodesŜ At the boundaryŜ w0 = S w0 ,Ŝ wN = S wN . Introducing coefficients matrix K ∈ R (N +1)×(N +1) with zero non-diagonal elements, the diagonal elements are given by In the numerical simulations, since fully implicit time discretizations allowing large time steps are preferred for solving long-time scale problems, thus we apply the implicit trapezoidal integration to (30) in the time direction. Discretizing (30) and boundary condition (8) with the CGM operators in the space direction, we obtain In (42),Î is a (N + 2) × (N + 2)-dimensional matrix withÎ 11 = 0,Î N +2,N +2 = 0,Î ii = 1, for i = 2, 3, · · · , N + 1, the matrices A, B and vector b(S n+1 w ) represent the boundary conditions. As a result of the flux boundary condition at x L and the Dirichlet boundary condition at x R , the only non-zero element in A is A N +2,N +2 = 1, and B differs fromB in the last three rows, 0 0 · · · · · · 0 0 0 1/8 −1/8 0 · · · · · · 0 0 0 −1/8 1/8 0 · · · · · · 0 0 0 The In order to solve (42) and (43) we apply the iteration method. By introducing the superscript l as an iteration counter, the algorithm for each time step is as follows: Remark: The application of MFD to partial differential equations constitutes an active filed of research Lipnikov et al. (2014). Formal analysis of MFD for the Richards equation can be achieved by combining the convergence results in Brezzi et al. (2005), with the equivalence between MFD and multipoint flux approximation (MPFA) established in Stephansen (2012) and the convergence proof of MPFA for Richards equation in Klausen et al. (2008). However, difficulties arise when establishing convergence for (30), because of the nonlinearity and the reformulation. In Section 4, we present numerical results to demonstrate the convergence of the method.    (1996) as well as the constants and Brooks-Corey models Brooks and Corey (1966) are listed in Table 1 and Table  2. DiCarlo observed that for the highest (q = 2.0e-3 [m s −1 ]) and lowest (q = 1.32e-07 [m s −1 ]) fluxes, the saturation profiles are monotonic with distance and no saturation overshoot is observed, while all of the intermediate fluxes exhibit saturation overshoots. In this section we will study the numerical behaviours of the three models presented in Section 2. In Eq. (7) when S w = 0 or S w = 1 the equation is degenerate. From Fig. 1 in DiCarlo (2007) we get the initial capillary pressure p 0 c ≈ 1600[Pa], using the Brooks-Corey capillary pressure model in Table 2 we can find for the imbibition process, when water saturation S w = 0.003, the Brooks-Corey capillary pressure p c (S w ) = 1566 [Pa]. So in the numerical simulations, we set S R w = 0.003 and S max w = 1 − 1.0e-03. The initial saturation is given by where S L w = 0.025, the initial phases pressure difference is p n − p w = P c (S w (x, 0)). The reason we set S L w > S R w is that, when S L w = 0.003, in Eq. (43) we have λ n (0.003)f (0.003) ≈ 6.7e-16, then the boundary saturation obtained by solving (30) will exceed 1; when S w is big enough, for example S L w = 0.025, we have λ n (0.025)f (0.025) ≈ 9.1e-13, the boundary saturation will not exceed 1, see Fig. 4(a). For the numerical simulations, S L w is small enough and will not influence the behaviour of the models. Before we carry out the numerical simulations, the parameters β and τ dr appear in Eq. (14) and Eq. (29) have to be decided. Here we choose β = 1.0e05, the ratio τ dr /τ im is set to be 0.2.
First, we test the accuracy of schemes (42) and (43). Since exact solutions of Eq. (12) are not known, the numerical solutions on fine grids are taken as reference solutions. For space and time accuracy tests, the reference grids are N = 1024, ∆t = 0.01 and N = 512, ∆t = 0.001, respectively. Setting S L w = 0.45, S R w = 0.3, q = 1.32e-04, τ = 4.0e03, fixing time or space step size, the L 2 errors and orders are obtained in Tables 3 and 4. Table 3 shows that in the space direction, schemes (42) and (43) are second order when applying to different models, and the L 2 errors of the three models are consistent. However, as a result of the hysteresis effects, in the time direction the L 2 errors increase when schemes (42) and (43) are applied to Model 2 and Model 3, also the convergence rates drop from two to about one. Setting S L w = 0.025, S R w = 0.003, q = 1.32e-04, τ = 4.0e03, the saturation and pressure profiles obtained by Model 1, 2 and 3 are presented in Fig. 2. The imbibition fronts obtained by the three models are similar, but the saturations and pressures at the plateaus and behind drainage fronts are different. For Model 1, the value of the plateau saturation is constant, behind the drainage front, oscillations appear and the phases  pressure difference follows the imbibition capillary pressure. For Model 2, the plateau saturations decrease a little from the imbibition front to the drainage front, behind the drainage front there is a slight oscillation in the saturation profile and the phase pressure difference moves to the imbibition capillary pressure. For Model 3, because of the hysteresis in the capillary pressure and the dynamic coefficient, no oscillation appears in the saturation profile. As a result, the pressure keeps constant and follows the drainage capillary pressure. The phases pressure difference-saturation curves are presented in Fig. 2(d). The result obtained by Model 3 shows similar behaviour as the measured data in Fig. 6 in DiCarlo (2007). Fig. 3 shows that schemes (42) and (43) preserve Eq. (9) with high accuracy for all three models. The evolutions of saturation and pressure at the left boundary are presented in Fig. 4. Fig. 4(a) shows the saturation obtained by Model 1 drops to the tail saturation after it reaches a high value, while the saturations of Model 2 and Model 3 keep the high values for a while. This phenomenon can be explained by the capillary pressure hysteresis in Model 2 and Model 3 as is shown in Fig. 4(b). From t = 0 to t ≈ 100 the hysteretic pressures in Model 2 and 3 increase, the pressure gradients keep the saturations stay at high values, when phases pressure differences reach the equilibrium drainage pressure, the pressure gradients vanish and the saturations move to the asymptotic tail values. The pressure curves obtained by Model 2 and Model 3 are also different. In Model 3, when pressure moves to the equilibrium drainage pressure, the hysteretic dynamic coefficient τ hyst also decreases from τ im to τ dr , thus keeps the saturation constant at the tail and the pressure stays at the equilibrium drainage pressure. The evolutions of the boundary saturation and pressure of Model 3 have good agreement with the observed and calculated profiles in Shiozawa and Fujimaki (2004). Before we apply different fluxes q to the three models, we have to know the end times for the simulations. In DiCarlo (2004), the end times for the experiments are not given, but the inner diameter of the tube is given as d = 1.27e-02 [m], then the total volume of water injected into each tube can be calculated by , where x k is the sample point in DiCarlo (2004). Thus we can calculate the end times using T end = V olume qπ(d/2) 2 [s]. To our best knowledge, the τ values are not known for the 20/30 sand, thus we have to first try different values of τ and then find the best match with the experiments. DiCarlo (2004) observed at the highest 2.0e-03 and lowest 1.32e-07 fluxes the saturation profiles are monotonic with distance and no saturation overshoot is observed. In order to find suitable values of τ for the highest and lowest fluxes, in Fig. 6(a) we plot τ as a function of q using a log-log diagram. Realizing the near log-log relationship between τ and q, we set τ = 40 for q = 2.0e-03, τ = 2.0e06 for q = 1.32e-07. The number of nodes used in space is N = 256, the values of τ (τ im ), τ dr , τ s , time steps as well as the end times are presented in Table 5. As can be seen the end times for  simulations are near to the calculated times except when q = 1.32e-07, we will explain this later. For all three models, when the flux q is 1.32e-03, the imbibition front moves quickly while the change in hysteretic capillary pressure is slow. In order to show the overshoot saturation phenomenon, we have to enlarge the interval to [0, 1]. In Fig. 5 we compare the numerical solutions with the experiments. Figs. 5(a), 5(b) and 5(c) show that all three models can obtain saturation overshoots for q = 1.32e-06, q = 1.32e-05, q = 1.32e-04, q = 1.32e-03. From Fig. 5(a) we can see oscillations behind the drainage fronts when q = 1.32e-06, q = 1.32e-05 and q = 1.32e-04, while Fig. 5(b) only present weaker oscillations for q = 1.32e-04 and 1.32e-05 and there is no oscillation behind the drainage fronts in Fig. 5(c). Althouth Table 5 shows that for q = 1.32e-04, 1.32e-05 and .32e-06, the τ dr values are slightly larger than τ s , no oscillation appears behind the drainage fronts. This may be cuased by the hystersis in the capillary pressure, because in Fig. 5(b) the oscillations are weaker than Fig. (5(a)) even without hysteresis in dynamic capillary coefficient.
In Fig. 5(d) we plot the relationship between phases pressure difference and saturation obtained by Model 3 for all fluxes. At the imbibition front, the phases pressure difference is smaller than the equilibrium imbibition pressure. Behind the front, the pressure-saturation follows the equilibrium drainage pressure while finally stops at the tail saturation. This figure shows similar pressure-saturation behaviour as the experiments presented in Fig. 6 in DiCarlo (2007). As can be seen, at the lowest flux, significant pressure overshoot can still be obtained by Model 3. This phenomenon was also observed in the experiment in DiCarlo (2007). Selker et al. (1992) shows that saturation overshoot is associated with pressure overshoot. Thus we guess even that at low flux, Model 3 can still produce saturation overshoot. Let the space interval be [0, 1] and T end = 96000, the saturation and pressure profiles computed by Model 3 are presented in Fig. 5(f). It shows at q = 1.32e-07 very small saturation overshoot appears.    Fig. 5(a), 5(b), 5(c) also show that the computed saturations at the left boundary differ from the experiments especially when q = 1.32e-07. We ascribe this to the limitation of the Brooks-Corey model. Assuming that after a long time, the saturation at the boundary reaches the equilibrium state, we set ∂Sw ∂t | x=0 = 0, ∂Sw ∂x | x=0 = 0. Then we obtain Using the parameters for imbibition process in Table 1 and the Brooks-Corey model in Table 2, from Eq. (47) we get . ( The saturations obtained by Eq. (48), the numerical simulations and experiments at x = 0 are presented in Table 6. The measured volumetric water saturation at q = 1.32e-07 is twice as high as the analytical one. Thus, the end time of the numerical simulation is much shorter than the calculated time.
Since Fig. 5(c) shows more realistic profiles than Fig. 5(a) and 5(b), we will focus on this model and apply more fluxes to test its effectiveness. In Fig. 6(a) we plot τ , T end , ∆t and q used in Fig. 5 with solid triangles. Then we use logarithmic interpolation to get the values of τ , T end and ∆t for intermediate fluxes, these parameters are plotted using open triangles. Fig. 6(b) plots the values of τ as functions of tip and tail saturations. For both imbibition and drainage, the values of τ increase as water saturations decrease. This trend seems agree with the measured data in Manthey et al. (2005); Das and Mirzaei (2012); Mirzaei and Das (2013).
In Fig. 6(c) the computed tip and tail saturations from Model 3 are compared with the measured data in DiCarlo (2004). As the flux increases, the tip saturation increases very fast for flux value in interval [1.0e-06, 1.0e-04], while the tip saturation increases slowly when flux q is above 1.32e-04. For tail saturations, both the experimental data and computed results follow the analytical curve given by Eq. (48) when flux is bigger than 1.0e-06. Fig. 6(d) plots the tip length versus flux obtained by Model 3. For flux values between 1.0e-05 and 1.0e-03, the tip length increases monotonically with the flux. This trend matches Fig. 12 in DiCarlo (2004). Fig. 6(e) presents the phases pressure differences at the imbibition front and the tail obtained by Model 3. For intermediate fluxes, the phases pressure differences at the tail follow the equilibrium drainage capillary pressure, while the phases pressure differences at the imbibition front are below the equilibrium imbibition capillary pressure as a result of the dynamic capillary pressure effect. DiCarlo (2007) defined the overshoot in capillary pressure and the overshoot of phases pressure difference is given as Overshoot(p n − p w ) = Tail(p n − p w ) − Front(p n − p w ).
Fig. 6(f) plots the phases pressure differences as well as the overshoots for different fluxes. The phases pressure differences at the imbibition front are higher at low flux values than at high flux values. At q ≈ 1e-05, the phases pressure difference at the imbibition front reaches a minimum while the pressure overshoot reaches a maximum.  (e) Cyan and magenta lines are the equilibrium imbibition and drainage capillary pressures obtained by the Brooks-Corey model. Cyan and magenta circles are phases pressure differences at the imbibition front and the tail. (f) Cyan and magenta circles denote phases pressure difference at the imbibition front and the tail, red circles are the overshoot.

Conclusion
In this study we applied the Castillo-Grone's mimetic operators and the implicit trapezoidal rule to solve two-phase flow models including dynamic capillary pressure (with constant and hysteretic coefficient) and capillary pressure hysteresis in porous media. Numerical simulations show that the second-order mimetic operators mimics the Green-Gauss-Stokes theorem with high accuracy for all three models. The hysteretic dynamic capillary pressure model with capillary pressure hysteresis produce realistic saturation overshoot and pressure overshoot phenomena as observed in DiCarlo (2004DiCarlo ( , 2007.