# Viscous stabilization of 2D drainage displacements with trapping

## Abstract

We investigate the stabilization mechanisms due to viscous forces in the invasion front during drainage displacement in two-dimensional porous media using a network simulator. We find that in horizontal displacement the capillary pressure difference between two different points along the front varies almost linearly as function of height separation in the direction of the displacement. The numerical result supports arguments taking into account the loopless displacement pattern where nonwetting fluid flow in separate strands (paths). As a consequence, we show that existing theories developed for viscous stabilization, are not compatible with drainage when loopless strands dominate the displacement process.

Immiscible displacement of one fluid by another fluid in porous media generates front structures and patterns ranging from compact to ramified and fractal . When a nonwetting fluid displaces a wetting fluid (drainage) at low injection rate, the nonwetting fluid generates a pattern of fractal dimension equal to the cluster formed by invasion percolation . The displacement is controlled solely by the capillary pressure, that is the pressure difference between the two fluids across a pore meniscus. At high injection rate and when the viscosity of the nonwetting fluid is higher or equal to the viscosity of the wetting fluid, the width of the displacement front stabilizes and a more compact pattern is generated 

The purpose of the present letter is to investigate the stabilization mechanisms of the front due to viscous forces.To study the stabilization mechanisms we consider two-dimensional (2D) horizontal drainage at different injection rates. Since the displacement is performed within the plane we neglect gravity. We present simulations where we have calculated the capillary pressure difference $`\mathrm{\Delta }P_c`$ between two different pore menisci along the front separated a height $`\mathrm{\Delta }h`$ in the direction of the displacement \[Fig. 1(a)\]. The simulations are based on a network model that properly describes the dynamics of the fluid-fluid displacement as well as the capillary and viscous pressure buildup . Simulations show that for a wide range of injection rates and different fluid viscosities $`\mathrm{\Delta }P_c`$ varies almost linearly with $`\mathrm{\Delta }h`$ (Figs. 2 and 3). Assuming a power law behavior $`\mathrm{\Delta }P_c\mathrm{\Delta }h^\kappa `$ we find $`\kappa =1.0\pm 0.1`$. This is a surprising result because the viscous force field that stabilizes the front, is non homogeneous due to trapping of wetting fluid behind the front and to the fractal behavior of the front structure.

Based on the observation that the displacement structures are characterized by loopless strands of nonwetting fluid \[Fig. 1(a)\], we also present arguments being supported by our numerical findings. We conjecture that the arguments might affect the behavior of the front width $`w_s`$ as function of the capillary number $`C_a`$. Here $`C_a`$ denotes the ratio between viscous and capillary forces and in the following $`C_aQ\mu _{nw}/\mathrm{\Sigma }\gamma `$, where $`Q`$ is the injection rate, $`\mathrm{\Sigma }`$ is the cross section of the inlet, and $`\mu _{nw}`$ is the viscosity of the nonwetting phase.

In the literature  there has been suggested slightly different scaling behavior of $`w_s`$ as function of $`C_a`$ and a general consensus has not yet been reached. However, none of them consider the evidence observed here that the displacement patterns are loopless and that nonwetting fluid only flows in strands to displace wetting fluid. As a consequence, we show that earlier proposed theories  can not be used to describe drainage when loopless nonwetting strands dominate the displacements.

Before we present the numerical results and the theoretical evidence, we briefly introduce the network model. The model porous medium consists of a square lattice of cylindrical tubes oriented at $`45^{}`$ to the longest side of the lattice \[Fig. 1(a)\]. Four tubes meet at each intersection where we put a node having no volume. The disorder is introduced by (1) assigning the tubes a radius $`r`$ chosen at random inside a defined interval or (2) moving the intersections a randomly chosen distance away from their initial positions. In (1) all tubes have equal length $`d`$ but different $`r`$. (2) results in a distorted square lattice giving the tubes different lengths. Here $`r=d/2\alpha `$ where $`\alpha `$ is the aspect ratio between the tube length and its radius.

The tubes are initially filled with a wetting fluid of viscosity $`\mu _w`$ and a nonwetting fluid of viscosity $`\mu _{nw}\mu _w`$, is injected at constant injection rate $`Q`$ along the bottom row (inlet). The viscosity ratio $`M`$ is defined as $`M\mu _{nw}/\mu _w`$. The wetting fluid is displaced and flows out along the top row (outlet). There are periodic boundary conditions in the orthogonal direction. The fluids are assumed immiscible, hence an interface (a meniscus) is located where the fluids meet in the tubes. The capillary pressure $`p_c`$ of a meniscus is given by $`p_c=(2\gamma /r)\left[1\mathrm{cos}(2\pi x/d)\right]`$. The first term is Young-Laplace law for a cylindrical tube when perfect wetting is assumed and in the second term $`x`$ is the position of the meniscus in the tube ($`0xd`$). Thus, with respect to the capillary pressure we treat the tubes as if they were hourglass shaped with effective radii following a smooth function. By letting $`p_c`$ vary as above, we include the effect of local readjustments of the menisci at pore level  which is important for the description of burst dynamics . The detailed modeling of $`p_c`$ costs computation time, but is necessary in order to properly simulate the capillary pressure behavior along the front.

