Keywords

1 Introduction

The scheme of IMplicit Pressure Explicit Saturation (IMPES) is a conditionally stable which is usually used to solve the two-phase flow in porous media. In the IMPES scheme, the pressure equation is treated implicitly while the saturation is updated explicitly. Hence it takes a very small time step size, in particular with heterogeneous porous media. The IMPES scheme has been improved in several versions (e.g. [1,2,3]). The temporal discretization scheme is considered an important factor that affects the efficiency of numerical reservoir simulators. The application of traditional single-scale temporal scheme is restricted by the rapid changes of the pressure and saturation with capillarity and concentrations if applicable. Applying time splitting strategies has a significant improvement compare to the traditional schemes. Time splitting method has been considered in a number of publications such as [4,5,6,7,8,9,10]. For example, in Refs. [9], an explicit sub-timing scheme is provided. On the other hand, an implicit time-stepping scheme have been proposed by Bhallamudi et al. [7] and Park et al. [8].

In the recent years, nanoparticles have been used in many engineering branches including petroleum applications such as enhanced oil recovery. A model of nanoparticles transport in porous media have been established by Ju and Fan [13] based on particle migration in porous media model [14]. El-Amin et al. [15,16,17] have investigated the problem of nanoparticles transport in porous media. In Ref. [17], they presented an extended model to include mixed relative permeabilities and a negative capillary pressure. Dimensional analysis of the problem of nanoparticles transport in porous media have been presented by El-Amin et al. [18]. The nanoparticles transport in anisotropic porous media have been studied numerically using the multipoint flux approximation by Salama et al. [19]. Also, a numerical simulation of drag reduction effects by hydrophobic nanoparticles adsorption method in water flooding processes has been presented by Chen et al. [20]. The problem of dynamic update of an anisotropic permeability field with nanoparticles transport in porous media has been considering by Chen et al. [21]. In Ref. [22], the authors presented a nonlinear iterative IMPES-IMC scheme to solve the governing system of the nanoparticles transport in subsurface.

In this paper, we propose a time-stepping IMPES-IMC scheme to solve the governing equations of the nanoparticles transport associated with two-phase flow in subsurface. The CCFD method has employed for the spatial discretization. The time-splitting scheme has applied together with the CFL condition to adaptive the time-steps sizes. Finally, some numerical experiments are provided.

2 Modeling and Mathematical Formulation