The volume flux $`q_{ij}`$ through a tube between the $`i`$th and the $`j`$th node is given by Washburn equation : $`q_{ij}=(\sigma _{ij}k_{ij}/\mu _{ij})(p_jp_ip_{c,ij})/d_{ij}`$. Here $`k_{ij}`$ is the permeability of the tube, $`\sigma _{ij}`$ is the average cross section of the tube, $`p_i`$ and $`p_j`$ is the pressures at node $`i`$ and $`j`$ respectively, and $`p_{c,ij}`$ is the sum of the capillary pressures of the menisci inside the tube. A tube partially filled with both liquids, is allowed to contain one or two menisci. Furthermore, $`\mu _{ij}`$ denotes the effective viscosity given by the sum of the volume fractions of each fluid inside the tube multiplied by their respective viscosities. Inserting the above equation for $`q_{ij}`$ into Kirchhoff equations at every node (volume flux conservation), $`_jq_{ij}=0`$, constitutes a set of linear equations which are to be solved for $`p_i`$. The set of equations is solved by using the Conjugate Gradient method with the constraint that $`Q`$ is held fixed. See Refs.  for details on the numerical scheme updating the menisci and solving $`p_i`$.

The front between the two phases is detected by running a Hoshen-Kopelman algorithm  on the lattice. The front width is defined as the standard deviation of the distances between each meniscus along the front and the average front position in the direction of the displacement. $`\mathrm{\Delta }P_c`$ as function of $`\mathrm{\Delta }h`$ is calculated by taking the mean of the capillary pressure differences between all pairs of menisci separated a height $`\mathrm{\Delta }h`$ along the front. The capillary pressure difference between a pair of menisci is calculated by taking the capillary pressure of the meniscus closest to the inlet minus the capillary pressure of the meniscus closest to the outlet \[Fig. 1(a)\].

Figure 2 shows $`\mathrm{\Delta }P_c`$ as function of $`\mathrm{\Delta }h`$ for simulations performed at three different $`C_a`$’s with $`M=100`$ or $`1`$. The simulations with $`M=100`$ were performed on a $`25\times 35`$ nodes lattice with $`\mu _{nw}=10`$ P, $`\mu _w=0.10`$ P, and $`\gamma =30`$ dyn/cm. The disorder was introduced by choosing the tube radii at random in the interval $`0.05dr_{ij}d`$. The tube length was $`d=0.1`$ cm. The simulations with $`M=1`$ were performed on a distorted lattice of $`40\times 60`$ nodes where $`0.02\text{cm}d_{ij}0.18`$ cm and $`r_{ij}=d_{ij}/2\alpha `$ with $`\alpha =1.25`$. Here $`\mu _{nw}=\mu _w=0.5`$ P. To obtain reliable average quantities we did 10–30 simulations at each $`C_a`$ with different sets of random $`r_{ij}`$ or $`d_{ij}`$.

From Fig. 2 we observe that $`\mathrm{\Delta }P_c`$ increases roughly linearly as function of $`\mathrm{\Delta }h`$. At lowest $`C_a`$ no clear stabilization of the front was observed due to the finite size of the system. At higher $`C_a`$ the viscous gradient stabilizes the front. The gradient causes the capillary pressure of the menisci closest to the inlet to exceed the capillary pressure of the menisci lying in the uppermost part. Thus, the menisci closest to the inlet will more easily penetrate a narrow tube compared to menisci further down stream. This will eventually stabilize the front.

To save computation time and thereby be able to study $`\mathrm{\Delta }P_c`$ on larger lattices in the small $`C_a`$ regime, we have generated bond invasion percolation (IP) patterns with trapping on lattices of $`200\times 300`$ nodes. The IP patterns were generated on the bonds in a square lattice with the bonds oriented diagonally at $`45^{}`$. Hence, the bonds correspond to the tubes in our network model. Each bond was assigned a random number $`f_{ij}`$ in the interval $`[0,1]`$. A small stabilizing gradient $`g=0.05`$ was applied, giving an occupation threshold $`t_{ij}`$ of every bond: $`t_{ij}=f_{ij}+gh_{ij}`$ . Here $`h_{ij}`$ denotes the height of the bond above the bottom row. The occupation of bonds started at the bottom row, and the next bond to be occupied was always the bond with the lowest threshold value from the set of empty bonds along the invasion front. The generated IP patterns are similar to the site-bond IP patterns in  and we assume they are statistical equal to structures that would have been obtained in a corresponding complete displacement simulation.