The mathematical model of the problem under consideration consists of water saturation, Darcy’s law, nanoparticles concentration, deposited nanoparticles concentration on the pore-wall, and entrapped nanoparticles concentration in the pore-throat. Moreover, the variations in both porosity and permeability due to the nanoparticles deposition/entrapment on/in the pores have been taken into consideration. In the following we introduce the governing equations briefly, (for details see Refs. [15,16,17,18,19, 21, 22]:

Momentum Conservation (Darcy’s Law):

$$\begin{aligned} \mathbf {u}_{\alpha } = -\frac{k_{\alpha }}{\mu _\alpha } \nabla \varPhi _{\alpha }, \quad \quad \alpha =w,n, \end{aligned}$$
(1)

where

$$ k_\alpha =k_{r\alpha } \mathbf {K}, \quad \varPhi _{\alpha } = p_{\alpha } + \rho _{\alpha } g \nabla z, \quad \alpha =w,n, $$

Mass Conservation (Saturation Equations):

$$\begin{aligned} \phi \frac{\partial s_{\alpha }}{\partial t} + \nabla \cdot \mathbf {u}_{\alpha } = q_{\alpha }, \quad \alpha =w,n, \end{aligned}$$
(2)

where \(s_w + s_n = 1\), \(\phi \) is the porosity, g is the gravitational acceleration, z is the depth. \(\mathbf {u}_\alpha \) is the velocity, \(\varPhi _{\alpha }\) is the pressure potential, \(p_{\alpha }\) is the pressure, \(\mu _\alpha \) is the viscosity, \(\rho _{\alpha }\) is the density, \(k_{r\alpha }\) is the relative permeability, \(q_{\alpha }\) is the external mass flow rate, \(s_\alpha \) is the saturation; all of the phase \(\alpha \). \(\mathbf {K}\) is the permeability tensor. w stands for the wetting phase (water), and n stands for the non-wetting phase (oil).

Providing the following definitions:

  • The capillary pressure: \(p_c \left( s_w\right) = p_n - p_w\).

  • The total velocity: \(\mathbf {u}_t = \mathbf {u}_w + \mathbf {u}_n\).

  • The flow fraction: \(f_w = \lambda _w/\lambda _t\).

  • The phase mobility: \(\lambda _{\alpha } = k_{r\alpha } / \mu _\alpha \).

  • The total mobility: \(\lambda _t\).

  • The capillary pressure potential: \(\varPhi _c = p_c + \left( \rho _n - \rho _w \right) g \nabla z\).

  • The total source mass transfer: \(q_t= q_w + q_n\).

After some mathematical manipulations and referring to Refs. [11], the pressure equation can be rewritten as,

$$\begin{aligned} \nabla \cdot \mathbf {u}_t = - \nabla \cdot \lambda _t \mathbf {K} \nabla {\varPhi }_w -\nabla \cdot \lambda _n \mathbf {K} \nabla {\varPhi }_c = q_t. \end{aligned}$$
(3)

Because this equation contents the capillary pressure which is a function of saturation, it will be coupled with the following saturation equation to calculate the pressure,

$$\begin{aligned} \phi \frac{\partial s_w}{\partial t} - q_w = - \nabla \cdot \lambda _w \mathbf {K} \nabla {\varPhi }_w. \end{aligned}$$
(4)

However, the saturation is updated using the following form,

$$\begin{aligned} \phi \frac{\partial s_w}{\partial t} - q_w = - \nabla \cdot (f_w \mathbf {u}_a). \end{aligned}$$
(5)

where \(\mathbf {u}_w=f_w \mathbf {u}_a\) and \(\mathbf {u}_a = - \lambda _t \mathbf {K} \nabla {\varPhi }_w\).

Nanoparticles Transport Model

Assuming that the nanoparticles exist only in the water phase of one size interval. So, The transport equation of the nanoparticles in the water phase is given as [12,13,14,15,16,17,18,19, 21, 22],

$$\begin{aligned} \phi \frac{\partial \left( s_w c \right) }{\partial t} + \nabla \cdot \left( \mathbf {u}_w c - D\nabla c \right) = R+Q_c, \end{aligned}$$
(6)

Nanoparticles Surface Deposition

The surface deposition is expressed by,

$$\begin{aligned} \frac{\partial c_{s1}}{\partial t} = \left\{ \begin{array} {cc} \gamma _d |\mathbf {u}_w| c, &{} \mathbf {u}_w \le u_r \\ \\ \gamma _d |\mathbf {u}_w| c - \gamma _e |\mathbf {u}_w - u_r| c_{s1}, &{} \mathbf {u}_w > u_r \end{array}\right. \end{aligned}$$
(7)

Nanoparticles Throat Entrapment

The rate of entrapment of the nanoparticles is,

$$\begin{aligned} \frac{\partial c_{s2}}{\partial t} = \gamma _{pt} |\mathbf {u}_w| c, \end{aligned}$$
(8)

where c is the nanoparticles concentrations. \(c_{s1}\) and \(c_{s2}\) are, respectively, the deposited nanoparticles concentration on the pore surface, and the entrapped nanoparticles concentration in pore throats. \(\tau \) is the tortuosity parameter. \(Q_c\) is the rate of change of particle volume belonging to a source/sink term. \(\gamma _d\) is the rate coefficients for surface retention of the nanoparticles. \(\gamma _e\) is the rate coefficients for entrainment of the nanoparticles. \(u_r\) is the critical velocity of the water phase. where \(\gamma _{pt}\) is the pore throat blocking constants. The diffusion-dispersion tensor is defined by,

$$\begin{aligned} D = \phi s_w \tau D_t, \quad D_t = D^\mathrm{Br} + D^\mathrm{disp} \end{aligned}$$
(9)

where \(D^\mathrm{Br}\) is the Brownian diffusion and \(D^\mathrm{disp}\) is the dispersion coefficient which is defined by [1],

$$\begin{aligned} \phi s_w \tau D^\mathrm{disp} = d_{t,w} |\mathbf {u}_w | \mathbf {I} + (d_{l,w} - d_{t,w}) \frac{\mathbf {u}_w \mathbf {u}_w^T}{|\mathbf {u}_w|} \end{aligned}$$
(10)

Thus,

$$\begin{aligned} D = (\phi s_w \tau D^\mathrm{Br} + d_{t,w} |\mathbf {u}_w | ) \mathbf {I} + (d_{l,w} - d_{t,w}) \frac{\mathbf {u}_w \mathbf {u}_w^T}{|\mathbf {u}_w|} \end{aligned}$$
(11)

where \(d_{l,w}\) and \(d_{t,w}\) are the longitudinal and transverse dispersion coefficients, respectively. R is the net rate of loss of nanoparticles which is defined by,

$$\begin{aligned} R = \frac{\partial c_{s1}}{\partial t} + \frac{\partial c_{s2}}{\partial t} \end{aligned}$$
(12)

Initial and Boundary Conditions

The initial conditions are,

$$\begin{aligned} s_w = s_w^0, \quad c = c_{s1} = c_{s2} = 0 \quad \mathrm{in} \quad \varOmega \quad \quad \mathrm{at} \quad t=0, \end{aligned}$$
(13)

The boundary conditions are given as,

$$\begin{aligned} p_w \left( \mathrm{or} \ p_n\right) = p^D \quad \mathrm{on} \quad \varGamma _D, \end{aligned}$$
(14)
$$\begin{aligned} \mathbf {u}_t \cdot \mathbf {n} = q^N, \quad s_w = S^N, \quad c = c_0, \quad c_{s1} = c_{s2} = 0 \quad \mathrm{on} \quad \varGamma _N. \end{aligned}$$
(15)

where \(\mathbf {n}\) is the outward unit normal vector to \(\partial \varOmega \), \(p^D\) is the pressure on \(\varGamma _D\) and \(q^N\) the imposed inflow rate on \(\varGamma _N\), respectively. \(\varOmega \) is the computational domain such that the boundary \(\partial \varOmega \) is Dirichlet \(\varGamma _D\) and/or Neumann \(\varGamma _N\) boundaries, i.e. \(\partial \varOmega = \varGamma _D \cup \varGamma _N\) and .

3 Multi-scale Time Splitting Method

The concept of time splitting method is to use different time step size for each equation has a time derivative. In the above-described method, the pressure is coupled with the saturation in each time-step. The time-step size for the pressure can be taken larger than the those of saturation and nanoparticle concentrations. So, for the pressure the total time interval [0, T] is divided into \(N_p\), time-steps as \(0=t_0< t_1< \cdots < t_{N_p=T}\). Thus, the time-step length assigned for pressure is, \(\varDelta t^k = t^{k+1} - t^k\). Since the saturation varies more rapidly than the pressure, we use a smaller time-step size for the saturation equation. That is, each interval, \((t^k,t^{k+1}]\), will be divided into \(N_{p,s}\) subintervals as \((t^k,t^{k+1}] = \cup _{l=0}^{N_{p,s}-1} (t^{k,l},t^{k,l+1}]\). On the other hand, as the concentration varies more rapidly than the pressure (and may be saturation), we also use a smaller time-step size for the concentration equations. Thus, we partition each subinterval \((t^{k,l},t^{k,l+1}]\) into \(N_{p,s,c}\) subsubintervals as \((t^{k,l},t^{k,l+1}] = \cup _{m=0}^{N_{p,s,c}-1} (t^{k,l,m},t^{k,l,m+1}]\). Therefore, the system of governing equations, (3), (4), (6), (7) and (8), is solved based on the adaptive time-splitting technique. The backward Euler time discretization is used for the equations of pressure, saturation, concentration and the two volume concentration. We linearized the capillary pressure function, \(\varPhi _c\), in terms of saturation using this formula,

$$\begin{aligned} \varPhi _c \left( s_w^{*} \right) \cong \varPhi _c \left( s_w^{k} \right) + {\varPhi }'_c \left( s_w^{k} \right) \left[ s_w^{k+1} - s_w^{k} \right] , \end{aligned}$$
(16)

where \({\varPhi }'_c\) is derivative of \(\varPhi _c\). The quantity, \([s_w^{k+1} - s_w^{k}]\), can be calculated from the saturation equation,

$$\begin{aligned} s_w^{k+1} - s_w^{k} = \frac{\varDelta t^k}{\phi \left( c_{s1}^{k}, c_{s2}^{k} \right) } \left[ q_w^{k+1} - \nabla \cdot \lambda _t \left( s_w^{k} \right) \mathbf {K} \left( c_{s1}^{k},c_{s2}^{k} \right) \nabla \varPhi _w^{k+1} \right] . \end{aligned}$$
(17)

In addition to the pressure equation,

$$\begin{aligned} - \nabla \cdot \lambda _t \left( s_w^{k} \right) \mathbf {K} \left( c_{s1}^{k},c_{s2}^{k} \right) \nabla \varPhi _w^{k+1} - \nabla \cdot \lambda _n \left( s_w^{k} \right) \mathbf {K} \left( c_{s1}^{k},c_{s2}^{k} \right) \nabla \varPhi _c \left( s_w^{*} \right) = q_t^{k+1}. \end{aligned}$$
(18)

Then, the above coupled system (16), (17) and (18) is solved implicitly to obtain the pressure potential. Therefore, the saturation is updated explicitly with using the upwind scheme for the convection term as,

$$\begin{aligned} \phi \left( c_{s1}^{k}, c_{s2}^{k} \right) \frac{s_w^{k,l+1} - s_w^{k,l}}{\varDelta t^l} + \nabla \cdot \left( f_w^{k} \mathbf {u}_a^{k+1} \right) = q_w^{k,l+1}. \end{aligned}$$
(19)

Therefore, the nanoparticles concentration, deposited nanoparticles concentration on the pore–walls and entrapped nanoparticles concentration in the pore–throats are computed implicitly as follow,

$$\begin{aligned}&\phi \left( c_{s1}^{k},c_{s2}^{k} \right) \frac{s_w^{k+1} c^{k,l,m+1} - s_w^k c^{k,l,m}}{\varDelta t^m} + \nabla \cdot \left\{ \mathbf {u}_w^{k+1} c^{k,l,m+1} - D\left( s_w^{k+1},\mathbf {u}_w^{k+1},c_{s1}^{k},c_{s2}^{k} \right) \right. \nonumber \\&\left. - R \left( \mathbf {u}_w^{k+1},c_{s1}^{k}\right) \right\} \nabla c^{k,l,m+1} = Q_c^{k,l,m+1}, \end{aligned}$$
(20)
$$\begin{aligned}&\frac{c^{k,l,m+1}_{s1} - c^{k,l,m}_{s1}}{\varDelta t^m} = \nonumber \\&\left\{ \begin{array} {cc} \gamma _d |\mathbf {u}_w^{k+1}| c^{k+1}, &{} \mathbf {u}_w^{k+1} \le u_r \\ \\ \gamma _d |\mathbf {u}_w^{k+1}| c^{k+1} - \gamma _e |\mathbf {u}_w^{k+1} - u_r| c^{k,l,m+1}_{s1}, &{} \mathbf {u}_w^{k,l,m+1}> u_r \end{array}\right. \end{aligned}$$
(21)

and,

$$\begin{aligned} \frac{c^{k,l,m+1}_{s2} - c^{k,l,m}_{s2}}{\varDelta t^m} = \gamma _{pt} |\mathbf {u}_w^{k+1}| c^{k,l,m+1} \end{aligned}$$
(22)

Finally, other parameters such as permeability, porosity, \(\lambda _w\),\(\lambda _n\),\(\lambda _t\) and \(f_w\) are updated each loop.

4 Adaptive Time-Stepping

In our algorithm, we have checked the Courant–Friedrichs–Lewy (CFL) condition to guarantee its satisfactory (i.e. CFL<1). In order to achieve this idea, we have to define the following CFLs,

$$\begin{aligned} \mathrm {CFL}_{s,x} = \frac{\mathbf {u}_x \varDelta t^{k,l} }{\varDelta x}, \quad \mathrm {CFL}_{s,y} = \frac{\mathbf {u}_y \varDelta t^{k,l} }{\varDelta y}, \end{aligned}$$
(23)

for saturation equation, and,

$$\begin{aligned} \mathrm {CFL}_{c,x} = \frac{\mathbf {u}_x \varDelta t^{k,l,m} }{\varDelta x}, \quad \mathrm {CFL}_{c,y} = \frac{\mathbf {u}_y \varDelta t^{k,l,m} }{\varDelta y}, \end{aligned}$$
(24)

and for concentration equations. It may be noted that the CFL depends on the ratio \(\varDelta t/\varDelta x\) which can be fixed at larger time-steps and larger mesh-size. Thus, when we use a lager domain (larger mesh size) then we can use larger time step size. In the code, the initial time-step for the saturation equation is taken as the pressure time-step, i.e., \(\varDelta t^{k,0} = \varDelta t^{k}\), and the initial time-step for the concentration equation is taken as the saturation time-step, i.e., \(\varDelta t^{k,l,0} = \varDelta t^{k,l}\). Then, we check if CFL\(_{s,x}>1\) or CFL\(_{s,y}>1\), the saturation time-step will be divided by 2 and the CFL\(_{s,x}\) and CFL\(_{s,y}\) will be recalculated. This procedure will be repeated until satisfying the condition CFL\(_{s,x}<1\) and CFL\(_{s,y}<1\), then the final adaptive saturation time-step will be obtained. Similarly, we check if the condition, CFL\(_{c,x}>1\) or CFL\(_{c,y}>1\) is satisfied, the concentration time-step will be divided by 2 and therefore, we recalculate both CFL\(_{c,x}\) and CFL\(_{c,y}\). We repeat this procedure to reach the condition CFL\(_{c,x}<1\) and CFL\(_{c,y}<1\), then we obtain the final adaptive concentration time-step.

5 Numerical Tests

In order to examine the performance of the current scheme, we introduce some numerical examples in this section. Firstly, we introduce the required physical parameters used in the computations. Then, we study the performance of the scheme by introducing some results for the adaptive time steps based on values of the corresponding Courant–Friedrichs–Lewy (CFL). Then we present some results for the distributions of water saturation and nanoparticles concentrations. In this study, given the normalized wetting phase saturation, \( S = (s_w - s_{wr})/(1-s_{nr} - s_{wr}), 0 \le S \le 1, \) the capillary pressure is defined as, \( p_c = - p_e \log S, \) and the relative permeabilities are defined as, \(k_{rw} = k^0_{rw} S^2\), \(k_{rn} = k^0_{rn} \left( 1 - S\right) ^2\), \(k^0_{rw} = k_{rw} \left( S = 1\right) \), \(k^0_{rw} = k_{rw} \left( S = 1\right) \). \(p_e\) is the capillary pressure parameter, \(s_{wr}\) is the irreducible water saturation and \(s_{nr}\) is the residual oil saturation after water flooding. The values and units of the physical parameters are inserted in Table 1.

In Table 2, some error estimates for water saturation, nanoparticles concentration, and deposited nanoparticles concentration, are provided for various values of the number of time steps k. The reference case (\(s^{n+1}_{w,r}=0.9565\), \(c^{n+1}_{r}=0.0360\), and \(c^{n+1}_{s1,r}=1.3233\times 10^{-7}\)) is calculated at the point (0.06,0.053) for non-adaptive dense mesh of 200 \(\times \) 150 for 0.6 \(\times \) 0.4 m, at 2000 time step. It can be seen from Table 2 that the error decreases as number of number of time steps increases. Table 3 shows the real cost time of running the adaptive scheme for different values of the number of time steps k. One may notice from this table that the real cost time decreases as the number of time steps k increases. It seems more time is needed for the adaptation process.

Table 1. Values of the physical parameters.
Table 2. Error estimates for various values of number of time steps k.
Table 3. Real cost time of running the adaptive scheme for different values of the number of time steps k.
Fig. 1.
figure 1

Adaptive time-step sizes, \(\varDelta t^l\) and \(\varDelta t^m\) against the number of steps of the outer loop k and the number of the inner loops l and m: Case 1 (k = 50).

We use a real permeability map of dimensions 120 \(\times \) 50, which is extracted from Ref. [23]. The permeabilities vary in a large scope and are they highly heterogeneous. We consider a domain of size 40m \(\times \)16m which is discretized into \(120\times 50\) uniform rectangles grids. The injection rate was 0.01 Pore-Volume-Injection (PVI) and continued the calculation until 0.5 PV. In this example, we choose the number of steps of the outer loop to be: (Case 1, k = 50, Fig. 1; Case 1, k = 100, Fig. 1; Case 1, k = 200, Fig. 3). In these figures, the adaptive time-step sizes, \(\varDelta t^l\) and \(\varDelta t^m\) are plotted against the number of steps of the outer loop k and the number of the inner loops l and m. It can be seen from these figures (Figs. 1, 2 and 3) that the behavior of adaptive \(\varDelta t^l\) and \(\varDelta t^m\) are very similar. This may be because the velocity is large and dominate the CFL. Also, both \(\varDelta t^l\) and \(\varDelta t^m\) start with large values then they gradually become smaller and smaller as k increases. On the other hand, for the first two cases when k = 50, 100, \(\varDelta t^l\) and \(\varDelta t^m\) are small when l and m are small, then they increase to reach a peak, then they are gradually decreasing. However, in the third case when k = 200, both \(\varDelta t^l\) and \(\varDelta t^m\) start with large values then they gradually become smaller and smaller as l and m, respectively, increases.

Fig. 2.
figure 2

Adaptive time-step sizes, \(\varDelta t^l\) and \(\varDelta t^m\) against the number of steps of the outer loop k and the number of the inner loops l and m: Case 2 (k = 100).

Fig. 3.
figure 3

Adaptive time-step sizes, \(\varDelta t^l\) and \(\varDelta t^m\) against the number of steps of the outer loop k and the number of the inner loops l and m: Case 3 (k = 200).

Fig. 4.
figure 4

Variations of saturation of the real heterogenous permeability case.

Fig. 5.
figure 5

Variations of nanoparticles concentration of the real heterogenous permeability case.

Fig. 6.
figure 6

Variations of deposited nanoparticles concentration of the real heterogenous permeability case.

Variations of saturation of the heterogenous permeability are shown in Fig. 4. It is noteworthy that the distribution for water saturation is discontinuous due to the high heterogeneity. Thus, for example, we may note higher water saturation at higher permeability regions. Similarly, one may observe the discontinuity in the nanoparticles concentration as illustrated in Fig. 5. The contrast of the nanoparticles concentration distribution arisen from the heterogeneity of the medium can also be noted here. Moreover, the contours of the deposited nanoparticles concentration are shown in Fig. 6. The behavior of the entrapped nanoparticles concentration \(c_{s2}\) is very similar to the behavior of \(c_{s1}\) because of the similarity in their governing equation.

6 Conclusions

In this paper, we introduce an efficient time-stepping scheme with adaptive time-step sizes based on the CFL calculation. Hence, we have calculated the CFL\(_{s,x}\), CFL\(_{s,y}\), CFL\(_{c,x}\) and CFL\(_{c,y}\) at each substep and check if the CFL condition is satisfied (i.e. CFL\(<1\)). We applied this scheme with the IMES-IMC scheme to simulate the problem of nanoparticles transport in two-phase flow in porous media. The model consists of five differential equations, namely, pressure, saturation, nanoparticles concentration, deposited nanoparticles concentration on the pore-walls, and entrapped nanoparticles concentration in pore–throats. The capillary pressure is linearized and used to couple the pressure and the saturation equations. Then, the saturation equation is solved explicitly to update the saturation at each step. Therefore, the nanoparticles concentration equation is treated implicitly. Finally, the deposited nanoparticles concentration on the pore-walls equation, and entrapped nanoparticles concentration in pore-throats equation are solved implicitly. The CCFD method was used to discretize the governing equations spatially. In order to show the efficiency of the proposed scheme, we presented some numerical experiments. The outer pressure time-step size was selected then the inner ones, namely, the saturation subtime-step and the concentration subtime-step were calculated and adaptive by the CFL condition. We presented three cases with different values of the outer pressure time-step. Moreover, distributions of water saturation, nanoparticles concentration, and deposited nanoparticles concentration on pore-wall are shown in graphs.