When the IP patterns became well developed with trapped (wetting) clusters of sizes between the bond length and the front width, the tubes in our network model were filled with nonwetting and wetting fluid according to occupied and empty bonds in the IP lattice. Moreover, the radii $`r_{ij}`$ of the tubes were mapped to the random numbers $`f_{ij}`$ of the bonds as $`r_{ij}=[0.05+0.95(1f_{ij})]d`$. Thus, $`0.05dr_{ij}d`$ and we set the tube length $`d=0.1\text{cm}`$. Note that $`r_{ij}`$ is mapped to $`1f_{ij}`$ because in our IP algorithm the next bond to be invaded is the one with the lowest threshold value, opposite to the network model, where the widest tube will be invaded first.

After the initiation of the tube network was completed, the network model was started and the simulations were run a limited number of time steps before it was stopped. The number of time steps where chosen sufficiently large to let the menisci along the front adjust according to the viscous pressure set up by the injection rate.

Totally, we generated four IP patterns with different sets of $`f_{ij}`$ and every pattern was loaded into the network model. The result of the calculated $`\mathrm{\Delta }P_c`$ versus $`\mathrm{\Delta }h`$ is shown in Fig. 3 for $`C_a=9.5\times 10^5`$ and $`M=100`$. If we assume a power law $`\mathrm{\Delta }P_c\mathrm{\Delta }h^\kappa `$, we find $`\kappa =1.0\pm 0.1`$. The slope of the straight line in Fig. 3 is 1.0. We have also calculated $`\mathrm{\Delta }P_c`$ for $`C_a=2\times 10^6`$ with $`M=1`$ and $`M=100`$ by using one of the generated IP patterns. The result of those simulations is consistent with Fig. 3.

Wilkinson  was the first to use percolation theory to deduce a power law between $`w_s`$ and $`C_a`$ when only viscous forces stabilize the front. In 3D, where trapping of wetting fluid is assumed to be of little importance, he suggested $`w_sC_{a}^{}{}_{}{}^{\alpha }`$ and $`\alpha =\nu /(1+t\beta +\nu )`$. Here $`t`$ is the conductivity exponent and $`\beta `$ is the order parameter exponent in percolation. Blunt et al.  used a similar approach, however, they found $`\alpha =\nu /(1+t+\nu )`$ in 3D. This is identical to the result of Lenormand  discussing limits of fractal patterns between capillary fingering and stable displacement in 2D porous media. Blunt et al. also deduced a scaling relation for the pressure drop $`\mathrm{\Delta }P_{nw}`$ across a height difference $`\mathrm{\Delta }h`$ in the nonwetting phase of the front and found $`\mathrm{\Delta }P_{nw}\mathrm{\Delta }h^{t/\nu +1}`$. Later on, Xu et al.  used the arguments of Gouyet et al.  and Wilkinson  to show that $`\mathrm{\Delta }P_{nw}\mathrm{\Delta }h^{t/\nu +d_\text{E}1\beta /\nu }`$, where $`d_\text{E}`$ is the Euclidean dimension of the space in which the front is embedded. They also argued that $`\mathrm{\Delta }P_c=\mathrm{\Delta }P_{nw}\mathrm{\Delta }P_w`$ where $`\mathrm{\Delta }P_w`$ denoting the pressure drop in the wetting phase of the front, is linearly dependent on $`\mathrm{\Delta }h`$ due to the compact displaced fluid \[see Fig. 1(a)\]. Thus, the result of Xu et al. would in 2D predict $`\mathrm{\Delta }P_c\mathrm{\Delta }h^{1.9}`$ where we have used $`t=1.3`$, $`\nu =4/3`$, $`\beta =5/36`$, and $`d_\text{E}=2`$. Our simulations give $`\mathrm{\Delta }P_c\mathrm{\Delta }h^\kappa `$ and $`\kappa =1.0\pm 0.1`$. Below we present an alternative view on the displacement pattern from those first suggested by Wilkinson. The alternative view is based upon the loopless nonwetting strands and is supported by our numerical result.

The simulated displacement patterns show that the nonwetting fluid contains no closed loops \[Fig. 1(a)\] because wetting fluid may be trapped in single tubes, due to volume conservation \[Fig. 1(b)\]. Because of fluid trapping in single tubes, the invading fluid flows in separate strands that cannot coalesce. We note that the definition in in Fig. 1(b) can be easily generalized to 3D , since increasing the coordination number of the lattice does not change the trapping rule. Therefore, we expect loopless patterns to develop in 3D lattices and our arguments that we present below should apply there too. We also note that trapping of wetting fluid is more difficult in real porous media due to a more complex topology of pores and throats there. Loopless IP patterns have earlier been observed in Refs. .

From Fig. 1(a) we may separate the displacement pattern into two parts: one consisting of the frontal region continuously covering new tubes, and the other consisting of the more static structure behind the front. The frontal region is supplied by nonwetting fluid through strands connecting the frontal region to the inlet. When the strands approach the frontal region they are more likely to split. Since we are dealing with a square lattice, a splitting strand may create either two or three new strands. As the strands proceed further into the frontal region they split again and again and eventually they cover the frontal region completely \[see Fig. 1(a)\].

On IP patterns without loops  the length $`l`$ of the minimum path between two points separated an Euclidean distance $`R`$ scales like $`lR^{D_s}`$ where $`D_s`$ is the fractal dimension of the shortest path. We assume that the displacement pattern of the frontal region for length less than the correlation length (in our case $`w_s`$) is statistically equal to IP patterns in . Therefore, the length of a strand in the frontal region is proportional to $`\mathrm{\Delta }h^{D_s}`$ when $`\mathrm{\Delta }h`$ is less than $`w_s`$. If we assume that on the average every tube in the lattice has same mobility ($`k_{ij}/\mu _{ij}`$), this causes the fluid pressure within a single strand to drop like $`\mathrm{\Delta }h^{D_s}`$ as long as the strand does not split. When the strand splits volume conservation causes the volume fluxes through the new strands to be less than the flux in the strand before it splits. Hence, following a path where strands split will cause the pressure to drop as $`\mathrm{\Delta }h^\kappa `$ where $`\kappa D_s`$.

From the above arguments we conclude that the pressure drop $`\mathrm{\Delta }P_{nw}`$, in the nonwetting phase of the frontal region (that is the strands) should scale as $`\mathrm{\Delta }P_{nw}\mathrm{\Delta }h^\kappa `$ where $`\kappa D_s`$. In 2D two different values for $`D_s`$ have been reported: $`D_s=1.22`$  and $`D_s=1.14`$ . Both values are consistent with our simulations finding $`\kappa =1.0\pm 0.1`$.

The evidence that $`\kappa 1.0`$ may influence the scaling of $`w_s`$ as function of $`C_a`$. At low $`C_a`$ simulations show that $`\mathrm{\Delta }\widehat{P}_cC_a\mathrm{\Delta }h^{1.0}`$ . Here $`\mathrm{\Delta }\widehat{P}_c`$ denotes the capillary pressure difference when the front is stationary. That means, $`\mathrm{\Delta }\widehat{P}_c`$ excludes situations where nonwetting fluid rapidly invades new tubes due to local instabilities. At sufficiently low $`C_a`$ the displacement can be mapped to percolation giving $`\mathrm{\Delta }\widehat{P}_cff_c\xi ^{1/\nu }`$ . Here $`f`$ is the occupation probability of the bonds, $`f_c`$ is the percolation threshold, and $`\xi w_s`$ is the correlation length. By combining the above relations, we obtain $`w_sC_{a}^{}{}_{}{}^{\alpha }`$ where $`\alpha =\nu /(1+\nu \kappa )`$. In 2D $`\nu =4/3`$ and inserting $`\kappa =1.0`$ gives $`\alpha 0.57`$. At high $`C_a`$ we expect a crossover to another type of behavior since it is not clear if the mapping to percolation  is valid there. We note that Wilkinson’s result  gives $`\alpha 0.38`$ in 2D.

In summary we conclude that $`\mathrm{\Delta }P_c\mathrm{\Delta }h^\kappa `$ where our simulations gives $`\kappa =1.0\pm 0.1`$. By describing the displacement structure in terms of loopless strands  we have argued that $`\kappa D_s`$, where $`D_s`$ is the fractal dimension of the shortest path between two points on IP patterns without loops. In 2D two values of $`D_s`$ has been reported ($`1.14`$  and $`1.22`$ ) and both are consistent with our numerical result $`\kappa 1.0`$. We conclude that earlier suggested theories  are not compatible in situations where a loopless pattern with nonwetting strands dominate the displacement. We have also shown that $`\alpha `$ in $`w_sC_{a}^{}{}_{}{}^{\alpha }`$, may be influenced by the evidence that $`\kappa D_s`$. Work is in progress to investigate our arguments in 3D and the effect of loops on $`\kappa `$.

The authors thank J. Feder and E. G. Flekkøy for valuable comments. The work is supported by the Norwegian Research Council (NFR) through a “SUP” program and we acknowledge them for a grant of computer time.